24小时热门版块排行榜

返回列表

zxqworld

新虫 (初入文坛)

应助: 0 (幼儿园)
帖子: 8
在线: 2.1小时
虫号: 2053599
注册: 2012-10-10
专业: 有机合成

[求助] 求大神帮忙看一下我写的积分法求取动力学参数的程式已有2人参与

function KineticsEst1_int
% 动力学参数辨识: 用积分法进行反应速率分析得到速率常数k和反应级数n
clear all
clc

global X
t =[30:30:120];
X1=[67.37 58.28 45.32 32.1];
X2=[7.88 17.88 17.45 26.43];
X3=[2.48 9.09 11.84 13.65];
X4=[1.58 0 0.94 0];
X5=[20.69 14.75 24.45 27.82];
% 非线性拟合
k0 = [0.00594 0.00391 0.00222 0.00810 0.00020 0.01071 ];
n0 =[1 1 1 1.43 1.38 1.26];
beta0=[k0; n0];
tspan = [30:120];
X0 = [67.37 7.88 2.48 1.58 20.69];
[beta,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@OptObjFunc,beta0,[],[],[],tspan,X0,X);
ci = nlparci(beta,resid,jacobian);

% 拟合效果图(实验与拟合的比较)
figure(1)
[t4plot X4plot] = ode45(@KineticsEqs,[tspan],X0,[],beta);
plot(t,X(:,1),'bo',t4plot,X4plot,'k-')
legend('Exp','Model')
xlabel('时间 min')
ylabel('收率C_C, %')

% 残差关于拟合值的残差图
[t Xc] = ode45(@KineticsEqs,tspan,X0,[],k,n);
figure
plot(Xc,resid,'*')
xlabel('收率拟合值（%）')
ylabel('残差R (%)')
refline(0,0)

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

% ------------------------------------------------------------------
function f = OptObjFunc(k,n,tspan,X0,X)
tspan=[30:30:120]
[t Xc] = ode45(@KineticsEqs,tspan,X0,[],k,n);
f1 = Xc(:,1) - X(:,1);
f2 = Xc(:,2) - X(:,2);
f3 = Xc(:,3) - X(:,3);
f4 = Xc(:,4) - X(:,4);
f5 = Xc(:,5) - X(:,5);
f = f1+ f2+ f3+ f4+ f5;

% ------------------------------------------------------------------
function dXdt = KineticsEqs(t,X,k,n)
dXdt = ...
[ ( -(k(1)+k(2)+k(3)+k(4))*X(1)^n(1) )
( k(2)*X(1)^n(2))
( k(1)*X(1)^n(1)+k(5)*X(4)^n(5)+k(6)*X(5)^n(6) )
( k(3)^X(1)^n(3)-k(5)*X(4)^n(5))
( k(4)*X(1)^n(4)-k(6)*X(5)^n(6))
];

程序错误好多，但是不知道怎么改，尤其是非线性拟合部分，关于矩阵beta，我是参照一本书上的一个例题写的，但是他只有一个方程，我有5个，还请各位大神帮忙看一下。。万分感谢！！

回复此楼