|
楼主 |
发表于 2007-7-21 19:29
|
显示全部楼层
本帖最后由 牛小贱 于 2014-6-29 13:14 编辑
AX=B
A是系数矩阵,B是列矩阵
附程序如下,望各位高手不吝赐教,小女子暂且先谢过:
- clear
- yy=1e-6;zz=1e-6;tt=1e-3;yita=4.6857e-7; R=3.3333e-5; seta=4.027e+6; B=0.5;rou=2450;
- N=1;
- A=5e-4;
- m=200;n=500;
- 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);
- D=939140;G=-935140;
- X=zeros(2*m+1,3*n)+eps;
- QX=zeros(2*m+1,3*n)+eps;
- U=ones(2*m+1,3*n)*eps;
- QU=ones(2*m+1,3*n)*eps;
- vy=zeros(2*m+1,3*n);
- vz=zeros(2*m+1,3*n);
- C=-[(N*vy)/(2*yy)]-yita/yy^2; E=[(N*vy)/(yy^2)]-yita/yy^2;
- F=[(N*(vz+R))/(2*zz)]+yita/zz^2;
- H=-[(N*(vz+R))/(2*zz)]+yita/zz^2;
- A1=zeros(2*m+1);
- for j=2:n-1
- i=m;
- while i*yy>800*(j*zz)^2
- i=i-1;
- end
-
- A1(1+m-i,1+m-i)=1;
- A1(1+m+i,1+m+i)=1;
- for k=2+m-i:m+i
- A1(k,k-1)=C(k,k-1);A1(k,k)=D;A1(k,k+1)=E(k,k+1);
- end
- AA=zeros(2*i+1);
- u=1;v=1;
- for h=1+m-i:1+m+i
- for g=1+m-i:1+m+i
- AA(u,v)=A1(h,g);
- v=v+1;
- end
- u=u+1;
- v=1;
- end
- Y=zeros(2*m+1,1);
- for f=2+m-i:m+i
- Y(f,1)=F(f,j-1)*X(f,j)+G*X(f,j)+H(f,j+1)*X(f,j+1)-seta*(B^2)*X(f,j)/rou;
- end
- Y(1+m-i,1)=2*eps/yy^2;
- Y(1+m+i,1)=2*eps/yy^2;
- Y1=zeros(2*i+1,1);
- u=1;
- for h=1+m-i:1+m+i
- Y1(u,1)=Y(h,1);
- u=u+1;
- end
- X1=zeros(2*i+1,1);
- X1=(AA+eps)\Y1;
- Xj=zeros(2*m+1,1);
- u=1;
- for h=1+m-i:1+m+i
- Xj(h,1)=X1(u,1);
- u=u+1;
- end
- X(:,j)=Xj;
- end
- X
复制代码
|
|