回 帖 发 新 帖 刷新版面

主题:[讨论]这是从一篇文章摘下的程序,我是一名maltlab初学者,因急需使用此程序,恳请帮助!

这是从一篇文章摘下的程序,我是一名maltlab初学者,因急需使用此程序,恳请帮助!

这有两个程序,上面的可以执行,下面的却无法运行,请大牛们帮忙看看,先谢了!


第一个程序
clear all
clc

global Cam
t=[0 4 8 12 16 20 24 28 32 36 40 44 48 52];
CAm=[0.006 0.047 0.511 2.070 3.170 3.944 4.283 4.236 4.163 3.965 3.804 3.681 3.645 3.585]'; %动力学数据

% 非线性拟合
beta0 = [0.5816 4.250]; 
tspan = [0 4 8 12 16 20 24 28 32 36 40 44 48 52];
CA0 = 0.006;
lb=[0.2 0.5]; ub=[0.6 4.0]
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@OptObjFunc5,beta0,lb,ub,[],tspan,CA0,CAm)
ci = nlparci(beta,resid,jacobian)

% 拟合效果图(实验与拟合的比较)
[t4plot CA4plot] = ode45(@KineticsEqs5,[tspan(1) tspan(end)],CA0,[],beta);
plot(t,CAm,'bo',t4plot,CA4plot,'k-')
legend('Exp','Model')
xlabel('时间t, s')
ylabel('浓度C_A, mol/L')

% 残差关于拟合值的残差图
[t CAc] = ode45(@KineticsEqs5,tspan,CA0,[],beta);
figure
plot(CAc,resid,'*')
xlabel('浓度拟合值(mol/L)')
ylabel('残差R (mol/L)')
refline(0,0)

% 参数辨识结果
fprintf('Estimated Parameters:\n')
fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tn = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2))

function f = OptObjFunc5(beta,tspan,CA0,CAm)
[t CAc] = ode45(@KineticsEqs5,tspan,CA0,[],beta);
f = CAc - CAm;

% ------------------------------------------------------------------
function dCAdt = KineticsEqs5(t,CA,beta)
dCAdt = beta(1)*(1-CA/beta(2))*CA; % k= beta(1), n= beta(2)


第二个程序
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)];

回复列表 (共1个回复)

沙发

第一个程序也不能运行啊

我来回复

您尚未登录,请登录后再回复。点此登录或注册