matlab问题

zzz.dan

我的代码如下：
4 C% W& r4 [8 [2 Vclose all;
1 f5 L: v& X$ K# B" Sclear all;
L=1; EI=5; T=2; rho=2; M=6;
0 L9 Q4 }! \9 S* O7 j5 }" C% Ynt=80000; nx=20; tf=15; Ttr=50;
dx=L/nx; dt=tf/nt;
+ F; X8 r, x6 H# zx = linspace(0,L,nx+1); t_tr = linspace(0,tf,Ttr+1);
w=zeros(nx+1,nt+1);
: G) w" q! C8 e* ]$ Lw3D_free=zeros(Ttr+1,nx+1);
" w) e2 W' f2 b2 W% ?d=zeros(nt+1,1);
- A6 t6 g# A0 Y; Z# @* S+ Bfor j=1:nt+1
( \$ N  T  m1 v6 T9 O& v    d(j) = 1 + sin(2*pi*(j-1)*dt) + cos(3*pi*(j-1)*dt) ;
% `- D" y0 X4 A( [+ wend
- v/ L, r4 c. q( q' r6 o& TBr=90;Bl=-90;k=100;k1=1;k2=4;r=1;m=1;
w=zeros(nx+1,nt+1); w3D_model=zeros(Ttr+1,nx+1);
2 [0 W7 ?  I6 ]u_model=zeros(nt+1,1);u3D_model=zeros(Ttr+1,1);
: J  s& y8 q% q) O; I5 Iua=zeros(nt+1,1); d_estimate=zeros(nt+1,1);
d1=zeros(nt+1,1);
dbar=max(abs(d));
d_estimate(1:2)= dbar;
for i=1:nx+1
$ B6 m6 z  g: C4 s  y8 X  w(i,1)=(i-1)*dx;
  w(i,2)=w(i,1);
8 I( r# d4 a6 B( t# p4 A. T  m" yend
" K5 W4 I0 Q( `" Bw(1,=0;w(2,=w(1,;
1 j' w  l0 r& v3 pw3D_model(1,:)=w(:,1);w3D_model(2,:)=w(:,2);
+ C. e# _7 E" _; L2 w) j. O% Nfor j=3:nt+1
9 T1 ]) C- h5 u    for i=3:nx-1
        wxx=(w(i+1,j-1)-2*w(i,j-1)+w(i-1,j-1))/dx^2;
1 ]7 c: C3 I& e+ e4 y6 {3 O        wxxxx=(w(i+2,j-1)-4*w(i+1,j-1)+6*w(i,j-1)-4*w(i-1,j-1)+w(i-2,j-1))/dx^4;
        w(i,j)=2*w(i,j-1)-w(i,j-2)+(T*wxx-EI*wxxxx)*dt^2/rho;
    end

6 o, J- j  K# F( G    wxl=(w(nx+1,j-1)-w(nx,j-1))/dx;
+ H  U, }' j3 l" Z. D) b$ S    wxxxl=(w(nx+1,j-1)-3*w(nx,j-1)+3*w(nx-1,j-1)-w(nx-2,j-1))/dx^3;
4 x; L% Z# A2 M5 r    wxxxtl=(w(nx+1,j-1)-3*w(nx,j-1)+3*w(nx-1,j-1)-w(nx-2,j-1)-w(nx+1,j-2)+3*w(nx,j-2)-3*w(nx-1,j-2)+w(nx-2,j-2))/(dt*dx^3);
    wxtl=(w(nx+1,j-1)-w(nx,j-1)-w(nx+1,j-2)+w(nx,j-2))/(dx*dt);
: @+ g) }% c2 @    wtl=(w(nx+1,j-1)-w(nx+1,j-2))/dt;
    ua(j)=wtl+k1*wxl-k2*wxxxl;
    d_estimate(j)= (abs(ua(j))+d_estimate(j-1))/(1+r*dt);
    u_model(j)= (-EI*wxxxl+T*wxl-k1*M*wxtl+k2*M*wxxxtl-k*ua(j)-sign(ua(j))*d_estimate(j))/m;
    dv=(u_model(j)-u_model(j-1))/dt;
9 h/ \  E1 L- B' E6 @, c    if  dv>0
        d0=-m*Br;
    else
! }9 n% ], h+ A2 w2 p+ {7 ^        d0=-m*Bl;
3 `. M( w8 X( L8 v/ s5 y    end
7 c* c# ~% S1 }2 N2 \    d1(j)=d0+d(j);
    w(nx+1,j)=2*w(nx+1,j-1) - w(nx+1,j-2) + (m*u_model(j-1)+d1(j-1)+EI*wxxxl-T*wxl)*dt^2/M;
# t9 [/ z( {2 n    w(nx,j)=(w(nx+1,j)+w(nx-1,j))/2;
    if mod(j-1,(nt/Ttr))==0
        w3D_model(1+(j-1)*Ttr/nt,:)=w(:,j);
* s" p& K. O- ~6 A0 Y        u3D_model(1+(j-1)*Ttr/nt)=u_model(j);
    end
end
/ \: i" a: P3 _7 _: Y+ ]figure(1);
suRF(x,t_tr,w3D_model);
1 d5 l, m4 Z6 Vtitle('Displacement of beam with robust boundary control');
5 X2 W* h1 J% l5 i2 Q: y% Wylabel('Time (s)','Fontsize',12);
xlabel('x(m)','Fontsize',12);
6 v/ B3 c; H' K* `9 H1 u9 ?1 D9 lzlabel('w(x,t)(m)','Fontsize',12);
figure(2);
plot(t_tr,u3D_model,'r')
xlabel('time');
ylabel('u');

5 ?7 ~: e3 a+ B! m) _! Y& S+ B% u; d
. g6 J# C( H: S! G以上我的代码没出问题，但图像出不来,求助!
! N! S0 E0 c8 X  k( }/ T

小小鲁班

w(nx+1,j)=2*w(nx+1,j-1) - w(nx+1,j-2) + (m*u_model(j-1)+d1(j-1)+EI*wxxxl-T*wxl)*dt^2/M;
这行代码近似看成w(nx+1,j)=2*w(nx+1,j-1)，每次迭代乘以二，当迭代超过1000多次就会出现inf，进而出现NAN
修改意见可以减少迭代次数，或者改用sym符号计算试试

Uifhjvv

来学习一下

勇往直前11

zaiyiaaaa