| 查看: 4752 | 回复: 12 | ||
| 【悬赏金币】回答本帖问题，作者Lucky 鑫将赠送您 5 个金币 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
[求助]
matlab非线性微分方程组参数拟合（回归）已有4人参与
|
||
|
请问，下图这样的方程，利用非线性回归方法估计参数k和n，该如何创立m文件呢？ function dXdt = KineticsEqs(t,X,k,beta) global k dXdt = ... [(-k(1)*X(1)^beta(13)*X(2)^beta(1)-k(2)*X(1)^beta(2)*X(2)^beta(14)-k(3)*X(1)^beta(3)*X(4)^beta(15)-... k(4)*X(1)^beta(16)*X(5)^beta(4)+k(10)*X(4)^beta(10)*X(1)^beta(22)+ k(11)*X(3)^beta(11)*X(1)^beta(23)+k(5)*X(2)^beta(17)*X(1)^n(5)) (k(1)*X(1)^beta(13)*X(2)^beta(1)+k(6)*X(2)^beta(6)*X(4)^beta(18)+k(7)*X(5)^beta(7)*X(2)^beta(19)-... k(12)*X(3)^beta(12)*X(2)^beta(24)-k(5)*X(2)^beta(5)*X(1)^beta(17)) (k(2)*X(1)^beta(2)*X(2)^beta(14)+k(12)*X(3)^beta(12)*X(2)^beta(24)-k(9)*X(3)^beta(9)*X(5)^beta(21)-k(11)*X(3)^beta(11)*X(1)^beta(23)) (k(3)*X(1)^beta(3)*X(4)^beta(15)-k(10)*X(4)^beta(10)*X(1)^beta(22)-k(6)*X(2)^beta(6)*X(4)^beta(18)+k(8)*X(5)^beta(8)*X(4)^beta(20)) (k(4)*X(1)^beta(16)*X(5)^beta(4)-k(7)*X(5)^beta(7)*X(2)^beta(19)-k(8)*X(5)^beta(8)*X(4)^beta(20)+k(9)*X(3)^beta(9)*X(5)^beta(21)) ]; function f = OptObjFunc(k,tspan,beta,x0,yexp) % 目标函数 [t beta Xsim] = ode45(@KineticsEqs,tspan,x0,[],k,beta); ysim(:,1) = Xsim(2:end,1); ysim(:,2) = Xsim(2:end,2); ysim(:,3) = Xsim(2:end,3); ysim(:,4) = Xsim(2:end,4); f = [ysim(:,1)-yexp(:,1); ysim(:,2)-yexp(:,2); ysim(:,3)-yexp(:,3);ysim(:,4)-yexp(:,4)]; tspan=[30 60 90 120]; x0=[67.37 7.88 2.48 21.17 1.1]; k0=[1 1 1 1 1 1 1 1 1 1 1 1]; beta0=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; KineticsData; yexp = Kinetics(:,2:6); [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@OptObjFunc,k0,[],[],[],tspan,beta,x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1)) fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2)) fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3)) fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4)) fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5)) fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6)) fprintf('\tk7= %.4f ± %.4f\n',k(7),ci(7,2)-k(7)) fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8)) fprintf('\tk9 = %.4f ± %.4f\n',k(9),ci(9,2)-k(9)) fprintf('\tk10 = %.4f ± %.4f\n',k(10),ci(10,2)-k(10)) fprintf('\tk11 = %.4f ± %.4f\n',k(11),ci(11,2)-k(11)) fprintf('\tk12 = %.4f ± %.4f\n',k(12),ci(12,2)-k(12)) %时间t RE KE ES OT AL Kinetics=... [ 30 67.37 7.88 2.48 21.17 1.10 60 58.28 13.88 11.09 12.19 2.16 90 45.32 17.65 11.84 22.06 3.13 120 32.10 26.43 13.65 23.98 3.84 ]  发自小木虫Android客户端 [ last edited by 独孤神宇 on 2018-11-9 at 21:25 ] [ last edited by 独孤神宇 on 2018-11-9 at 21:25 ] [ Last edited by 独孤神宇 on 2018-11-9 at 21:27 ] |
回复此楼» 猜你喜欢
2024年度“中德博士后”项目启动申请，截止时间6月30号
已经有0人回复
-
新泽西理工招实验凝聚态物理博士
已经有2人回复
-
物理学I论文润色/翻译怎么收费?
已经有152人回复
-
24年博士招生
已经有10人回复
-
真空封石英管 北京
已经有5人回复
-
在NVIDIA GPU上运行Quantum Espresso
已经有0人回复
-
微信公众号“电介质Dielectrics”，报道了， 聊城大学王晓春教授团队优异面外压电材料
已经有0人回复
-
求助将cif文件转化为jade
已经有0人回复
-
金属氧化物薄膜空气退火和真空退火对表面粗糙度的影响
已经有0人回复
9楼2018-12-28 13:52:16
5