24小时热门版块排行榜    

Znn3bq.jpeg
查看: 1840  |  回复: 5
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

Mr_杰克

新虫 (初入文坛)

[求助] [求助] 1stopt 微分方程代码求助 已有2人参与

这两天在用1stopt研究一个动力学微分方程的拟合求动力学参数,看起来应该是比较简单的方程,但是研究了很长时间,自己也根据网上能搜到的例子编了两个程序跑了下,但是结果相差比较大,我也不知道程序编的对不对,请懂的大神指点一二,跪谢了!我使用的是试用版的1stopt 5.0

方程是这样的:
dx/dt=-p1*x;
dy/dt=p1*x-p2*y;
初始值:t=0,此时x=3.5,y=0
测量的数据是t 和 y。具体数据请见程序中。p1和p2应为正。

程序1:(dx/dt=-p1*x   这个方程,我先积分成了x=x0*exp(-p1*t),也不知道对不对?)
//y'=dy/dt=3.5*p1*exp(-p1*t)-p2*y
Constant V =0; //Initial Values
Constant dt=1;
Parameter p(2);
Variable t, y;
StartProgram [Pascal];
Procedure MainModel;
var i,t: integer;
    y0 : double;
Begin
    y0 := V;
    for i := 0 to DataLength - 1 do begin
        y0 := (3.5*p1*exp(-p1*t)-p2*y0)*dt+y0;
        y := y0;

    end;
End;
EndProgram;
Data;
//0     0
10     0.010515
14     0.023065
17     0.037836
20     0.062683
22     0.063288
24     0.064385
26     0.063359
28     0.062942
30     0.061384
32     0.060358
35     0.065102
45     0.059101
60     0.043515

这个程序的输出报告和图:[求助] 1stopt 微分方程代码求助
                                          [求助] 1stopt 微分方程代码求助-1

关于程序1,我的主要疑惑是:1.不知道这个程序编的对不对?
                                                2.关于 dt 这个数值,在程序中的作用,我比较担心的是,我的t是不连续取值的,dt会不会使得t在程序中变成连续取值了?从图上看我的这一担心就更明显,因为横坐标只取到了12



程序2:
InitialODEValue t=0,x=3.5,y=0;
Parameters p1,p2;
Variable t[Output],y[Output];
ODEFunction x'=-p1*x;
            y'=p1*x-p2*y;

Data;
//0              0
10                0.010514645
14                0.023065468
17                0.037836003
20                0.062682935
22                0.06328842
24                0.064385481
26                0.063358784
28                0.062942392
30                0.061383892
32                0.060358109
35                0.065102196
45                0.059100765
60                0.043515354

这个程序的输出报告和图:     [求助] 1stopt 微分方程代码求助-2
                                               [求助] 1stopt 微分方程代码求助-3

关于程序2,我的主要疑惑是:1.不知道这个程序编的对不对?
                                                2.关于ODEFunction 语句的使用是否正确?ODEFunction语句对微分函数的定义默认是怎么样的?在函数中我都没定义t,他是怎么计算的?也涉及到t的取值问题,ODEFunction是不是自动把t按照1,2,3,4,5……取值了?


总结一下问题:1.请高手指导我这个动力学微分方程应该如何编程,能给代码就最好了;(授人以鱼,但是是最主要的。。。)
                        2.请指出我自己编写程序的问题,以及关于两个程序中t的取值,是否连续?还是间歇,如何确定(授人以渔)     
                        3.关于ODEFunction的定义,关于程序1中dt的定义。


谢谢!!!
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +1
CODE:
InitialODEValue t=0,x=3.5,y=0;
Parameters p1,p2;
Variable t,y;
ODEFunction x'=-p1*x;
            y'=p1*x-p2*y;
Data;
//0              0
10                0.010514645
14                0.023065468
17                0.037836003
20                0.062682935
22                0.06328842
24                0.064385481
26                0.063358784
28                0.062942392
30                0.061383892
32                0.060358109
35                0.065102196
45                0.059100765
60                0.043515354

结果稳定唯一,但效果不是很好:

Root of Mean Square Error (RMSE):0.0119983163868176
Sum of Square Error:0.00187147474953628
Correlation Coef. (R): 0.77730265287911
R-Square: 0.604199414172902
Determination Coef. (DC): 0.514455569529589
F-Statistic: 22.6549813151864

Parameter                  Best Estimate
--------------------        -------------
p1        0.00131507520277041
p2        0.0675323910740055
[求助] 1stopt 微分方程代码求助-4
c2.jpg

4楼2015-04-08 11:32:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 6 个回答

月只蓝

主管区长 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +1
去看看这位大神关于 微分方程拟合的回帖吧:http://muchong.com/bbs/space.php?uid=291104&view=post
MATLAB、MS小问题、普通问题请发帖求助!时间精力有限,恕不接受无偿私信求助。
2楼2015-04-08 10:40:32
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Mr_杰克

新虫 (初入文坛)

引用回帖:
2楼: Originally posted by 月只蓝 at 2015-04-08 10:40:32
去看看这位大神关于 微分方程拟合的回帖吧:http://muchong.com/bbs/space.php?uid=291104&view=post

谢谢版大关注!您是使用Matlab的是吧?

嗯,他的帖子也看了不少,但还是有不少疑惑,我已经站短他了,希望大神本人能现身说法吧。谢谢!
3楼2015-04-08 10:53:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Mr_杰克

新虫 (初入文坛)

引用回帖:
4楼: Originally posted by dingd at 2015-04-08 11:32:00
InitialODEValue t=0,x=3.5,y=0;
Parameters p1,p2;
Variable t,y;
ODEFunction x'=-p1*x;
            y'=p1*x-p2*y;
Data;
//0              0
10                0.010514645
14                0.023 ...

谢谢!跟我的程序2一样,结果也一样看来我编的是对的?那能麻烦大神解答下程序1为什么跑出来的结果不一样吗?
5楼2015-04-08 12:05:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考博] 求博导|生物质基多孔碳/超级电容方向,已有相关成果,寻能源材料/碳材料方向老师 +3 猪猪人Zzz 2026-04-12 3/150 2026-04-17 19:10 by 阳阳阳^_^
[教师之家] 山东双非院校考核超级无底线,领导幸灾乐祸,教师遭殃恐 +4 qut2026 2026-04-11 8/400 2026-04-17 16:10 by 会飞的猪157
[考研] 0854求调剂 +21 门路摸摸 2026-04-15 25/1250 2026-04-17 15:45 by qzxyhcsy
[论文投稿] 有没有接收比较快的sci期刊呀,最好在一个月之内的,研三孩子求毕业 20+4 之护着 2026-04-16 5/250 2026-04-17 10:02 by bobvan
[考研] 297,工科调剂? +4 河南农业大学-能 2026-04-14 4/200 2026-04-16 22:52 by wulijun2012
[考研] 一志愿沪9,生物学326求调剂 +9 刘墨墨 2026-04-15 9/450 2026-04-16 17:14 by 崔崔崔cccc
[基金申请] RY:中国产出的科学垃圾论文,绝对数量和比例都世界第一 +7 zju2000 2026-04-14 18/900 2026-04-16 11:36 by 欢乐颂叶蓁
[考研] 0854调剂 +13 长弓傲 2026-04-12 16/800 2026-04-15 13:45 by fenglj492
[考研] 310求调剂 +16 666真好 2026-04-11 18/900 2026-04-15 13:28 by 黑科技矿业
[考研] 297工科调剂? +14 河南农业大学-能 2026-04-13 15/750 2026-04-15 13:25 by 黑科技矿业
[考研] 366求调剂 +11 不知名的小卅 2026-04-11 11/550 2026-04-14 15:50 by zs92450
[考研] 考研调剂 +13 长弓傲 2026-04-13 14/700 2026-04-14 14:44 by zs92450
[考研] 食品与营养(0955)271求调剂 +15 升格阿达 2026-04-12 16/800 2026-04-14 13:18 by 浮若_安生
[考研] 085600材料与化工329分求调剂 +24 叶zilin 2026-04-13 25/1250 2026-04-14 09:20 by 试管破裂
[考研] 考研求调剂 +12 子木呐 2026-04-12 13/650 2026-04-14 01:19 by 王珺璞
[考研] 2026硕士调剂_能动_河南农业大学 +4 河南农业大学-能 2026-04-12 4/200 2026-04-13 22:01 by bljnqdcc
[考研] +10 李多米lee. 2026-04-12 11/550 2026-04-12 22:58 by yuyin1233
[考研] 一志愿浙大生物325分求调剂 +9 zysheng 2026-04-12 9/450 2026-04-12 22:31 by yuyin1233
[考研] 344 材料专业 求调剂211 无地域要求 +8 hualkop 2026-04-11 8/400 2026-04-12 22:24 by fqwang
[考研] 359求调剂 +5 胃痉挛累了 2026-04-11 5/250 2026-04-11 19:55 by lbsjt
信息提示
请填处理意见