为您推荐:
精华内容
最热下载
问答
  • 5星
    1.44MB m0_52957036 2020-09-14 08:52:51
  • 5星
    198KB weixin_48549810 2021-01-06 14:45:40
  • 5星
    3KB qq_46753404 2021-03-15 22:29:06
  • 5星
    11KB weixin_42696333 2021-09-10 20:59:17
  • 5星
    3KB weixin_42696333 2021-09-10 21:17:12
  • 5星
    567KB SKCQTGZX 2021-09-03 15:03:25
  • 原问题地址:matlab 分式拟合,matlab 微分方程组的参数拟合_催眠神兔的博客-CSDN博客 微分方程式: x'=dx/dt=a*0.0321*(b-x)-d*x-dy/dt, y'= dy/dt=0.25*p1*exp(-p1*t)*x , 四个待求参数:a、b、d、p1 t、x、y...

    原问题地址:matlab 分式拟合,matlab 微分方程组的参数拟合_催眠神兔的博客-CSDN博客

    微分方程式:

    x'=dx/dt=a*0.0321*(b-x)-d*x-dy/dt,

    y'= dy/dt=0.25*p1*exp(-p1*t)*x ,

    四个待求参数:a、b、d、p1

    t、x、y数据见下面:

    //0 0 0  //这是初值

    0,0,0,

    0.1,0.486966799,0.048018378,

    0.167,1.6657,0.05823,

    0.2,0.860306078,0.060834243,

    0.3,1.156255213,0.064254733,

    0.4,1.390856542,0.065167644,

    0.5,1.67518,0.06638,

    0.6,1.724247244,0.065476325,

    0.7,1.841108525,0.065493681,

    0.8,1.93374543,0.065498314,

    0.9,2.007179471,0.06549955,

    1,1.92438,0.05641,

    1.1,2.111536196,0.065499968,

    1.2,2.148115682,0.065499991,

    1.3,2.177112544,0.065499998,

    1.4,2.200098596,0.065499999,

    1.5,2.218319829,0.0655,

    1.6,2.232763952,0.0655,

    1.7,2.244213928,0.0655,

    1.8,2.253290418,0.0655,

    1.9,2.260485427,0.0655,

    2,2.16156,0.07359,

    2.1,2.270710216,0.0655,

    2.2,2.274294246,0.0655,

    2.3,2.277135335,0.0655,

    2.4,2.27938749,0.0655,

    2.5,2.281172792,0.0655,

    2.6,2.282588016,0.0655,

    2.7,2.283709875,0.0655,

    2.8,2.284599183,0.0655,

    2.9,2.285304144,0.0655,

    3,2.44976,0.15449,

    3.1,2.286305961,0.0655,

    3.2,2.286657121,0.0655,

    3.3,2.286935489,0.0655,

    3.4,2.287156153,0.0655,

    3.5,2.287331076,0.0655,

    3.6,2.287469738,0.0655,

    3.7,2.287579657,0.0655,

    3.8,2.287666791,0.0655,

    3.9,2.287735862,0.0655,

    4,2.287790616,0.0655,

    4.1,2.287834019,0.0655,

    4.2,2.287868426,0.0655,

    4.3,2.2878957,0.0655,

    4.4,2.287917321,0.0655,

    4.5,2.287934459,0.0655,

    4.6,2.287948045,0.0655,

    4.7,2.287958815,0.0655,

    4.8,2.287967352,0.0655

    分析:

    微分方程中

    dx=a*0.0321*(b-x)-d*x-dy=0.0321*(a*b-a*x)-d*x-dy=0.0321*a*b-0.0321*a*x-d*x-dy=0.0321*ab-(0.0321*a-d)*x-dy=0.0321*ab-ad*x-dy

    等效于:
    dx=0.0321*ab-ad*x-dy

    故拟合参数a、b、d、p1中,前三个参数a、b、d不是独立的,如果进行拟合,将有无穷多最优解,故为过拟合现象。因而简化为三个拟合参数ab、ad、p1。

    Lu脚本代码(使用OpenLu+lugslmath计算,可从www.forcal.net下载):

    !!!using["luopt","math"]; //使用命名空间
    f(t,x,y,dx,dy, params :: ab,ad,p1)=
    {
        dy=0.25*p1*exp(-p1*t)*x,
        //dx=a*0.0321*(b-x)-d*x-dy, //dx=0.0321*(a*b-a*x)-d*x-dy=0.0321*a*b-0.0321*a*x-d*x-dy=0.0321*ab-(0.0321*a-d)*x-dy=0.0321*ab-ad*x-dy
        dx=0.0321*ab-ad*x-dy,
        0 //必须返回0
    };
    目标函数(_ab,_ad,_p1 : i,s,tyz : tyArray,tA,max, ab,ad,p1)=
    {
        ab=_ab, ad=_ad, p1=_p1, //传递优化变量
        //最后一个参数100表示gsl_ode函数在计算时,最多循环计算100次,这样可以提高速度
        tyz=gsl_ode[@f,nil,0.0,tA,ra1(0,0), 1e-6, 1e-6, gsl_rk8pd, 1e-6,100],
        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{ //存放实验数据t、x、y
            "0,0,0,
    0.1,0.486966799,0.048018378,
    .........省略数据
    4.8,2.287967352,0.0655"
        },
        len[tyArray,0,&max], tA=tyArray(all:0), //用len函数取矩阵的行数,tA取矩阵的列
        Opt1[@目标函数] //Opt1函数全局优化
        //Opt[@目标函数] //也可以使用此语句
    };

    结果(ab,ad,p1,最小值):

    196.3535183967741         2.775743871883714         20.18977085052178         0.9084940395557356

    以下Lu代码可绘制图形:

    !!!using["luopt","math","win"]; //使用命名空间
    f(t,x,y,dx,dy, params :: ab,ad,p1)=
    {
        dy=0.25*p1*exp(-p1*t)*x,
        //dx=a*0.0321*(b-x)-d*x-dy, //dx=0.0321*(a*b-a*x)-d*x-dy=0.0321*a*b-0.0321*a*x-d*x-dy=0.0321*ab-(0.0321*a-d)*x-dy=0.0321*ab-ad*x-dy
        dx=0.0321*ab-ad*x-dy,
        0 //必须返回0
    };
    目标函数(_ab,_ad,_p1 : i,s,tyz : tyArray,tA,max, ab,ad,p1)=
    {
        ab=_ab, ad=_ad, p1=_p1, //传递优化变量
        //最后一个参数100表示gsl_ode函数在计算时,最多循环计算100次,这样可以提高速度
        tyz=gsl_ode[@f,nil,0.0,tA,ra1(0,0), 1e-6, 1e-6, gsl_rk8pd, 1e-6,100],
        i=0, s=0, while{++i<max, s=s+[tyz(i,1)-tyArray(i,1)]^2+[tyz(i,2)-tyArray(i,2)]^2},
        s
    };
    init(main:tyz:tyArray,tA,max,ab,ad,p1)=
    {
        tyArray=matrix{ //存放实验数据t、x、y
            "0,0,0,
    0.1,0.486966799,0.048018378,
    
    ... ...省略数据
    
    4.8,2.287967352,0.0655"
        },
    
        len[tyArray,0,&max], tA=tyArray(all:0), //用len函数取矩阵的行数,tA取矩阵的列
        ab=0.0, ad=0.0, p1=0.0, Opt1[@目标函数, optstarter,&ab,&ad,&p1,0],   //Opt1函数全局优化
        tyz=gsl_ode[@f,nil,0.0,tA,ra1(0,0), 1e-6, 1e-6, 2, 1e-6,100],  //获得最终结果
        outa[tyz],                  //outa输出结果
        luShareX2[tyArray, tyz]     //绘制共享X轴视图
    };
    ChartWnd[@init];                //显示窗口并初始化

    图形:

    展开全文
    forcal 2021-07-31 04:56:26
  • time CLA CHMF-0.292699813 296.95 371.150 0 0.10420537730 0.004714244 0.09771414460 0.011371834 0.09162074190 0.017292886 ...

    time        CLA        CHMF

    -0.292699813        296.95        371.15

    0        0        0.104205377

    30        0.004714244        0.097714144

    60        0.011371834        0.091620741

    90        0.017292886        0.083395946

    120        0.024604931        0.076880763

    180        0.034326394        0.067098618

    240        0.043770544        0.057635083

    300        0.050232723        0.049835334

    420        0.0630261        0.035480465

    540        0.074273358        0.024649385

    -1.187859399        296.95        371.15

    0        0.000998896        0.104359628

    7        0.010982214        0.093943616

    15        0.021126176        0.082156616

    30        0.037770911        0.064592594

    45        0.049576229        0.050176205

    60        0.062465132        0.038429704

    90        0.070549574        0.025632208

    120        0.080625062        0.014099057

    150        0.086861135        0.008531942

    180        0.092640851        0.003482618

    -0.306744318        296.95        371.15

    0        0        1.064590325

    30        0.035820956        0.99951873

    60        0.091231251        0.937544756

    90        0.134419728        0.863743863

    120        0.175771334        0.812009778

    180        0.261121043        0.714875158

    240        0.322707477        0.612848744

    300        0.374201765        0.547687408

    360        0.43844506        0.459826754

    480        0.542643307        0.350285652

    -1.051427816        296.95        371.15

    0        0.005972939        1.009605954

    15        0.111613907        0.807385517

    30        0.239243562        0.653177373

    45        0.323494645        0.551727619

    60        0.401448109        0.416047561

    90        0.515684189        0.306615783

    120        0.611196296        0.16394207

    150        0.637330955        0.128679898

    180        0.657593008        0.075611085

    240        0.736999768        0.019444307

    -1.134063515        297.9        414.45

    0.883333333        0        0.103346707

    1.883333333        0.001512522        0.100892698

    2.883333333        0.01367608        0.088755571

    3.883333333        0.020184345        0.082041377

    4.883333333        0.044262628        0.054390229

    5.883333333        0.056989435        0.041914107

    6.883333333        0.077149736        0.020978875

    7.883333333        0.085874714        0.011700555

    8.883333333        0.090409481        0.002783166

    10.88333333        0.097292737        0

    -0.316500027        297.9        414.45

    0        0        1.070970111

    2        0.003948991        1.056056057

    3        0.020186072        1.037126132

    5        0.072206172        0.953098035

    7        0.193384568        0.788352458

    10        0.331724523        0.630124042

    13        0.46896552        0.434878839

    16        0.553992395        0.340920015

    20        0.64828259        0.221513537

    25        0.720649283        0.129803699

    -1.134691646        297.9        414.45

    0.883333333        0.000951425        1.015091849

    1.883333333        0.022369218        0.964249856

    2.883333333        0.109778689        0.841350168

    3.883333333        0.200071993        0.735440236

    4.883333333        0.265436739        0.647223865

    5.883333333        0.534716911        0.339015497

    7.883333333        0.746720695        0.087667964

    9.883333333        0.80490024        0.026375542

    11.88333333        0.814507102        0.003078599

    13.88333333        0.844622073        0

    -0.069553776        297.9        454.05

    0        0        0.105069261

    0.883333333        0        0.1033609

    2.883333333        0.001240168        0.100456984

    5.883333333        0.031816596        0.059602151

    7.883333333        0.041520146        0.046256377

    8.883333333        0.055889344        0.029139086

    9.883333333        0.057040816        0.028828737

    11.88333333        0.071843951        0.01235586

    15.88333333        0.077975291        0.001686074

    -0.132438303        297.9        454.05

    0        0        0.107494764

    0.883333333        0        0.107741493

    4.883333333        0.034722983        0.064167986

    6.883333333        0.069460981        0.023443174

    7.883333333        0.073060273        0.009201692

    8.883333333        0.07989269        0.007082479

    9.883333333        0.081267505        0.002823071

    10.88333333        0.087028096        0.001622815

    12.88333333        0.086598397        0

    -0.069742593        297.9        454.05

    0        0        1.045091436

    0.883333333        0        1.019821604

    2.883333333        0.027121965        0.980370094

    5.883333333        0.237319106        0.612502376

    7.883333333        0.327882492        0.463633654

    8.883333333        0.45866343        0.262292609

    9.883333333        0.476993327        0.245222172

    11.88333333        0.544647928        0.113480966

    15.88333333        0.599978991        0.03028414

    -0.132277891        297.9        454.05

    0        0        1.073088085

    1        0        1.042241346

    4.883333333        0.246301959        0.661045394

    6.883333333        0.517594712        0.290417203

    7.883333333        0.559102723        0.22683689

    8.883333333        0.662520515        0.082134805

    9.883333333        0.701897526        0.056834349

    10.88333333        0.723673938        0.009283515

    12.88333333        0.73410445        0.00186463

    -0.301540898        297.9        414.45

    0.883333333        0        0.101558639

    1.883333333        0.00563862        0.096383724

    3.883333333        0.017255774        0.084089401

    5.883333333        0.026288775        0.072388068

    7.883333333        0.039693626        0.058040049

    10.88333333        0.055995293        0.04294237

    13.88333333        0.06437459        0.032954514

    16.88333333        0.076142889        0.021339873

    20.88333333        0.083887249        0.010987306

    25.88333333        0.09228189        0

    -0.529462294        296.95000         371.15000

    0        0        0.104396538

    30        0.01086087        0.092633826

    60        0.025903593        0.077374404

    90        0.035707982        0.066417356

    120        0.046876147        0.055103546

    180        0.059555879        0.03988887

    240        0.071565642        0.027640828

    300        0.076601728        0.021842596

    420        0.088010322        0.010088388

    540        0.095707003        0.003371728

    -2.07036895        296.95000         371.15000

    0        0.005354629        0.101113608

    5        0.015596272        0.089273904

    10        0.031459676        0.070240021

    20        0.048850215        0.050642749

    30        0.064007378        0.034673716

    45        0.080303313        0.017800588

    60        0.083870308        0.014048253

    80        0.091867922        0.00226478

    100        0.096076813        0.000809425

    120        0.094754675        0

    -0.531088765        296.95000         371.15000

    0        0        1.085663867

    30        0.092227585        0.959687156

    60        0.196777828        0.815955815

    90        0.26771054        0.678294715

    120        0.339669208        0.584243991

    180        0.44493791        0.388930384

    240        0.546276325        0.31232313

    300        0.61758637        0.234488071

    360        0.656355053        0.157608349

    480        0.705466486        0.074621884

    -1.96250305        296.95000         371.15000

    0        0.009942844        0.943763649

    7        0.103997441        0.792109828

    15        0.26409945        0.602246938

    30        0.422141533        0.384870483

    45        0.559939068        0.204035091

    60        0.625774343        0.110517741

    90        0.704614433        0.03486973

    120        0.730085821        0.011077271

    180        0.738694061        0

    300        0.725429428        0

    -0.543478572        297.90000         414.45000

    0.883333333        0        0.103371092

    1.883333333        0        0.101492209

    2.883333333        0.010860825        0.0913963

    3.883333333        0.018569505        0.083314629

    4.883333333        0.031600628        0.069524342

    5.883333333        0.038345432        0.061930036

    6.883333333        0.05561085        0.044439081

    8.883333333        0.068330924        0.029366472

    10.88333333        0.079149507        0.017914315

    13.88333333        0.094041084        0.001441428

    -2.047061922        297.90000         414.45000

    0.883333333        0        0.103068737

    1.383333333        0.000594026        0.101344014

    1.883333333        0.004078897        0.099061111

    2.383333333        0.010968284        0.091356357

    2.883333333        0.025168999        0.075723863

    3.383333333        0.034499333        0.067241358

    3.883333333        0.056895681        0.044080416

    4.883333333        0.076428187        0.02052236

    5.883333333        0.079637269        0.01927162

    6.883333333        0.100845748        0

    -0.544306924        297.90000         414.45000

    0.883333333        0        1.009874204

    1.883333333        0        1.009745063

    2.883333333        0        1.010440284

    3.883333333        0.011090946        0.994311593

    4.883333333        0.047614676        0.926359288

    6.883333333        0.048051427        0.929046366

    8.883333333        0.127466089        0.817776336

    10.88333333        0.23769185        0.685748319

    12.88333333        0.300362249        0.620264623

    15.88333333        0.474651916        0.397967837

    -2.067390913        297.90000         414.45000

    0.883333333        0.003733023        0.971632292

    1.383333333        0.006034621        0.962671549

    1.883333333        0.039335525        0.916364421

    2.383333333        0.072289459        0.887838939

    2.883333333        0.077965623        0.874929815

    3.883333333        0.41942945        0.461344221

    4.883333333        0.68586173        0.136323751

    5.883333333        0.769533999        0.044770056

    6.883333333        0.807748957        0.002235684

    7.883333333        0.792673928        0

    展开全文
    weixin_36242516 2021-04-18 06:23:57
  • 5星
    1KB weixin_42683394 2021-09-29 03:36:33
  • Lu 微分方程参数优化(拟合) 例子1数学模型: dy/dt=k*y*z+0.095*b*z dz/dt=-b*z-0.222*z 实验数据(ti; yi): ti yi0.012.3291015630.683.8512927831.14.500939361.636.7491722472.079.1120628722.679....

    欢迎访问 Lu程序设计

    Lu 微分方程参数优化(拟合)

    例子1 数学模型:

          dy/dt=k*y*z+0.095*b*z
          dz/dt=-b*z-0.222*z

        实验数据(ti; yi):

     ti        yi

    0.01    2.329101563
    0.68    3.851292783
    1.1     4.50093936
    1.63    6.749172247
    2.07    9.112062872
    2.67    9.691657366
    3.09    11.16928013
    3.64    10.91451823
    4.65    16.44279204
    5.1     18.29615885
    5.58    21.63989025
    6.11    25.78611421
    6.63    26.34282633
    7.06    26.50581287
    7.62    27.63951823
    8.66    35.02757626
    9.04    35.5562035
    9.63    36.10393415

        已知z在t=0.01时的初值为5000,求k和b的最优值?

        Lu代码1:

    !!!using["luopt","math"]; //使用命名空间
    f(t,y,z,dy,dz, params::k,b)=
    {
        dy=k*y*z+0.095*b*z,
        dz=-b*z-0.222*z,
        0 
    //必须返回0
    };
    目标函数(_k,_b : i,s,tyz : tyArray,tA,max,k,b)=
    {
        k=_k,b=_b,
     //传递优化变量,函数f中要用到k,b
        //最后一个参数100表示gsl_ode函数在计算时,最多循环计算100次,这样可以提高速度
        tyz=gsl_ode[@f,nil,0.0,tA,ra1(2.329101563,5000), 1e-6, 1e-6, 2, 1e-6,100],
        i=0, s=0, while{++i<max, s=s+[tyz(i,1)-tyArray(i,1)]^2},
        s
    };
    main(::tyArray,tA,max)=
    {
        tyArray=matrix{
     //存放实验数据ti,yi
            "0.01 2.329101563
            0.68 3.851292783
            1.1 4.50093936
            1.63 6.749172247
            2.07 9.112062872
            2.67 9.691657366
            3.09 11.16928013
            3.64 10.91451823
            4.65 16.44279204
            5.1 18.29615885
            5.58 21.63989025
            6.11 25.78611421
            6.63 26.34282633
            7.06 26.50581287
            7.62 27.63951823
            8.66 35.02757626
            9.04 35.5562035
            9.63 36.10393415"
        },
        len[tyArray,0,&max], tA=tyArray(all:0),
     //用len函数取矩阵的行数,tA取矩阵的列
        Opt1[@目标函数] //Opt1函数全局优化
        //Opt[@目标函数] //也可以使用此语句
    };

        结果(前面的数为最优参数,最后一个数是极小值。下同):忽略Lu警告信息

    1.408312606232919e-004 -1.465812213944895e-004 24.29124283367381

        Lu代码2:优化(拟合)并绘图

    !!!using["luopt","math","win"]; //使用命名空间
    f(t,y,z,dy,dz,params::k,b)=
    {
        dy=k*y*z+0.095*b*z,
        dz=-b*z-0.222*z,
        0
    };
    目标函数(_k,_b : i,s,tyz : tyArray,tA,max,k,b)=
    {
        k=_k,b=_b,       //传递优化变量,函数f中要用到k,b
        tyz=gsl_ode[@f,nil,0.0,tA,ra1(2.329101563,5000), 1e-6, 1e-6, 2, 1e-6,100],
        i=0, s=0, while{++i<max, s=s+[tyz(i,1)-tyArray(i,1)]^2},
        s
    };
    init(main:tyz:tyArray,tA,max,k,b)=
    {
        tyArray=matrix{            //存放实验数据ti,yi
            "0.01 2.329101563
            0.68 3.851292783
            1.1  4.50093936
            1.63 6.749172247
            2.07 9.112062872
            2.67 9.691657366
            3.09 11.16928013
            3.64 10.91451823
            4.65 16.44279204
            5.1  18.29615885
            5.58 21.63989025
            6.11 25.78611421
            6.63 26.34282633
            7.06 26.50581287
            7.62 27.63951823
            8.66 35.02757626
            9.04 35.5562035
            9.63 36.10393415"
        },
        len[tyArray,0,&max], tA=tyArray(all:0), //用len函数取矩阵的行数,tA取矩阵的列
        k=0.0, b=0.0, Opt1[@目标函数, optstarter,&k,&b,0],   //Opt1函数全局优化
        tyz=gsl_ode[@f,nil,0.0,tyArray(all:0),ra1(2.329101563,5000), 1e-6, 1e-6, 2, 1e-6,100],  //获得最终结果
        outa[tyz],                  //outa输出结果
        luShareX2[tyArray, tyz]     //绘制共享X轴视图
    };
    ChartWnd[@init];                //显示窗口并初始化

        结果:

     

    例子2 数学模型:需要拟合的微分方程组为:

            df1/dt=a1*(x-f2)^a2,
            df2/dt=a3*df1/dt-a4*f2^a5

        其中:x=0.3。拟合参数为a1,a2,a3,a4,a5。试验获取的数据为(t,f1);t=0时,f2=0.15。

     t       f1
    0   , 0        ,
    10  , 0.002174 ,
    20  , 0.002663 ,
    30  , 0.002934 ,
    40  , 0.003113 ,
    50  , 0.003248 ,
    60  , 0.003354 ,
    70  , 0.003442 ,
    80  , 0.003514 ,
    90  , 0.003578 ,
    100 , 0.003635 ,
    110 , 0.003686 ,
    120 , 0.00373 ,
    130 , 0.003774 ,
    140 , 0.003813 ,
    150 , 0.003847 ,
    160 , 0.003882 ,
    170 , 0.003913 ,
    180 , 0.003942 ,
    190 , 0.00397  ,
    200 , 0.003996 ,
    210 , 0.00402  ,
    220 , 0.004044 ,
    230 , 0.004067 ,
    240 , 0.004087 ,
    250 , 0.004107 ,
    260 , 0.004126 ,
    270 , 0.004143 ,
    280 , 0.004159 ,
    290 , 0.004174 ,
    300 , 0.00419  ,
    310 , 0.004207 ,
    320 , 0.00422  ,
    330 , 0.004235 ,
    340 , 0.004248 ,
    350 , 0.004263 ,
    360 , 0.004276 ,
    370 , 0.004287 ,
    380 , 0.004301 ,
    390 , 0.00431  ,
    400 , 0.004323 ,
    410 , 0.004333 ,
    420 , 0.004344 ,
    430 , 0.004352 ,
    440 , 0.004362 ,
    450 , 0.004375 ,
    460 , 0.004383 ,
    470 , 0.00439  ,
    480 , 0.004399 ,
    490 , 0.004408 ,
    500 , 0.004416 ,
    510 , 0.004424 ,
    520 , 0.00443  ,
    530 , 0.00444  ,
    540 , 0.004449 ,
    550 , 0.004453 ,
    560 , 0.004458 ,
    570 , 0.004468 ,
    580 , 0.004474 ,
    590 , 0.004483 ,
    600 , 0.004488 ,
    610 , 0.004494 ,
    620 , 0.004499 ,
    630 , 0.004505 ,
    640 , 0.00451  ,
    650 , 0.004516 ,
    660 , 0.004522 ,
    670 , 0.004527 ,
    680 , 0.004532 ,
    690 , 0.004537 ,
    700 , 0.004541 ,
    710 , 0.004548 ,
    720 , 0.004552 ,
    730 , 0.004557 ,
    740 , 0.004561 ,
    750 , 0.004566 ,
    760 , 0.004569 ,
    770 , 0.004575 ,
    780 , 0.004579 ,
    790 , 0.004584 ,
    800 , 0.004589 ,
    810 , 0.004594 ,
    820 , 0.004596 ,
    830 , 0.004601 ,
    840 , 0.004604 ,
    850 , 0.004609 ,
    860 , 0.004611 ,
    870 , 0.004614 ,
    880 , 0.00462  ,
    890 , 0.004622 ,
    900 , 0.004626 ,
    910 , 0.00463  ,
    920 , 0.004632 ,
    930 , 0.004636 ,
    940 , 0.004638 ,
    950 , 0.004641 ,
    960 , 0.004647 ,
    970 , 0.004649 ,
    980 , 0.004653 ,
    990 , 0.004655 ,
    1000, 0.004658

        Lu代码1:

    !!!using["luopt","math"];
    f(t,f1,f2,df1,df2,pp : x : a1,a2,a3,a4,a5)=
    {
        x=0.3,
        df1=a1*(x-f2)^a2,
        df2=a3*df1-a4*f2^a5,
        0
    };
    J(_a1,_a2,_a3,_a4,_a5 : tf,i,s : tArray,tA,max,a1,a2,a3,a4,a5)=
    {
        a1=_a1, a2=_a2, a3=_a3, a4=_a4, a5=_a5,
        tf=gsl_ode[@f,nil,0.0,tA,
    ra1(0,0.15), 1e-6, 1e-6, 2, 1e-6,100],
        i=0, s=0, while{++i<max, s=s+[tf(i,1)-tArray(i,1)]^2},
        s
    };
    main(::tArray,tA,max)=
    {
        max=101,
        tArray=matrix[max,2].SetArray{
            0 ,  0 ,
            10 , 0.002174 ,
            20 , 0.002663 ,
            30 , 0.002934 ,
            40 , 0.003113 ,
            50 , 0.003248 ,
            60 , 0.003354 ,
            70 , 0.003442 ,
            80 , 0.003514 ,
            90 , 0.003578 ,
            100 , 0.003635 ,
            110 , 0.003686 ,
            120 , 0.00373 ,
            130 , 0.003774 ,
            140 , 0.003813 ,
            150 , 0.003847 ,
            160 , 0.003882 ,
            170 , 0.003913 ,
            180 , 0.003942 ,
            190 , 0.00397 ,
            200 , 0.003996 ,
            210 , 0.00402 ,
            220 , 0.004044 ,
            230 , 0.004067 ,
            240 , 0.004087 ,
            250 , 0.004107 ,
            260 , 0.004126 ,
            270 , 0.004143 ,
            280 , 0.004159 ,
            290 , 0.004174 ,
            300 , 0.00419 ,
            310 , 0.004207 ,
            320 , 0.00422 ,
            330 , 0.004235 ,
            340 , 0.004248 ,
            350 , 0.004263 ,
            360 , 0.004276 ,
            370 , 0.004287 ,
            380 , 0.004301 ,
            390 , 0.00431 ,
            400 , 0.004323 ,
            410 , 0.004333 ,
            420 , 0.004344 ,
            430 , 0.004352 ,
            440 , 0.004362 ,
            450 , 0.004375 ,
            460 , 0.004383 ,
            470 , 0.00439 ,
            480 , 0.004399 ,
            490 , 0.004408 ,
            500 , 0.004416 ,
            510 , 0.004424 ,
            520 , 0.00443 ,
            530 , 0.00444 ,
            540 , 0.004449 ,
            550 , 0.004453 ,
            560 , 0.004458 ,
            570 , 0.004468 ,
            580 , 0.004474 ,
            590 , 0.004483 ,
            600 , 0.004488 ,
            610 , 0.004494 ,
            620 , 0.004499 ,
            630 , 0.004505 ,
            640 , 0.00451 ,
            650 , 0.004516 ,
            660 , 0.004522 ,
            670 , 0.004527 ,
            680 , 0.004532 ,
            690 , 0.004537 ,
            700 , 0.004541 ,
            710 , 0.004548 ,
            720 , 0.004552 ,
            730 , 0.004557 ,
            740 , 0.004561 ,
            750 , 0.004566 ,
            760 , 0.004569 ,
            770 , 0.004575 ,
            780 , 0.004579 ,
            790 , 0.004584 ,
            800 , 0.004589 ,
            810 , 0.004594 ,
            820 , 0.004596 ,
            830 , 0.004601 ,
            840 , 0.004604 ,
            850 , 0.004609 ,
            860 , 0.004611 ,
            870 , 0.004614 ,
            880 , 0.00462 ,
            890 , 0.004622 ,
            900 , 0.004626 ,
            910 , 0.00463 ,
            920 , 0.004632 ,
            930 , 0.004636 ,
            940 , 0.004638 ,
            950 , 0.004641 ,
            960 , 0.004647 ,
            970 , 0.004649 ,
            980 , 0.004653 ,
            990 , 0.004655 ,
            1000, 0.004658
        },
        tA=tArray(all,0),
        Opt1[@J] 
        //Opt1函数全局优化
        //Opt[@目标函数] //也可以使用此语句

    };

        结果(没有多次运行,或许还有更优解):

    17082.05250017119 8.719683388872223 18.73041824823782 -2.443865629726595e-004 2.783118159118608 1.426677323284488e-009

        Lu代码2:优化(拟合)并绘图

    !!!using["luopt","math","win"];
    f(t,f1,f2,df1,df2,pp : x : a1,a2,a3,a4,a5)=
    {
        x=0.3,
        df1=a1*(x-f2)^a2,
        df2=a3*df1-a4*f2^a5,
        0
    };
    J(_a1,_a2,_a3,_a4,_a5 : tf,i,s : tArray,tA,max,a1,a2,a3,a4,a5)=
    {
        a1=_a1, a2=_a2, a3=_a3, a4=_a4, a5=_a5,
        tf=gsl_ode[@f,nil,0.0,tA,
    ra1(0,0.15), 1e-6, 1e-6, 2, 1e-6,100],
        i=0, s=0, while{++i<max, s=s+[tf(i,1)-tArray(i,1)]^2},
        s
    };
    init(main:tf:tArray,max,tA,a1,a2,a3,a4,a5)=
    {
        max=101,
        tArray=matrix[max,2].SetArray{
            0 ,  0 ,
            10 , 0.002174 ,
            20 , 0.002663 ,
            30 , 0.002934 ,
            40 , 0.003113 ,
            50 , 0.003248 ,
            60 , 0.003354 ,
            70 , 0.003442 ,
            80 , 0.003514 ,
            90 , 0.003578 ,
            100 , 0.003635 ,
            110 , 0.003686 ,
            120 , 0.00373 ,
            130 , 0.003774 ,
            140 , 0.003813 ,
            150 , 0.003847 ,
            160 , 0.003882 ,
            170 , 0.003913 ,
            180 , 0.003942 ,
            190 , 0.00397 ,
            200 , 0.003996 ,
            210 , 0.00402 ,
            220 , 0.004044 ,
            230 , 0.004067 ,
            240 , 0.004087 ,
            250 , 0.004107 ,
            260 , 0.004126 ,
            270 , 0.004143 ,
            280 , 0.004159 ,
            290 , 0.004174 ,
            300 , 0.00419 ,
            310 , 0.004207 ,
            320 , 0.00422 ,
            330 , 0.004235 ,
            340 , 0.004248 ,
            350 , 0.004263 ,
            360 , 0.004276 ,
            370 , 0.004287 ,
            380 , 0.004301 ,
            390 , 0.00431 ,
            400 , 0.004323 ,
            410 , 0.004333 ,
            420 , 0.004344 ,
            430 , 0.004352 ,
            440 , 0.004362 ,
            450 , 0.004375 ,
            460 , 0.004383 ,
            470 , 0.00439 ,
            480 , 0.004399 ,
            490 , 0.004408 ,
            500 , 0.004416 ,
            510 , 0.004424 ,
            520 , 0.00443 ,
            530 , 0.00444 ,
            540 , 0.004449 ,
            550 , 0.004453 ,
            560 , 0.004458 ,
            570 , 0.004468 ,
            580 , 0.004474 ,
            590 , 0.004483 ,
            600 , 0.004488 ,
            610 , 0.004494 ,
            620 , 0.004499 ,
            630 , 0.004505 ,
            640 , 0.00451 ,
            650 , 0.004516 ,
            660 , 0.004522 ,
            670 , 0.004527 ,
            680 , 0.004532 ,
            690 , 0.004537 ,
            700 , 0.004541 ,
            710 , 0.004548 ,
            720 , 0.004552 ,
            730 , 0.004557 ,
            740 , 0.004561 ,
            750 , 0.004566 ,
            760 , 0.004569 ,
            770 , 0.004575 ,
            780 , 0.004579 ,
            790 , 0.004584 ,
            800 , 0.004589 ,
            810 , 0.004594 ,
            820 , 0.004596 ,
            830 , 0.004601 ,
            840 , 0.004604 ,
            850 , 0.004609 ,
            860 , 0.004611 ,
            870 , 0.004614 ,
            880 , 0.00462 ,
            890 , 0.004622 ,
            900 , 0.004626 ,
            910 , 0.00463 ,
            920 , 0.004632 ,
            930 , 0.004636 ,
            940 , 0.004638 ,
            950 , 0.004641 ,
            960 , 0.004647 ,
            970 , 0.004649 ,
            980 , 0.004653 ,
            990 , 0.004655 ,
            1000, 0.004658
        },
        tA=tArray
    (all:0),
        a1=0.0, a2=0.0, a3=0.0, a4=0.0, a5=0.0,
        Opt1[@J, optstarter,&a1,&a2,&a3,&a4,&a5,0],
      //Opt1函数全局优化
        tf=gsl_ode[@f,nil,0.0,tA,ra1(0,0.15), 1e-6, 1e-6, 2, 1e-6,100],  //获得最终结果
        outa[tf],                   //输出结果
        luShareX2[tArray, tf]       //绘制共享X轴视图
    };
    ChartWnd[@init];               
     //显示窗口并初始化

        结果:

     

    例子3 数学模型:

    有关于y1,y2,y3,y4的微分方程:

    k1 = 3.5*10^4;
    k2 = 1.0*10^3;
    k3 = a*ka+b*kb+c*kc;

    r1 = k1*y1*y2;
    r2 = k2*y1*y3;
    r3 = k3*y1*y4;

    dy1/dt = -r1-r2-r3;
    dy2/dt = -r1;
    dy3/dt = -r2+r1;
    dy4/dt = -r3+r2;


    其中ka,kb,kc为拟合参数;常数a、b、c为已知数据;k1,k2,k3,r1,r2,r3为中间变量;t=0时,y1=30.0*10^-6; y2=10.2*10^-6; y3=4.1*10^-6; y4=6.7。

    有4组数据:

    当[a b c]=[100 5 20]时
    t=[0, 10, 20, 30, 40]
    y1=[30.0*10^-6, 17.5*10^-6, 15*10^-6, 14*10^-6, 13.5*10^-6 ]

    当[a b c]=[120 15 25]时
    t=[0 10 20 30 40]
    y1=[30.0*10^-6, 12.5*10^-6, 10*10^-6, 8*10^-6, 7*10^-6 ]

    当[a b c]=[130 25 30]时
    t=[0 10 20 30 40]
    y1=[30.0*10^-6, 11.5*10^-6, 9*10^-6, 7.5*10^-6, 6*10^-6 ]

    当[a b c]=[140 30 35]时
    t=[0 10 20 30 40]
    y1=[30.0*10^-6, 10.5*10^-6, 8.5*10^-6, 6.5*10^-6, 5*10^-6 ]

    求ka,kb,kc?

    分析:观察式子 k3 = a*ka+b*kb+c*kc ,当[a b c]实验数据少于3组时,将有无穷组最优的ka,kb,kc。现在有4组数据,故只有一组最优解。

    Lu代码:

    !!!using["luopt","math"];
    f(t,y1,y2,y3,y4,z1,z2,z3,z4,pp : k1,k2,k3,r1,r2,r3 : a,b,c,ka,kb,kc)=
    {
        k1 = 3.5e4,
        k2 = 1.0e3,
        k3 = a*ka+b*kb+c*kc,

        r1 = k1*y1*y2,
        r2 = k2*y1*y3,
        r3 = k3*y1*y4,

        z1 = -r1-r2-r3,
        z2 = -r1,
        z3 = -r2+r1,
        z4 = -r3+r2,

        0
    };
    目标函数(kka,kkb,kkc : i,s,yy : tArray,y11Array,y12Array,y13Array,y14Array,max,a,b,c,ka,kb,kc)=
    {
        ka=kka, kb=kkb, kc=kkc,
        s=0,

        a=100, b=5, c=20, 
    //[a b c]=[100 5 20]
        yy=gsl_ode[@f,nil,0.0,tArray,ra1(30.0e-6, 10.2e-6, 4.1e-6, 6.7), 1e-6, 1e-6, 2, 1e-6,100],
        i=0, while{++i<max, s=s+[yy(i,1)-y11Array(i)]^2},

        a=120, b=15, c=25, 
    //[a b c]=[120, 15, 25]
        yy=gsl_ode[@f,nil,0.0,tArray,ra1(30.0e-6, 10.2e-6, 4.1e-6, 6.7), 1e-6, 1e-6, 2, 1e-6,100],
        i=0, while{++i<max, s=s+[yy(i,1)-y12Array(i)]^2},

        a=130, b=25, c=30,
     //[a b c]=[130, 25, 30]
        yy=gsl_ode[@f,nil,0.0,tArray,ra1(30.0e-6, 10.2e-6, 4.1e-6, 6.7), 1e-6, 1e-6, 2, 1e-6,100],
        i=0, while{++i<max, s=s+[yy(i,1)-y13Array(i)]^2},

        a=140, b=30, c=35,
     //[a b c]=[140, 30, 35]
        yy=gsl_ode[@f,nil,0.0,tArray,ra1(30.0e-6, 10.2e-6, 4.1e-6, 6.7), 1e-6, 1e-6, 2, 1e-6,100],
        i=0, while{++i<max, s=s+[yy(i,1)-y14Array(i)]^2},

        sqrt[s/(max*4-4)]
    };
    main(::tArray,y11Array,y12Array,y13Array,y14Array,max)=
    {
        max=5,
        tArray=ra1[0, 10, 20, 30, 40],
        y11Array=ra1[30.0e-6, 17.5e-6, 15e-6, 14e-6, 13.5e-6],
        y12Array=ra1[30.0e-6, 12.5e-6, 10e-6, 8e-6, 7e-6],
        y13Array=ra1[30.0e-6, 11.5e-6, 9e-6, 7.5e-6, 6e-6],
        y14Array=ra1[30.0e-6, 10.5e-6, 8.5e-6, 6.5e-6, 5e-6],
        Opt1[@目标函数, optrange,-1e5,1e5,-1e5,1e5,-1e5,1e5] 
    //Opt1函数优化
        //Opt[@目标函数, optrange,-1e5,1e5,-1e5,1e5,-1e5,1e5] //也可以使用此语句
    };

    结果(ka,kb,kc,目标函数终值):

    9.15291371378818e-005 3.357752241197708e-004 -5.37879598667895e-004 1.178551159055296e-006

    验证(t,实验值,拟合值):

    当[a b c]=[100 5 20]时

    10. 1.75e-005 1.73918e-005
    20. 1.5e-005  1.54895e-005
    30. 1.4e-005  1.40236e-005
    40. 1.35e-005 1.28548e-005

    当[a b c]=[120 15 25]时

    10. 1.25e-005 1.45484e-005
    20. 1.e-005   1.09127e-005
    30. 8.e-006   8.29305e-006
    40. 7.e-006   6.35421e-006

    当[a b c]=[130 25 30]时

    10. 1.15e-005 1.29864e-005
    20. 9.e-006   8.73647e-006
    30. 7.5e-006  5.94617e-006
    40. 6.e-006   4.07251e-006

    当[a b c]=[140 30 35]时

    10. 1.05e-005 1.30756e-005
    20. 8.5e-006  8.85422e-006
    30. 6.5e-006  6.06631e-006
    40. 5.e-006   4.18281e-006

    例子4 数学模型:

    微分方程组为:

    dn1/dt = a*n1*n2-b*n1;
    dn2/dt = c*n1*n2;


    其中a,b,c为拟合参数;t=0时,n1=n0; n2=n0。

    数据:

    t    1    2     3     4     5     6

    n1   1    3.5   2.8   2.4   1.5   0.9

    分析:由于n0未知,故共有4个拟合参数:n0,a,b,c。

    Lu代码:

    !!!using["luopt","math"];
    f(t,n1,n2,dn1,dn2,pp ::a,b,c)=
    {
        dn1=a*n1*n2-b*n1,
        dn2=c*n1*n2,
        0
    };
    目标函数(n0,aa,bb,cc : i,s,nn : tArray,n1Array,max,a,b,c)=
    {
        a=aa, b=bb, c=cc,
        nn=gsl_ode[@f,nil,0.0,tArray,
    ra1(n0, n0), 1e-6, 1e-6, 2, 1e-6,100],
        i=0, s=0, while{++i<max, s=s+[nn(i,1)-n1Array(i)]^2},
        sqrt[s/(max-1)]
    };
    main(::tArray,n1Array,max)=
    {
        max=7,      
    //补充数据t=0,n1取任意值(本例n1=0),故共有7组数据
        tArray=ra1{0, 1, 2, 3, 4, 5, 6},
        n1Array=ra1{0, 1, 3.5, 2.8, 2.4, 1.5, 0.9},
        Opt1[@目标函数, optwaysimdeep, optwayconfra, optrange,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10] 
    //Opt1函数优化
        //Opt[@目标函数, optwaysimdeep, optwayconfra, optrange,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10] //也可以使用此语句
    };

    结果(n0,a,b,c,误差):

    4.842109347483411e-002 74.58161099198458 0.3525934220396638 -0.7028769944300635 0.1446429325392884

    验证(t,实验值,拟合值):

    1. 1.   1.00662
    2. 3.5  3.47867
    3. 2.8  2.96124
    4. 2.4  2.12753
    5. 1.5  1.50219
    6. 0.9  1.05743

    例子5 动力学方程参数估计:

    微分方程组为:dC/dt=-k0*exp{-Ea/[R*(273.15+T)]}*(A^m)*(B^n)*(C^o)

    单次实验T,A,B都是固定的,而C随时间变化并测量出来,所以T,A,B可以说是可变常数,R=8.314,拟合参数为k0,Ea,m,n,o,实验数据如下:

    T(℃)  A  B(MPa)  t/min  C

    260 0.05 4 0   3759.7
    260 0.05 4 10  2772.45
    260 0.05 4 20  2100
    260 0.05 4 30  1679.03
    260 0.05 4 40  1356.12
    260 0.05 4 50  1122.72
    260 0.05 4 60  818.6
    260 0.05 4 70  643.28
    260 0.05 4 80  485.26
    260 0.05 4 90  361.16
    260 0.05 4 100 290.17

    260 0.1  4 0   3759.7
    260 0.1  4 10  2450
    260 0.1  4 20  1803.37
    260 0.1  4 30  1246.87
    260 0.1  4 40  895
    260 0.1  4 50  496.14
    260 0.1  4 60  410.13
    260 0.1  4 70  179.72
    260 0.1  4 80  142.38
    260 0.1  4 90  94.73
    260 0.1  4 100 93.51

    260 0.15 4 0   3759.7
    260 0.15 4 10  2300
    260 0.15 4 20  1424.68
    260 0.15 4 30  884.85
    260 0.15 4 40  546.36
    260 0.15 4 50  336.83
    260 0.15 4 60  203.33
    260 0.15 4 70  117.97
    260 0.15 4 80  66.95
    260 0.15 4 90  54.49
    260 0.15 4 100 49.41

    260 0.2  4 0   3759.7
    260 0.2  4 10  2150
    260 0.2  4 20  1200
    260 0.2  4 30  650
    260 0.2  4 40  380
    260 0.2  4 50  230
    260 0.2  4 60  130
    260 0.2  4 70  80
    260 0.2  4 80  65
    260 0.2  4 90  50
    260 0.2  4 100 45

    200 0.1  4 0   3759.7
    200 0.1  4 10  2893.56
    200 0.1  4 20  2037.91
    200 0.1  4 30  1319.34
    200 0.1  4 40  1005.13
    200 0.1  4 50  773.12
    200 0.1  4 60  562.77
    200 0.1  4 70  384.27
    200 0.1  4 80  267.1
    200 0.1  4 90  171.38
    200 0.1  4 100 157.51

    220 0.1  4 0   3759.7
    220 0.1  4 10  2600
    220 0.1  4 20  1691.77
    220 0.1  4 30  1126.16
    220 0.1  4 40  902.03
    220 0.1  4 50  628.67
    220 0.1  4 60  371.19
    220 0.1  4 70  240.16
    220 0.1  4 80  139.84
    220 0.1  4 90  116.87
    220 0.1  4 100 83.8

    240 0.1  4 0   3759.7
    240 0.1  4 10  2350.51
    240 0.1  4 20  1450
    240 0.1  4 30  1000
    240 0.1  4 40  708.77
    240 0.1  4 50  375.68
    240 0.1  4 60  271.72
    240 0.1  4 70  175.88
    240 0.1  4 80  123.3
    240 0.1  4 90  90
    240 0.1  4 100 66.69

    260 0.1  4 0   3759.7
    260 0.1  4 10  2080
    260 0.1  4 20  1200
    260 0.1  4 30  884.85
    260 0.1  4 40  546.36
    260 0.1  4 50  336.83
    260 0.1  4 60  203.33
    260 0.1  4 70  117.97
    260 0.1  4 80  66.95
    260 0.1  4 90  54.49
    260 0.1  4 100 49.41

    280 0.1  4 0   3759.7
    280 0.1  4 10  1935.79
    280 0.1  4 20  1007.44
    280 0.1  4 30  758.08
    280 0.1  4 40  500
    280 0.1  4 50  315.13
    280 0.1  4 60  190.44
    280 0.1  4 70  100
    280 0.1  4 80  50
    280 0.1  4 90  40
    280 0.1  4 100 30

    260 0.1  4 0   3759.7
    260 0.1  4 10  2080
    260 0.1  4 20  1200
    260 0.1  4 30  749.04
    260 0.1  4 40  500
    260 0.1  4 50  280
    260 0.1  4 60  203.33
    260 0.1  4 70  117.97
    260 0.1  4 80  66.95
    260 0.1  4 90  54.49
    260 0.1  4 100 49.41

    260 0.1  3.5 0   3759.7
    260 0.1  3.5 10  2244.58
    260 0.1  3.5 20  1484.05
    260 0.1  3.5 30  749.04
    260 0.1  3.5 40  547.23
    260 0.1  3.5 50  331.18
    260 0.1  3.5 60  267.75
    260 0.1  3.5 70  198.36
    260 0.1  3.5 80  130.51
    260 0.1  3.5 90  100.47
    260 0.1  3.5 100 93.91

    260 0.1  3 0   3759.7
    260 0.1  3 10  2483.29
    260 0.1  3 20  1669.27
    260 0.1  3 30  1162.42
    260 0.1  3 40  827.86
    260 0.1  3 50  517.85
    260 0.1  3 60  357.71
    260 0.1  3 70  274.21
    260 0.1  3 80  185.53
    260 0.1  3 90  129.69
    260 0.1  3 100 88.35

    260 0.1  2.5 0   3759.7
    260 0.1  2.5 10  2527.4
    260 0.1  2.5 20  1857.92
    260 0.1  2.5 30  1219.77
    260 0.1  2.5 40  850.2
    260 0.1  2.5 50  628.67
    260 0.1  2.5 60  515.01
    260 0.1  2.5 70  359.75
    260 0.1  2.5 80  265.34
    260 0.1  2.5 90  172.36
    260 0.1  2.5 100 149.93

    260 0.1  2 0   3759.7
    260 0.1  2 10  2722.63
    260 0.1  2 20  2089.13
    260 0.1  2 30  1485.29
    260 0.1  2 40  1062.05
    260 0.1  2 50  831.48
    260 0.1  2 60  766.41
    260 0.1  2 70  547.1
    260 0.1  2 80  461.41
    260 0.1  2 90  347.28
    260 0.1  2 100 244.06

    分析:一共有14组实验数据,另外注意到时间序列都一样[0,10,20,30,40,50,60,70,80,90,100],这样可以简化代码。

    Lu代码:

    !!!using["luopt","math"];        //使用命名空间
    f(t,C,dC,pp::k0,Ea,m,n,o,T,A,B)= dC=-k0*exp{-Ea/[8.314*(273.15+T)]}*(A^m)*(B^n)*(C^o), 0;
    g(kk0,EEa,mm,nn,oo : i,j,s,tCC : tT,tA,tB,tTime,tC,max,tmax,k0,Ea,m,n,o,T,A,B)=
    {
        k0=kk0, Ea=EEa, m=mm, n=nn, o=oo,   //传递优化变量
        s=0, i=0,
        while{i<max,
            T=tT[i], A=tA[i], B=tB[i],
            tCC=gsl_ode[@f,nil,0.0,tTime,ra1(tC(0,i)), 1e-6, 1e-6, 2, 1e-6,100],
            j=0, while{++j<tmax, s=s+[tCC(j,1)-tC(j,i)]^2},
            i++
        },
        s
    };
    main(::tT,tA,tB,tTime,tC,max,tmax)=
    {
        max=14,                    //实验数据组数
        tT=ra1[260,260,260,260,200,220,240,260,280,260,260,260,260,260],
        tA=ra1[0.05,0.1,0.15,0.2,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1],
        tB=ra1[4,4,4,4,4,4,4,4,4,4,3.5,3,2.5,2],
        tmax=11,                   //时间序列
        tTime=ra1[0,10,20,30,40,50,60,70,80,90,100],
        tC=matrix[tmax,max].SetArray{ //存放实验数据Ci
    "3759.7   3759.7    3759.7    3759.7  3759.7    3759.7    3759.7    3759.7    3759.7    3759.7    3759.7    3759.7    3759.7    3759.7
    2772.45   2450      2300      2150    2893.56   2600      2350.51   2080      1935.79   2080      2244.58   2483.29   2527.4    2722.63
    2100      1803.37   1424.68   1200    2037.91   1691.77   1450      1200      1007.44   1200      1484.05   1669.27   1857.92   2089.13
    1679.03   1246.87   884.85    650     1319.34   1126.16   1000      884.85    758.08    749.04    749.04    1162.42   1219.77   1485.29
    1356.12   895       546.36    380     1005.13   902.03    708.77    546.36    500       500       547.23    827.86    850.2     1062.05
    1122.72   496.14    336.83    230     773.12    628.67    375.68    336.83    315.13    280       331.18    517.85    628.67    831.48
    818.6     410.13    203.33    130     562.77    371.19    271.72    203.33    190.44    203.33    267.75    357.71    515.01    766.41
    643.28    179.72    117.97    80      384.27    240.16    175.88    117.97    100       117.97    198.36    274.21    359.75    547.1
    485.26    142.38    66.95     65      267.1     139.84    123.3     66.95     50        66.95     130.51    185.53    265.34    461.41
    361.16    94.73     54.49     50      171.38    116.87    90        54.49     40        54.49     100.47    129.69    172.36    347.28
    290.17    93.51     49.41     45      157.51    83.8      66.69     49.41     30        49.41     93.91     88.35     149.93    244.06"
        },
        Opt1[@g, optwaysimdeep, optwayconfra]   //Opt1函数全局优化 //Opt[@g, optwaysimdeep, optwayconfra] //也可以使用此语句	
    };

    结果(k0,Ea,m,n,o,误差):

    0.5911010179120556 11993.11965328726 0.5943806176877708 0.5842721259729485 1.094210900092406 2087691.254199586

    例子6 求解微分方程组的参数:

    求解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

    代码:

    !!!using["luopt","math"]; //使用命名空间
    f(t,CB,CH,dCB,dCH,pp::K1,K2,a,m)=
    {
        dCB=-K1*CB^a*(K1/K2*CB^a)^m,
        dCH=K2*(K1/K2*CB^a)^m,
        0
    };
    目标函数(_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=gsl_ode[@f,nil,0.0,tA,ra1(10.9237,0), 1e-6, 1e-6, 2, 1e-6,100],
        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{ 
    //存放实验数据
    "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取矩阵的列
        Opt1[@目标函数,optwaysimdeep, optwayconfra] //Opt1函数全局优化
        //Opt[@目标函数, optwaysimdeep, optwayconfra] //也可以使用此语句
    };

    结果(K1,K2,a,m,目标函数值):

    2.007221403771439e-002 3.477931047132456e-002 0.7598653691173753 -1.205432420176749 18.90800296978626

    绘图:

    !!!using("win","math");
    f(t,CB,CH,dCB,dCH,pp :: K1,K2,a,m)=
    {
        dCB=-K1*CB^a*(K1/K2*CB^a)^m,
        dCH=K2*(K1/K2*CB^a)^m,
        0
    };
    init(x,tx::K1,K2,a,m)=
        x=matrix[
    "
    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
    "
        ],
        tx=x(all:0),
        K1=2.007221403771439e-002, K2=3.477931047132456e-002, a=0.7598653691173753, m=-1.205432420176749,
        luShareX2(x, gsl_ode[@f,nil,0.0,tx,ra1(10.9237,0), 1e-6, 1e-6, 2, 1e-6,100]);
    ChartWnd[@init];
       //显示窗口并初始化

    分析:本例所求解的应是最优解,但图形显示数据与曲线不一致,故怀疑数据与模型不匹配。

    例子7 常微分方程组的参数拟合:

    根据实验数据拟合出细胞对水的渗透参数lp和对保护剂的渗透系数ps。
    微分方程组是这样的:
    da/dtt=lp*r*temp*(0.289*(v0-vb)/a-0.289+b/a-4.809)*4.83579*c^0.6667
    db/dt=ps*(4.809-b/a)*4.83579*c^0.6667
    dc/dt=da/dt+vcpa*db/dt
    其中,a,b,c是因变量,t是自变量,其他的r,temp,v0,vb,vcpa均为常量,如下:r=0.08206,temp=298,v0=0.9402*10^(-12),vb=0.2642*10^(-12),vcpa=63.75*10^(-6)
    实验只能测出a和c值,b值无法测出。

    实验数据:

    t  a      c
    0  0.719  1
    21 0.2503 0.5314
    29 0.1412 0.4222
    38 0.1893 0.4703
    46 0.2268 0.5078
    54 0.2552 0.5362
    71 0.2818 0.5628
    87 0.3122 0.5932
    104 0.3434 0.6244
    120 0.3724 0.6534
    137 0.3957 0.6767
    153 0.4068 0.6878
    170 0.428 0.709
    195 0.4514 0.7324
    228 0.4744 0.7554
    261 0.4977 0.7787
    294 0.5154 0.7964
    335 0.5344 0.8154
    377 0.5424 0.8234
    426 0.5709 0.8519
    476 0.587 0.868
    526 0.6022 0.8832

    分析:因b值无法测出,故追加b的初值_b为优化参数。下面代码中本来没有[tyz(i,0)-tac(i,0)]^2,但发现所求解存在异常,故添加此代码。

    代码:

    !!!using["luopt","math"]; //使用命名空间
    初值(::r,temp,v0,vb,vcpa)= r=0.08206,temp=298,v0=0.9402*10^(-12),vb=0.2642*10^(-12),vcpa=63.75*10^(-6);
    f(t,a,b,c,da,db,dc,pp :: r,temp,v0,vb,vcpa,lp,ps)=
    {
        da=lp*r*temp*(0.289*(v0-vb)/a-0.289+b/a-4.809)*4.83579*c^0.6667,
        db=ps*(4.809-b/a)*4.83579*c^0.6667,
        dc=da+vcpa*db,
        0
    };
    目标函数(_lp,_ps,_b : i,s,tyz : tac,t,max,lp,ps)=
    {
        lp=_lp,ps=_ps, 
    //传递优化变量,函数f中要用到lp,ps
        tyz=gsl_ode[@f,nil,0.0,t,ra1(0.719,_b, 1), 1e-6, 1e-6, 2, 1e-6,100],
        i=0, s=0, while{++i<max, s=s+[tyz(i,0)-tac(i,0)]^2+[tyz(i,1)-tac(i,1)]^2+[tyz(i,3)-tac(i,2)]^2},
        s
    };
    main(::tac,t,max)=
    {
        tac=matrix{ 
    //存放实验数据
    "0 0.719 1
    21 0.2503 0.5314
    29 0.1412 0.4222
    38 0.1893 0.4703
    46 0.2268 0.5078
    54 0.2552 0.5362
    71 0.2818 0.5628
    87 0.3122 0.5932
    104 0.3434 0.6244
    120 0.3724 0.6534
    137 0.3957 0.6767
    153 0.4068 0.6878
    170 0.428 0.709
    195 0.4514 0.7324
    228 0.4744 0.7554
    261 0.4977 0.7787
    294 0.5154 0.7964
    335 0.5344 0.8154
    377 0.5424 0.8234
    426 0.5709 0.8519
    476 0.587 0.868
    526 0.6022 0.8832"
    },
        len[tac,0,&max], t=tac(all:0), 
    //用len函数取矩阵的行数,tA取矩阵的列
        Opt1[@目标函数] //Opt1函数全局优化
        //Opt[@目标函数] //也可以使用此语句
    };

    结果(lp,ps,_b,目标函数值):

    2.473962346617309e-003 -3.626422822558594e-003 1.308637154884405 0.1199412216587657

    绘图:

    !!!using("win","math");
    初值(::r,temp,v0,vb,vcpa)= r=0.08206,temp=298,v0=0.9402*10^(-12),vb=0.2642*10^(-12),vcpa=63.75*10^(-6);
    f(t,a,b,c,da,db,dc,pp :: r,temp,v0,vb,vcpa,lp,ps)=
    {
        da=lp*r*temp*(0.289*(v0-vb)/a-0.289+b/a-4.809)*4.83579*c^0.6667,
        db=ps*(4.809-b/a)*4.83579*c^0.6667,
        dc=da+vcpa*db,
        0
    };
    init(x,y,tx,max::lp,ps)=
        x=matrix[
    "
    0 0.719 1
    21 0.2503 0.5314
    29 0.1412 0.4222
    38 0.1893 0.4703
    46 0.2268 0.5078
    54 0.2552 0.5362
    71 0.2818 0.5628
    87 0.3122 0.5932
    104 0.3434 0.6244
    120 0.3724 0.6534
    137 0.3957 0.6767
    153 0.4068 0.6878
    170 0.428 0.709
    195 0.4514 0.7324
    228 0.4744 0.7554
    261 0.4977 0.7787
    294 0.5154 0.7964
    335 0.5344 0.8154
    377 0.5424 0.8234
    426 0.5709 0.8519
    476 0.587 0.868
    526 0.6022 0.8832
    "
        ],
        cwAttach[typeSplit], cwResizePlots(2,2,2),
        tx=x(all:0).reshape(), max=len[tx],
        cwAddCurve{tx, x(all:1).reshape(), max, 0},
        cwAddCurve{tx, x(all:2).reshape(), max, 1},
        lp=2.473962346617309e-003,ps=-3.626422822558594e-003,
        y=gsl_ode[@f,nil,0.0,tx,ra1(0.719,1.308637154884405, 1), 1e-6, 1e-6, 2, 1e-6,100],
        cwAddCurve{tx, subg(y,all:1).reshape(), max, 0},
        cwAddCurve{tx, subg(y,all:2).reshape(), max, 2},
        cwAddCurve{tx, subg(y,all:3).reshape(), max, 1};
    ChartWnd[@init]; 
     //显示窗口并初始化

    例子8 常微分方程组的参数拟合:

    dz/dt=1/tv*(((lp*t)^av)*(z^(-av-1))-((lp*t)^(-2*av))*(z^(2*av-1)))
    ds/dt=(ue*(ae-1)*(lp*t)^(ae-2)+ue*(2*ae+1)*(lp*t)^(-2*ae-2)+uv*(av-1)*(lp*t)^(av-2)*z^(-av)+uv*(2*av+1)*(lp*t)^(-2*av-2)*z^(2*av))*lp-(uv*av*(lp*t)^(av-1)*z^(-av-1)+2*uv*av*(lp*t)^(-2*av-
    1)*z^(2*av-1))*dz/dt

    ue,ae,uv,av,tv为待拟合参数
    z为中间变量
    lp=0.00354
    t-s为实验数据:
    t=282.5,361.0,439.4,517.9,596.4,674.8,753.3,831.8,910.2,988.7;
    s=0,29099,36228,40626,44290,47686,50964,54181,57363,60521;

    分析:因z值无法测出,故追加z的初值_z为优化参数。

    代码:

    !!!using["luopt","math"]; //使用命名空间
    f(t,z,s,dz,ds,pp : lp : ue,ae,uv,av,tv)=
    {
        lp=0.00354,
        dz=1/tv*(((lp*t)^av)*(z^(-av-1))-((lp*t)^(-2*av))*(z^(2*av-1))),
        ds=(ue*(ae-1)*(lp*t)^(ae-2)+ue*(2*ae+1)*(lp*t)^(-2*ae-2)+uv*(av-1)*(lp*t)^(av-2)*z^(-av)+uv*(2*av+1)*(lp*t)^(-2*av-2)*z^(2*av))*lp-(uv*av*(lp*t)^(av-1)*z^(-av-1)+2*uv*av*(lp*t)^(-2*av-1)*z^(2*av-1))*dz,
        0
    };
    目标函数(_ue,_ae,_uv,_av,_tv,_z : i,ss,ts : t,s,max,ue,ae,uv,av,tv)=
    {
        ue=_ue,ae=_ae,uv=_uv,av=_av,tv=_tv, 
    //传递优化变量,函数f中要用到ue,ae,uv,av,tv
        ts=gsl_ode[@f,nil,0.0,t,ra1(_z,0), 1e-6, 1e-6, 2, 1e-6,100],
        i=0, ss=0, while{++i<max, ss=ss+[ts(i,2)-s(i)]^2},
        sqrt[ss/max]
    };
    main(::t,s,max)=
    {
        t=ra1[282.5,361.0,439.4,517.9,596.4,674.8,753.3,831.8,910.2,988.7],
        s=ra1[0,29099,36228,40626,44290,47686,50964,54181,57363,60521],
        max=len[t],
        Opt1[@目标函数] 
    //Opt1函数全局优化
        //Opt[@目标函数] //也可以使用此语句
    };

    结果(ue,ae,uv,av,tv,_z,目标函数值):

    -9114964.001609007 -0.4453565389059405 7335590868.341645 -2.754919559106283e-004 -10485633.9903882 2.654977892066841 244.5944885704022

    -249.2203422970038 -2.000969463745669 1822318656671169 -5.117661642056397e-013 -30625.54305602086 2.426847712155262e-010 151.3043733914658(此组解不能确定是否有效)

    绘图:

    !!!using("win","math");
    f(t,z,s,dz,ds,pp : lp : ue,ae,uv,av,tv)=
    {
        lp=0.00354,
        dz=1/tv*(((lp*t)^av)*(z^(-av-1))-((lp*t)^(-2*av))*(z^(2*av-1))),
        ds=(ue*(ae-1)*(lp*t)^(ae-2)+ue*(2*ae+1)*(lp*t)^(-2*ae-2)+uv*(av-1)*(lp*t)^(av-2)*z^(-av)+uv*(2*av+1)*(lp*t)^(-2*av-2)*z^(2*av))*lp-(uv*av*(lp*t)^(av-1)*z^(-av-1)+2*uv*av*(lp*t)^(-2*av-1)*z^(2*av-1))*dz,
        0
    };
    init(x,t,tx,max,y::ue,ae,uv,av,tv)=
        x = matrix[
    "
    282.5,361.0,439.4,517.9,596.4,674.8,753.3,831.8,910.2,988.7
    0,29099,36228,40626,44290,47686,50964,54181,57363,60521
    "
        ],
        t=ra1[282.5,361.0,439.4,517.9,596.4,674.8,753.3,831.8,910.2,988.7],
        ue=-9114964.001609007, ae=-0.4453565389059405, uv=7335590868.341645, av=-2.754919559106283e-004, tv=-10485633.9903882,
        luShareX2(x', gsl_ode[@f,nil,0.0,t,ra1(2.654977892066841,0), 1e-6, 1e-6, 2, 1e-6,100]), outa[gsl_ode[@f,nil,0.0,t,ra1(2.654977892066841,0), 1e-6, 1e-6, 2, 1e-6,100]];
    ChartWnd[@init];
     //显示窗口并初始化


    版权所有© Lu程序设计 2002-2021,保留所有权利
    E-mail: forcal@sina.com
      QQ:630715621
    最近更新: 2021年07月29日

    展开全文
    forcal 2021-07-30 19:51:44
  • 1stOpt计算:Elapsed Time (Hr:Min:Sec:Msec): 00:00:19:284Root of Mean Square Error (RMSE): 1.18378166169196Sum of Square Error: 7.00669511279087Correlation Coef. (R): 0.999185213441573R-Square: 0....

    1stOpt计算:

    Elapsed Time (Hr:Min:Sec:Msec): 00:00:19:284

    Root of Mean Square Error (RMSE): 1.18378166169196

    Sum of Square Error: 7.00669511279087

    Correlation Coef. (R): 0.999185213441573

    R-Square: 0.998371090760282

    Adjusted R-Square: 0.996742181520565

    Determination Coef. (DC): 0.998162335813469

    F-Statistic: 557.971518394492

    Parameter                  Best Estimate

    --------------------        -------------

    r        -6.91097295660219E-14

    xm        13056.0823655475

    a        -48.8514044361061

    ====== Output Results ======

    File: Data File-1

    No        t        Target x        Calculated x

    1        2011        7898.8        7896.87227419657

    2        2012        7919.98        7920.46086793675

    3        2013        7939.49        7941.10309493816

    4        2014        7960.06        7959.45928858892

    5        2015        7976.3        7975.98950863453,

    展开全文
    weixin_31881743 2021-04-20 04:06:05
  • weixin_35344136 2021-04-18 14:38:07
  • weixin_42303078 2021-04-18 14:38:10
  • shengzimao 2020-10-03 21:34:44
  • chimchim12138 2020-06-24 10:15:50
  • weixin_30616309 2021-04-18 14:37:41
  • weixin_39622217 2021-04-18 14:38:15
  • weixin_31889919 2021-04-18 14:38:23
  • 6.22MB qq_41934573 2021-08-26 22:42:44
  • weixin_31032795 2021-04-18 08:58:08
  • wlfc 2021-08-11 15:15:10
  • weixin_36337403 2021-04-21 04:28:54
  • weixin_42309722 2021-04-20 06:53:28
  • wlfc 2021-08-20 20:06:16
  • weixin_33056525 2021-04-21 04:29:52
  • weixin_35319894 2021-04-18 08:58:04
  • weixin_35695434 2020-12-23 12:25:09
  • weixin_35141129 2021-04-22 11:37:30
  • weixin_28837357 2021-05-08 15:47:22
  • wlfc 2021-08-10 10:52:42
  • weixin_33277564 2021-04-23 20:30:15
  • weixin_34791959 2021-04-21 15:37:00
  • weixin_33013787 2021-04-21 15:36:58
  • wlfc 2021-08-09 15:25:10
  • weixin_39851809 2021-04-20 08:45:41
  • weixin_38425472 2020-06-19 12:28:29

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 7,445
精华内容 2,978
关键字:

复杂微分方程参数拟合