24小时热门版块排行榜

返回列表

当前主题已经存档。

chshy1234

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 384.3
散金: 18
帖子: 52
在线: 227.1小时
虫号: 615611
注册: 2008-10-02
性别: GG
专业: 化学工艺

[交流] 【求助】，MATLAB程序错在哪了？？？

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> chemical_kinetic at 17
[t,x]=ode45(@myfun,[t0:0.1:t1],x0,optimset('display','iter'),A);

运行下面程序出现上述错误：各位高手帮忙分析一下错在哪了？？？
function chemical_kinetic
%计算反应物随浓度的变化关系
clear all
clc
%输入初值
A0=1.50;
B0=0;C0=0;D0=0;
k1=1.00;k2=0.50;k3=0.30;k4=0.10;
x0=[B0 C0 D0];             %赋浓度初值
t0=0;t1=20;                %赋时间起始值
for f=0.5:0.1:1.0
switch f
      case 0.5
%计算反应速率方程中的系数矩阵
A=[-k1*f*A0 k2 -k3*f*A0 0;k1*f*A0 -k2 -k3*f*A0 0 ;0 0 k3*f*A0 -k4;0 0 0 k4];
%调用ode45函数计算常微分方程组
[t,x]=ode45(@myfun,[t0:0.1:t1],x0,optimset('display','iter'),A);
%作反应物浓度随时间的变化图
figure(1)
plot(t,x(:,1),'k-')
xlabel('t')
ylabel('B浓度')
title('B的浓度变化曲线(f=0.5)')
figure(2)
plot(t,x(:,2),'k-')
xlabel('t')
ylabel('C浓度')
title('C的浓度变化曲线(f=0.5)')
figure(3)
plot(t,x(:,3),'k-')
xlabel('t')
ylabel('D浓度')
title('D的浓度变化曲线(f=0.5)')
fprintf('\n f=0.5\n')
fprintf('\n B=%.4g\n',x(length(x),1))
fprintf('\n C=%.4g\n',x(length(x),2))
fprintf('\n D=%.4g\n',x(length(x),3))
   case 0.75
A=[-k1*f*A0 k2 -k3*f*A0 0;k1*f*A0 -k2 -k3*f*A0 0 ;0 0 k3*f*A0 -k4;0 0 0 k4];
[t,x]=ode45(@myfun,[t0:0.1:t1],x0,optimset('display','iter'),A);
figure(4)
plot(t,x(:,1),'k-')
xlabel('t')
ylabel('B浓度')
title('B的浓度变化曲线(f=1.6)')
figure(5)
plot(t,x(:,2),'k-')
xlabel('t')
ylabel('C浓度')
title('C的浓度变化曲线(f=1.6)')
figure(6)
plot(t,x(:,3),'k-')
xlabel('t')
ylabel('D浓度')
title('D的浓度变化曲线(f=1.6)')
fprintf('\n f=1.6\n')
fprintf('\n B=%.4g\n',x(length(x),1))
fprintf('\n C=%.4g\n',x(length(x),2))
fprintf('\n D=%.4g\n',x(length(x),3))
   case 1.0
A=[-k1*f*A0 k2 -k3*f*A0 0;k1*f*A0 -k2 -k3*f*A0 0 ;0 0 k3*f*A0 -k4;0 0 0 k4];
[t,x]=ode45(@myfun,[t0:0.1:t1],x0,optimset('display','iter'),A);
figure(7)
plot(t,x(:,1),'k-')
xlabel('t')
ylabel('B浓度')
title('B的浓度变化曲线(f=2.3)')
figure(8)
plot(t,x(:,2),'k-')
xlabel('t')
ylabel('C浓度')
title('C的浓度变化曲线(f=2.3)')
figure(9)
plot(t,x(:,3),'k-')
xlabel('t')
ylabel('D浓度')
title('D的浓度变化曲线(f=2.3)')
fprintf('\n f=2.3\n')
fprintf('\n B=%.4g\n',x(length(x),1))
fprintf('\n C=%.4g\n',x(length(x),2))
fprintf('\n D=%.4g\n',x(length(x),3))
end
end
%------------------
function y=myfun(t,x,A)
y=A*[1;x(1);x(1);x(2)];

[ Last edited by sunxiao on 2009-3-9 at 08:30 ]

回复此楼