6个微分方程系统

我想用matlab解6个微分方程组。我创建了一组6次微分方程如命名为Untitled.m6个微分方程系统

function ydot=Untitled(t,y) 
ydot = zeros(6,1); 
%y(1)=A 
%y(2)=B 
%y(3)=C 
%y(4)=D 
%y(5)=P 
%y(6)=T; 

A=0.50265; 
k11=(333/106.7)*1.15*1000*exp(-59660/(8.314*960)); 
k31=(333/40)*73.6*exp(-47820/(8.314*960)); 
k32=(333/14.4)*1.79*exp(-30950/(8.314*960)); 
k21=(106.7/40)*426*exp(-68830/(8.314*960)); 
k22=(106.7/14.4)*0.000599*exp(-57740/(8.314*960)); 
Pcat=1450; 
g=9.81; 
%phi=exp(-(59100*exp(-67210/(8.314*y(6))))*((0.50265*1450*33)/143.64)); 
H11=393000; 
H31=795000; 
H32=1200000; 
H21=1150000; 
H22=151000; 
E=1-((285.765*17.56)/((6.1*1450)+(17.56*285.765))); 
Fcat=143.64; 
Cpcat=1087; 
%Cp=1000*(y(1)*3.3+y(2)*3.3+y(3)*3.3+y(4)*1.087); 
F=19.95; 




ydot(1)= -(A*(1-E)*Pcat*(k11+k31+k32)*(exp(-(59100*exp(-67210/(8.314*y(6))))*((0.50265*1450*33)/143.64)))*y(1)*y(1))/F; 
ydot(2)= (A*(1-E)*Pcat*(k11*y(1)*y(1)-(k21+k22)*y(2))*(exp(-(59100*exp(-67210/(8.314*y(6))))*((0.50265*1450*33)/143.64))))/F; 
ydot(3)= (A*(1-E)*Pcat*(k31*y(1)*y(1)+ k21*y(2))*(exp(-(59100*exp(-67210/(8.314*y(6))))*((0.50265*1450*33)/143.64))))/F; 
ydot(4)= (A*(1-E)*Pcat*(k32*y(1)*y(1)+k22*y(2))*(exp(-(59100*exp(-67210/(8.314*y(6))))*((0.50265*1450*33)/143.64))))/F; 
ydot(5)= -(Pcat*g*(1-E)); 
%dydz(6)= ((1-E)*Pcat*A*(((k11+k31+k32)*phi*y(1)*y(1)*-H1)+ ((k11*y(1)*y(1)-(k21+k22)*y(2))*phi*-H2)+((k31*y(1)*y(1)+ k21*y(2))*phi*H3)+((k32*y(1)*y(1)+k22*y(2))*phi*H4)))/(F*Cp+Fcat*Cpcat) 
ydot(6) = (((exp(-(59100*exp(-67210/(8.314*y(6))))*((0.50265*1450*33)/143.64)))*(1-E)*Pcat*A)*(y(1)*y(1)*((k11*H11)+(k31*H31)+(k32*H32))+y(2)*((k21*H21)+(k22*H22)))/((F*(1000*(y(1)*3.3+y(2)*3.3+y(3)*3.3+y(4)*1.087)))+(Fcat*Cpcat)));

然后我已创建的另一个文件使用ODE45求解器求解方程的函数米文件如下

function ode() 

options = odeset('RelTol',1e-1,'AbsTol',[1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 ]); 
Y0=[1.0;0.0;0.0;0.0;180000.0;959.81]; 
zspan=0:0.1:33; 

[z,y]= ode45(@Untitled,zspan,Y0,options); 

figure 
hold on 
plot(z,y(:,5));

的代码正在编译，但不能求解微分方程组，只能在变量的初始值处给出直线图。例如，y（1）的初始值是1，所以我在y = 1时得到一个图。

谁能告诉我什么错误

来源

2014-02-12 user3277001

Untitled.m给出了一个警告：“输入参数z可能未被使用，但后面的输入参数被使用”;所以我猜这是不正确得到差异方程式...但我不知道该怎么办...请帮助 – user3277001

这只是意味着't'没有被用于计算衍生物。我测试了'Untitled（0，Y0）'，并得到了'ydot = [0; 0; 0; 0; -5148; 0]'，因此它看起来应该是除了'y（5）'之外的直线。检查你的衍生物！ – David

看来你衍生物严重缩放，或只是简单的不正确。

的ydot的第五元素是

-(Pcat*g*(1-E));

仅包含常数，因此是恒定。因此，这将导致一条负斜率的直线。

的ydot另外5个元素包含术语如

exp(-59100*exp(-67210/(...))*(...))

即几乎肯定会下溢，这意味着，该值过小在双精度数中表达，并因此将等于0。本也将产生直线，这次是水平的并且等于它们的初始值。

价值exp(-745.1332...)（等于4.940656458412465e-324）是您能够去的最低价。而且这个值已经非常棘手，因为它已经变质了（exp(-724) < realmin）。

所以要么你搞砸你的方程，要么你需要重新规模的方程，以便在解空间内的上/下溢概率减少到零。

来源

2014-02-12 09:15:15

6个微分方程系统

回答

相关问题