Lu 微分方程参数优化（拟合）

欢迎访问 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