24小时热门版块排行榜

返回列表

泡花alkali

新虫 (小有名气)

应助: 0 (幼儿园)
金币: 820.7
散金: 40
帖子: 106
在线: 49.6小时
虫号: 2633230
注册: 2013-09-03
专业: 高分子合成化学

[求助] 关于自由基聚合参数拟合的一个matlab代码，求大神看看错哪了？

我在跑一篇文献上给出的动力学模型，但是最后结果不太对，怀疑是其中的参数有问题，于是我打算用文献中的实验数据去拟合参数，模型中的物质有6种，文献中数据有浓度值的有一种（即单体转化率），参数一共有9个，所以我打算固定其中8个参数不变，去拟合1个参数kp，代码如下，但是跑不出来，求大神改改。
function KineticsMX180
clear all;clc
global t kp kd ki kt ka km a b f kp0 tspan
t=[0 7 10 20 30];
CAm=[4.45 3.231 2.759 1.936 1.317];
CA0=[4.45 0 0 0 0 0]';
kd=2.37e-4;                   % 初级自由基生成速率常数,1/min
ki=41.22;                      % 链引发速率常数,L/(mol·min)
kt=2.83e11;                   % 链终止速率常数,L/(mol·min)
ka=18.03;                      % 初级自由基迁入聚合物相速率,1/min
km=42.83;                      % 链转移速率常数,L/(mol·min)
a=0.4; b=0.4; f=0.8;
kp0=[4.25e8];
%beta0=[kd ki kt ka km a b f kp];
tspan=t;
[kp,resnorm,resid,exitflag,output,lambda,jacobian]=lsqnonlin(@OptObjFunc,kp0,[],[],[],tspan,CA0,CAm);
ci=nlparci(kp,resid,jacobian);
[t4plot CA4plot]=ode23(@KineticsEqs,[tspan(1) tspan(end)],CA0,[],kp);
[t CAc]=ode23(@KineticsEqs,tspan,CA0,[],kp);

function f=OptObjFunc(kp,tspan,CA0,CAm)
global t kp kd ki kt ka km a b f kp0 tspan
[t CAc]=ode23(@KineticsEqs,tspan,CA0,[],kp);
CAc
f=(CAm'-CAc);
function dCAdt=KineticsEqs(t,CA,kp)
global t kp kd ki kt ka km a b f kp0 tspan
kd=2.37e-4;                   % 初级自由基生成速率常数,1/min
ki=41.22;                      % 链引发速率常数,L/(mol·min)
kt=2.83e11;                   % 链终止速率常数,L/(mol·min)
ka=18.03;                      % 初级自由基迁入聚合物相速率,1/min
km=42.83;                      % 链转移速率常数,L/(mol·min)
a=0.4; b=0.4; f=0.8;
CA
c1=CA(1,: );
c2=CA(2,: );
c3=CA(3,: );                                  % C=k2*c1+k3*c2+k4*c3+k5*c4 Co=c6+c7+c8+c9
c4=CA(4,: );
c5=CA(5,: );
c6=CA(6,: );
dc1dt=-kd*c1^a*c6^b;
dc2dt=f*kd*c1^a*c6^b-ka*c2;
dc3dt=ka*c2-ki*c3*c5-kt*c3*c4;
dc4dt=ki*c3*c5-kt*c3*c4;
dc5dt=kp*c5*c4;
dc6dt=-3*kd*c1^a*c6^b;
dCAdt=[dc1dt;dc2dt;dc3dt;dc4dt;dc5dt;dc6dt];

其中c1~c6的：后加了空格是因为它会自动变成微笑表情。。

回复此楼