2014-02-05 232 views
1

我想用matlab求解6个微分方程组。我创建了一组6次微分方程如命名为Untitled.m在matlab中用ODE45求解6个微分方程组

function ydot=Untitled(z,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(-al*tc) 
al=59100*exp(-67210/(8.314*T)) 
tc=(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)*phi*y(1)*y(1))/F 
ydot(2)= (A*(1-E)*Pcat*(k11*y(1)*y(1)-(k21+k22)*y(2))*phi)/F 
ydot(3)= (A*(1-E)*Pcat*(k31*y(1)*y(1)+ k21*y(2))*phi)/F 
ydot(4)= (A*(1-E)*Pcat*(k32*y(1)*y(1)+k22*y(2))*phi)/F 
ydot(5)= -(Pcat*g*(1-E)) 
ydot(6) = ((phi*(1-E)*Pcat*A)*(y(1)*y(1)((k11*H11)+(k31*H31)+(k32*H32))+y(2)((k21*H21)+ (k22*H22)))/((F*Cp)+(Fcat*Cpcat))) 


%UNTITLED Summary of this function goes here 
% Detailed explanation goes here 

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

function main 
options = odeset('RelTol',1e-6); %,'AbsTol',[1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 ]); 
Y0=[1.0;0.0;0.0;0.0;180000.0;959.81] 
zspan=0:0.5:33; 
[z,y]= ode45(@Untitled,zspan,Y0,options); 
figure 
hold on 
plot(z,y(:,1)); 
plot(z,y(:,2),':'); 

但我收到错误

??? Error using ==> feval 
Error: File: Untitled.m Line: 41 Column: 37 
()-indexing must appear last in an index expression. 
Error in ==> odearguments at 110 
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0. 
Error in ==> ode45 at 173 
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ... 
Error in ==> ode at 6 
[z,y]= ode45(@Untitled,zspan,Y0,options); 

任何人都可以请帮我这个,因为它是真的对我很重要,我是新来的MATLAB

thanx提前:)

回答

2

在代码中搜索)(,它在语法上不允许。可能你错过了一个*之间。

+0

你说得对,它在'ydot(6)'行 – am304