| 查看: 419 | 回复: 1 | ||
| 【悬赏金币】回答本帖问题，作者臣本布衣将赠送您 150 个金币 | ||
[求助]
求助matlab进行反应动力学拟合，大神帮忙看看代码哪里有问题已有1人参与
|
||
|
动力学方程如下： (dC_A)/dt=-(k_1+k_2)*K_A*C_A/(1+K_A*C_A+K_B*C_B+K_C*C_C+K_D*C_D+K_E*C_E) (dC_B)/dt=k_1*K_A*C_A/(1+K_A*C_A+K_B*C_B+K_C*C_C+K_D*C_D+K_E*C_E) (dC_C)/dt=(k_2*K_A*C_A-(k_3+k_4)*K_C*C_C+k_5*K_E*C_E)/(1+K_A*C_A+K_B*C_B+K_C*C_C+K_D*C_D+K_E*C_E) (dC_D)/dt=k_3*K_C*C_C/(1+K_A*C_A+K_B*C_B+K_C*C_C+K_D*C_D+K_E*C_E) (dC_E)/dt=(k_4*K_C*C_C-k_5*K_E*C_E)/(1+K_A*C_A+K_B*C_B+K_C*C_C+K_D*C_D+K_E*C_E) 用matlab语言的ode45解算子求上述微分方程组，以各物质浓度计算值和实验值的残差平方和为目标函数，用lsqnonline解算子求解上述最优化问题。 代码如下： function piadatfit clc; clear all; close all; %己烯异构化反应动力学数据拟合％输入实验数据 tspan=[0 30 36 45 60 90 180]; cexp=[0.84691 0 0 0 0; 0.65164 0 0.10305 0.02491 0.00485; 0.65409 0 0.11407 0.02995 0.00681; 0.60309 0 0.14508 0.03510 0.00949; 0.08397 0 0.65655 0.04303 0.01346; 0.05952 0.00752 0.61224 0.05047 0.01814; 0 0.05346 0.53471 0.06414 0.02753]; c0=[0.84691 0 0 0 0]; %输入参数初值及范围 k0=[0.1 2 0.05 0.05 0.05 0.3 0.3 0.3 0.1 0.1]; lb=[0 0 0 0 0 0 0 0 0 0]; ub=[10 10 10 10 10 10 10 10 10 10]; %非线性最小二乘拟合 options.algorithm = 'levenberg-marquardt'; [k,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@objpia,k0,lb,ub,options,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-',tspan,cexp(:,3),'g*',tplot,cplot(:,3),'g-',tspan,cexp(:,4),'rs',tplot,cplot(:,4),'r-',tspan,cexp(:,5),'cd',tplot,cplot(:,5),'c-') xlabel('time(min)') ylabel('concentration(mol/l)') legend('caexp','cacal','cbexp','cbcal','ccexp','cccal','cdexp','cdcal','ceexp','cecal') function f=objpia(k0,cexp,tspan,c0) %最小二乘拟合目标函数 [t,c]=ode45(@piakin,tspan,c0,[],k0); f1=c(:,1)-cexp(:,1); f2=c(:,2)-cexp(:,2); f=[f1;f2]; function dcdt=piakin(t,c,k) %反应动力学方程 r1=k(1)*k(6)*c(1)/(1+k(6)*c(1)+k(7)*c(2)+k(8)*c(3)+k(9)*c(4)+k(10)*c(5)); r2=k(2)*k(6)*c(1)/(1+k(6)*c(1)+k(7)*c(2)+k(8)*c(3)+k(9)*c(4)+k(10)*c(5)); r3=k(3)*k(8)*c(3)/(1+k(6)*c(1)+k(7)*c(2)+k(8)*c(3)+k(9)*c(4)+k(10)*c(5)); r4=k(4)*k(8)*c(3)/(1+k(6)*c(1)+k(7)*c(2)+k(8)*c(3)+k(9)*c(4)+k(10)*c(5)); r5=k(5)*k(10)*c(5)/(1+k(6)*c(1)+k(7)*c(2)+k(8)*c(3)+k(9)*c(4)+k(10)*c(5)); dcdt(1)=-r1-r2; dcdt(2)=r1; dcdt(3)=r2-r3-r4+r5; dcdt(4)=r3; dcdt(5)=r4-r5; dcdt=dcdt'; matlab小白，程序是模仿教材写的，应该存在不少问题，请大神帮忙修改一下！ |
回复此楼» 猜你喜欢
氧化硅薄膜玻璃 求购
已经有4人回复
-
vasp极化强度
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有141人回复
-
vasp计算的一些问题
已经有0人回复
-
国科温州研究院 项目聘用有意“二战”的优秀本科毕业生
已经有31人回复
-
SCAPS-1D学习交流
已经有2人回复
-
深圳大学跟哈工大联合博士
已经有8人回复
-
Journal of alloys and compounds期刊
已经有2人回复
-
请教一个概率问题
已经有0人回复
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15356.8
- 散金: 101
- 红花: 234
- 帖子: 3409
- 在线: 1222.7小时
- 虫号: 291104
- 注册: 2006-10-28
【答案】应助回帖
感谢参与，应助指数 +1
|
微分方程拟合问题，1stOpt试试，更简单、方便： Sum Squared Error (SSE): 0.257211769123463 Root of Mean Square Error (RMSE): 0.0925944147205909 Correlation Coef. (R): 0.867645804757432 R-Square: 0.752809242513172 Parameter Best Estimate --------- ------------- k1 0.00075219124176959 k2 0.0093680714188618 k3 0.0883840107805217 k4 9.99999999961713 k5 9.99999999989503 k6 9.99999999999999 k7 4.61483024965692E-14 k8 0.0268675102514994 k9 5.40027468998149E-19 k10 0.307240694154773 |
2楼2022-08-13 12:58:37
5