Lu 微分方程参数优化（拟合）
例子1 数学模型：
dy/dt=k*y*z+0.095*b*z
dz/dt=-b*z-0.222*z
实验数据（ti; yi）：
ti yi
0.012.329101563
0.683.851292783
1.14.50093936
1.636.749172247
2.079.112062872
2.679.691657366
3.0911.16928013
3.6410.91451823
4.6516.44279204
5.118.29615885
5.5821.63989025
6.1125.78611421
6.6326.34282633
7.0626.50581287
7.6227.63951823
8.6635.02757626
9.0435.5562035
9.6336.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, gsl_rk4, 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, gsl_rk4, 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, gsl_rk4, 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, gsl_rk4, 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, gsl_rk4, 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, gsl_rk4, 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, gsl_rk4, 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, gsl_rk4, 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, gsl_rk4, 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, gsl_rk4, 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, gsl_rk4, 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, gsl_rk4, 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, 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{ //存放实验数据
"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, gsl_rk4, 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, gsl_rk4, 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, gsl_rk4, 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, gsl_rk4, 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
-310.2286982374551 -1.932248685024008 962425592.5302644 -1.280249791209473e-006 -22306228.22878688 6.553909044211868e-008 144.0493473573013（此组解不能确定是否有效）
绘图：
!!!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, gsl_rk4, 1e-6,100]), outa[gsl_ode[@f,nil,0.0,t,ra1(2.654977892066841,0), 1e-6, 1e-6, gsl_rk4, 1e-6,100]];
ChartWnd[@init]; //显示窗口并初始化
例子9 隐函数式微分方程参数拟合：
根据两个式子（一个隐函数式微分方程，一个代数方程）拟合7个未知参数fac1~fac7。式中lambv'=d(lambv)/d(lamb)，lambv为中间变量。
lambv'=((1.0/3/fac7)*(fac1*((lamb/lambv')^fac2-(lambv'/lamb)^(0.5*fac2))+fac3*((lamb/lambv')^fac4-(lambv'/lamb)^(0.5*fac4))+fac5*((lamb/lambv')^fac6-(lambv'/lamb)^(0.5*fac6))));
p=(pa+fac1*((lamb/lambv)^fac2/lamb-(lambv/lamb)^(fac2*0.5)/lamb)+fac3*((lamb/lambv)^fac4/lamb-(lambv/lamb)^(fac4*0.5)/lamb)+fac5*((lamb/lambv)^fac6/lamb-(lambv/lamb)^(fac6*0.5)/lamb));
实验数据：
|
t |
lamb |
pa |
p |
|
0 |
0.959926667 |
-0.218715994 |
-0.060966222 |
|
0.25 |
0.963063333 |
-0.200840559 |
-0.046836278 |
|
0.5 |
0.96612 |
-0.183553472 |
-0.035042256 |
|
0.75 |
0.969016667 |
-0.167289626 |
-0.024438988 |
|
1 |
0.97171 |
-0.152268899 |
-0.01706765 |
|
1.25 |
0.97409 |
-0.139075651 |
-0.010939712 |
|
1.5 |
0.976123333 |
-0.127862697 |
-0.006795087 |
|
1.75 |
0.977786667 |
-0.118729759 |
-0.002919842 |
|
2 |
0.979000001 |
-0.112089924 |
-0.001981819 |
|
2.25 |
0.979743333 |
-0.108031322 |
-0.001163676 |
|
2.5 |
0.98001 |
-0.106577018 |
-0.003491456 |
|
2.75 |
0.979766667 |
-0.107904033 |
-0.009723697 |
分析：未知中间变量lambv用解方程法求解。
Lu代码：
!!!using["luopt","math"]; //使用命名空间
g(lambv :: fac1,fac2,fac3,fac4,fac5,fac6,fac7,lamb)= ((1.0/3/fac7)*(fac1*((lamb/lambv)^fac2-(lambv/lamb)^(0.5*fac2))+fac3*((lamb/lambv)^fac4-(lambv/lamb)^(0.5*fac4))+fac5*((lamb/lambv)^fac6-(lambv/lamb)^(0.5*fac6))))-lambv;
f(t,_lamb,lambv,pp : x : fac1,fac2,fac3,fac4,fac5,fac6,fac7,lamb) =
lamb=_lamb, x=lamb,
if{pqrt[@g,&x,1e-6]<=0, return(1)}, //连分式法解方程，返回迭代次数，迭代次数<=0，求解失败
lambv=x,
0; //成功返回0
pp(pa,p,lamb,lambv :: fac1,fac2,fac3,fac4,fac5,fac6,fac7) =
(pa+fac1*((lamb/lambv)^fac2/lamb-(lambv/lamb)^(fac2*0.5)/lamb)+fac3*((lamb/lambv)^fac4/lamb-(lambv/lamb)^(fac4*0.5)/lamb)+fac5*((lamb/lambv)^fac6/lamb-(lambv/lamb)^(fac6*0.5)/lamb)) - p;
目标函数(_fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7 : i,s,tf : t_A, t_lamb, t_pa, t_p, max, fac1,fac2,fac3,fac4,fac5,fac6,fac7,la)=
{
fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7, //传递优化变量
la[0]=t_lamb[0,0],
//最后一个参数50表示gsl_ode函数在计算时,最多循环计算50次,这样可以提高速度
tf=gsl_ode[@f, nil, 0.0, t_A, la, 1e-6, 1e-6, gsl_rk4, 1e-6,50],
i=-1, s=0, while{++i<max,
s=s+pp[t_pa(i,0),t_p(i,0),t_lamb(i,0),tf(i,1)]^2.0
},
s
};
main(: tArray : t_A, t_lamb, t_pa, t_p, max, la)=
{
tArray=matrix{ //存放实验数据//t,lamb,pa,p
"0 0.959926667 -0.218715994 -0.060966222
0.25 0.963063333 -0.200840559 -0.046836278
0.5 0.96612 -0.183553472 -0.035042256
0.75 0.969016667 -0.167289626 -0.024438988
1 0.97171 -0.152268899 -0.01706765
1.25 0.97409 -0.139075651 -0.010939712
1.5 0.976123333 -0.127862697 -0.006795087
1.75 0.977786667 -0.118729759 -0.002919842
2 0.979000001 -0.112089924 -0.001981819
2.25 0.979743333 -0.108031322 -0.001163676
2.5 0.98001 -0.106577018 -0.003491456
2.75 0.979766667 -0.107904033 -0.009723697"
},
len[tArray,0,&max], t_A=tArray(all:0), t_lamb=tArray(all:1), t_pa=tArray(all:2), t_p=tArray(all:3), //用len函数取矩阵的行数，t_A等取矩阵的列
la=ra1[0], //预先申请数组
Opt1[@目标函数, optwaysimdeep, optwayconfra] //Opt1函数全局优化
};
结果;
-0.2713835810685396 7.198887084780003 -0.1092823587604948 0.5705253143288862 -0.2717119875995815 -3.59893525610079 -0.2229088929623526 2.501917969994917e-002
绘图代码：
!!!using["luopt","math","win"]; //使用命名空间
g(lambv :: fac1,fac2,fac3,fac4,fac5,fac6,fac7,lamb)= ((1.0/3/fac7)*(fac1*((lamb/lambv)^fac2-(lambv/lamb)^(0.5*fac2))+fac3*((lamb/lambv)^fac4-(lambv/lamb)^(0.5*fac4))+fac5*((lamb/lambv)^fac6-(lambv/lamb)^(0.5*fac6))))-lambv;
f(t,_lamb,lambv,pp : x : fac1,fac2,fac3,fac4,fac5,fac6,fac7,lamb) =
lamb=_lamb, x=lamb,
if{pqrt[@g,&x,1e-6]<=0, return(1)}, //连分式法解方程，返回迭代次数，迭代次数<=0，求解失败
lambv=x,
0; //成功返回0
pp(pa,lamb,lambv :: fac1,fac2,fac3,fac4,fac5,fac6,fac7) =
(pa+fac1*((lamb/lambv)^fac2/lamb-(lambv/lamb)^(fac2*0.5)/lamb)+fac3*((lamb/lambv)^fac4/lamb-(lambv/lamb)^(fac4*0.5)/lamb)+fac5*((lamb/lambv)^fac6/lamb-(lambv/lamb)^(fac6*0.5)/lamb));
set(_fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7 :: fac1,fac2,fac3,fac4,fac5,fac6,fac7)=
{
fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7 //传递优化变量
};
init(main : tArray, tf, i, k, t_pp,t_lambv : t_A, t_lamb, t_pa, t_p, max)=
{
tArray=matrix{ //存放实验数据//t,lamb,pa,p
"0 0.959926667 -0.218715994 -0.060966222
0.25 0.963063333 -0.200840559 -0.046836278
0.5 0.96612 -0.183553472 -0.035042256
0.75 0.969016667 -0.167289626 -0.024438988
1 0.97171 -0.152268899 -0.01706765
1.25 0.97409 -0.139075651 -0.010939712
1.5 0.976123333 -0.127862697 -0.006795087
1.75 0.977786667 -0.118729759 -0.002919842
2 0.979000001 -0.112089924 -0.001981819
2.25 0.979743333 -0.108031322 -0.001163676
2.5 0.98001 -0.106577018 -0.003491456
2.75 0.979766667 -0.107904033 -0.009723697"
},
len[tArray,0,&max], t_A=tArray(all:0), t_lamb=tArray(all:1), t_pa=tArray(all:2), t_p=tArray(all:3), //用len函数取矩阵的行数，t_A等取矩阵的列
set[-0.2713835810685396,7.198887084780003, -0.1092823587604948, 0.5705253143288862, -0.2717119875995815, -3.59893525610079, -0.2229088929623526 ],
t_pp=new[real_s,max], t_lambv=new[real_s,max],
//最后一个参数50表示gsl_ode函数在计算时,最多循环计算50次,这样可以提高速度
tf=gsl_ode[@f, nil, 0.0, t_A, ra1(t_lamb[0,0]), 1e-6, 1e-6, gsl_rk4, 1e-6,50],
i=-1, while{++i<max, t_pp[i]=pp[t_pa(i,0),t_lamb(i,0),tf(i,1)], t_lambv[i]=tf(i,1)},
t_A=reshape(t_A),
o[t_lambv, t_pp],
cwAttach[typeSplit], cwResizePlots(1,2,2), //左二右二分裂
k=cwAddCurve{t_A, t_p.reshape(), max, 0}, //给0子图添加曲线
cwSetScatter(k,0), //设置绘制点
cwSetDataLineSize(5, k, 0), //设置点的大小
cwAddCurve{t_A, t_pp, max, 0}, //给0子图添加曲线
cwAddCurve{t_A, t_lamb.reshape(), max, 1}, //给1子图添加曲线
cwAddCurve{t_A, t_pa.reshape(), max, 2}, //给2子图添加曲线
cwAddCurve{t_A, t_lambv, max, 3} //给3子图添加曲线
};
ChartWnd[@init];
图形：
例子10 含有未知中间变量同时需要传递其他中间变量的微分方程参数拟合
拟合参数：fac1～fac7
未知中间变量：lambv
需传递中间变量：lamb, pa
常量：A=0.02,w=2*pi*0.1;
lambdif=A*w*cos(w*t);
lambvdif=((1.0/3/fac7)*(fac1*((lamb/lambv)^fac2-(lambv/lamb)^(0.5*fac2))+fac3*((lamb/lambv)^fac4-(lambv/lamb)^(0.5*fac4))+fac5*((lamb/lambv)^fac6-(lambv/lamb)^(0.5*fac6))));
微分方程：p'=pa+fac1*((fac2-1)*lamb^(fac2-2)*lambv^(-fac2)*lambdif-fac2*lamb^(fac2-1)*lambv^(-fac2-1)*lambvdif
+(0.5*fac2+1)*lamb^(-0.5*fac2-2)*lambv^(0.5*fac2)*lambdif-0.5*fac2*lamb^(-0.5*fac2-1)*lambv^(0.5*fac2-1)*lambvdif)
+fac3*((fac4-1)*lamb^(fac4-2)*lambv^(-fac4)*lambdif-fac4*lamb^(fac4-1)*lambv^(-fac4-1)*lambvdif+(0.5*fac4+1)*lamb^(-0.5*fac4-2)*lambv^(0.5*fac4)*lambdif
-0.5*fac4*lamb^(-0.5*fac4-1)*lambv^(0.5*fac4-1)*lambvdif)+fac5*((fac6-1)*lamb^(fac6-2)*lambv^(-fac6)*lambdif
-fac6*lamb^(fac6-1)*lambv^(-fac6-1)*lambvdif+(0.5*fac6+1)*lamb^(-0.5*fac6-2)*lambv^(0.5*fac6)*lambdif
-0.5*fac6*lamb^(-0.5*fac6-1)*lambv^(0.5*fac6-1)*lambvdif);
数据：t,lamb,pa,p 同例子9
分析：
1、微分方程求解时要传递中间变量：gsl支持的Lu扩展数学库中求解微分方程的函数gsl_ode中提供了该功能。
进一步讨论：微分方程求解传递中间变量时，按时间t的改变量进行插值，可提高求解精度。例如，进行线性插值：lamb=lamb0+(lamb1-lamb0)*[(t-t0)/(t1-t0)], pa=pa0+(pa1-pa0)*[(t-t0)/(t1-t0)]。
2、lambv为未知中间变量，使用幂级数逼近lambv进行拟合的方法，需增加k0,k1,k2,k3四个拟合参数；当然拟合参数多时精度高，但耗时长。
Lu代码：
!!!using["luopt","math"]; //使用命名空间
f(t,p,dp, i : lamb,pa, A,w,lambdif,lambvdif,lambv: t_A, t_lamb, t_pa, fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4) =
{
lamb=t_lamb[i]+(t_lamb[i+1]-t_lamb[i])*[(t-t_A[i])/(t_A[i+1]-t_A[i])], pa=t_pa[i]+(t_pa[i+1]-t_pa[i])*[(t-t_A[i])/(t_A[i+1]-t_A[i])], //线性插值计算lamb,pa
lambv=k0+k1*lamb+k2*lamb*lamb+k3*lamb^3+k4*lamb^4, //lambv=f(lamb,k0,k1,k2,k3,k4)，使用幂级数逼近
A=0.02,w=2*pi*0.1,
lambdif=A*w*cos(w*t),
lambvdif=((1.0/3/fac7)*(fac1*((lamb/lambv)^fac2-(lambv/lamb)^(0.5*fac2))+fac3*((lamb/lambv)^fac4-(lambv/lamb)^(0.5*fac4))+fac5*((lamb/lambv)^fac6-(lambv/lamb)^(0.5*fac6)))),
dp=pa+fac1*((fac2-1)*lamb^(fac2-2)*lambv^(-fac2)*lambdif-fac2*lamb^(fac2-1)*lambv^(-fac2-1)*lambvdif
+(0.5*fac2+1)*lamb^(-0.5*fac2-2)*lambv^(0.5*fac2)*lambdif-0.5*fac2*lamb^(-0.5*fac2-1)*lambv^(0.5*fac2-1)*lambvdif)
+fac3*((fac4-1)*lamb^(fac4-2)*lambv^(-fac4)*lambdif-fac4*lamb^(fac4-1)*lambv^(-fac4-1)*lambvdif+(0.5*fac4+1)*lamb^(-0.5*fac4-2)*lambv^(0.5*fac4)*lambdif
-0.5*fac4*lamb^(-0.5*fac4-1)*lambv^(0.5*fac4-1)*lambvdif)+fac5*((fac6-1)*lamb^(fac6-2)*lambv^(-fac6)*lambdif
-fac6*lamb^(fac6-1)*lambv^(-fac6-1)*lambvdif+(0.5*fac6+1)*lamb^(-0.5*fac6-2)*lambv^(0.5*fac6)*lambdif
-0.5*fac6*lamb^(-0.5*fac6-1)*lambv^(0.5*fac6-1)*lambvdif),
0 //必须返回0
};
目标函数(_fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7,_k0,_k1,_k2,_k3,_k4 : tf : t_A, t_p, fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4)=
{
fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7, k0=_k0, k1=_k1, k2=_k2, k3=_k3, k4=_k4, //传递优化变量
//最后一个参数50表示gsl_ode函数在计算时,最多循环计算50次,这样可以提高速度
tf=gsl_ode[@f, nil, nil, t_A, ra1(-0.060966222), 1e-6, 1e-6, gsl_rk4, 1e-6,50],
sum{[tf(all:1).reshape()-t_p].^2.0} //代码矢量化加速，特别适合大数据量
};
main(: tArray : t_A, t_lamb, t_pa, t_p)=
{
tArray=matrix{ //存放实验数据//t,lamb,pa,p
"0 0.959926667 -0.218715994 -0.060966222
0.25 0.963063333 -0.200840559 -0.046836278
0.5 0.96612 -0.183553472 -0.035042256
0.75 0.969016667 -0.167289626 -0.024438988
1 0.97171 -0.152268899 -0.01706765
1.25 0.97409 -0.139075651 -0.010939712
1.5 0.976123333 -0.127862697 -0.006795087
1.75 0.977786667 -0.118729759 -0.002919842
2 0.979000001 -0.112089924 -0.001981819
2.25 0.979743333 -0.108031322 -0.001163676
2.5 0.98001 -0.106577018 -0.003491456
2.75 0.979766667 -0.107904033 -0.009723697"
},
t_A=tArray(all:0).reshape(), t_lamb=tArray(all:1).reshape(), t_pa=tArray(all:2).reshape(), t_p=tArray(all:3).reshape(), //t_A等取矩阵的列，并转为一维数组
Opt1[@目标函数, optwaysimdeep, optwayconfra] //Opt1函数全局优化
};
结果（fac1～fac7，k0,k1,k2,k3,k4，最小值）：
-8.801359745756983e-004 -6.149102784564371e-009 -1.11091297833596e-018 -4.917163751655897 0.4722303344423368 0.7444565347550796 6.428499397806392e-002 31250.06922251955 -3733.557030362606 15223.1661928006 -2.760725760727475e-008 -45252.4633856661 1.835645989327923e-006
绘图代码：
!!!using["luopt","math","win"]; //使用命名空间
f(t,p,dp, i : lamb,pa, A,w,lambdif,lambvdif,lambv: t_A, t_lamb, t_pa, fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4) =
{
lamb=t_lamb[i]+(t_lamb[i+1]-t_lamb[i])*[(t-t_A[i])/(t_A[i+1]-t_A[i])], pa=t_pa[i]+(t_pa[i+1]-t_pa[i])*[(t-t_A[i])/(t_A[i+1]-t_A[i])], //线性插值计算lamb,pa
lambv=k0+k1*lamb+k2*lamb*lamb+k3*lamb^3+k4*lamb^4, //lambv=f(lamb,k0,k1,k2,k3,k4)，使用幂级数逼近
A=0.02,w=2*pi*0.1,
lambdif=A*w*cos(w*t),
lambvdif=((1.0/3/fac7)*(fac1*((lamb/lambv)^fac2-(lambv/lamb)^(0.5*fac2))+fac3*((lamb/lambv)^fac4-(lambv/lamb)^(0.5*fac4))+fac5*((lamb/lambv)^fac6-(lambv/lamb)^(0.5*fac6)))),
dp=pa+fac1*((fac2-1)*lamb^(fac2-2)*lambv^(-fac2)*lambdif-fac2*lamb^(fac2-1)*lambv^(-fac2-1)*lambvdif
+(0.5*fac2+1)*lamb^(-0.5*fac2-2)*lambv^(0.5*fac2)*lambdif-0.5*fac2*lamb^(-0.5*fac2-1)*lambv^(0.5*fac2-1)*lambvdif)
+fac3*((fac4-1)*lamb^(fac4-2)*lambv^(-fac4)*lambdif-fac4*lamb^(fac4-1)*lambv^(-fac4-1)*lambvdif+(0.5*fac4+1)*lamb^(-0.5*fac4-2)*lambv^(0.5*fac4)*lambdif
-0.5*fac4*lamb^(-0.5*fac4-1)*lambv^(0.5*fac4-1)*lambvdif)+fac5*((fac6-1)*lamb^(fac6-2)*lambv^(-fac6)*lambdif
-fac6*lamb^(fac6-1)*lambv^(-fac6-1)*lambvdif+(0.5*fac6+1)*lamb^(-0.5*fac6-2)*lambv^(0.5*fac6)*lambdif
-0.5*fac6*lamb^(-0.5*fac6-1)*lambv^(0.5*fac6-1)*lambvdif),
0 //必须返回0
};
set(_fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7,_k0,_k1,_k2,_k3,_k4 :: fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4)=
{
fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7, k0=_k0, k1=_k1, k2=_k2, k3=_k3, k4=_k4 //传递优化变量
};
init(init: tArray,t_lambv,i,tf,k,lamb : t_A, t_lamb, t_pa, t_p, max, fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4)=
{
tArray=matrix{ //存放实验数据//t,lamb,pa,p
"0 0.959926667 -0.218715994 -0.060966222
... ... 数据省略
2.75 0.979766667 -0.107904033 -0.009723697"
},
len[tArray,0,&max], t_A=tArray(all:0).reshape(), t_lamb=tArray(all:1).reshape(), t_pa=tArray(all:2).reshape(), t_p=tArray(all:3).reshape(),
set[-8.801359745756983e-004, -6.149102784564371e-009, -1.11091297833596e-018, -4.917163751655897, 0.4722303344423368, 0.7444565347550796, 6.428499397806392e-002, 31250.06922251955, -3733.557030362606, 15223.1661928006, -2.760725760727475e-008, -45252.4633856661 ],
t_lambv=new[real_s,max],
i=-1, while{++i<max, lamb=t_lamb[i], t_lambv[i]=k0+k1*lamb+k2*lamb*lamb+k3*lamb^3+k4*lamb^4},
tf=gsl_ode[@f, nil, nil, t_A, ra1(-0.060966222), 1e-6, 1e-6, gsl_rk4, 1e-6,50],
cwAttach[typeSplit], cwResizePlots(1,2,2), //左二右二分裂
k=cwAddCurve{t_A, t_p, max, 0}, //给0子图添加曲线
cwSetScatter(k,0), //设置绘制点
cwSetDataLineSize(5, k, 0), //设置点的大小
cwAddCurve{t_A, tf(all:1).reshape(), max, 0}, //给0子图添加曲线
cwAddCurve{t_A, t_lamb, max, 1}, //给1子图添加曲线
cwAddCurve{t_A, t_pa, max, 2}, //给2子图添加曲线
cwAddCurve{t_A, t_lambv, max, 3} //给3子图添加曲线
};
ChartWnd[@init];
图形：
版权所有© Lu程序设计
2002-2021，保留所有权利
E-mail: forcal@sina.com
QQ：630715621
最近更新:
2021年08月21日