24小时热门版块排行榜    

查看: 2074  |  回复: 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的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +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的回帖

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的回帖

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的回帖

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 ...

另外再请教一个问题,1stopt能否像Origin非线性拟合时那样给出拟合的参数的标准误?
6楼2015-04-08 12:07:39
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 Mr_杰克 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 国自然面上复盘~欢迎讨论 (金币+15) +15 晴天加油 2026-08-26 16/800 2026-08-29 18:28 by symmetry
[论文投稿] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +3 gy1nBQXYQJqL 2026-08-29 4/200 2026-08-29 18:15 by 4FFAWE8HcgUD
[基金申请] 有没有仍没收到信息的 +5 德尚中行 2026-08-27 6/300 2026-08-29 18:10 by manplx
[公派出国] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +4 ASdOkHsho7FD 2026-08-28 6/300 2026-08-29 17:30 by 4FFAWE8HcgUD
[考博] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +5 ASdOkHsho7FD 2026-08-28 8/400 2026-08-29 17:22 by 4FFAWE8HcgUD
[基金申请] 为什么到现在没收到通知? +4 tannykie 2026-08-29 4/200 2026-08-29 17:14 by lmz0216
[硕博家园] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +4 ASdOkHsho7FD 2026-08-28 10/500 2026-08-29 16:49 by zICmwzsBXjbN
[公派出国] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 ASdOkHsho7FD 2026-08-28 4/200 2026-08-29 14:02 by jCd0dEvKHShX
[教师之家] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +4 ASdOkHsho7FD 2026-08-28 5/250 2026-08-29 11:49 by jCd0dEvKHShX
[基金申请] 中青基了要发朋友圈吗? +3 349506619 2026-08-28 3/150 2026-08-29 10:49 by 木风9012
[基金申请] 2026年叶企孙基金 +4 bud_bud 2026-08-27 7/350 2026-08-29 07:23 by foolishmani
[基金申请] 基金不中,共勉 +11 eulota 2026-08-26 11/550 2026-08-28 14:22 by 火星超人xi
[基金申请] 怎么查啊 +6 huang1991js 2026-08-26 6/300 2026-08-28 08:42 by winsaint
[基金申请] 面上合作单位盖章 +5 ssyjh 2026-08-27 5/250 2026-08-27 20:50 by gdfollow
[基金申请] 我不理解! +15 Edward_pc 2026-08-26 23/1150 2026-08-26 20:34 by zzuzxg
[基金申请] 为什么国自然不能直接公布 +4 bjdxyxy 2026-08-26 4/200 2026-08-26 13:12 by qingmu1201
[基金申请] 国际合作可查了,中了面上 (EPI+1)(金币+50) +18 Ldrop2023 2026-08-26 18/900 2026-08-26 11:15 by cmrandy
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
[基金申请] 今天务委会开完了,明天出结果吗 +19 angus9576 2026-08-25 23/1150 2026-08-26 10:03 by zp519
信息提示
请填处理意见