Matlab讨论区 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

ahu250 发表于 2009-10-10 15:01

求教高人指点优化程序

最近对一个含有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=;
b1=1:n-1; b2=1:n-2 ; b=;
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=;
b1=1:n-1; b2=1:n-2 ; b=;
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=;
b1=1:n-1; b2=1:n-2 ; b=;
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=;
b1=1:n-1; b2=1:n-2 ; b=;
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,:)=;

%完成对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=;
b1=1:n-1; b2=1:n-2 ; b=;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=;
b1=1:n-1; b2=1:n-2 ; b=;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=;
b1=1:n-1; b2=1:n-2 ; b=;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=;
b1=1:n-1; b2=1:n-2 ; b=;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

ChaChing 发表于 2009-10-10 17:35

楼主认为别人要看完并懂你的程序需要多久时间? 若是LZ会看吗?
所以建议楼主最好简化下问题, 不然要别人参与机会较少些!

ahu250 发表于 2009-10-10 22:44

:@L 主要把偏微分方程离散化成一介齐线性常微分方程组后求解，解得形似一般是y=exp(At),遇到的问题就是离散节点数一旦增加，矩阵A相应也增加的很快，exp（At）计算量就越大，越消耗时间，运行程序就出不了结果，我该怎么办，有办法解决没，尽管已经使用了稀疏矩阵进行简化但好像还是没什么效果

liushuiwuxin 发表于 2010-8-5 20:30

"偏微分方程离散化成一介齐线性常微分方程组"想请问楼主这个是用什么方法可以得到呀？可否介绍一本书看看？谢谢啦！

页: [1]

声振论坛's Archiver

求教高人指点优化程序