matlab问题

zzz.dan

我的代码如下：
8 C9 H$ M: `2 C: B/ K2 r7 Dclose all;
clear all;
L=1; EI=5; T=2; rho=2; M=6;
nt=80000; nx=20; tf=15; Ttr=50;
dx=L/nx; dt=tf/nt;
5 C* W. s- ~$ S4 A( ^x = linspace(0,L,nx+1); t_tr = linspace(0,tf,Ttr+1);
  y6 @) l! J6 T+ @+ j( bw=zeros(nx+1,nt+1);
7 k) q( y0 a7 Y/ c/ Z4 Ew3D_free=zeros(Ttr+1,nx+1);
5 ]1 i1 S3 j1 W. ad=zeros(nt+1,1);
for j=1:nt+1
' `) B4 u0 B! f' d. q/ a    d(j) = 1 + sin(2*pi*(j-1)*dt) + cos(3*pi*(j-1)*dt) ;
5 [5 R) \0 \8 P' }' ^$ U- m  o- b) E2 @end
Br=90;Bl=-90;k=100;k1=1;k2=4;r=1;m=1;
" x' S' @: v+ |0 V" [! aw=zeros(nx+1,nt+1); w3D_model=zeros(Ttr+1,nx+1);
u_model=zeros(nt+1,1);u3D_model=zeros(Ttr+1,1);
+ @$ Y1 }: B; y$ V; }ua=zeros(nt+1,1); d_estimate=zeros(nt+1,1);
; G, d* _3 w3 A, H0 l4 e. Bd1=zeros(nt+1,1);
dbar=max(abs(d));
d_estimate(1:2)= dbar;
for i=1:nx+1
- G- e5 k3 i8 Z2 P. r/ m/ W3 d. h0 H  w(i,1)=(i-1)*dx;
  w(i,2)=w(i,1);
end
w(1,=0;w(2,=w(1,;
. n5 P7 |2 }, Q% o0 V/ Dw3D_model(1,:)=w(:,1);w3D_model(2,:)=w(:,2);
for j=3:nt+1
    for i=3:nx-1
; u' m, A1 s, f3 b/ `/ q$ w        wxx=(w(i+1,j-1)-2*w(i,j-1)+w(i-1,j-1))/dx^2;
        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;
2 r# i  ~3 i5 H! X+ p    end
- n  R/ `' o  Z3 c  |% g
* t: U5 U8 |9 O7 L, s    wxl=(w(nx+1,j-1)-w(nx,j-1))/dx;
9 I6 j6 a$ o  C* a    wxxxl=(w(nx+1,j-1)-3*w(nx,j-1)+3*w(nx-1,j-1)-w(nx-2,j-1))/dx^3;
    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);
7 ^, F1 v1 V# p( Q/ j    wtl=(w(nx+1,j-1)-w(nx+1,j-2))/dt;
9 v/ \, N0 O7 I7 b    ua(j)=wtl+k1*wxl-k2*wxxxl;
/ j3 K+ a4 o- I" S    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;
    if  dv>0
        d0=-m*Br;
3 w; L) p4 O5 R2 E! W    else
8 ~& L' ], E+ S& P# v        d0=-m*Bl;
    end
    d1(j)=d0+d(j);
6 f/ v- s7 }0 L) {4 u. |2 r+ `    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;
) o- W/ I5 F: d8 H7 s' N    w(nx,j)=(w(nx+1,j)+w(nx-1,j))/2;
    if mod(j-1,(nt/Ttr))==0
. ?5 ?& Z" b; e! q. C        w3D_model(1+(j-1)*Ttr/nt,:)=w(:,j);
        u3D_model(1+(j-1)*Ttr/nt)=u_model(j);
' E9 N" k# _6 a' W1 O3 h  o4 f6 U    end
; i+ F8 s" h4 _! ]. fend
figure(1);
1 h# t0 p1 }/ V4 {suRF(x,t_tr,w3D_model);
6 u4 o$ R$ |8 a+ c4 d2 v2 }title('Displacement of beam with robust boundary control');
ylabel('Time (s)','Fontsize',12);
xlabel('x(m)','Fontsize',12);
7 l& x. Y8 |. k7 J8 R2 szlabel('w(x,t)(m)','Fontsize',12);
- X( G: }# R) s  J1 |figure(2);
plot(t_tr,u3D_model,'r')
  r; s. g& |1 o8 N5 @; C& y* {xlabel('time');
; @: z9 Z+ s2 c: C+ eylabel('u');
1 [$ ^  U: A& S" c
5 a. y+ K$ O( `& u7 g3 A) S
以上我的代码没出问题，但图像出不来,求助!
. \/ v/ v+ W: |* L4 D8 k& k

zaiyiaaaa

勇往直前11

Uifhjvv

来学习一下

小小鲁班

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
3 U5 h0 d2 B- x" R, k6 p9 b修改意见可以减少迭代次数，或者改用sym符号计算试试