马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
最近对一个含有15个偏微分的方程组进行数值离散化后,编程画图,n是离散的节点数目,n取较小值比如20,30的时候,用不了一分钟就能画出结果图,但是n作为节点当n取100或者更多时运行一直不出结果,最多运行一整天也没结果,:'( 也没提示出错,希望哪位高手能指点一下,有无改进的地方还是哪里出问题,我电脑配置2G内存,主频2.0的,下面是主要程序,希望哪位高人能指点我一下,感激不尽:handshake
比如运行huahuatutu(100,5,50,0.1,0.01,1000,200,5)
主要程序如下:
function huahuatutu(n,delta,N,lambda0,min,max,tt,cc)
tic
lambda1=lambda0;lambda2=lambda0;lambda3=lambda0;lambdap=0.05;mur=30;mu2=90;mu4=1;mum=12;muM=4;
;
%区间间隔大小:1/delta;单位(小时)
%对A11进行赋值,为了节省储存空间,把A看作稀疏矩阵进行赋值,命令sparse(a,b,c)中,a,b,c分别表示原矩阵中非零元素所对应的行标,列标
%,以及非零元素自身
a1=1:n-1 ; a2=2:n-1 ; a=[a1 a2];
b1=1:n-1 ; b2=1:n-2 ; b=[b1 b2];
lambda(1:N)=min ;lambda(N+1:n-2)=max;a11=[-(lambda0+lambdap+delta+lambda),-(lambda0+lambdap+max),delta*ones(1,n-2)];
A11=sparse(a,b,a11);
%完成了对A11矩阵的赋值
%定义一个首行为1,其余元素为0的矩阵
BB=spalloc(n-1,n-1,n-1);
BB(1,:)=1;
% spalloc(m,n,s)表示创建一个m行n列的稀疏矩阵,其中预留s个非零元素的储存空间为赋值准备
A12=spalloc(n-1,n-1,1);
A13=spalloc(n-1,n-1,1);
A14=spalloc(n-1,n-1,1);
a15=ones(1,n-1)*mur; A15=sparse(a1,b1,a15);
A16=spalloc(n-1,n-1,1);
A17=spalloc(n-1,n-1,1);
A18=spalloc(n-1,n-1,1);
A19=BB*mu2;
A1A=spalloc(n-1,n-1,1);
A1B=spalloc(n-1,n-1,1);
A1C=spalloc(n-1,n-1,1);
A1D=spalloc(n-1,n-1,1);
A1E=BB*muM;
A1F=BB*mu4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A21=eye(n-1,n-1)*lambda0;
%开始对A22赋值
a1=1:n-1 ; a2=2:n-1 ; a=[a1 a2];
b1=1:n-1 ; b2=1:n-2 ; b=[b1 b2];
a22=[-(lambda1+lambdap+delta+lambda),-(lambda1+lambdap+max),delta*ones(1,n-2)];
A22=sparse(a,b,a22);
%完成对A22的赋值
A23=spalloc(n-1,n-1,1);
A24=spalloc(n-1,n-1,1);
A25=spalloc(n-1,n-1,1);
a26=ones(1,n-1)*mur; A26=sparse(a1,b1,a26);
A27=spalloc(n-1,n-1,1);
A28=spalloc(n-1,n-1,1);
A29=spalloc(n-1,n-1,1);
A2A=BB*mu2;
A2B=spalloc(n-1,n-1,1);
A2C=spalloc(n-1,n-1,1);
A2D=BB*mum;
A2E=spalloc(n-1,n-1,1);
A2F=spalloc(n-1,n-1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A31=spalloc(n-1,n-1,1);
a32=ones(1,n-1)*lambda1; A32=sparse(a1,b1,a32);
%开始对A33赋值
a1=1:n-1 ; a2=2:n-1 ; a=[a1 a2];
b1=1:n-1 ; b2=1:n-2 ; b=[b1 b2];
a33=[-(lambda2+lambdap+delta+lambda),-(lambda2+lambdap+max),delta*ones(1,n-2)];
A33=sparse(a,b,a33);
%完成对A33的赋值
A34=spalloc(n-1,n-1,1);
A35=spalloc(n-1,n-1,1);
A36=spalloc(n-1,n-1,1);
a37=ones(1,n-1)*mur; A37=sparse(a1,b1,a37);
A38=spalloc(n-1,n-1,1);
A39=spalloc(n-1,n-1,1);
A3A=spalloc(n-1,n-1,1);
A3B=spalloc(n-1,n-1,1);
A3C=spalloc(n-1,n-1,1);
A3D=spalloc(n-1,n-1,1);
A3E=spalloc(n-1,n-1,1);
A3F=spalloc(n-1,n-1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A41=spalloc(n-1,n-1,1);
A42=spalloc(n-1,n-1,1);
a43=ones(1,n-1)*lambda2; A43=sparse(a1,b1,a43);
%开始对A44赋值
a1=1:n-1 ; a2=2:n-1 ; a=[a1 a2];
b1=1:n-1 ; b2=1:n-2 ; b=[b1 b2];
a44=[-(lambda3+lambdap+delta+lambda),-(lambda3+lambdap+max),delta*ones(1,n-2)];
A44=sparse(a,b,a44);
%完成对A44的赋值
A45=spalloc(n-1,n-1,1);
A46=spalloc(n-1,n-1,1);
A47=spalloc(n-1,n-1,1);
a48=ones(1,n-1)*mur; A48=sparse(a1,b1,a48);
A49=spalloc(n-1,n-1,1);
A39=spalloc(n-1,n-1,1);
A4A=spalloc(n-1,n-1,1);
A4B=spalloc(n-1,n-1,1);
A4C=spalloc(n-1,n-1,1);
A4D=spalloc(n-1,n-1,1);
A4E=spalloc(n-1,n-1,1);
A4F=spalloc(n-1,n-1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a51=ones(1,n-1)*lambdap; A51=sparse(a1,b1,a51);
A52=spalloc(n-1,n-1,1);
A53=spalloc(n-1,n-1,1);
A54=spalloc(n-1,n-1,1);
a55=ones(1,n-1)*(-mur); A55=sparse(a1,b1,a55);
A56=spalloc(n-1,n-1,1);
A57=spalloc(n-1,n-1,1);
A58=spalloc(n-1,n-1,1);
A59=spalloc(n-1,n-1,1);
A5A=spalloc(n-1,n-1,1);
A5B=spalloc(n-1,n-1,1);
A5C=spalloc(n-1,n-1,1);
A5D=spalloc(n-1,n-1,1);
A5E=spalloc(n-1,n-1,1);
A5F=spalloc(n-1,n-1,1);
A61=spalloc(n-1,n-1,1);
a62=ones(1,n-1)*lambdap; A62=sparse(a1,b1,a62);
A63=spalloc(n-1,n-1,1);
A64=spalloc(n-1,n-1,1);
A65=spalloc(n-1,n-1,1);
A66=A55;
A67=spalloc(n-1,n-1,1);
A68=spalloc(n-1,n-1,1);
A69=spalloc(n-1,n-1,1);
A6A=spalloc(n-1,n-1,1);
A6B=spalloc(n-1,n-1,1);
A6C=spalloc(n-1,n-1,1);
A6D=spalloc(n-1,n-1,1);
A6E=spalloc(n-1,n-1,1);
A6F=spalloc(n-1,n-1,1);
A71=spalloc(n-1,n-1,1);
A72=spalloc(n-1,n-1,1);
A73=A62;
A74=spalloc(n-1,n-1,1);
A75=spalloc(n-1,n-1,1);
A76=spalloc(n-1,n-1,1);
A77=A55;
A78=spalloc(n-1,n-1,1);
A79=spalloc(n-1,n-1,1);
A7A=spalloc(n-1,n-1,1);
A7B=spalloc(n-1,n-1,1);
A7C=spalloc(n-1,n-1,1);
A7D=spalloc(n-1,n-1,1);
A7E=spalloc(n-1,n-1,1);
A7F=spalloc(n-1,n-1,1);
A81=spalloc(n-1,n-1,1);
A82=spalloc(n-1,n-1,1);
A83=spalloc(n-1,n-1,1);
A84=A62;
A85=spalloc(n-1,n-1,1);
A86=spalloc(n-1,n-1,1);
A87=spalloc(n-1,n-1,1);
A88=A55;
A89=spalloc(n-1,n-1,1);
A8A=spalloc(n-1,n-1,1);
A8B=spalloc(n-1,n-1,1);
A8C=spalloc(n-1,n-1,1);
A8D=spalloc(n-1,n-1,1);
A8E=spalloc(n-1,n-1,1);
A8F=spalloc(n-1,n-1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%开始对A91赋值
A91=spalloc(n-1,n-1,n-1);A91(1,:)=[lambda,max];
%完成对A91的赋值
A92=spalloc(n-1,n-1,1);
A93=spalloc(n-1,n-1,1);
A94=spalloc(n-1,n-1,1);
A95=spalloc(n-1,n-1,1);
A96=spalloc(n-1,n-1,1);
A97=spalloc(n-1,n-1,1);
A98=spalloc(n-1,n-1,1);
%开始对A99赋值
a1=1:n-1 ; a2=2:n-1 ; a=[a1 a2];
b1=1:n-1 ; b2=1:n-2 ; b=[b1 b2]; a99=[-(mu2+delta)*ones(1,n-2),-mu2,delta*ones(1,n-2)];
A99=sparse(a,b,a99);
%完成对A99的赋值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%A10,10; A11.11 A12.12 相同
A10A=A99;
A11B=A99;
A12C=A99;
%开始对A1313赋值
a1=1:n-1 ; a2=2:n-1 ; a=[a1 a2];
b1=1:n-1 ; b2=1:n-2 ; b=[b1 b2]; a13d=[-(mum+delta)*ones(1,n-2),-mum,delta*ones(1,n-2)];
A13D=sparse(a,b,a13d);
%完成对A1313的赋值
%开始对A1414赋值
a1=1:n-1 ; a2=2:n-1 ; a=[a1 a2];
b1=1:n-1 ; b2=1:n-2 ; b=[b1 b2]; a14e=[-(muM+delta)*ones(1,n-2),-muM,delta*ones(1,n-2)];
A14E=sparse(a,b,a14e);
%完成对A1414的赋值
%开始对A1515赋值
a1=1:n-1 ; a2=2:n-1 ; a=[a1 a2];
b1=1:n-1 ; b2=1:n-2 ; b=[b1 b2]; a15f=[-(mu4+delta)*ones(1,n-2),-mu4,delta*ones(1,n-2)];
A15F=sparse(a,b,a15f);
%完成对A1515的赋值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A9A=spalloc(n-1,n-1,1);
A9B=spalloc(n-1,n-1,1);
A9C=spalloc(n-1,n-1,1);
A9D=spalloc(n-1,n-1,1);
A9E=spalloc(n-1,n-1,1);
A9F=spalloc(n-1,n-1,1);
A101=spalloc(n-1,n-1,1);
%开始对A102赋值
A102=A91;
%完成对A102的赋值
A103=spalloc(n-1,n-1,1);
A104=spalloc(n-1,n-1,1);
A105=spalloc(n-1,n-1,1);
A106=spalloc(n-1,n-1,1);
A107=spalloc(n-1,n-1,1);
A108=spalloc(n-1,n-1,1);
A109=spalloc(n-1,n-1,1);
A10B=spalloc(n-1,n-1,1);
A10C=spalloc(n-1,n-1,1);
A10D=spalloc(n-1,n-1,1);
A10E=spalloc(n-1,n-1,1);
A10F=spalloc(n-1,n-1,1);
A111=spalloc(n-1,n-1,1);
A112=spalloc(n-1,n-1,1);
%开始对A113赋值
A113=A91;
%完成对A113的赋值
A114=spalloc(n-1,n-1,1);
A115=spalloc(n-1,n-1,1);
A116=spalloc(n-1,n-1,1);
A117=spalloc(n-1,n-1,1);
A118=spalloc(n-1,n-1,1);
A119=spalloc(n-1,n-1,1);
A11A=spalloc(n-1,n-1,1);
A11C=spalloc(n-1,n-1,1);
A11D=spalloc(n-1,n-1,1);
A11E=spalloc(n-1,n-1,1);
A11F=spalloc(n-1,n-1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%
A121=spalloc(n-1,n-1,1);
A122=spalloc(n-1,n-1,1);
A123=spalloc(n-1,n-1,1);
%开始对A124赋值
A124=A91;
%完成对A124的赋值
A125=spalloc(n-1,n-1,1);
A126=spalloc(n-1,n-1,1);
A127=spalloc(n-1,n-1,1);
A128=spalloc(n-1,n-1,1);
A129=spalloc(n-1,n-1,1);
A12A=spalloc(n-1,n-1,1);
A12B=spalloc(n-1,n-1,1);
A12D=spalloc(n-1,n-1,1);
A12E=spalloc(n-1,n-1,1);
A12F=spalloc(n-1,n-1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%
A131=spalloc(n-1,n-1,1);
A132=spalloc(n-1,n-1,1);
A133=spalloc(n-1,n-1,1);
A134=spalloc(n-1,n-1,1);
A135=spalloc(n-1,n-1,1);
A136=spalloc(n-1,n-1,1);
A137=spalloc(n-1,n-1,1);
A138=spalloc(n-1,n-1,1);
A139=spalloc(n-1,n-1,1);
A13A=spalloc(n-1,n-1,1);
A13B=BB*mu2;
A13C=spalloc(n-1,n-1,1);
A13E=spalloc(n-1,n-1,1);
A13F=spalloc(n-1,n-1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%
A141=spalloc(n-1,n-1,1);
A142=spalloc(n-1,n-1,1);
A143=spalloc(n-1,n-1,1);
A144=spalloc(n-1,n-1,1);
A145=spalloc(n-1,n-1,1);
A146=spalloc(n-1,n-1,1);
A147=spalloc(n-1,n-1,1);
A148=spalloc(n-1,n-1,1);
A149=spalloc(n-1,n-1,1);
A14A=spalloc(n-1,n-1,1);
A14B=spalloc(n-1,n-1,1);
A14C=BB*mu2;
A14D=spalloc(n-1,n-1,1);
A14F=spalloc(n-1,n-1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%
A151=spalloc(n-1,n-1,1);
A152=spalloc(n-1,n-1,1);
A153=spalloc(n-1,n-1,1);
A154=BB*lambda3;
A155=spalloc(n-1,n-1,1);
A156=spalloc(n-1,n-1,1);
A157=spalloc(n-1,n-1,1);
A158=spalloc(n-1,n-1,1);
A159=spalloc(n-1,n-1,1);
A15A=spalloc(n-1,n-1,1);
A15B=spalloc(n-1,n-1,1);
A15C=spalloc(n-1,n-1,1);
A15D=spalloc(n-1,n-1,1);
A15E=spalloc(n-1,n-1,1);
A=[A11,A12,A13,A14,A15,A16,A17,A18,A19,A1A,A1B,A1C,A1D,A1E,A1F;
A21,A22,A23,A24,A25,A26,A27,A28,A29,A2A,A2B,A2C,A2D,A2E,A2F;
A31,A32,A33,A34,A35,A36,A37,A38,A39,A3A,A3B,A3C,A3D,A3E,A3F;
A41,A42,A43,A44,A45,A46,A47,A48,A49,A4A,A4B,A4C,A4D,A4E,A4F;
A51,A52,A53,A54,A55,A56,A57,A58,A59,A5A,A5B,A5C,A5D,A5E,A5F;
A61,A62,A63,A64,A65,A66,A67,A68,A69,A6A,A6B,A6C,A6D,A6E,A6F;
A71,A72,A73,A74,A75,A76,A77,A78,A79,A7A,A7B,A7C,A7D,A7E,A7F;
A81,A82,A83,A84,A85,A86,A87,A88,A89,A8A,A8B,A8C,A8D,A8E,A8F;
A91,A92,A93,A94,A95,A96,A97,A98,A99,A9A,A9B,A9C,A9D,A9E,A9F;
A101,A102,A103,A104,A105,A106,A107,A108,A109,A10A,A10B,A10C,A10D,A10E,A10F;
A111,A112,A113,A114,A115,A116,A117,A118,A119,A11A,A11B,A11C,A11D,A11E,A11F;
A121,A122,A123,A124,A125,A126,A127,A128,A129,A12A,A12B,A12C,A12D,A12E,A12F;
A131,A132,A133,A134,A135,A136,A137,A138,A139,A13A,A13B,A13C,A13D,A13E,A13F;
A141,A142,A143,A144,A145,A146,A147,A148,A149,A14A,A14B,A14C,A14D,A14E,A14F;
A151,A152,A153,A154,A155,A156,A157,A158,A159,A15A,A15B,A15C,A15D,A15E,A15F];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p0=spalloc((n-1)*15,1,1);
p0(1,1)=1;
for t=1:tt
solution=expm(A*t/cc)*p0;
avail(t)=sum(solution(1:4*(n-1)))/sum(solution);
end
plot(avail)
toc |