| 查看: 751 | 回复: 1 | ||
[交流]
【求助】解刚性微分方程组 已有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)); |
回复此楼
信彼南山
木虫 (著名写手)
- 应助: 33 (小学生)
- 金币: 4142.9
- 散金: 1221
- 红花: 16
- 帖子: 1178
- 在线: 233.5小时
- 虫号: 1133529
- 注册: 2010-10-27
- 专业: 导航、制导与传感技术
2楼2011-04-07 00:32:21
