欢迎访问 Lu程序设计

Lu 解方程、优化(拟合)例子

例子1 解方程组:

    其中:θv=0.05, qv=0.05, Da=0.3, δ=0.5, Z=10, Le=2.0, Q=2.0, σ=0.15, Rf=1.0。自变量为U,Rv,θf。

    分析:该方程组较难求解,故使用优化方法;由于包含了大量积分运算,故耗时较长。

    Lu代码:

!!!using["math","luopt"];
常量定义(::θv,Q,qv,Da,δ,Le,σ,Z,Rf)= θv=0.05, Q=2.0, qv=0.05, Da=0.3, δ=0.5, Le=2.0, σ=0.15, Z=10.0, Rf=1.0;
I(a:x:Da,U,Rv)= x=Da/U+U, exp[Da/U*(a-Rv)+a*U]*[a*a/x-2*a/x^2+2/x^3];
I1(s::U)= s^(-2)*exp(-U*s);
I2(s::U,Rv)= s^(-2)*exp(-U*s)*[I(s)-I(Rv)];
I3(s::U,Le)= s^(-2)*exp(-U*Le*s);
F1(UU,RRv,θθf::U,Rv,θf,θv,Rf,qv,Da,δ)= U=UU,Rv=RRv,θf=θθf, θv*gsl_qagiu(@I1,Rf)/ gsl_qagiu(@I1,Rv)-{qv*Da*δ* gsl_qag(@I2,Rf,Rv)};
F2(UU,RRv,θθf:x,y:U,Rv,θf,θv,Rf,Q,qv,Da,δ,Le)= U=UU,Rv=RRv,θf=θθf, x=Rf^(-2)*exp(-U*Rf),y= gsl_qagiu(@I1,Rf),
  θf*x / y-Q*x+qv*Da*δ*x*{I(Rv)-I(0)+ gsl_qag(@I2,Rf,Rv)/y} - 1/Le* Rf^(-2)*exp(-U*Rf*Le)/ gsl_qagiu(@I3,Rf);
F3(UU,RRv,θθf::U,Rv,θf,θv,Rf,Le,σ,Z)= U=UU,Rv=RRv,θf=θθf, 1/Le* Rf^(-2)*exp(-U*Rf*Le)/ gsl_qagiu(@I3,Rf) - {[σ+(1-σ)*θf]^2*exp{Z/2*(θf-1)/[σ+(1-σ)*θf]}};
F(UU,RRv,θθf)= F1(UU,RRv,θθf)^2+F2(UU,RRv,θθf)^2+F3(UU,RRv,θθf)^2;
Opt1[@F]   //Opt1函数全局优化

    结果(U,Rv,θf,误差):

0.8643894660384006  -9.986043390866924e-007  1.079488506612208  9.920352296080973e-031

例子2 数学模型:

    其中x=240;t=0的时候,f1=0,f2=71。拟合参数a1、a2和a3。实验数据(t,f1)在文件data1.txt中:

0.       0.
3E-4     0
8E-4     0
0.0011   0
0.002    9.09091E-4
0.0033   9.09091E-4
0.0058   9.09091E-4
0.0106   9.09091E-4
0.0197   9.09091E-4
0.0378   9.09091E-4
0.0736   0.00182
0.1447   0.00364
0.287    0.00545
0.5708   0.00818
1.1395   0.01182
1.1417   0.01182
1.1431   0.01182
1.2108   0.01182
1.2147   0.01182
1.2467   0.01182
1.2517   0.01182
1.2589   0.01182
1.2597   0.01182
1.4667   0.01273
1.4672   0.01273
1.8067   0.01364
1.807    0.01364
1.8089   0.01364
1.8097   0.01364
1.8106   0.01364
1.8128   0.01364
1.8164   0.01364
1.8195   0.01364
1.8214   0.01364
1.8233   0.01364
1.8256   0.01364
1.8311   0.01364
1.832    0.01364
1.8322   0.01364
1.837    0.01364
1.837    0.01364
2.7097   0.01545
2.71     0.01545
2.7125   0.01545
2.7136   0.01545
2.7442   0.01545
2.7445   0.01545
2.7453   0.01545
2.7461   0.01545
2.8128   1.636e-002
3.2728   1.636e-002
3.392    1.727e-002
3.7397   1.727e-002
3.8586   1.727e-002
4.0272   1.727e-002
4.4492   1.818e-002
4.545    1.818e-002
4.7931   1.818e-002
5.0847   1.909e-002
5.2003   1.909e-002
5.4625   1.909e-002
5.6206   1.909e-002
6.1597   2.e-002
6.2464   2.e-002
6.4933   2.e-002
6.7872   2.e-002
6.8247   2.e-002
6.9442   2.e-002
7.015    2.e-002
7.1361   2.091e-002
7.2964   2.091e-002
7.3925   2.091e-002
7.4745   2.091e-002
7.5781   2.091e-002
7.7561   2.091e-002
7.8542   2.091e-002
8.0386   2.091e-002
8.2875   2.091e-002
8.355    2.091e-002
8.3939   2.091e-002
8.4356   2.091e-002
8.567    2.091e-002
8.6008   2.182e-002
9.4197   2.182e-002
9.4689   2.182e-002
10.3439  2.273e-002
10.6614  2.273e-002
10.7181  2.273e-002
10.7706  2.273e-002
10.9411  2.273e-002
11.0047  2.273e-002
11.3175  2.273e-002
11.4417  2.273e-002
11.477   2.273e-002
12.5083  0.02364
13.5047  0.02455
14.5117  2.455e-002
14.5506  2.455e-002
14.6125  2.455e-002
14.785   2.455e-002
15.4081  2.545e-002
15.5286  2.545e-002
15.9611  2.545e-002
16.0197  2.545e-002
16.0678  2.545e-002
16.127   2.545e-002
16.1697  2.545e-002
16.2211  2.545e-002

    Lu代码1:

!!!using["luopt","math","win","sys"];
f(t,f1,f2,df1,df2,p : x : a1,a2,a3)=
{
    x=240.0,
    df1=sinh(a1*x*(1-f2)),
    df2=a2*df1*(1-f2/a3)/x,
    0
};
J(aa1,aa2,aa3 : tf,s : max,tA,f1A,a1,a2,a3)=
{
    a1=aa1, a2=aa2, a3=aa3,
    tf=gsl_ode[@f,nil,0.0,tA,ra1(0,71), 1e-6, 1e-6, gsl_rk4, 1e-6,100],
    s=sum{[tf(all:1)-f1A].^2.0},            //代码矢量化加速
    sqrt[s/max]                             //均方差作为目标函数值
};
main(::tArray,max,tA,f1A)=
{
    //因数据较多,故将数据保存在文件data1.txt中,程序运行时用函数GetFileName选择该文件,并用函数readtxt读取该文件。
    tArray=matrix[GetFileName().readtxt()],  //读取文件data1.txt中的数据并转换为一个矩阵
    len[tArray,0,&max], tA=tArray(all:0), f1A=tArray(all:1),  //用len函数取矩阵的行数,tA和f1A取矩阵的列
    Opt1[@J]                                //Opt1函数全局优化
};

    结果(前面的数为最优参数,最后一个数是极小值。下同):

-2.459730744985535e-006  67875.73531627061  1.346955517639491  2.783775439907169e-004

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

!!!using["luopt","math","win","sys"];
f(t,f1,f2,df1,df2,p : x : a1,a2,a3)=
{
    x=240.0,
    df1=sinh(a1*x*(1-f2)),
    df2=a2*df1*(1-f2/a3)/x,
    0
};
J(aa1,aa2,aa3 : tf,s : max,tA,f1A,a1,a2,a3)=
{
    a1=aa1, a2=aa2, a3=aa3,
    tf=gsl_ode[@f,nil,0.0,tA,ra1(0,71), 1e-6, 1e-6, gsl_rk4, 1e-6,100],
    s=sum{[tf(all:1)-f1A].^2.0},            //代码矢量化加速
    sqrt[s/max]                             //均方差作为目标函数值
};
init(main:tf:tArray,max,tA,f1A,a1,a2,a3)=
{
    //因数据较多,故将数据保存在文件data1.txt中,程序运行时用函数GetFileName选择该文件,并用函数readtxt读取该文件。
    tArray=matrix[GetFileName().readtxt()],  //读取文件data1.txt中的数据并转换为一个矩阵
    len[tArray,0,&max], tA=tArray(all:0), f1A=tArray(all:1),    //用len函数取矩阵的行数,tA和f1A取矩阵的列
    a1=0.0, a2=0.0, a3=0.0,
    Opt1[@J, optstarter,&a1,&a2,&a3,0],       //Opt1函数全局优化
    tf=gsl_ode[@f,nil,0.0,tA,ra1(0,71)],       //获得最终结果
    outa[tf],                                //输出结果
    luShareX2[tArray, tf]                    //绘制共享X轴视图
};
ChartWnd[@init];                             //显示窗口并初始化

    结果:

例子3 微分方程过拟合及处理

微分方程式:

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脚本代码:

!!!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_rk4, 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_rk4, 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, gsl_rk4, 1e-6,100],
//获得最终结果
    outa[tyz],
//outa输出结果
    luShareX2[tyArray, tyz]
//绘制共享X轴视图
};
ChartWnd[@init];
//显示窗口并初始化

图形:

例子4 缺少部分参数的微分方程拟合-本例耗时较长

微分方程组如下:

dx/dt=a*x-b*x*y
dy/dt=-c*y+d*x*y

数据如下:
t  x(t)  y(t)
11 45.79 41.40
12 53.03 38.90
13 64.05 36.78
14 75.40 36.04
15 90.36 33.78
16 107.14 35.40
17 127.79 34.68
18 150.77 36.61
19 179.65 37.71
20 211.82 41.98
21 249.91 45.72
22 291.31 53.10
23 334.95 65.44
24 380.67 83.00
25 420.28 108.74
26 445.56 150.01
27 447.63 205.61
28 414.04 281.60
29 347.04 364.56
30 265.33 440.30
31 187.57 489.68
32 128.00 512.95
33 85.25 510.01
34 57.17 491.06
35 39.96 462.22
36 29.22 430.15
37 22.30 396.95
38 16.52 364.87
39 14.41 333.16
40 11.58 304.97
41 10.41 277.73
42 10.17 253.16
43 7.86 229.66
44 9.23 209.53
45 8.22 190.07
46 8.76 173.58
47 7.90 156.40
48 8.38 143.05
49 9.53 130.75
50 9.33 117.49
51 9.72 108.16
52 10.55 98.08
53 13.05 88.91
54 13.58 82.28
55 16.31 75.42
56 17.75 69.58
57 20.11 62.58
58 23.98 59.22
59 28.51 54.91
60 31.61 49.79
61 37.13 45.94
62 45.06 43.41
63 53.40 41.30
64 62.39 40.28
65 72.89 37.71
66 86.92 36.58
67 103.32 36.98
68 121.70 36.65
69 144.86 37.87
70 171.92 39.63
71 202.51 42.97
72 237.69 46.95
73 276.77 54.93
74 319.76 64.61
75 362.05 81.28
76 400.11 105.50
77 427.79 143.03
78 434.56 192.45
79 410.31 260.84
80 354.18 339.39
81 278.49 413.79
82 203.72 466.94
83 141.06 494.72
84 95.08 499.37
85 66.76 484.58
86 45.41 460.63
87 33.13 429.79
88 25.89 398.77
89 20.51 366.49
90 17.11 336.56
91 12.69 306.39
92 11.76 279.53
93 11.22 254.95
94 10.29 233.50
95 8.82 212.74
96 9.51 193.61
97 8.69 175.01
98 9.53 160.59
99 8.68 146.12
100 10.82 131.85


求参数a,b,c,d以及x,y的初值。

分析:本例中x和y的初值未知,追加为拟合参数,故拟合参数有6个:a,b,c,d,x0,y0。另外,在数据第一行补充 0  0  0,与t=0时初值相对应,但拟合时初值 0  0 被忽略。

Lu代码:

!!!using["luopt","math"]; //使用命名空间
f(t,x,y,dx,dy, params :: a,b,c,d)=
{
    dx = a*x-b*x*y,
    dy = -c*y+d*x*y,
    0
//必须返回0
};
目标函数(_a,_b,_c,_d,x0,y0 : i,s,ty : tyArray,tA,max, a,b,c,d)=
{
    a=_a, b=_b, c=_c, d=_d,
//传递优化变量
    //最后一个参数50表示gsl_ode函数在计算时,最多循环计算50次,这样可以提高速度

    ty=gsl_ode[@f,nil,0.0,tA,ra1(x0, y0), 1e-6, 1e-6, gsl_rk4, 1e-6,50],
    i=0, s=0, while{++i<max, s=s+[ty(i,1)-tyArray(i,1)]^2+[ty(i,2)-tyArray(i,2)]^2},
    s
};
main(::tyArray,tA,max)=
{
    tyArray=matrix{
//存放实验数据ti,xi,yi
        "0 0 0
 //补充t=0时的数据,拟合时后面的两个0被忽略
        11 45.79 41.40
        ... ...省略数据
        100 10.82 131.85"
    },
    len[tyArray,0,&max], tA=tyArray(all:0),
//用len函数取矩阵的行数,tA取矩阵的列
    Opt1[@目标函数, optwaygrow]
//Opt1函数全局优化
};

结果:

0.2146380490718624 1.207150888688551e-003 0.1032935041696845 9.484393405386229e-004 10.66928857150351 104.905115896753 1283.613374197572

绘图代码:

!!!using["luopt","math","win"]; //使用命名空间
f(t,x,y,dx,dy, params :: a,b,c,d)=
{
    dx = a*x-b*x*y,
    dy = -c*y+d*x*y,
    0
//必须返回0
};
init(main:ty:tyArray,tA,max, a,b,c,d)=
{
    tyArray=matrix{
//存放实验数据ti,xi,yi
        "0 0 0
        11 45.79 41.40
        ... ...省略数据
        100 10.82 131.85"
    },
    len[tyArray,0,&max], tA=tyArray(all:0),
//用len函数取矩阵的行数,tA取矩阵的列
    a=0.2146380490718624, b=1.207150888688551e-003, c=0.1032935041696845, d=9.484393405386229e-004,
    ty=gsl_ode[@f,nil,0.0,tA,ra1(10.66928857150351, 104.905115896753), 1e-6, 1e-6, gsl_rk4, 1e-6,50],
//获得最终结果
    outa[ty],
//outa输出结果
    tyArray[0,1]=ty[0,1], tyArray[0,2]=ty[0,2],
//补充数据
    luShareX2[tyArray, ty]
//绘制共享X轴视图
};
ChartWnd[@init];
//显示窗口并初始化

图形:

例子5 缺少部分参数的微分方程拟合-本例耗时较长

已知4个微分方程如下:

dS/dt = − αSI − pS
dE/dt = αSI − βE − εE + qR
dI/dt = βE − θI
dR/dt = pS +εE+ θI − qR

已知的参数是:

t=[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15]
I=[11714 21352 21789 22403 16308 17017 17146 33229 60907 38875 33437 47282 51749 40530 31617]

怎么才可以求出方程中的其他数据:S,E,R ,α,p,β,ε,q,θ这些参数的值?

分析:S,E,R缺少初值数据,故作为拟合参数与α,p,β,ε,q,θ一起进行拟合。

Lu代码:

!!!using["luopt","math"]; //使用命名空间
f(t,S,E,I,R,dS,dE,dI,dR, params :: α,p,β,ε,q,θ)=
{
    dS = -α*S*I - p*S,
    dE = α*S*I - β*E - ε*E + q*R,
    dI = β*E - θ*I,
    dR = p*S +ε*E+ θ*I - q*R,
    0
//必须返回0
};
目标函数(_α,_p,_β,_ε,_q,_θ,S,E,R : i,s,ty : tA,IA,max, α,p,β,ε,q,θ)=
{
    α=_α, p=_p, β=_β, ε=_ε, q=_q, θ=_θ,
//传递优化变量
    //最后一个参数50表示gsl_ode函数在计算时,最多循环计算50次,这样可以提高速度

    ty=gsl_ode[@f,nil,0.0,tA,ra1(S, E, 11714, R), 1e-6, 1e-6, gsl_rk4, 1e-6,50],
    i=0, s=0, while{++i<max, s=s+[ty(i,3)-IA(i)]^2},
    sqrt[s/(max-1)]
};
main(::tA,IA,max)=
{
    tA=ra1[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15],
    IA=ra1[11714, 21352, 21789, 22403, 16308, 17017, 17146, 33229, 60907, 38875, 33437, 47282, 51749, 40530, 31617],
    max=len[tA],
    Opt1[@目标函数, optwaycom, optwaygrow]
//Opt1函数全局优化
   
//Opt1[@目标函数, optwaycom, optwaynewton] //Opt1函数全局优化
   
//Opt1[@目标函数, optwaycom, optwaylm] //Opt1函数全局优化
};

结果(未多次运行,可能还有更优解):

2.482748838436231e-004 -10.85224138992894 2.202455853536803e-005 -11.93998584165819 13.2304263331211 0.9250513971698489 -6.614778299707771e-010 1442590253.096928 -1370808306.615013 2359.28301212148

1.769857225652821e-004 -7.627690405609265 0.6030462597114438 -6.852326341185071 9.282758472827229 -1.517826328917237 -2.166343265010053e-008 2589.497835198801 -11201.81806513364 2136.246710504934

绘图代码:

!!!using["luopt","math","win"]; //使用命名空间
f(t,S,E,I,R,dS,dE,dI,dR, params :: α,p,β,ε,q,θ)=
{
    dS = -α*S*I - p*S,
    dE = α*S*I - β*E - ε*E + q*R,
    dI = β*E - θ*I,
    dR = p*S +ε*E+ θ*I - q*R,
    0
//必须返回0
};
init(main:k:tA,IA,max,ty,α,p,β,ε,q,θ)=
{
    tA=ra1[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15],
    IA=ra1[11714, 21352, 21789, 22403, 16308, 17017, 17146, 33229, 60907, 38875, 33437, 47282, 51749, 40530, 31617],
    max=len[tA],
    α=2.482748838436231e-004, p=-10.85224138992894, β=2.202455853536803e-005, ε=-11.93998584165819, q=13.2304263331211, θ=0.92505139716984892,
    ty=gsl_ode[@f,nil,0.0,tA,ra1(-6.614778299707771e-010, 1442590253.096928, 11714, -1370808306.615013), 1e-6, 1e-6, gsl_rk4, 1e-6,50],
    o[ty],
    cwAttach[typeSplit], cwResizePlots(1,2,2),
//左二右二分裂
    k=cwAddCurve{tA, IA, max, 0},
//给0子图添加曲线
    cwSetScatter(k,0),       
//设置绘制点
    cwSetDataLineSize(5, k, 0),
  //设置点的大小
    cwAddCurve{tA, ty[all:3].reshape(), max, 0},
//给0子图添加曲线
    cwAddCurve{tA, ty[all:1].reshape(), max, 1},
//给1子图添加曲线
    cwAddCurve{tA, ty[all:2].reshape(), max, 2},
//给2子图添加曲线
    cwAddCurve{tA, ty[all:4].reshape(), max, 3}
//给3子图添加曲线
};
ChartWnd[@init];

图形:

例子6 复数拟合:利用cole-cole模型拟合介电常数

    e=(a+(b-a)/(1+(2*pi*f*c)^2))-1.0i*(d/(2*pi*f*(1/(36*pi)*1e-9))+((b-a)*2*pi*f*c)/(1+(2*pi*f*c)^2))

拟合参数:a,b,c,d

数据:复数 e = e1 + e2 * i (i为虚数单位)

f     e1[realPart] e2[imagPart]
20000000 78.6416 0.0101
20718771.93 78.5189 0.0451
21437543.86 78.5453 0.0435
22156315.78 78.587 0.023
22875087.71 78.6542 0.1911
23593859.64 78.6054 0.1855
24312631.57 78.6201 0.1524
25031403.49 78.5626 0.089
25750175.42 78.5386 0.068
26468947.35 78.5829 0.1993
27187719.28 78.6149 0.1347
27906491.2 78.6495 0.2324
28625263.13 78.6226 0.1863
29344035.06 78.5616 0.1389
30062806.99 78.5146 0.0626
30781578.91 78.5795 0.0778
31500350.84 78.5842 0.15
32219122.77 78.5785 0.1944
32937894.7 78.6011 0.1988
33656666.62 78.5833 0.2089
34375438.55 78.5298 0.1492
35610843.56 78.5337 0.1835
36846248.57 78.5676 0.2353
38081653.59 78.5411 0.2025
39317058.6 78.5219 0.209
40552463.61 78.5208 0.2245
41787868.62 78.5198 0.2761
43023273.63 78.5312 0.2217
44258678.64 78.5029 0.2126
45494083.66 78.5344 0.2328
46729488.67 78.5319 0.224
47964893.68 78.5077 0.2183
49200298.69 78.5345 0.2317
50435703.7 78.4988 0.2495
51671108.71 78.5091 0.2274
52906513.72 78.4915 0.2267
54141918.74 78.4901 0.249
55377323.75 78.5024 0.2342
56612728.76 78.4836 0.2522
57848133.77 78.4873 0.2532
59083538.78 78.4978 0.2615
61206918.23 78.4948 0.271
63330297.69 78.4762 0.2821
65453677.14 78.479 0.28
67577056.59 78.4771 0.2963
69700436.05 78.4804 0.2938
71823815.5 78.4822 0.3015
73947194.95 78.486 0.3105
76070574.4 78.482 0.319
78193953.86 78.4746 0.3237
80317333.31 78.4806 0.3389
82440712.76 78.4747 0.3411
84564092.22 78.4703 0.3535
86687471.67 78.4737 0.3644
88810851.12 78.4658 0.3889
90934230.58 78.4736 0.3944
93057610.03 78.482 0.4014
95180989.48 78.4542 0.4141
97304368.93 78.4759 0.416
99427748.39 78.4894 0.4195
101551127.8 78.491 0.4132
105200732.8 78.4802 0.4335
108850337.8 78.4632 0.4587
112499942.8 78.4881 0.4779
116149547.8 78.4726 0.5065
119799152.8 78.4908 0.5164
123448757.8 78.4856 0.5224
127098362.8 78.492 0.5305
130747967.8 78.4948 0.5082
134397572.8 78.5064 0.5592
138047177.8 78.4726 0.5863
141696782.8 78.478 0.5951
145346387.8 78.4849 0.5928
148995992.8 78.4968 0.6001
152645597.8 78.4783 0.5992
156295202.8 78.4736 0.6006
159944807.8 78.4737 0.6351
163594412.8 78.4618 0.652
167244017.8 78.4684 0.6446
170893622.8 78.4699 0.6806
174543227.7 78.4661 0.686
180816066.4 78.479 0.7161
187088905 78.4725 0.7234
193361743.6 78.4698 0.7509
199634582.2 78.4618 0.7711
205907420.8 78.4567 0.7959
212180259.4 78.457 0.8197
218453098 78.4535 0.8434
224725936.6 78.4542 0.8728
230998775.3 78.4548 0.8929
237271613.9 78.4428 0.9265
243544452.5 78.4528 0.9455
249817291.1 78.4494 0.9737
256090129.7 78.4509 0.9944
262362968.3 78.451 1.0089
268635806.9 78.4428 1.0312
274908645.5 78.446 1.0494
281181484.2 78.4415 1.0718
287454322.8 78.4424 1.101
293727161.4 78.4476 1.123
300000000 78.4461 1.1537
310781578.9 78.4525 1.1918
321563157.8 78.4471 1.2336
332344736.7 78.4527 1.2755
343126315.7 78.4449 1.3167
353907894.6 78.4407 1.3568
364689473.5 78.4349 1.4022
375471052.4 78.4317 1.4425
386252631.3 78.4325 1.4809
397034210.2 78.4263 1.5171
407815789.1 78.4243 1.5621
418597368.1 78.4242 1.6133
429378947 78.4334 1.6376
440160525.9 78.4304 1.6946
450942104.8 78.4303 1.7297
461723683.7 78.4214 1.7644
472505262.6 78.4189 1.8087
483286841.5 78.4175 1.8453
494068420.4 78.4149 1.8888
504849999.4 78.4115 1.9302
515631578.3 78.4059 1.9722
534162653.4 78.3997 2.0431
552693728.6 78.3949 2.1094
571224803.8 78.3907 2.1788
589755879 78.391 2.2449
608286954.1 78.3883 2.3201
626818029.3 78.3852 2.3927
645349104.5 78.38 2.4643
663880179.7 78.3711 2.5338
682411254.8 78.3642 2.6052
700942330 78.3575 2.6762
719473405.2 78.3524 2.7449
738004480.3 78.3515 2.8165
756535555.5 78.3502 2.8895
775066630.7 78.3433 2.9579
793597705.9 78.3361 3.0293
812128781 78.3288 3.0984
830659856.2 78.3219 3.1682
849190931.4 78.3161 3.2349
867722006.5 78.3078 3.304
886253081.7 78.3039 3.3728
918103773.5 78.2928 3.4979
949954465.3 78.2817 3.625
981805157.1 78.2674 3.7441
1013655849 78.2551 3.8599
1045506541 78.2436 3.9782
1077357232 78.2297 4.1022
1109207924 78.2163 4.229
1141058616 78.2042 4.3487
1172909308 78.1894 4.4612
1204760000 78.1752 4.5803
1236610691 78.1562 4.7061
1268461383 78.1413 4.832
1300312075 78.1275 4.947
1332162767 78.1113 5.0617
1364013459 78.0944 5.1861
1395864150 78.0747 5.3089
1427714842 78.0588 5.4322
1459565534 78.0435 5.5502
1491416226 78.0232 5.6666
1523266918 78.0036 5.788
1578010993 77.9694 5.9989
1632755067 77.9384 6.1999
1687499142 77.9009 6.4073
1742243217 77.8615 6.6087
1796987292 77.8266 6.8135
1851731367 77.7883 7.0213
1906475442 77.744 7.2225
1961219517 77.7061 7.4315
2015963592 77.6628 7.6281
2070707667 77.6176 7.8345
2125451742 77.5775 8.0377
2180195817 77.531 8.2393
2234939892 77.4807 8.4404
2289683967 77.4391 8.6434
2344428042 77.3859 8.841
2399172116 77.3332 9.0447
2453916191 77.2872 9.2392
2508660266 77.23 9.4387
2563404341 77.1786 9.6421
2618148416 77.1261 9.8366
2712240995 77.0262 10.1777
2806333575 76.9288 10.5137
2900426154 76.8265 10.8515
2994518733 76.7183 11.1866
3088611312 76.6111 11.5187
3182703891 76.4991 11.8533
3276796471 76.3844 12.1827
3370889050 76.2637 12.5094
3464981629 76.1481 12.8349
3559074208 76.0205 13.1668
3653166787 75.8923 13.4861
3747259366 75.7675 13.8109
3841351946 75.632 14.1324
3935444525 75.4994 14.4493
4029537104 75.3629 14.7689
4123629683 75.2199 15.0811
4217722262 75.0831 15.3966
4311814842 74.9318 15.7106
4405907421 74.7874 16.018
4500000000 74.6375 16.3285

Lu 代码:

!!!using["luopt","math","win"]; //使用命名空间

g(a,b,c,d,f)= (a+(b-a)/(1+(2*pi*f*c)^2))-1.0i*(d/(2*pi*f*(1/(36*pi)*1e-9))+((b-a)*2*pi*f*c)/(1+(2*pi*f*c)^2));

目标函数(a,b,c,d : i,s, e11,e22 : max, f, e1, e2)=
{
    i=-1, s=0, while{++i<max, toreal[g(a,b,c,d,f(i,0)), &e11,&e22], s=s+[e1(i,0)-e11]^2+[e2(i,0)-e22]^2 },
    s
};

main(: tArray : max, f, e1, e2)=
{
    tArray=matrix{
//存放实验数据 //f e1[realPart] e2[imagPart]
    "20000000 78.6416 0.0101
   
// 省略数据
    4500000000 74.6375 16.3285"
    },
    len[tArray,0,&max], f=tArray(all:0), e1=tArray(all:1),e2=tArray(all:2),
//用len函数取矩阵的行数,f等取矩阵的列
    Opt1[@目标函数]
//Opt1函数全局优化
};

结果(a,b,c,d,最小值):

5.918198522529668 78.489752592936 -8.397205875077787e-012 -6.746825318924285e-005 0.4566408877616998

绘图代码:

!!!using["luopt","math","win"]; //使用命名空间

g(a,b,c,d,f)= (a+(b-a)/(1+(2*pi*f*c)^2))-1.0i*(d/(2*pi*f*(1/(36*pi)*1e-9))+((b-a)*2*pi*f*c)/(1+(2*pi*f*c)^2));

ge(a,b,c,d : i,s,e1,e2 : max, f, e12)=
{
    i=-1, while{++i<max, toreal[g(a,b,c,d,f(i,0)), &e1, &e2], e12[i,0]=f(i,0), e12[i,1]=e1, e12[i,2]=e2}
};

init(main : tArray : max, f, e12)=
{
    tArray=matrix{
//存放实验数据 //f e1[realPart] e2[imagPart]
    "20000000 78.6416 0.0101
    20718771.93 78.5189 0.0451
    //省略数据
    4500000000 74.6375 16.3285"
    },
    len[tArray,0,&max], f=tArray(all:0),
//用len函数取矩阵的行数,f等取矩阵的列
    e12=new[real_s,max,3],
    ge[ 5.918198522529668 , 78.489752592936 , -8.397205875077787e-012 , -6.746825318924285e-005],
    luShareX2[tArray, e12]
//绘制共享X轴视图
};
ChartWnd[@init];
//显示窗口并初始化

图形:

例子7 复数拟合:用一阶德拜方程拟合介电常数(复数)

介电常数在每个频点的数据为复数,现在知道部分频点的介电常数,需要用一阶德拜方程来拟合。
一阶德拜方程: y=(a+(b-a)/(1+(2*pi*f*c)^2))-j*(d/(2*pi*f*(1/(36*pi)*10^-9))+((b-a)*2*pi*f*c)/(1+(2*pi*f*c)^2));其中f为频率,a b c d 为需要拟合的系数。

数据:频率f  介电常数y实部  介电常数y虚部

3000000 13.5 17.85
4000000 10.61 15.28
5000000 8.8 13.35
6000000 7.58 11.86
7000000 6.7 10.67
8000000 6.04 9.72
9000000 5.53 8.93
10000000 5.12 8.27
15000000 3.92 6.07
20000000 3.35 4.83
25000000 3 4.02
30000000 2.78 3.46
40000000 2.52 2.72
50000000 2.37 2.25
60000000 2.26 1.92
70000000 2.18 1.67
80000000 2.11 1.47
90000000 2.048 1.3
100000000 2.01 1.17

Lu代码:

!!!using["luopt","math","win"]; //使用命名空间

g(a,b,c,d,f)= (a+(b-a)/(1+(2*pi*f*c)^2))-1.0i*(d/(2*pi*f*(1/(36*pi)*1e-9))+((b-a)*2*pi*f*c)/(1+(2*pi*f*c)^2));

目标函数(a,b,c,d : i,s, e11,e22 : max, f, e1, e2)=
{
    i=-1, s=0, while{++i<max, toreal[g(a,b,c,d,f(i,0)), &e11,&e22], s=s+[e1(i,0)-e11]^2+[e2(i,0)-e22]^2 },
    s
};

main(: tArray : max, f, e1, e2)=
{
    tArray=matrix{
//存放实验数据 //f e1[realPart] e2[imagPart]
    "3000000 13.5 17.85
    ... ...数据省略
    100000000 2.01 1.17
    "
    },
    len[tArray,0,&max], f=tArray(all:0), e1=tArray(all:1),e2=tArray(all:2),
//用len函数取矩阵的行数,f等取矩阵的列
    Opt1[@目标函数]
//Opt1函数全局优化
};

结果(a,b,c,d,最小值):

2.393920870707986 18.24067405097755 -3.758382749525807e-008 -1.610640545623026e-003 3.654009362095743

绘图代码:

!!!using["luopt","math","win"]; //使用命名空间

g(a,b,c,d,f)= (a+(b-a)/(1+(2*pi*f*c)^2))-1.0i*(d/(2*pi*f*(1/(36*pi)*1e-9))+((b-a)*2*pi*f*c)/(1+(2*pi*f*c)^2));

ge(a,b,c,d : i,s,e1,e2 : max, f, e12)=
{
    i=-1, while{++i<max, toreal[g(a,b,c,d,f(i,0)), &e1, &e2], e12[i,0]=f(i,0), e12[i,1]=e1, e12[i,2]=e2}
};

init(main : tArray : max, f, e12)=
{
    tArray=matrix{ //存放实验数据 //f e1[realPart] e2[imagPart]
    "3000000 13.5 17.85
    ... ... 数据省略
    100000000 2.01 1.17"
    },
    len[tArray,0,&max], f=tArray(all:0),
//用len函数取矩阵的行数,f等取矩阵的列
    e12=new[real_s,max,3],
    ge[ 2.393920870707986, 18.24067405097755, -3.758382749525807e-008, -1.610640545623026e-003],
    luShareX2[tArray, e12]
//绘制共享X轴视图
};
ChartWnd[@init];
//显示窗口并初始化

图形;

例子8 拟合包含复数的函数的参数

公式: y = yibxilong_infite+((yibxilong_s-yibxilong_infite)/((1+i*Omega*tao)^bet_a));

拟合参数:yibxilong_infite,yibxilong_s,tao,bet_a -- 当作实数时数据与模型匹配度低,故本例均当作复数,匹配度有所提高

数据:Omega y[realPart] y[imagPart]

2 4.37523 0.62119
2.08 4.33845 0.62521
2.16 4.34997 0.63053
2.24 4.31506 0.65416
2.32 4.31287 0.64956
2.4 4.2921 0.6785
2.48 4.27592 0.67085
2.56 4.27271 0.68946
2.64 4.24622 0.69245
2.72 4.24252 0.69161
2.8 4.2235 0.70592
2.88 4.2131 0.69066
2.96 4.20274 0.7088
3.04 4.17474 0.69278
3.12 4.18037 0.69368
3.2 4.15111 0.68694
3.28 4.15204 0.67542
3.36 4.13015 0.68071
3.44 4.11932 0.65229
3.52 4.11557 0.66384
3.6 4.09273 0.64346
3.68 4.10344 0.64157
3.76 4.08152 0.63953
3.84 4.09644 0.63258
3.92 4.07714 0.64613
4 4.07915 0.62682
4.08 4.08866 0.64295
4.16 4.07995 0.64287
4.24 4.09006 0.65073
4.32 4.07778 0.66372
4.4 4.08659 0.65938
4.48 4.09006 0.69362
4.56 4.08228 0.68802
4.64 4.10068 0.71379
4.72 4.0723 0.72193
4.8 4.09413 0.73359
4.88 4.0731 0.75905
4.96 4.0816 0.74975
5.04 4.08084 0.78647
5.12 4.06507 0.77138
5.2 4.0818 0.79578
5.28 4.0533 0.79477
5.36 4.08002 0.7944
5.44 4.0605 0.81608
5.52 4.06823 0.79508
5.6 4.06843 0.82243
5.68 4.05722 0.79832
5.76 4.0792 0.81293
5.84 4.0543 0.81464
5.92 4.07444 0.80338
6 4.05518 0.8267
6.08 4.06451 0.80252
6.16 4.07049 0.84089
6.24 4.04613 0.8307
6.32 4.06599 0.84915
6.4 4.02474 0.86132
6.48 4.04207 0.85187
6.56 4.01568 0.8869
6.64 4.00774 0.861
6.72 4.01495 0.90078
6.8 3.98407 0.88764
6.88 4.01325 0.90488
6.96 3.97776 0.9195
7.04 4.00142 0.90944
7.12 3.99027 0.94878
7.2 3.99063 0.9286
7.28 4.01148 0.96782
7.36 3.98202 0.96614
7.44 4.02715 0.9764
7.52 3.99927 1.00541
7.6 4.0286 0.9848
7.68 4.02259 1.03537
7.76 4.0183 1.00616
7.84 4.04591 1.04756
7.92 4.00783 1.03853
8 4.04809 1.04458
8.08 4.00583 1.06889
8.16 4.02548 1.04709
8.24 4.00542 1.08586
8.32 3.97856 1.05657
8.4 3.99808 1.0767
8.48 3.94704 1.08118
8.56 3.97794 1.066
8.64 3.93227 1.11335
8.72 3.9337 1.07041
8.8 3.93626 1.13698
8.88 3.88795 1.09824
8.96 3.93594 1.15048
9.04 3.86228 1.15887
9.12 3.91482 1.16699
9.2 3.8589 1.22328
9.28 3.88337 1.18775
9.36 3.88369 1.2631
9.44 3.85848 1.23161
9.52 3.90621 1.2869
9.6 3.85492 1.29122
9.68 3.90851 1.28457
9.76 3.8738 1.3413
9.84 3.89225 1.29721
9.92 3.89533 1.37074
10 3.86391 1.32386
10.08 3.91378 1.35101
10.16 3.84783 1.35357
10.24 3.89741 1.33465
10.32 3.84317 1.3746
10.4 3.84082 1.31404
10.48 3.84381 1.36793
10.56 3.79267 1.32282
10.64 3.83601 1.34366
10.72 3.75153 1.35447
10.8 3.78495 1.32402
10.88 3.73595 1.37583
10.96 3.72531 1.32322
11.04 3.73699 1.39732
11.12 3.66998 1.37475
11.2 3.71147 1.40152
11.28 3.64639 1.43016
11.36 3.6951 1.40445
11.44 3.66645 1.49943
11.52 3.64921 1.45173
11.6 3.67278 1.53109
11.68 3.62043 1.50571
11.76 3.70031 1.54225
11.84 3.63168 1.57811
11.92 3.67956 1.54815
12 3.64643 1.61486
12.08 3.65567 1.54865
12.16 3.68399 1.61469
12.24 3.62965 1.58142
12.32 3.69052 1.59474
12.4 3.61013 1.60768
12.48 3.65312 1.56458
12.56 3.60661 1.61913
12.64 3.58948 1.55235
12.72 3.60285 1.6151
12.8 3.52316 1.58248
12.88 3.55328 1.60553
12.96 3.46259 1.6284
13.04 3.4799 1.59773
13.12 3.42929 1.67722
13.2 3.38338 1.63751
13.28 3.37557 1.70799
13.36 3.29446 1.68937
13.44 3.32208 1.7159
13.52 3.23213 1.76005
13.6 3.23184 1.72701
13.68 3.18662 1.79163
13.76 3.1535 1.73895
13.84 3.16797 1.78983
13.92 3.10565 1.77051
14 3.1318 1.75787
14.08 3.07723 1.77302
14.16 3.11093 1.70777
14.24 3.12093 1.74051
14.32 3.10783 1.68152
14.4 3.15819 1.68976
14.48 3.11168 1.66027
14.56 3.18533 1.61469
14.64 3.16197 1.64149
14.72 3.21006 1.57209
14.8 3.22436 1.6219
14.88 3.21211 1.55072
14.96 3.25935 1.57393
15.04 3.20934 1.55437
15.12 3.26857 1.52734
15.2 3.2262 1.55872
15.28 3.23214 1.51635
15.36 3.23624 1.56765
15.44 3.18766 1.52781
15.52 3.21679 1.5537
15.6 3.15515 1.5535
15.68 3.18251 1.53893
15.76 3.13452 1.58305
15.84 3.12766 1.55072
15.92 3.11925 1.6005
16 3.07438 1.58123
16.08 3.09502 1.59905
16.16 3.04198 1.61886
16.24 3.0635 1.62175
16.32 3.03214 1.68653
16.4 3.01135 1.66188
16.48 3.01369 1.71296
16.56 2.97412 1.71076
16.64 3.0108 1.75704
16.72 2.9509 1.79014
16.8 2.94444 1.79493
16.88 2.91727 1.84866
16.96 2.92702 1.82523
17.04 2.93201 1.88693
17.12 2.87968 1.89585
17.2 2.90914 1.91016
17.28 2.86698 1.89892
17.36 2.89035 1.86992
17.44 2.87119 1.90592
17.52 2.86254 1.87056
17.6 2.85872 1.88602
17.68 2.82335 1.8428
17.76 2.86058 1.8239
17.84 2.81505 1.82252
17.92 2.82796 1.7813

Lu脚本代码:

!!!using["luopt","math","win"]; //使用命名空间

y(yibxilong_infite,yibxilong_s,tao,bet_a,Omega)= yibxilong_infite+((yibxilong_s-yibxilong_infite)/((1+1i*Omega*tao)^bet_a));

目标函数(yibxilong_infite,y1,yibxilong_s,y2,tao,t,bet_a,b : i,s, e11,e22 : max, f, e1, e2)=
{
    i=-1, s=0, while{++i<max, toreal[y(yibxilong_infite$y1,yibxilong_s$y2,tao$t,bet_a$b,f(i,0)), &e11,&e22], s=s+[e1(i,0)-e11]^2+[e2(i,0)-e22]^2 },
    s
};

main(: tArray : max, f, e1, e2)=
{
    tArray=matrix{ //存放实验数据 //f e1[realPart] e2[imagPart]
    "2 4.37523 0.62119
    2.08 4.33845 0.62521
    。。。省略数据
    17.92 2.82796 1.7813"
    },
    len[tArray,0,&max], f=tArray(all:0), e1=tArray(all:1),e2=tArray(all:2), //用len函数取矩阵的行数,f等取矩阵的列
    Opt1[@目标函数, optwaysimdeep, optwayconfra] //Opt1函数全局优化
};

结果:

3.65747430967335 1.6268180322264 4.318672716624971 0.3945758454948024 -8.820890091794004e-004 8.261459468079664e-002 -0.7118355007262194 2.904330877799051e-003 1.837635420785162

绘图代码:

!!!using["luopt","math","win"]; //使用命名空间

y(yibxilong_infite,yibxilong_s,tao,bet_a,Omega)= yibxilong_infite+((yibxilong_s-yibxilong_infite)/((1+1i*Omega*tao)^bet_a));

init(main : tArray,ty,i,e11,e22 : max, f, e1, e2)=
{
    tArray=matrix{ //存放实验数据 //f e1[realPart] e2[imagPart]
    "2 4.37523 0.62119
    ... ...省略数据
    17.92 2.82796 1.7813"
    },
    len[tArray,0,&max], f=tArray(all:0), e1=tArray(all:1),e2=tArray(all:2),
//用len函数取矩阵的行数,f等取矩阵的列
    ty=copy[tArray],
    i=-1, while{++i<max,
        toreal[y(3.657474304181978 $ 1.626818030899468 , 4.318672818363945 $ 0.3945758491533642 , -8.820892851511673e-004 $ 8.261459473645558e-002 , -0.711835511823157 $ 2.904319481003881e-003 ,f(i,0)), &e11,&e22],
        ty(i,1)=e11, ty(i,2)=e22
    },
    luShareX2[tArray, ty]
//绘制共享X轴视图
};
ChartWnd[@init];
//显示窗口并初始化

图形:

例子9 带参数(中间变量)传递的微分方程拟合

微分方程:y'=p1*y+p2*x

拟合参数:p1, p2

数据: t   y   x

0 1.25664E-07 0
0.005 0.114139113 0.156918191
0.01 0.224758442 0.31286893
0.015 0.324559128 0.466890728
0.02 0.41208094 0.618033989
0.025 0.487772126 0.765366865
0.03 0.552888138 0.907980999
0.035 0.608846185 1.044997129
0.04 0.656951948 1.175570505
0.045 0.698318143 1.298896097
0.05 0.733862018 1.414213562
0.055 0.764329017 1.520811931
0.06 0.790321332 1.618033989
0.065 0.812323879 1.705280329
0.07 0.83072574 1.782013048
0.075 0.845837104 1.847759065
0.08 0.857902342 1.902113033
0.085 0.867109907 1.944739841
0.09 0.873599694 1.975376681
0.095 0.87746832 1.993834667
0.1 0.878772672 2

分析:微分方程计算时要传递参数x(中间变量),这对Lu脚本来说不是问题,gsl支持的Lu扩展数学库中求解微分方程的函数gsl_ode中提供了该功能。

进一步讨论:微分方程求解传递中间变量时,按时间t的改变量进行插值,可提高求解精度。例如,进行线性插值,需传递一个数组 x4[t0,t1,x0,x1],这样 x=x0+(x1-x0)*[(t-t0)/(t1-t0)] 。

Lu脚本代码1:传递当前中间变量,gsl_ode单步计算

!!!using["luopt","math","win"]; //使用命名空间
f(t,y,dy,x4 : t0,t1,x0,x1,x : p1,p2)=
{
    in[x4, 0, &t0,&t1,&x0,&x1],
//取出数据,x4即从gsl_ode传过来的参数
    x=x0+(x1-x0)*[(t-t0)/(t1-t0)],
//线性插值计算x
    dy=p1*y+p2*x,
    0 //成功返回0
};
目标函数(_p1,_p2 : i,s,y0,tf : tA,yA,xA,max,x4,t2,y1, p1,p2)=
{
    p1=_p1, p2=_p2,
//传递优化变量
    y0=1.25664E-07,
//初值
    i=0, s=0, while{++i<max,
        SetArray{x4 : tA[i-1,0], tA[i,0], xA[i-1,0], xA[i,0]},
//准备传递数据
        t2[0]=tA[i-1,0], t2[1]=tA[i,0], y1[0]=y0,
        //最后一个参数50表示gsl_ode函数在计算时,最多循环计算50次,这样可以提高速度
        tf=gsl_ode[@f, nil, x4, t2, y1, 1e-6, 1e-6, gsl_rk4, 1e-6,50],
        y0=tf(1,1), s=s+[y0-yA(i,0)]^2.0
    },
    s
};
main(: tyArray : tA,yA,xA,max,x4,t2,y1)=
{
    tyArray=matrix{
//存放实验数据
    "0 1.25664E-07 0
    0.005 0.114139113 0.156918191
    ... ... 省略数据
    0.1 0.878772672 2
    "
    },
    len[tyArray,0,&max], tA=tyArray(all:0), yA=tyArray(all:1), xA=tyArray(all:2),
//用len函数取矩阵的行数,tA取矩阵的列
    x4=ra1[0,0,0,0], t2=ra1[0,0], y1=ra1[0],
//准备数组
    Opt1[@目标函数]
//Opt1函数全局优化
};

结果:

-32997.3853793671 15830.42381251212 0.1177179718140771

Lu脚本代码2:使用 nil 参数,传递当前中间变量索引,gsl_ode整体计算

!!!using["luopt","math","win"]; //使用命名空间
f(t,y,dy,i : x : tA,xA,p1,p2)=
{
    x=xA[i]+(xA[i+1]-xA[i])*[(t-tA[i])/(tA[i+1]-tA[i])],
//线性插值计算x
    dy=p1*y+p2*x,
    0
//成功返回0
};
目标函数(_p1,_p2 : i,s,tf : tA,yA,xA,max, p1,p2)=
{
    p1=_p1, p2=_p2,
//传递优化变量
    //最后一个参数50表示gsl_ode函数在计算时,最多循环计算50次,这样可以提高速度
    tf=gsl_ode[@f, nil, nil, tA, ra1(1.25664E-07), 1e-6, 1e-6, gsl_rk4, 1e-6,50],
    i=0, s=0, while{++i<max, s=s+[tf(i,1)-yA(i)]^2.0 },
    s
};
main(: tyArray : tA,yA,xA,max)=
{
    tyArray=matrix{
//存放实验数据
    "0 1.25664E-07 0
    0.005 0.114139113 0.156918191
    ... ... 数据省略
    0.1 0.878772672 2
    "
    },
    len[tyArray,0,&max], tA=tyArray(all:0).reshape(), yA=tyArray(all:1).reshape(), xA=tyArray(all:2).reshape(),
//用len函数取矩阵的行数,tA取矩阵的列
    Opt1[@目标函数]
//Opt1函数全局优化
};

结果:

-32995.38534216554 15829.22215148295 0.1177179963080844

因本模型与数据匹配不好,故本例不提供绘图,本例的意义在于微分方程拟合时如何传递中间参数。


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