Subscript indices must either be real positive integers or logicals.m=20;n=20;t=0.01;l=0.1;h=0.1;r=1000;w=0.5;s=0.001;dumax=0;dvmax=0;dpmax=0;dx=l/(m+1);dy=l/(n+1);ou=zeros(20);ou(1,:)=1;ov=zeros(20);op=zeros(20);op(:,:)=10315;for i=2:m-1for j=2:n-1o

Subscript indices must either be real positive integers or logicals.m=20;n=20;t=0.01;l=0.1;h=0.1;r=1000;w=0.5;s=0.001;dumax=0;dvmax=0;dpmax=0;dx=l/(m+1);dy=l/(n+1);ou=zeros(20);ou(1,:)=1;ov=zeros(20);op=zeros(20);op(:,:)=10315;for i=2:m-1for j=2:n-1o
Subscript indices must either be real positive integers or logicals.
m=20;
n=20;
t=0.01;
l=0.1;
h=0.1;
r=1000;
w=0.5;
s=0.001;
dumax=0;
dvmax=0;
dpmax=0;
dx=l/(m+1);
dy=l/(n+1);
ou=zeros(20);
ou(1,:)=1;
ov=zeros(20);
op=zeros(20);
op(:,:)=10315;
for i=2:m-1
for j=2:n-1
oA1(i,j)=(ou(i+1,j)-2*ou(i,j)+ou(i-1,j))/dx^2;
oA2(i,j)=(ou(i,j+1)-2*ou(i,j)+ou(i,j-1))*dx^2;
oA3(i,j)=(ou(i+1,j)^2-ou(i,j)^2)/dx;
oA4(i,j)=(ou(i,j+1)*ov(i,j+1)-ou(i,j)*ov(i,j))/dy;
oA(i,j)=1/r*(oA1(i,j)+oA2(i,j))-oA3(i,j)-oA4(i,j);
su(i,j)=t*(oA(i,j)-(t/dx)*(op(i+1,j)-op(i,j))+ou(i,j));
end
end
for i=2:m-1
for j=2:n-1
oB1(i,j)=(ov(i+1,j)-2*ov(i,j)+ou(i-1,j))/dx^2;
oB2(i,j)=(ov(i,j+1)-2*ov(i,j)+ou(i,j-1))*dx^2;
oB3(i,j)=(ov(i+1,j)^2-ov(i,j)^2)/dx;
oB4(i,j)=(ov(i,j+1)*ou(i,j+1)-ov(i,j)*ou(i,j))/dy;
oB(i,j)=1/r*(oB1(i,j)+oB2(i,j))-oB3(i,j)-oB4(i,j);
sv(i,j)=t*(oB(i,j)-(t/dx)*(op(i,j+1)-op(i,j))+ov(i,j));
end
end
for i=2:m-1
for j=2:n-1
a=2*t*(1/dx^2+1/dy^2);
b=t/dx^2;
c=t/dy^2;
d(i,j)=1/dx*(ou(i,j)-ou(i-1,j))+1/dy*(ov(i,j)-ov(i,j-1));
p1(i,j)=(b*op(i+1,j)+op(i-1,j)+c*op(i,j+1)+c*op(i,j-1)-d(i,j))/a;
p1(i,j)=w*p1(i,j)+(1-w)*p1(i,j);
end
end
for i=2:m-1
for j=2:n-1
nu(i,j)=ou(i,j)-t*(p1(i+1,j)-p1(i,j))/dx;
nv(i,j)=ov(i,j)-t*(p1(i,j+1)-p1(i,j))/dy;
np(i,j)=op(i,j)+p1(i,j);
end
end
for i=2:m
for j=2:n
du=abs(nu(i,j)-ou(i,j));
dv=abs(nv(i,j)-ov(i,j));
dp=abs(np(i,j)-op(i,j));
if du>dumax
dumax=du;
end
if dv>dvmax
dvmax=dv;
end
if dp>dpmax
dpmax=dp;
end
end
end
su=nu;
sv=nv;
sp=np;
if dumax

Subscript indices must either be real positive integers or logicals.m=20;n=20;t=0.01;l=0.1;h=0.1;r=1000;w=0.5;s=0.001;dumax=0;dvmax=0;dpmax=0;dx=l/(m+1);dy=l/(n+1);ou=zeros(20);ou(1,:)=1;ov=zeros(20);op=zeros(20);op(:,:)=10315;for i=2:m-1for j=2:n-1o
你第一个程序中第5行“a=d1/(4*na);b=d1/(4*nb);c=d1/(4*nc);d=d1/(4*nd); ”这里你把d的值改了变成了74.93（而你for循环的递增