matlab问题

zzz.dan · 发表于 2020-9-27 13:49

EDA365欢迎您登录！

您需要登录才可以下载或查看，没有帐号？注册

x

我的代码如下：
close 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;
x = linspace(0,L,nx+1); t_tr = linspace(0,tf,Ttr+1);
w=zeros(nx+1,nt+1);
w3D_free=zeros(Ttr+1,nx+1);
d=zeros(nt+1,1);
for j=1:nt+1
d(j) = 1 + sin(2*pi*(j-1)*dt) + cos(3*pi*(j-1)*dt) ;
end
Br=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);
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));
d_estimate(1:2)= dbar;
for i=1:nx+1
w(i,1)=(i-1)*dx;
w(i,2)=w(i,1);
end
w(1,

=0;w(2,

=w(1,

;
w3D_model(1,:)=w(:,1);w3D_model(2,:)=w(:,2);
for j=3:nt+1
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

wxl=(w(nx+1,j-1)-w(nx,j-1))/dx;
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);
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;
if  dv>0
      d0=-m*Br;
else
      d0=-m*Bl;
end
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;
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);
      u3D_model(1+(j-1)*Ttr/nt)=u_model(j);
end
end
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')
xlabel('time');
ylabel('u');

以上我的代码没出问题，但图像出不来,求助!

小小鲁班 · 发表于 2020-9-27 14:50

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 · 发表于 2020-9-27 15:28

来学习一下

勇往直前11 · 发表于 2020-9-28 11:08

zaiyiaaaa · 发表于 2020-9-28 14:00

帐号		自动登录	找回密码
密码			注册

matlab问题

EDA365欢迎您登录！

浏览过的版块

推荐内容 /1