24小时热门版块排行榜

返回列表

nadou

新虫 (初入文坛)

应助: 0 (幼儿园)
帖子: 1
在线: 37分钟
虫号: 1228257
注册: 2011-03-10

[交流] 【求助】解刚性微分方程组已有1人参与

我要解一个44维刚性微分方程组，使用ode23s，但计算结果明显不对，在此求助大家！微分方程组是关于化学反应动力学方面的
代码如下：
A=zeros(5,1);
E=zeros(5,1);
tspan=[0 2820];
y0=zeros(1,44);
y0(1,44)=0.003425;
[t,y]=ode23s('rateequation', tspan, y0,[]);
function ydot=rateequation(t,y)
A1=3.98E+12;
E1=281.4E+3;
A2=1.58E+10;
E2=49.14E+3;
A3=1.00E+14;
E3=130.2E+3;
A4=3.16E+7;
E4=57.54E+3;
A5=1.00E+10;
E5=0;
T=653.15;
R=8.314;                                     %J/mol/K
k1=A1*exp(E1/R/T);                   %求解动力学速率常数
k2=A2*exp(E2/R/T);
k3=A3*exp(E3/R/T);
k4=A4*exp(E4/R/T);
k5=A5*exp(E5/R/T);
ydot=zeros(44,1);
ydot(1)=k3*y(42);                      %所求微分方程组，44维，y（）为某一组分的浓度
ydot(2)=k3*y(43);
ydot(3)=k3*y(43);
ydot(4)=k3*y(42);
ydot(5)=k3*y(41);
ydot(6)=k3*y(40);
ydot(7)=k3*y(39);
ydot(8)=28*k4*y(44)*y(24)+2*k5*(y(15)*y(19)+2*y(16)*y(18)+2*y(17)*y(17));
ydot(9)=14*k4*y(44)*y(26)+2*k5*(y(15)*y(21)+2*y(16)*y(19)+2*y(17)*y(18));
ydot(10)=28*k4*y(44)*y(28)+2*k5*(y(15)*y(23)+2*y(16)*y(21)+2*y(17)*y(19)+2*y(18)*y(18));
ydot(11)=28*k4*y(44)*y(30)+2*k5*(y(15)*y(25)+2*y(16)*y(23)+2*y(17)*y(21)+2*y(18)*y(19));
ydot(12)=28*k4*y(44)*y(32)+2*k5*(y(15)*y(27)+2*y(16)*y(25)+2*y(17)*y(23)+2*y(18)*y(21)+2*y(19)*y(19));
ydot(13)=28*k4*y(44)*y(34)+2*k5*(y(15)*y(29)+2*y(16)*y(27)+2*y(17)*y(25)+2*y(18)*y(23)+2*y(19)*y(21));
ydot(14)=28*k4*y(44)*y(36)+2*k5*(y(15)*y(31)+2*y(16)*y(29)+2*y(17)*y(27)+2*y(18)*y(25)+2*y(19)*y(23)+2*y(21)*y(21));
ydot(15)=k1*y(44)-2*k5*y(15)*(y(19)+y(21)+y(23)+y(25)+y(27)+y(29)+y(31))+k3*y(39);
ydot(16)=k1*y(44)-4*k5*y(16)*(y(18)+y(19)+y(21)+y(23)+y(25)+y(27)+y(29))+k3*y(37);
ydot(17)=k1*y(44)-4*k5*y(17)*(y(17)+y(18)+y(19)+y(21)+y(23)+y(25)+y(27))+k3*y(41);
ydot(18)=k1*y(44)-4*k5*y(18)*(y(16)+y(17)+y(18)+y(19)+y(21)+y(23)+y(25))+k3*y(42);
ydot(19)=k1*y(44)-2*k5*y(19)*(y(15)+2*y(16)+2*y(17)+2*y(18)+2*y(19)+2*y(21)+2*y(23))+k3*y(43)-4*k2*y(19);
ydot(20)=4*k2*y(19)-28*k4*y(20)*y(44);
ydot(21)=k1*y(44)-2*k5*y(21)*(y(15)+2*y(16)+2*y(17)+2*y(18)+2*y(19)+2*y(21))+k3*y(43)-4*k2*y(21);
ydot(22)=4*k2*y(21)-28*k4*y(22)*y(44);
ydot(23)=k1*y(44)-2*k5*y(23)*(y(15)+2*y(16)+2*y(17)+2*y(18)+2*y(19))+k3*y(42)-4*k2*y(23);
ydot(24)=4*k2*y(23)-28*k4*y(24)*y(44);
ydot(25)=k1*y(44)-2*k5*y(25)*(y(15)+2*y(16)+2*y(17)+2*y(18))+k3*y(41)-4*k2*y(25);
ydot(26)=4*k2*y(25)-28*k4*y(26)*y(44);
ydot(27)=k1*y(44)-2*k5*y(27)*(y(15)+2*y(16)+2*y(17))+k3*y(40)-4*k2*y(27);
ydot(28)=4*k2*y(27)-28*k4*y(28)*y(44);
ydot(29)=k1*y(44)-2*k5*y(29)*(y(15)+2*y(16))+k3*y(39)-4*k2*y(29);
ydot(30)=4*k2*y(29)-28*k4*y(30)*y(44);
ydot(31)=k1*y(44)-2*k5*y(15)*y(31)+k3*y(38)-4*k2*y(31);
ydot(32)=4*k2*y(31)-28*k4*y(32)*y(44);
ydot(33)=k1*y(44)+k3*y(37)-4*k2*y(33);
ydot(34)=4*k2*y(33)-28*k4*y(34)*y(44);
ydot(35)=k1*y(44)-4*k2*y(35);
ydot(36)=4*k2*y(35)-28*k4*y(36)*y(44);
ydot(37)=2*y(44)*k4*(2*y(24)+y(26)+2*y(28)+2*y(30)+2*y(32)+2*y(34)+2*y(36))-k3*y(37);
ydot(38)=2*y(44)*k4*(2*y(24)+y(26)+2*y(28)+2*y(30)+2*y(32)+2*y(34)+2*y(36))-k3*y(38);
ydot(39)=2*y(44)*k4*(2*y(24)+y(26)+2*y(28)+2*y(30)+2*y(32)+2*y(34)+2*y(36))-k3*y(39);
ydot(40)=2*y(44)*k4*(2*y(24)+y(26)+2*y(28)+2*y(30)+2*y(32)+2*y(34)+2*y(36))-k3*y(40);
ydot(41)=2*y(44)*k4*(2*y(24)+y(26)+2*y(28)+2*y(30)+2*y(32)+2*y(34)+2*y(36))-k3*y(41);
ydot(42)=2*y(44)*k4*(2*y(24)+y(26)+2*y(28)+2*y(30)+2*y(32)+2*y(34)+2*y(36))-k3*y(42);
ydot(43)=2*y(44)*k4*(2*y(24)+y(26)+2*y(28)+2*y(30)+2*y(32)+2*y(34)+2*y(36))-k3*y(43);
ydot(44)=-7*k1*y(44)-14*y(44)*k4*(2*y(24)+y(26)+2*y(28)+2*y(30)+2*y(32)+2*y(34)+2*y(36));

回复此楼