matlab问题

zzz.dan

我的代码如下：
0 D( H; ^) u1 F$ X5 z( R4 mclose all;
clear all;
* P$ v, w0 n" A; d' T$ q& p- Z% P' UL=1; EI=5; T=2; rho=2; M=6;
nt=80000; nx=20; tf=15; Ttr=50;
6 o3 D% U; m4 I/ f  Fdx=L/nx; dt=tf/nt;
x = linspace(0,L,nx+1); t_tr = linspace(0,tf,Ttr+1);
w=zeros(nx+1,nt+1);
% e  A3 f" Q4 f6 k2 u4 [) cw3D_free=zeros(Ttr+1,nx+1);
d=zeros(nt+1,1);
for j=1:nt+1
0 t! m7 }# i% F+ ?1 p    d(j) = 1 + sin(2*pi*(j-1)*dt) + cos(3*pi*(j-1)*dt) ;
+ H/ g8 B  ^6 H6 j: W1 `end
Br=90;Bl=-90;k=100;k1=1;k2=4;r=1;m=1;
# \: P  P0 K- fw=zeros(nx+1,nt+1); w3D_model=zeros(Ttr+1,nx+1);
u_model=zeros(nt+1,1);u3D_model=zeros(Ttr+1,1);
ua=zeros(nt+1,1); d_estimate=zeros(nt+1,1);
d1=zeros(nt+1,1);
dbar=max(abs(d));
; L1 Z0 w* G4 P1 td_estimate(1:2)= dbar;
for i=1:nx+1
& R( ?5 e: X  A# ~6 h  w(i,1)=(i-1)*dx;
, h/ b- U/ E' {0 `4 m+ X) }  w(i,2)=w(i,1);
end
w(1,=0;w(2,=w(1,;
( X( `, R' }8 \# `8 O0 gw3D_model(1,:)=w(:,1);w3D_model(2,:)=w(:,2);
8 d6 N' X4 f' X2 Zfor j=3:nt+1
! _. v: d" I: o" Z    for i=3:nx-1
        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;
    end

8 u7 E6 w1 }, k, x7 a    wxl=(w(nx+1,j-1)-w(nx,j-1))/dx;
( X4 d; C5 W% q$ X% ?2 p    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 S4 w7 Z# ]% f) k    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);
0 T6 E1 J6 H/ i9 Z  i" y5 @8 p    wtl=(w(nx+1,j-1)-w(nx+1,j-2))/dt;
! E% D- }% q. U9 W3 @) N    ua(j)=wtl+k1*wxl-k2*wxxxl;
0 R$ P3 w- ~1 G: n& R" O    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;
$ ^/ y) w1 T' l4 P3 a    dv=(u_model(j)-u_model(j-1))/dt;
    if  dv>0
        d0=-m*Br;
* ]- u- x* i& |$ A# w    else
        d0=-m*Bl;
    end
! F/ k& j8 q; G, }  P    d1(j)=d0+d(j);
' k; K0 r  u' l    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;
1 C1 }7 \9 z; F  f    w(nx,j)=(w(nx+1,j)+w(nx-1,j))/2;
. F8 I8 r; d: c1 o    if mod(j-1,(nt/Ttr))==0
        w3D_model(1+(j-1)*Ttr/nt,:)=w(:,j);
        u3D_model(1+(j-1)*Ttr/nt)=u_model(j);
    end
7 p( a" A" l( hend
figure(1);
suRF(x,t_tr,w3D_model);
title('Displacement of beam with robust boundary control');
ylabel('Time (s)','Fontsize',12);
xlabel('x(m)','Fontsize',12);
zlabel('w(x,t)(m)','Fontsize',12);
figure(2);
plot(t_tr,u3D_model,'r')
& ?8 x# S- [; w& `5 Fxlabel('time');
ylabel('u');

1 f5 j: i- Q* Y5 O' d6 O
以上我的代码没出问题，但图像出不来,求助!
% Q3 z: g: U$ r" e' R: O. s

小小鲁班

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;
2 z! S3 t5 B4 O这行代码近似看成w(nx+1,j)=2*w(nx+1,j-1)，每次迭代乘以二，当迭代超过1000多次就会出现inf，进而出现NAN
修改意见可以减少迭代次数，或者改用sym符号计算试试

Uifhjvv

来学习一下

勇往直前11

zaiyiaaaa