可以帮我看看程序写的对不对吗

问题遇到的现象和发生背景

可以帮我看看程序写的对不对吗，能运行不报错，就是数据nan

问题相关代码，请勿粘贴截图

ddvx(1,j)=(vx(1,j)-2vx(1,j-1)+vx(1,j-2))/T^2;%vx(0,t)对t求2阶导
ddvy(1,j)=(vy(1,j)-2vy(1,j-1)+vy(1,j-2))/T^2;%vy(0,t)对t求2阶导
ddvz(1,j)=(vz(1,j)-2*vz(1,j-1)+vz(1,j-2))/T^2;%vz(0,t)对t求2阶导

%qxa是q_x(L,t)对L求一阶偏导 dqx是q_x(L,t)对t求一阶偏导
%qxaa是q_x(L,t)对L求二阶偏导 ddqx是q_x(L,t)对t求二阶偏导
%既跟a有关又跟t有关%
for i=3:nx-2
vxa=(vx(i,j-1)-vx(i-1,j-1))/dx;%vx_a(a,t)一阶导
vya=(vy(i,j-1)-vy(i-1,j-1))/dx;%vx_a(a,t)
vza=(vz(i,j-1)-vz(i-1,j-1))/dx;%vx_a(a,t)
vxaa=(vx(i+1,j-1)-2vx(i,j-1)+vx(i-1,j-1))/dx^2;%二阶导vx_aa(a,t)
vyaa=(vy(i+1,j-1)-2vy(i,j-1)+vy(i-1,j-1))/dx^2;%二阶导vx_aa(a,t)
vzaa=(vz(i+1,j-1)-2*vz(i,j-1)+vz(i-1,j-1))/dx^2;%二阶导vx_aa(a,t)

  dvx(i,j-1)=(vx(i,j-1)-vx(i,j-2))/T;%vx对t求一阶导
  dvy(i,j-1)=(vy(i,j-1)-vy(i,j-2))/T;%
  dvz(i,j-1)=(vz(i,j-1)-vz(i,j-2))/T;%
  
  fx=(1.2+0.2*sin(0.6*pi*j)+0.4*sin(0.4*pi*j)+0.6*sin(0.2*pi*j))*i;%fx(L,t)
  %(1)式：qx(i,j)
  qx(i,j)=T^2*(-ddvx(1,j-1)+(fx-kesei*dvx(i,j-1)+TTx(i,1)*vxa+TT(i,1)*vxaa
          +2*wx(i,1)*(vxa^3+vxa*vya^2)+2*w(i,1)*(3*vxa^2*vxaa+2*vxa*vya*vyaa+vxaa*vya^2)
          +EA*(1.5*vxa^2*vxaa+vza*vxaa+vzaa*vxa+vxa*vya*vyaa+0.5*vya^2*vxaa))/rho(i,1))+2*qx(i,j-1)-qx(i,j-2);

% qxa(i,j-1)=vxa; %vx_a(L,t)=qx_a(L,t)
% qxaa(i,j-1)=vxaa;
dqx(i,j-1)=(qx(i,j-1)-qx(i,j-2))/T;%qx对t求一阶导
vx(i,j)=vx(1,j)+qx(i,j);%v和q之间的关系式

  fy=(1.2+0.2*sin(0.6*pi*j)+0.4*sin(0.4*pi*j)+0.6*sin(0.2*pi*j))*i;%fy(L,t)

%(2)式
qy(i,j)=T^2*(-ddvy(1,j-1)+(fy-keseidvy(i,j-1)+TTx(i,1)vya+TT(i,1)vyaa
+2wx(i,1)(vya^3+vyavxa^2)+2w(i,1)(3vya^2vyaa+2vyavxavxaa+vyaavxa^2)
+EA*(1.5vya^2vyaa+vzavyaa+vzaavya+vyavxavxaa+0.5vxa^2vyaa))/rho(i,1))+2*qy(i,j-1)-qy(i,j-2);
% qya(i,j-1)=vya;
% qyaa(i,j-1)=vyaa;
dqy(i,j-1)=(qy(i,j-1)-qy(i,j-2))/T;%qy对t求一阶导
vy(i,j)=vy(1,j)+qy(i,j);%v和q之间的关系式

  fz=(0.6+0.2*sin(0.6*pi*j)+0.4*sin(0.4*pi*j)+0.6*sin(0.2*pi*j))/15*i;%fz(a,t)
  %(3)式
  qz(i,j)=T^2*(-ddvz(1,j-1)+(fz-kesei*dvz(i,j-1)+EA*(vzaa+vxa*vxaa+vya*vyaa))/rho(i,1))+2*qz(i,j-1)-qz(i,j-2);

% qza(i,j-1)=vza;
% qzaa(i,j-1)=vzaa;
dqz(i,j-1)=(qz(i,j-1)-qz(i,j-2))/T;%qz对t求一阶导
vz(i,j)=vz(1,j)+qz(i,j);%v和q之间的关系式

end


dvx(nx-1,j)=(vx(nx-1,j)-vx(nx-1,j-1))/T;
dvy(nx-1,j)=(vy(nx-1,j)-vy(nx-1,j-1))/T;
dvz(nx-1,j)=(vz(nx-1,j)-vz(nx-1,j-1))/T;

vxa(nx-1,j-1)=(vx(nx-1,j-1)-vx(nx-2,j-1))/dx;%差分法求对p的一阶导数
vxaa(nx-1,j-1)=(vx(i,j-1)-2*vx(i-1,j-1)+vx(i-2,j-1))/dx^2;
vya(nx-1,j-1)=(vy(nx-1,j-1)-vy(nx-2,j-1))/dx;%差分法求对p的一阶导数
vyaa(nx-1,j-1)=(vy(i,j-1)-2*vy(i-1,j-1)+vy(i-2,j-1))/dx^2;
vza(nx-1,j-1)=(vz(nx-1,j-1)-vz(nx-2,j-1))/dx;%差分法求对p的一阶导数
vzaa(nx-1,j-1)=(vz(i,j-1)-2*vz(i-1,j-1)+vz(i-2,j-1))/dx^2;

  qx(nx-1,j)=T^2*(-ddvx(1,j-1)+(fx-kesei*dvx(nx-1,j)+TTx(nx-1,1)*vxa(nx-1,j-1)...
             +TT(nx-1,1)*vxaa(nx-1,j-1)+2*wx(nx-1,1)*(vxa(nx-1,j-1)^3+vxa(nx-1,j-1)*vya(nx-1,j-1)^2)
             +2*w(nx-1,1)*(3*vxa(nx-1,j-1)^2*vxaa(nx-1,j-1)+2*vxa(nx-1,j-1)*vya(nx-1,j-1)*vyaa(nx-1,j-1)+vxaa(nx-1,j-1)*vya(nx-1,j-1)^2)
             +EA*(1.5*vxa(nx-1,j-1)^2*vxaa(nx-1,j-1)+vza(nx-1,j-1)*vxaa(nx-1,j-1)+vzaa(nx-1,j-1)*vxa(nx-1,j-1)
             +vxa(nx-1,j-1)*vya(nx-1,j-1)*vyaa(nx-1,j-1)+0.5*vya(nx-1,j-1)^2*vxaa(nx-1,j-1)))/rho(nx-1,1))+2*qx(nx-1,j-1)-qx(nx-1,j-2);%(1)

% qxa(nx-1,j-1)=vxa(nx-1,j-1);
% qxaa(nx-1,j-1)=vxaa(nx-1,j-1);
dqx(nx-1,j)=(qx(nx-1,j)-qx(nx-1,j-1))/T;%qx对t求一阶导
vx(nx-1,j)=vx(1,j)+qx(nx-1,j);%v和q之间的关系式
omega32=omega31;
omega31=omega3(j);
snx6=snx5;%s(r,t2)=s(r,t2)
snx5=vz(nx,j);%nx=r s(r,t)

运行结果及报错内容

img

我的解答思路和尝试过的方法

img