function piadatfit2
clear all;clc
tspan=[0 10 20 30 40 50 60];
ca0=1;
caexp=ca0*(1-[ 0.0;0.35;0.55;0.70;0.80;0.85;0.88]);
cbexp=ca0-caexp;
cexp=[caexp cbexp];
k0=[5 5];
c0=[1 0];
LB=[0 0];
UB=[+inf +inf];
[k,resnorm,residual]=lsqnonlin(@objpia,k0,LB,UB,[],cexp,tspan,c0)
[tplot cplot]=ode45(@piakin,tspan,c0,[],k);
plot(tspan,cexp(:,1),'bx',tplot,cplot(:,1),'b-',tspan,cexp(:,2),'ko',tplot,cplot(:,2),'k-')
function f=objpia(k,cexp,tspan,c0)
[t c]=ode45(@piakin,tspan,c0,[],k);
f1=c(:,1)-cexp(:,1);
f2=c(:,2)-cexp(:,2);
f=[f1;f2];
function dcdt=piakin(t,c,k)
dc1dt=-k(1)*c(1)^k(2);
dc2dt=k(1)*c(1)^k(2);
dcdt=[dc1dt dc2dt]';
我是假设(-ra)=k(1)*ca^(2)来做的,结果收敛还不错,拟合出k和n后积分可以求得ca与t的关系,或者直接看图也可,图片插入不了,你可以自己运行一下,下面是结果;
k =
0.0439 1.1570
resnorm =
7.5403e-004
residual =
0
0.0040
-0.0107
0.0020
0.0120
0.0016
-0.0097
0
-0.0040
0.0107
-0.0020
-0.0120
-0.0016
0.0097,