微分方程组计算速度的制约因素

lc622503

请问各位大虾一阶微分方程组

制约其计算速度的最大障碍是什么啊

我一个方程组就是改变了反应器的体积计算起来所用时间就大不一样了阿


先谢谢了

lc622503

function asm21
clear all,clc
y0=[23.8,35.7,25.6,0.77,5.17,29.75,107.1,23.8,50,30,1,1,23.8,35.7,25.6,0.77,5.17,29.75,107.1,23.8,50,30,1,1]
[t,y]=ode23(@f,[0:1:20],y0)
%---------------------
function dydt=f(t,y)
% parameters
par(1)=3;par(2)=0.06;par(3)=0.1;par(4)=0.2;par(5)=0.5;
par(6)=0.1;par(7)=6;par(8)=3;par(9)=0.8;par(10)=0.4;
par(11)=0.2;par(12)=4;par(13)=20;par(14)=4;par(15)=0.5;
par(16)=0.05;par(17)=0.01;par(18)=0.1;par(19)=3;par(20)=1.5;
par(21)=1;par(22)=0.2;par(23)=0.2;par(24)=0.2;par(25)=0.2;
par(26)=4;par(26)=0.05;par(28)=0.2;par(29)=0.01;par(30)=0.1;
par(31)=0.01;par(32)=0.34;par(33)=0.02;par(34)=0.01;par(35)=1;
par(36)=0.15;par(37)=0.5;par(38)=1;par(39)=0.5;par(40)=0.01;
%influent concentrations
i1=23.8; i2=35.7; i3=25.6; i4=0.77; i5=5.17
i6=29.75; i7=107.1; i8=23.8;i9=0; i10=0
i11=0; i12=0 ;oxygen1=0;oxygen2=3
%flow rates and volumes
qi=10.2;q1=16.1;v1=0.55;q2=1.6;v2=1.66
qr=5.9;qw=0.223;q0=10
% processes
proc1=par(1)*(oxygen1/(par(4)+oxygen1))*((y(7)/y(8))/(par(6)+y(7)/y(8)))*y(8)
proc2=par(1)*par(2)*(par(4)/(par(4)+oxygen1))*(y(4)/(par(5)+y(4)))*((y(7)/y(8))/(par(6)+y(7)/y(8)))*y(8)
proc3=par(1)*par(3)*(par(4)/(par(4)+oxygen1))*(par(5)/(par(5)+y(4)))*((y(7)/y(8))/(par(6)+y(7)/y(8)))*y(8)
proc4=par(7)*(oxygen1/(par(11)+oxygen1))*(y(2)/(par(12)+y(2)))*(y(2)/(y(1)+y(2)))*(y(3)/(par(16)+y(3)))...
       *(y(5)/(par(17)+y(5)))*y(8)
proc5=par(7)*(oxygen1/(par(11)+oxygen1))*(y(1)/(par(14)+y(1)))*(y(1)/(y(1)+y(2)))*(y(3)/(par(16)+y(3)))...
       *(y(5)/(par(17)+y(5)))*y(8)
proc6=par(7)*par(9)*(par(11)/(par(11)+oxygen1))*(y(2)/(par(12)+y(2)))*(y(2)/(y(1)+y(2)))*(y(3)/(par(16)+y(3)))...
       *(y(4)/(par(15)+y(4)))*(y(5)/(par(17)+y(5)))*y(8)
proc7=par(7)*par(9)*(par(11)/(par(11)+oxygen1))*(y(1)/(par(14)+y(1)))*(y(1)/(y(1)+y(2)))*(y(3)/(par(16)+y(3)))...
       *(y(4)/(par(15)+y(4)))*(y(5)/(par(17)+y(5)))*y(8)
proc8=par(8)*(par(11)/(par(11)+oxygen1))*(par(15)/(par(15)+y(4)))*(y(2)/(par(13)+y(2)))*y(8)
proc9=par(10)*y(8)
proc10=par(19)*(y(1)/(par(26)+y(1)))*((y(11)/y(10))/(par(31)+y(11)/y(10)))*y(10)
proc11=par(20)*(oxygen1/(par(25)+oxygen1))*(y(5)/(par(28)+y(5)))*((y(12)/y(10))/(par(34)+y(12)/y(10)))...
    *((par(32)-(y(11)/y(10)))/(par(33)+par(32)-(y(11)/y(10))))*y(10)
proc12=par(21)*(oxygen1/(par(25)+oxygen1))*(y(5)/(par(29)+y(5)))*((y(12)/y(10))/(par(34)+y(12)/y(10)))...
         *(y(3)/(par(27)+y(3)))*y(10)
proc13=par(22)*y(10)
proc14=par(23)*y(11)  
proc15=par(24)*y(12)
proc16=par(35)*(oxygen1/(par(37)+oxygen1))*(y(3)/(par(38)+y(3)))*(y(5)/(par(40)+y(5)))*y(9)
proc17=par(36)*y(9)
proc18=par(1)*(oxygen2/(par(4)+oxygen2))*((y(19)/y(20))/(par(6)+y(19)/y(20)))*y(20)
proc19=par(1)*par(2)*(par(4)/(par(4)+oxygen2))*(y(16)/(par(5)+y(16)))*((y(19)/y(20))/(par(6)+y(19)/y(20)))*y(20)
proc20=par(1)*par(3)*(par(4)/(par(4)+oxygen2))*(par(5)/(par(5)+y(16)))*((y(19)/y(20))/(par(6)+y(19)/y(20)))*y(20)
proc21=par(7)*(oxygen2/(par(11)+oxygen2))*(y(14)/(par(12)+y(14)))*(y(14)/(y(13)+y(14)))*(y(15)/(par(16)+y(15)))...
       *(y(17)/(par(17)+y(17)))*y(20)
proc22=par(7)*(oxygen2/(par(11)+oxygen2))*(y(13)/(par(14)+y(13)))*(y(13)/(y(13)+y(14)))*(y(15)/(par(16)+y(15)))...
       *(y(17)/(par(17)+y(17)))*y(20)
proc23=par(7)*par(9)*(par(21)/(par(11)+oxygen2))*(y(14)/(par(12)+y(14)))*(y(14)/(y(13)+y(14)))*(y(15)/(par(16)+y(15)))...
       *(y(16)/(par(15)+y(16)))*(y(17)/(par(17)+y(17)))*y(20)
proc24=par(7)*par(9)*(par(21)/(par(11)+oxygen2))*(y(13)/(par(14)+y(13)))*(y(13)/(y(13)+y(14)))*(y(15)/(par(16)+y(15)))...
       *(y(16)/(par(15)+y(16)))*(y(17)/(par(17)+y(17)))*y(20)
proc25=par(8)*(par(11)/(par(11)+oxygen2))*(par(15)/(par(15)+y(16)))*(y(14)/(par(13)+y(14)))*y(20)
proc26=par(10)*y(20)
proc27=par(19)*(y(13)/(par(26)+y(13)))*((y(23)/y(22))/(par(31)+y(23)/y(22)))*y(22)
proc28=par(20)*(oxygen2/(par(25)+oxygen2))*(y(17)/(par(28)+y(17)))*((y(24)/y(22))/(par(34)+y(24)/y(22)))...
    *((par(32)-(y(23)/y(22)))/(par(33)+par(32)-(y(23)/y(22))))*y(22)
proc29=par(21)*(oxygen2/(par(25)+oxygen2))*(y(17)/(par(29)+y(17)))*((y(24)/y(22))/(par(34)+y(24)/y(22)))...
         *(y(15)/(par(27)+y(15)))*y(22)
proc30=par(22)*y(22)
proc31=par(23)*y(23)      
proc32=par(24)*y(24)
proc33=par(35)*(oxygen2/(par(37)+oxygen2))*(y(15)/(par(38)+y(15)))*(y(17)/(par(40)+y(17)))*y(21)
proc34=par(36)*y(21)
%f
f1=(qi/v1)*i1+(qr/v1)*y(13)-(q1/v1)*y(1)-1.59*(proc5+proc7)+proc8-proc10+proc15
f2=(qi/v1)*i2+(qr/v1)*y(14)-(q1/v1)*y(2)+proc1+proc2+proc3-1.59*(proc4+proc6)-proc8
f3=(qi/v1)*i3+(qr/v1)*y(15)-(q1/v1)*y(3)+0.01*(proc1+proc2+proc3)-0.022*(proc4+proc6)-0.07*(proc5+proc7+proc12)...
    +0.03*proc8+0.031*(proc9+proc13+proc17)-4.24*proc16
f4=(qi/v1)*i4+(qr/v1)*y(16)-(q1/v1)*y(4)-0.21*(proc6+proc7)+4.17*proc16
f5=(qi/v1)*i5+(qr/v1)*y(17)-(q1/v1)*y(5)-0.004*(proc4+proc6)-0.02*(proc5+proc7+proc12+proc16)...
    +0.01*(proc8+proc9+proc13+proc17)+0.4*proc10-proc11+proc14
f6=(qi/v1)*i6+(qr/v1)*3*y(18)-(q1/v1)*y(6)+0.1*(proc9+proc13+proc17)
f7=(qi/v1)*i7+(qr/v1)*3*y(19)-(q1/v1)*y(7)-(proc1+proc2+proc3)+0.9*(proc9+proc13+proc17)
f8=(qi/v1)*i8+(qr/v1)*3*y(20)-(q1/v1)*y(8)+proc4+proc5+proc6+proc7-proc9
f9=(qi/v1)*i9+(qr/v1)*3*y(21)-(q1/v1)*y(9)+proc16-proc17
f10=(qi/v1)*i10+(qr/v1)*3*y(22)-(q1/v1)*y(10)+proc12-proc13
f11=(qi/v1)*i11+(qr/v1)*3*y(23)-(q1/v1)*y(11)-0.4*proc10+proc11-proc14
f12=(qi/v1)*i12+(qr/v1)*3*y(24)-(q1/v1)*y(12)+proc10-0.2*proc11-1.6*proc12-proc15
f13=(q1/v2)*y(1)-(q0+qr+qw/v2)*y(13)-1.59*(proc22+proc24)+proc25-proc27+proc32
f14=(q1/v2)*y(2)-(q0+qr+qw/v2)*y(14)+proc18+proc19+proc20-1.59*(proc21+proc23)-proc25
f15=(q1/v2)*y(3)-(q0+qr+qw/v2)*y(15)+0.01*(proc18+proc19+proc20)-0.022*(proc21+proc23)-0.07*(proc22+proc24+proc29)...
    +0.03*proc25+0.031*(proc26+proc30+proc34)-4.24*proc33
f16=(q1/v2)*y(4)-(q0+qr+qw/v2)*y(16)-0.21*(proc23+proc24)+4.17*proc33
f17=(q1/v2)*y(5)-(q0+qr+qw/v2)*y(17)-0.004*(proc21+proc23)-0.02*(proc22+proc24+proc29+proc33)...
    +0.01*(proc25+proc26+proc30+proc34)+0.4*proc27-proc28+proc31
f18=(q1/v2)*y(6)-2.5*((qr+qw)/v2)*y(18)+0.1*(proc26+proc30+proc34)
f19=(q1/v2)*y(7)-2.5*((qr+qw)/v2)*y(19)-(proc18+proc19+proc20)+0.9*(proc26+proc30+proc34)
f20=(q1/v2)*y(8)-2.5*((qr+qw)/v2)*y(20)+proc21+proc22+proc23+proc24-proc26
f21=(q1/v2)*y(9)-2.5*((qr+qw)/v2)*y(21)+proc33-proc34
f22=(q1/v2)*y(10)-2.5*((qr+qw)/v2)*y(22)+proc29-proc30
f23=(q1/v2)*y(11)-2.5*((qr+qw)/v2)*y(23)-0.4*proc27+proc28-proc31
f24=(q1/v2)*y(12)-2.5*((qr+qw)/v2)*y(24)+proc27-0.2*proc28-1.6*proc29-proc32
dydt=[f1;f2;f3;f4;f5;f6;f7;f8;f9;f10;f11;f12;f13;f14;f15;f16;f17;f18;f19;f20;f21;f22;f23;f24]
以上的程序我只是把f18到f24七个方程中七个2.5由原来的3改变而来
其余的地方没动
但是方程就从原来三分钟得到答案变成1个小时还得不到答案阿
还有就是
%flow rates and volumes
qi=10.2;q1=16.1;v1=0.55;q2=1.6;v2=1.66
qr=5.9;qw=0.223;q0=10
改变这里的流量也遇到同样的问题啊
模拟中必须改变这些参数
所以这些问题亟待解决阿
希望各位多多指教阿