|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
我调试的时候发现结果是:矩阵索引超出大小限制.检查好多遍也没有发现问题出在哪里.拜托各位高手指教啊!具体程序如下:
yy=0.001;zz=0.001;tt=0.01;yita=1; R=0.01;
N=1.20;A=0.05;
m=236;n=234;
X=zeros(2*m+1,3*n);
vy=zeros(2*m+1,3*n);
vz=zeros(2*m+1,3*n);
C=zeros(2*m+1,3*n); E=zeros(2*m+1,3*n);
F=zeros(2*m+1,3*n); H=zeros(2*m+1,3*n);
for t=1:10000
C=-[(N*vy)/(2*yy)]-yita/yy^2;D=1/tt+2*yita/yy^2;E=[(N*vy)/(yy^2)]-yita/yy^2;
F=[(N*(vz+R))/(2*zz)]+yita/zz^2;
G=1/tt-(2*yita)/zz^2;
H=-[(N*(vz+R))/(2*zz)]+yita/zz^2;
A1=zeros(2*m+1);
for j=1:n
i=m;
while i*yy>sqrt(A)-sqrt(A-j*zz)
i=i-1;
end
A1(1+m-i,j)=1;
A1(1+m+i,j)=1;
for k=2+m-i:m+i
A1(k,k-1)=C(k,k-1);A1(k,k)=D(k,k);A1(k,k+1)=E(k,k+1);
end
Xj=X(:,j);
Y=zeros(2*m+1,1);
for f=2+m-i:m+i
Y(f,1)=F(i,j-1)*X(i,j-1)+G(i,j)*X(i,j)+H(i,j+1)*X(i,j+1)-seta*(B^2)*X(i,j);
end
Xj=Y/A1;
end
Xj= (这边我想输出Xj的结果,怎么弄啊?)
end
X= (这边我想输出X的结果) |
|