24小时热门版块排行榜    

查看: 2709  |  回复: 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
回复此楼

» 猜你喜欢

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

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

lu_yu_lan

新虫 (初入文坛)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
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
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 7 个回答

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★
感谢参与,应助指数 +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
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lu_yu_lan

新虫 (初入文坛)

【答案】应助回帖

参考: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
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lu_yu_lan

新虫 (初入文坛)

引用回帖:
2楼: Originally posted by dingd at 2016-09-06 23:02:57
1stOpt试了下,或许还有更好的结果:

均方差(RMSE): 0.681320279985844
残差平方和(SSE): 15.7827090132796
相关系数(R): 0.971338834557611
相关系数之平方(R^2): 0.943499131519738
修正R平方(Adj. R^2):  ...

带入参数后为什么是如下结果,不知哪里有什么问题?

            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
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +3 YLHlRHNwYkce 2026-09-02 6/300 2026-09-03 15:23 by T0rGB46095mJ
[找工作] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +3 YLHlRHNwYkce 2026-09-02 9/450 2026-09-03 15:21 by T0rGB46095mJ
[硕博家园] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +3 rFEsUBKXRll0 2026-09-02 5/250 2026-09-03 15:10 by T0rGB46095mJ
[教师之家] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +3 ero8OE6tv9cu 2026-09-02 8/400 2026-09-03 14:48 by T0rGB46095mJ
[考研] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +3 ero8OE6tv9cu 2026-09-02 9/450 2026-09-03 14:48 by T0rGB46095mJ
[基金申请] 国社科又开始会评了,不知道这次命运如何 +10 雨打竹帘 2026-08-30 14/700 2026-09-03 13:17 by qsd10086
[教师之家] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 YLHlRHNwYkce 2026-09-02 8/400 2026-09-03 12:03 by T0rGB46095mJ
[基金申请] 麻烦专家们看看评委们的意见(F口面上) +8 gdd2018 2026-08-28 13/650 2026-09-03 06:22 by 羽毛枫f
[公派出国] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +3 ero8OE6tv9cu 2026-09-02 4/200 2026-09-03 03:14 by rM1TE0WVDIIY
[基金申请] 科研人应该花精力去思考如何解决问题,而不是去凝练问题 +7 瞬息宇宙 2026-09-01 14/700 2026-09-02 17:33 by ma0526
[基金申请] 基金系统什么内容也没有 30+4 winsaint 2026-08-27 10/500 2026-09-02 11:35 by 大不刘6
[基金申请] 学科评审组评审是指会评吗? +5 瞬息宇宙 2026-08-31 5/250 2026-09-02 10:13 by 雨冰共舞
[基金申请] 面上函评意见出来了,像什么等级? 20+4 Tsingking1 2026-08-27 17/850 2026-09-01 19:51 by 超级无敌华子
[论文投稿] 小白求助 投论文要求的highlights应该如何写 5+3 l1963982152 2026-08-29 4/200 2026-09-01 09:04 by 北京莱茵编辑
[基金申请] 面上意见出来了 +12 黄鸟于飞Chao 2026-08-29 23/1150 2026-08-31 18:57 by 黄鸟于飞Chao
[基金申请] 能否申诉? +7 echo8914667 2026-08-30 8/400 2026-08-31 17:00 by yihongxu
[基金申请] 中青基了要发朋友圈吗? +7 349506619 2026-08-28 7/350 2026-08-31 13:39 by 冼亮淀粉酶
[考博] 找导师 +6 yuanjiabao 2026-08-29 7/350 2026-08-30 14:40 by 生科新手
[基金申请] 我就是申请一个面上项目而已,这评审意见是按照杰青的条件评的吧? +6 gouxfjh 2026-08-28 11/550 2026-08-30 07:57 by gouxfjh
[基金申请] 国自然评审意见 +13 wangmingqi 2026-08-28 19/950 2026-08-29 10:22 by Poppy1104
信息提示
请填处理意见