[lkyi7104]

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


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