新人求助 非线性方程组迭代问题

[shuiyudao]

rt 最近在做一个算例 控制方程组是非线性的 有有限差分后 循环迭代求解 最后求补出来 显示NaN
程序如下：
gt=0.1;gb=0.1;q=sym0.2;ratio=2;
ht=0.15;hb=0.15;eg=20;dx=ratio/100;die=0;diedai=input('迭代步=');
t(102)=0;w(102)=0;p(102)=0;tol=0.0001;
for i=1:51
p(i)=q;
end
for i=1:diedai
die=die+1;
wr(1:102)=w(1:102);
tr(1:102)=t(1:102);
for j=2:101
aa(j)=wr(j+1)-wr(j-1);
aa(j)=aa(j)/2/dx;
bb(j)=tr(j+1)-tr(j-1);
bb(j)=bb(j)/2/dx;
zx(j)=aa(j)/sqrt(aa(j)*aa(j)+1);
yx(j)=1/sqrt(aa(j)*aa(j)+1);
w(j)=0.5*wr(j+1)+0.5*wr(j-1)+dx*dx/2/(tr(j)*yx(j)^3+gt+gb)*(p(j)+zx(j)*bb(j)-wr(j));
t(j)=0.5*tr(j+1)+0.5*tr(j-1);
t(j)=t(j)+dx*dx/2*zx(j)*yx(j)/(tr(j)*yx(j)^3+gt+gb)*(p(j)+zx(j)*bb(j)-wr(j))*bb(j);
t(j)=t(j)-dx*dx/2/yx(j)*(gt/ht/ht+gb/hb/hb)*(sqrt(tr(j)*tr(j)/eg/eg+1-aa(j)*aa(j))-1);
end
w(1)=w(3);
w(102)=w(100);
t(1)=t(3);
t(102)=-t(100);
pa=0;
pb=0;
for j=2:102
if w(j)~=0
ca=abs(w(j)-wr(j))/abs(w(j));
if double(ca)>tol
pa=pa+1;
end
elseif wr(j)==0&w(j)==0
pa=pa;
elseif wr(j)~=0&w(j)==0
pa=pa+1;
end
end
for j=2:102
if t(j)~=0
ca=abs(t(j)-tr(j))/abs(t(j));
if double(ca)>tol
pa=pa+1;
end
elseif tr(j)==0&t(j)==0
pb=pb;
elseif tr(j)~=0&t(j)==0
pb=pb+1;
end
end
if pa==0&pb==0
die
break
end
end