-
原问题地址：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 -
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 微分方程参数优化（拟合）
例子1 数学模型：
dy/dt=
k*y*z+0.095*b*z
dz/dt=-b*z-0.222*z实验数据（ti; yi）：
ti yi
0.012.3291015630.683.8512927831.14.500939361.636.7491722472.079.1120628722.679.6916573663.0911.169280133.6410.914518234.6516.442792045.118.296158855.5821.639890256.1125.786114216.6326.342826337.0626.505812877.6227.639518238.6635.027576269.0435.55620359.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, 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.004658Lu代码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
-
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
