| 查看: 2533 | 回复: 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 |
回复此楼» 猜你喜欢
求标准粉末衍射卡号 ICDD 01-076-1802
已经有0人回复
-
新西兰Robinson研究所招收全奖PhD
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有59人回复
-
石墨烯转移--二氧化硅衬底石墨烯
已经有0人回复
-
笼目材料中量子自旋液体基态的证据
已经有0人回复
-
数学教学论硕士可以读数学物理博士吗？
已经有0人回复
-
德国亥姆霍兹Hereon中心汉堡分部招镁合金腐蚀裂变SCC课题方向2026公派博士生
已经有4人回复
-
澳门大学 应用物理及材料工程研究院 潘晖教授课题组诚招博士后
已经有11人回复
-
求助NH4V4O10晶体的CIF文件
已经有0人回复
-
英国全奖博士招聘-深度学习与量子物理
已经有0人回复
-
间接带隙半导体有效质量求助
已经有0人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 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