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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » [求助] 1stopt 微分方程代码求助
51/1返回列表
查看: 1780  |  回复: 5
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

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的定义。

谢谢！！！
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2015-04-08 10:28:48
已阅   回复此楼   关注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的回帖
查看全部 6 个回答
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭