应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜     石溪大学接受考研调剂申请>

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 程序代码 » 用1stopt拟合常微分方程参数
【调剂】北京石油化工学院2024年16个专业接受调剂
71/1返回列表
查看: 2067  |  回复: 6
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

plkolili

铁虫 (初入文坛)

[求助] 用1stopt拟合常微分方程参数已有1人参与

求助各位大佬，帮忙用1stopt高版本跑下面的代码，感激不尽！！！
title "canshunihe";
constant yield=0.229, k=60.328, mu=2.282,b=2.854, q=3.933, n=1.4,alpha=73.408, gamma0=5.751, beta=2.188, phai=0.141,e=1000;
//initialodevalue t=[0,8000],y1=[0,], y2=[0,], y3=[8,],y4=[0.229,],y5=[0.001,];
variable  t,y3,y4;
parameter  c1,c2,c3;
//odeoptions = [s=1,a=4,p=50]; //s: step size = 1, a: algorithm (1: euler, 2: rk2, 3: rk3, 4: rk4), p: population size = 5;
plot y4,y5;
odefunction y1'=((y4/(1-y5)-y2-yield)/c1)^c2*y3^(-mu);
            y2'=b*(q-y2)*y1';
            y3'=alpha*y3^(-gamma0)+beta*y1'*y3^(-phai);
            y4'=e*(1e-3-y1');
            y5'=y5^c3;
            
// t   y4
rowdata;
0,100,200,300,400,500,600,700,800,900,1000,1500,2000,3000,4000,5500,6500,7500,8000;
8,8.1,8.2,8.3,8.4,8.5,8.6,8.7,8.8,8.9,9.0,9.1,9.2,9.3,9.4,9.5,9.6,9.7,9.8;
0,2.63,2.95,2.99,3.52,3.81,4.05,4.25,4.43,4.6,4.72,5.07,5.28,5.2,5.01,4.65,4.2,3.55,1.04;
回复此楼

» 猜你喜欢
黑发不知勤学早，白首方悔读书迟。
1楼2019-01-13 18:23:19
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与，应助指数 +1
参考下，微分方程拟合比较费时间，不保证是最优解：

均方差(RMSE): 0.532681607689156
残差平方和(SSR): 10.2149890261309
相关系数(R): 0.932921869702266
相关系数之平方(R^2): 0.870343214968771
修正R平方(Adj. R^2): 0.845683647429779
确定系数(DC): 0.718147888636154
F统计(F-Statistic): 28.6980923726451

参数                  最佳估算
--------------------        -------------
c1        -3.77485429494E-5
c2        -0.341995932934685
c3        1.1560796736969
y1 初值         4.64067042228588
y2 初值         2.77131321062938
y5 初值         2.15540575432628E-14
赞一下 回复此楼
2楼2019-01-13 19:55:40
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
plkolili: 金币+10, ★★★★★最佳答案, 感谢百忙之中，抽出时间帮我计算。 2019-01-13 21:05:38
下面效果更好些：
均方差(RMSE): 0.42661365444081
残差平方和(SSR): 6.55197156559236
相关系数(R): 0.980897479289642
相关系数之平方(R^2): 0.962159864876774
修正R平方(Adj. R^2): 0.954431386315503
确定系数(DC): 0.778949690316621
F统计(F-Statistic): 22.6179953241421

参数                  最佳估算
--------------------        -------------
c1        1.19861624726491
c2        32.5275008605798
c3        2.24973236410188
y1 初值         -1.56792343709876
y2 初值         -0.697815532746156
y5 初值         0.00313654380605316
赞一下 回复此楼
3楼2019-01-13 20:13:47
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

plkolili

铁虫 (初入文坛)
感谢dingd大神的回复。
我想要得到的结果是考虑断裂的应力应变曲线图，类似下图所示。走势应该是先上升后下降。我用您得到的参数，在matlab计算的时候，有些不符合实际的曲线。
可否让我看下拟合的结果图，再次感谢。
用1stopt拟合常微分方程参数
结果图.jpg
赞一下 回复此楼
黑发不知勤学早，白首方悔读书迟。
4楼2019-01-13 21:04:36
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

plkolili

铁虫 (初入文坛)
引用回帖:
2楼: Originally posted by dingd at 2019-01-13 19:55:40
参考下，微分方程拟合比较费时间，不保证是最优解：

均方差(RMSE): 0.532681607689156
残差平方和(SSR): 10.2149890261309
相关系数(R): 0.932921869702266
相关系数之平方(R^2): 0.870343214968771
修正R平 ...

还有就是，我仅仅知道y1和y2的初值为0，y5的初值为0.001.这样的话计算的值可能就会更符合实际，但是我又不知道怎么设定他们的初值。
赞一下 回复此楼
黑发不知勤学早，白首方悔读书迟。
5楼2019-01-13 21:14:58
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

plkolili

铁虫 (初入文坛)
能在帮我把下面的代码算一下么？
title "canshunihe";
constant yield=0.229, k=60.328, mu=2.282,b=2.854, q=3.933, n=1.4,alpha=73.408, gamma0=5.751, beta=2.188, phai=0.141,e=1000;
initialodevalue t=0,y1=0, y2=0, y3=8,y4=0.229,y5=0.001;
variable  t,y3,y4;
parameter  c1,c2,c3;
//odeoptions = [s=1,a=4,p=50]; //s: step size = 1, a: algorithm (1: euler, 2: rk2, 3: rk3, 4: rk4), p: population size = 5;
plot y4,y5;
odefunction y1'=((y4/(1-y5)-y2-yield)/c1)^c2*y3^(-mu);
            y2'=b*(q-y2)*y1';
            y3'=alpha*y3^(-gamma0)+beta*y1'*y3^(-phai);
            y4'=e*(1e-3-y1');
            y5'=y5^c3;

// t   y4
rowdata;
t=0,100,200,300,400,500,600,700,800,900,1000,1500,2000,3000,4000,5500,6500,7500,8000;
y3=8,8.1,8.2,8.3,8.4,8.5,8.6,8.7,8.8,8.9,9.0,9.1,9.2,9.3,9.4,9.5,9.6,9.7,9.8;
y4=0.229,2.63,2.95,2.99,3.52,3.81,4.05,4.25,4.43,4.6,4.72,5.07,5.28,5.2,5.01,4.65,4.2,3.55,1.04;
赞一下 回复此楼
黑发不知勤学早，白首方悔读书迟。
6楼2019-01-14 10:57:07
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

慢跑的蜗牛wc

铁虫 (初入文坛)
Parameters k1,k2,k3,k4,k5,k6,k7;
InitialODEValue t=0,a=0,b=0,c=0,d=0,e=0.0000074;
Variable t,a,b,c,d,e；
ODEFunction e'=-k1*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5-k2*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5-k3*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5;
                       a'=6*k1*3051*e*172/0.000001343*(1-1.13*Ln(e/0.0000074))^0.5+2*k2*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5/0.000001343+k4*b*3051-k5*a*d-3*k6*a*b+3*k7*c*3051;
                       b'=4*k2*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5/0.000001343+2*k3*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5/0.000001343-k4*b*3051+k5*a*d-k6*a*b+k7*c*3051;
                       c'=k3*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5/0.000001343+k6*a*b-k7*c*3051;
                       d'=4*k1*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5/0.000001343+k3*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5/0.000001343+k4*3051*b-k5*d*a;
Data;
//t                 H2 a                           CO  b                        CH4 c                  CO2 d
60                5.317534543        2.282574774        0.076404295        9.720719792
300               16.53258605        0.921333731        0.460228388        14.65292426
600                21.35578457        0.748749046        0.707176158        17.85455396
1800        30.09578114        0.572388504        1.333230207        22.92109313
2400        32.48226677        0.445459631        1.986111389        26.69658144
                               

有没有好心人帮忙运行一下这个程序呀，谢谢啦~
赞一下 回复此楼
7楼2021-07-13 16:38:32
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 plkolili 的主题更新
71/1返回列表
普通表情 龙 兔 虎 猫 论文翻译-EnPapers美助理教授团队-SCI英文润色
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[找工作] 普通院校药学硕士，做合成的，感觉找不到工作 +10 pom戴墨镜 2024-04-24 17/850 2024-04-26 23:23 by ZZZemmm
[考博] 真的好想读博！ +15 wangzhe_bs 2024-04-22 22/1100 2024-04-26 22:11 by 小木雄子
[基金申请] 基金开始函评了吗？ +16 wych1103 2024-04-25 16/800 2024-04-26 21:32 by 淀粉搬运工
[考研] 没学上 +6 季向阳 2024-04-26 12/600 2024-04-26 21:06 by 季向阳
[论文投稿] 研二光催化6月底四篇二区什么水平 5+5 wjtab 2024-04-22 15/750 2024-04-26 19:25 by wjtab
[有机交流] 环肽的合成 +3 徐来不惊 2024-04-25 5/250 2024-04-26 16:56 by 徐来不惊
[考研] 学硕专硕 +5 小蜗牛* 2024-04-26 5/250 2024-04-26 16:43 by 鱼翔浅底1
[考研] 0854－0855调剂 +8 shangannum1 2024-04-21 12/600 2024-04-26 16:42 by yz仔
[基金申请] 两类问题算是白选了~ +7 jurkat.1640 2024-04-23 12/600 2024-04-26 14:39 by lucky_my2010
[教师之家] 博士论文被抄袭 +25 和尚敲小木鱼 2024-04-22 42/2100 2024-04-26 13:55 by ZHONGWU_U
[考研] 381求调剂 +4 小刺猬987654321 2024-04-25 6/300 2024-04-26 10:57 by czl12138
[基金申请] "颜宁:基础研究应顶天立地"能做到基础研究同时顶天立地的才是牛人 +5 zju2000 2024-04-24 5/250 2024-04-26 09:36 by LittleBush
[考博] 取博导收留 5+4 zzb777888 2024-04-20 10/500 2024-04-26 08:52 by polymerfriend
[论文投稿] 一直找不到审稿人 +5 lizhengke06 2024-04-21 6/300 2024-04-25 14:01 by chongdong
[基金申请] 国社科项目，你们学校都限额申报吗？ +7 屡战屡败 2024-04-21 10/500 2024-04-25 12:10 by 屡战屡败
[博后之家] 南京大学-广州大学联合招聘博士后 欢迎广大优秀人才！！！ +4 黑魔变身啾 2024-04-20 12/600 2024-04-25 11:18 by dodonaomi
[考博] 24年 申博 化学/材料 一作6篇sci +9 wangyp123 2024-04-23 11/550 2024-04-24 19:01 by bangbangbiu
[论文投稿] 期刊推荐 20+4 木颜尘ip 2024-04-22 7/350 2024-04-24 10:06 by bobvan
[考博] 研二光催化6月底4篇2区 +7 wjtab 2024-04-22 11/550 2024-04-23 06:59 by byron2012
[高分子] 请问UV灯是365nm的，那么选光引发剂的波长选多少的？要完全一致吗？ +4 engledd2004 2024-04-21 4/200 2024-04-22 16:08 by wangcz23
信息提示
关闭
请填处理意见
关闭