当前位置: 首页 > 计算模拟 >用1stopt拟合常微分方程参数

用1stopt拟合常微分方程参数

作者 plkolili
来源: 小木虫 300 6 举报帖子
+关注

求助各位大佬,帮忙用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;
 返回小木虫查看更多

今日热帖
  • 精华评论
  • dingd

    参考下,微分方程拟合比较费时间,不保证是最优解:

    均方差(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

  • dingd

    下面效果更好些:
    均方差(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

  • plkolili

    感谢dingd大神的回复。
    我想要得到的结果是考虑断裂的应力应变曲线图,类似下图所示。走势应该是先上升后下降。我用您得到的参数,在matlab计算的时候,有些不符合实际的曲线。
    可否让我看下拟合的结果图,再次感谢。
    用1stopt拟合常微分方程参数
    结果图.jpg

  • 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.这样的话计算的值可能就会更符合实际,但是我又不知道怎么设定他们的初值。

  • 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;

  • 慢跑的蜗牛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
                                   

    有没有好心人帮忙运行一下这个程序呀,谢谢啦~,

猜你喜欢