| 查看: 2696 | 回复: 6 | |||
1135725495铁杆木虫 (著名写手)
|
[求助]
求解微分方程组的参数 已有2人参与
|
|
求解K1，K2，a, m 方程组为： CB'=-K1*CB^a*(k1/k2*CB^a)^m CH'=K2*(K1/K2*CB^a)^m t CB CH 60 10.9237 0 90 10.7462 0 120 10.4357 0.03778 135 10.1695 0.0432 150 9.7481 0.1203 165 9.2346 0.21242 180 8.6613 0.34579 195 8.0058 0.56225 210 7.2423 0.83487 225 6.4188 1.12793 240 5.5353 1.38079 255 4.5768 1.869 270 4.0146 2.5 285 3.5703 3.01 300 3.1158 3.54452 330 2.4438 4.312 360 1.9878 4.70402 390 1.6668 4.8548 |
回复此楼» 猜你喜欢
PRB投稿，Acknowledgment sent to author状态十天了
已经有3人回复
-
PVA溶解
已经有2人回复
-
物理学I论文润色/翻译怎么收费?
已经有115人回复
-
基金申请
已经有44人回复
-
CSC与新西兰维多利亚大学PhD奖学金项目
已经有0人回复
-
新西兰Robinson研究所 招聘CSC公派访问人员
已经有0人回复
-
帮我的英语口语老师找学生
已经有0人回复
-
什么时候开奖？
已经有13人回复
-
散金币祈福
已经有94人回复
-
青基已中
已经有1人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- matlab拟合微分方程组中的参数
已经有5人回复
- 关于matlab求解常微分方程或者常微分方程组
已经有4人回复
- matlab 如何求解多元二阶微分方程组
已经有7人回复
- 微分方程ode45求解，最小二乘法优化微分方程参数，程序运行求助
已经有5人回复
- 隐式的微分方程组求解问题。
已经有7人回复
- assempde求解二维椭圆型偏微分方程组参数设置的问题
已经有0人回复
- matlab非线性最优化求解，微分方程组的拟合，参数估计
已经有2人回复
- Matlab常微分方程组求解问题
已经有4人回复
- 微分方程组参数拟合的问题
已经有11人回复
- matlab微分方程组参数拟合，以周为单位求解，汇总后以年为单位进行数值比较
已经有4人回复
- 求解偏微分方程组pdepe函数参数设置问题
已经有0人回复
- 求解四阶微分方程组
已经有14人回复
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
【答案】应助回帖
★ ★ ★ ★ ★
感谢参与，应助指数 +1
1135725495: 金币+5, ★有帮助 2016-09-11 17:37:52
感谢参与，应助指数 +1
1135725495: 金币+5, ★有帮助 2016-09-11 17:37:52
|
1stOpt试了下，或许还有更好的结果： 均方差(RMSE): 0.681320279985844 残差平方和(SSE): 15.7827090132796 相关系数(R): 0.971338834557611 相关系数之平方(R^2): 0.943499131519738 修正R平方(Adj. R^2): 0.930539243440337 确定系数(DC): 0.890053744503595 F统计(F-Statistic): 87.6923836038994 参数 最佳估算 -------------------- ------------- k1 0.00556462413275392 a -0.471720992027294 k2 0.000955248259076988 m 3.39014801069688 |
2楼2016-09-06 23:02:57
【答案】应助回帖
|
参考：http://www.forcal.net/yyhz/luwffcnh.htm 代码： !!!using["IMSL","luopt","math"]; //使用命名空间 f(t,CB,CH,dCB,dCH::K1,K2,a,m)= { dCB=-K1*CB^a*(K1/K2*CB^a)^m, dCH=K2*(K1/K2*CB^a)^m }; 目标函数(_K1,_K2,_a,_m : i,s,tyz : tyArray,tA,max,K1,K2,a,m)= { K1=_K1, K2=_K2, a=_a, m=_m, //传递优化变量，函数f中要用到K1,K2,a,m tyz=ode[@f,tA,ra1(10.9237,0)], i=0, s=0, while{++i<max, s=s+[tyz(i,1)-tyArray(i,1)]^2+[tyz(i,2)-tyArray(i,2)]^2}, s }; main(::tyArray,tA,max)= { tyArray=matrix{ //存放实验数据ti,yi "60 10.9237 0 90 10.7462 0 120 10.4357 0.03778 135 10.1695 0.0432 150 9.7481 0.1203 165 9.2346 0.21242 180 8.6613 0.34579 195 8.0058 0.56225 210 7.2423 0.83487 225 6.4188 1.12793 240 5.5353 1.38079 255 4.5768 1.869 270 4.0146 2.5 285 3.5703 3.01 300 3.1158 3.54452 330 2.4438 4.312 360 1.9878 4.70402 390 1.6668 4.8548" }, len[tyArray,0,&max], tA=tyArray(all:0), //用len函数取矩阵的行数，tA取矩阵的列 ClearImslErr(), //清空IMSL错误输出 ERSET(0,0,0), //关闭IMSL所有警告 Opt[@目标函数,optwaysimdeep, optwayconfra], //Opt函数全局优化 ERSET(0,2,2), ERSET(0,1,0) //恢复IMSL警告 }; 结果（K1,K2,a,m，目标函数值）： 2.007221403771439e-002 3.477931047132456e-002 0.7598653691173753 -1.205432420176749 18.90800296978626 |
3楼2016-09-16 17:31:34
|
带入参数后为什么是如下结果，不知哪里有什么问题？ 60. 10.9237 0. 90. 10.4374 0.255129 120. 9.89908 0.531102 135. 9.60554 0.678648 150. 9.29215 0.833842 165. 8.95519 0.997972 180. 8.5897 1.17273 195. 8.18881 1.36041 210. 7.74256 1.56429 225. 7.23555 1.78928 240. 6.64186 2.04341 255. 5.91159 2.34155 270. 4.92268 2.71809 285. 3.116 3.31463 300. -1.#IND -1.#IND 330. -1.#IND -1.#IND 360. -1.#IND -1.#IND 390. -1.#IND -1.#IND |
4楼2016-09-17 08:11:19
【答案】应助回帖
|
如果怀疑t=60时初值的准确性，可将t=60时的初值也作为拟合参数。 代码： !!!using["IMSL","luopt","math"]; //使用命名空间 f(t,CB,CH,dCB,dCH::K1,K2,a,m)= { dCB=-K1*CB^a*(K1/K2*CB^a)^m, dCH=K2*(K1/K2*CB^a)^m }; 目标函数(_K1,_K2,_a,_m,CB60,CH60 : tyz : tyArray,tA,K1,K2,a,m)= { K1=_K1, K2=_K2, a=_a, m=_m, //传递优化变量，函数f中要用到K1,K2,a,m tyz=ode[@f,tA,ra1(CB60 ,CH60)], //补充CB60,CH60为优化参数 sum[(tyArray-tyz).^2.0,0] }; main(::tyArray,tA)= { tyArray=matrix{ //存放实验数据ti,yi "60 10.9237 0 90 10.7462 0 120 10.4357 0.03778 135 10.1695 0.0432 150 9.7481 0.1203 165 9.2346 0.21242 180 8.6613 0.34579 195 8.0058 0.56225 210 7.2423 0.83487 225 6.4188 1.12793 240 5.5353 1.38079 255 4.5768 1.869 270 4.0146 2.5 285 3.5703 3.01 300 3.1158 3.54452 330 2.4438 4.312 360 1.9878 4.70402 390 1.6668 4.8548" }, tA=tyArray(all:0), //用tA取矩阵的列 ClearImslErr(), //清空IMSL错误输出 ERSET(0,0,0), //关闭IMSL所有警告 Opt[@目标函数,optwaysimdeep, optwayconfra,optrange, -1e5,1e5, -1e5,1e5, -1e5,1e5, -1e5,1e5, 10.0,15.0, 0.0,0.03778], //Opt函数全局优化 ERSET(0,2,2), ERSET(0,1,0) //恢复IMSL警告 }; 结果（K1,K2,a,m，CB60,CH60，目标函数值）： 1.323582077985482e-002 3.022295480133575e-002 0.9970827404206801 -0.8620485362713961 12.34254172021993 2.767575346307145e-008 13.27488182003417 |
5楼2016-09-17 10:07:08
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
1135725495: 金币+10, ★有帮助 2016-09-25 09:16:27
1135725495: 金币+10, ★有帮助 2016-09-25 09:16:27
|
用matlab验证了一下1stopt的结果： function dy=cbch(t,x) K1 = 0.00556462413275392; K2 = 0.000955248259076988; a = -0.471720992027294; m = 3.39014801069688; dy=zeros(2,1); dy(1)=-K1*x(1)^a*(K1/K2*x(1)^a)^m; dy(2)=K2*(K1/K2*x(1)^a)^m; >> [T,X]=ode45('cbch',[60, 90, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 330, 360 , 390],[10.9237, 0]); [T,X] 警告: 在 t=2.898865e+02 处失败。在时间 t 处，若不将步长降至允许的最小值(9.094947e-13)以下，积分公差要求无法满足。 > In ode45 (line 308) ans = 60.0000 10.9237 0 90.0000 10.4374 0.2551 120.0000 9.8991 0.5311 135.0000 9.6055 0.6786 150.0000 9.2921 0.8338 165.0000 8.9552 0.9980 180.0000 8.5897 1.1727 195.0000 8.1888 1.3604 210.0000 7.7426 1.5643 225.0000 7.2355 1.7893 240.0000 6.6419 2.0434 255.0000 5.9068 2.3427 270.0000 4.9277 2.7169 285.0000 3.1174 3.3143 300, 330, 360 , 390 这几个值没有给出 |
6楼2016-09-21 05:48:12
【答案】应助回帖
|
绘图： 3#应是最优解，但图形显示数据与曲线不一致，故怀疑数据与模型不匹配。 数据与模型匹配时参考：微分方程组参数拟合问题求助：http://muchong.com/bbs/viewthread.php?tid=10232397&fpage=1  图像.png |
7楼2016-09-24 05:02:19
30
