| 查看: 215 | 回复: 1 | |||
| 当前主题已经存档。 | |||
chshy1234银虫 (小有名气)
|
[交流]
【求助】，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 ] |
回复此楼» 猜你喜欢
面上本子正文33页，违规吗？会被低分嘛？
已经有12人回复
-
Sci. Bull. 悲剧经验
已经有7人回复
-
提交了我也来说说感想
已经有5人回复
-
今年审到国自然15份，谈谈感受
已经有20人回复
-
国自然评分
已经有4人回复
-
评审有感
已经有22人回复
-
反应很差，大量原料没有反应
已经有6人回复
-
如果工作了想读博，可以边工作边读全日制嘛？
已经有5人回复
-
如何实现卤原子转化
已经有8人回复
-
河北省自然科学基金
已经有3人回复
zajage
金虫 (著名写手)
- 应助: 7 (幼儿园)
- 金币: 4648.5
- 散金: 1
- 红花: 1
- 帖子: 2078
- 在线: 160.1小时
- 虫号: 244937
- 注册: 2006-04-22
- 性别: GG
- 专业: 药剂学
2楼2009-02-17 19:21:46
30