用matlab求动力学参数k1-k7,金币没有,太寒碜请见谅
希望木虫大神帮忙用Matlab解解方程,可以酬金支付,
t/h C1 C2 C3 C4 C5
0 100 0 0 0 0
0.5 18.7 48.9 19.3 9.7 0.2
1 6.5 31.3 34.7 19.6 2.6
1.5 2.7 15.3 41.2 28.7 5.7
2 2.1 6.9 40.9 34.5 9.3
3 1.7 2.2 37.2 38.4 17.1
4 1.5 1.8 31.7 38.2 24.8
5 1.2 1.7 26.1 36.3 30.5
6 0.2 0.5 19.4 33.4 36.7
微分方程组
dC1dt = -k1*C1-k7*C1;
dC2dt = k1*C1-k2*C2-k5*C2;
dC3dt = k2*C2-k3*C3-k6*C6;
dC4dt = k3*C3-k4*C4;
dC5dt = k4*C4;
dC6dt = k5*C2+k6*C3+k7*C1;
我网上依葫芦画瓢填的报错结果出不来
代码:
function odes_fit
format long
clear all
clc
k0 = [0 0 0 0 0 0];
lb = -[1 1 1 1 1 1]*1e9;
ub = [1 1 1 1 1 1]*1e9;
data=...
[0 100 0 0 0 0 0
0.5 18.7 48.9 19.3 9.7 0.2 0.0;
1 6.5 31.3 34.7 19.6 2.6 0.0;
1.5 2.7 15.3 41.2 28.7 5.7 0.05;
2 2.1 6.9 40.9 34.5 9.3 0.1;
3 1.7 2.2 37.2 38.4 17.1 0.12;
4 1.5 1.8 31.7 38.2 24.8 0.16;
5 1.2 1.7 26.1 36.3 30.5 0.22;
6 0.2 0.5 19.4 33.4 36.7 0.3;
];
x0=data(1,2:end);
tspan = [data(:,1)'];
yexp = [data(2:end,2) data(2:end,3) data(2:end,4) data(2:end,5) data(2:end,6) data(2:end,7)];
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.9f \n',k(1))
fprintf('\tk2 = %.9f \n',k(2))
fprintf('\tk3 = %.9f \n',k(3))
fprintf('\tk4 = %.9f \n',k(4))
fprintf('\tk5 = %.9f \n',k(5))
fprintf('\tk5 = %.9f \n',k(6))
fprintf('\tk5 = %.9f \n',k(7))
figure(1)
ts=0(max(tspan)-min(tspan))/100):max(tspan);
[ts ys] = ode45(@KineticsEqs,ts,x0,[],k);
yy = [data(:,2) data(:,3) data(:,4) data(:,5) data(:,6) data(:,7)];
figure(1)
plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo');
figure(2)
plot(ts,ys(:,2),'r',tspan,yy(:,2),'ro');
figure(3)
plot(ts,ys(:,3),'k',tspan,yy(:,3),'ko');
figure(4)
plot(ts,ys(:,4),'g',tspan,yy(:,4),'go');
figure(5)
plot(ts,ys(:,5),'m',tspan,yy(:,5),'mo');
figure(6)
plot(ts,ys(:,6),'h',tspan,yy(:,6),'ho');
figure(7)
plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo',ts,ys(:,2),'r',tspan,yy(:,2),'ro',ts,ys(:,3),'k',tspan,yy(:,3),'ko',ts,ys(:,4),'g',tspan,yy(:,4),'go',ts,ys(:,5),'m',tspan,yy(:,5),'mo',ts,ys(:,5),'h',tspan,yy(:,5),'ho'),
legend('C1的计算值','C1的实验值','C2的计算值','C2的实验值','C3的计算值','C3的实验值','C4的计算值','C4的实验值','C5的计算值','C5的实验值','C6的计算值','C6的实验值','Location','best');
function f = ObjFunc(k,tspan,x0,yexp) % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);
Xsim1=Xsim(:,1);
Xsim2=Xsim(:,2);
Xsim3=Xsim(:,3);
Xsim4=Xsim(:,4);
Xsim5=Xsim(:,5);
Xsim6=Xsim(:,6);
ysim(:,1) = Xsim1(2:end);
ysim(:,2) = Xsim2(2:end);
ysim(:,3) = Xsim3(2:end);
ysim(:,4) = Xsim4(2:end);
ysim(:,5) = Xsim5(2:end);
ysim(:,6) = Xsim6(2:end);
f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3)) (ysim(:,5)-yexp(:,5)) (ysim(:,6)-yexp(:,6)) ...
(ysim(:,6)-yexp(:,6))];
function dCdt = KineticsEqs(t,C,k) % ODE模型方程
C1=C(1);C2=C(2);C3=C(3);C4=C(4);
k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);k7=k(7);
dC1dt = -k1*C1-k7*C1;
dC2dt = k1*C1-k2*C2-k5*C2;
dC3dt = k2*C2-k3*C3-k6*C6;
dC4dt = k3*C3-k4*C4;
dC5dt = k4*C4;
dC6dt = k5*C2+k6*C3+k7*C1;
dCdt = [dC1dt;dC2dt;dC3dt;dC4dt;dC5dt;dC6dt;]; 返回小木虫查看更多
如方程没有错误的话,其k值为
k1=2.9355; k2=1.318; k3=0.44525; k4=0.2476; k5=2.2315e-14; k6=2.2214e-14; k7=0.47562
C1的拟合精度R^2=0.99696
C2的拟合精度R^2=0.98746
C3的拟合精度R^2=0.74181
C4的拟合精度R^2=0.71926
C5的拟合精度R^2=0.94598,
请问楼主问题解决了吗?方便交流交流吗?