Labfans是一个针对大学生、工程师和科研工作者的技术社区。 论坛首页 | 联系我们(Contact Us)
MATLAB爱好者论坛-LabFans.com
返回   MATLAB爱好者论坛-LabFans.com > 工程数学软件 > MATLAB论坛
MATLAB论坛 一切MATLAB相关问题在此讨论。
回复
 
主题工具 显示模式
旧 2009-03-29, 17:51   #1
lkyi7104
初级会员
 
注册日期: 2009-03-12
年龄: 54
帖子: 3
声望力: 0
lkyi7104 正向着好的方向发展
默认 [求助]发酵动力学matlab程序不执行

请大家帮忙看一下,发酵动力学matlab程序不执行
程序如下:
k0=[0.5818 20.0957 0.028907]; %μmxa、Yx/s、beta参数初值
x0=[0.006 36.001 0.006]; %x0,s0,p0
t=[0 4 8 12 16 20 24 28 32 36 40 44 48 52];
tspan=[0 4 8 12 16 20 24 28 32 36 40 44 48 52];
%yexp实验数据[xl x2 x3]
yexp=[[0.006 0.047 0.511 2.070 3.170 3.944 4.283 4.236 4.1633 .965 3.804 3.681 3.645 3.585];
[36.001 35.788 34.735 32.215 27.187 20.593 10.802 5.850 3.129 2.251 1.850 1.499 1.111 0.960];
[0.006 0.007 0.008 0.009 0.030 0.139 0.605 1.033 1.621 2.251 2.544 2.840 3.249 3.300]]’:
lb=[0.2 0.05 0.01]; ub=[0.6 1 0.031]; %参数下、上限
%
[k, resnorm,resid,exitflag,output,lambda, jacobian]=…
lsnqonlin(@ObjFunc,k0,lb,ub,[],x0,yexp)
ci=nlparci(k, resid, jacoban)
%
yl=[0.006 0.047 0.511 2.070 3.170 3.944 4.283 4.236 4.1633 .965 3.804 3.681 3.645 3.585]';
y2=[36.001 35.788 34.735 32.215 27.187 20.593 10.802 5.850 3.129 2.251 1.850 1.499 1.111 0.960]';
y3=[0.006 0.007 0.008 0.009 0.030 0.139 0.605 1.033 1.621 2.251 2.544 2.840 3.249 3.300]';
[t4plot x4plot]=ode45(@kineticseqs,[tspan(1) tspan(end)],x0,[],k);
plot(t1,yl,'bo', tl,y2', 'g*',tl,y3,'r*', t4plot,x4Plot,'k-'),
legend('Exp-biomass','Exp-Na-citrate','Exp-产物活性单位','Model'),
xlabel('time(h)'),ylabel('菌体浓度,Na-citrate(g/l),产物活性单位(u)')
%
fprintf('Estimated Parameter\n')
fprintf('\tk=%.4f±%.4f\n', k(1), ci(1,2)-k(1))
fprintf('\tk=%.4f±%.4f\n', k(2), ci(2,2)-k(2))
fprintf('\tk=%.4f±%.4f\n', k(1), ci(1,2)-k(3))
%
[t x]=ode45(@kineticseqs,tspan,x0,[],k);
figure,plot(x,resid,'*')
xlabel('菌体、Na-citrate浓度拟合值(g/l)、产物活性单位(U)'),ylbael('残差R(g/l)')
M文件:
function f=ObjFunc (k,x0,yexp)
tspan=[0 4 8 12 16 20 24 28 32 36 40 44 48 52]';
[t x]=ode45(@kineticseqs,tspan,x0,[],k);
y(:,1)=x(:,l); y(:,2)=x(:,2); y(:,3)=x(:,3); %?
fl=y(:,l)-yexp(:,1); f2=y(:,2)-yexp(:,2); f3=y(:,3)-yexp(:,3);
f=[f1, f2, f3];
M文件:
function dxdt=kineticseqs (t,x,k) %模型方程
dxdt=[k(l)*x(l)*(1.0-x(l)/4.25044)
-k(l)/k(2)*x(l)*(1.0-x(l)/4.25044)
k(3)*x(l)];
lkyi7104 当前离线   回复时引用此帖
回复

主题工具
显示模式

发帖规则
不可以发表新主题
不可以发表回复
不可以上传附件
不可以编辑自己的帖子

启用 BB 代码
论坛启用 表情符号
论坛启用 [IMG] 代码
论坛禁用 HTML 代码


相似的主题
主题 主题作者 版面 回复 最后发表
如何实现gui窗口嵌套 beulah MATLAB论坛 1 2008-12-28 17:57


所有时间均为北京时间。现在的时间是 17:20


Powered by vBulletin
版权所有 ©2000 - 2025,Jelsoft Enterprises Ltd.