| 查看: 4748 | 回复: 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 ] |
回复此楼» 猜你喜欢
零点能
已经有2人回复
-
科研求助
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有184人回复
-
【25申博】博导推荐求助
已经有6人回复
-
青岛大学李山东教授招生
已经有2人回复
-
新泽西理工大学Prof. Junjie Yang招实验凝聚态物理博士
已经有1人回复
-
新泽西理工大学Prof. Junjie Yang招实验凝聚态物理博士
已经有2人回复
-
2024年度“中德博士后”项目启动申请，截止时间6月30号
已经有0人回复
-
新泽西理工招实验凝聚态物理博士
已经有2人回复
-
24年博士招生
已经有4人回复
独孤神宇
版主 (知名作家)
- 应助: 489 (硕士)
- 贵宾: 0.008
- 金币: 29653.8
- 散金: 802
- 红花: 121
- 沙发: 1
- 帖子: 5364
- 在线: 850.6小时
- 虫号: 3522474
- 注册: 2014-11-06
- 性别: GG
- 专业: 机械动力学
- 管辖: 计算模拟
5楼2018-11-08 19:45:39
2楼2018-11-08 14:40:24
3楼2018-11-08 15:08:53
4楼2018-11-08 17:31:49
hzlhm
至尊木虫 (著名写手)
- 应助: 374 (硕士)
- 金币: 16279.7
- 红花: 52
- 帖子: 2586
- 在线: 572.8小时
- 虫号: 1517335
- 注册: 2011-11-30
- 性别: GG
- 专业: 常微分方程与动力系统
6楼2018-11-09 19:45:02
7楼2018-11-10 15:23:09
8楼2018-11-26 20:59:27
9楼2018-12-28 13:52:16
强强腹黑大人
木虫 (正式写手)
- 应助: 1 (幼儿园)
- 金币: 5053.1
- 散金: 533
- 红花: 1
- 帖子: 369
- 在线: 100.9小时
- 虫号: 3303102
- 注册: 2014-07-02
- 性别: GG
- 专业: 安全科学与工程
10楼2019-01-03 21:29:58
5