• 5星
1.44MB m0_52957036 2020-09-14 08:52:51
• 5星
198KB weixin_48549810 2021-01-06 14:45:40
• 5星
3KB qq_46753404 2021-03-15 22:29:06
• 5星
11KB weixin_42696333 2021-09-10 20:59:17
• 5星
3KB weixin_42696333 2021-09-10 21:17:12
• 5星
567KB SKCQTGZX 2021-09-03 15:03:25
• 微分方程参数拟合-过拟合现象及处理 微分方程 过拟合 优化

原问题地址：matlab 分式拟合,matlab 微分方程组的参数拟合_催眠神兔的博客-CSDN博客 微分方程式： x'=dx/dt=a*0.0321*(b-x)-d*x-dy/dt, y'= dy/dt=0.25*p1*exp(-p1*t)*x , 四个待求参数：a、b、d、p1 t、x、y...
原问题地址：matlab 分式拟合,matlab 微分方程组的参数拟合_催眠神兔的博客-CSDN博客
微分方程式：
x'=dx/dt=a*0.0321*(b-x)-d*x-dy/dt,
y'= dy/dt=0.25*p1*exp(-p1*t)*x ,
四个待求参数：a、b、d、p1
t、x、y数据见下面：
//0 0 0  //这是初值
0,0,0,
0.1,0.486966799,0.048018378,
0.167,1.6657,0.05823,
0.2,0.860306078,0.060834243,
0.3,1.156255213,0.064254733,
0.4,1.390856542,0.065167644,
0.5,1.67518,0.06638,
0.6,1.724247244,0.065476325,
0.7,1.841108525,0.065493681,
0.8,1.93374543,0.065498314,
0.9,2.007179471,0.06549955,
1,1.92438,0.05641,
1.1,2.111536196,0.065499968,
1.2,2.148115682,0.065499991,
1.3,2.177112544,0.065499998,
1.4,2.200098596,0.065499999,
1.5,2.218319829,0.0655,
1.6,2.232763952,0.0655,
1.7,2.244213928,0.0655,
1.8,2.253290418,0.0655,
1.9,2.260485427,0.0655,
2,2.16156,0.07359,
2.1,2.270710216,0.0655,
2.2,2.274294246,0.0655,
2.3,2.277135335,0.0655,
2.4,2.27938749,0.0655,
2.5,2.281172792,0.0655,
2.6,2.282588016,0.0655,
2.7,2.283709875,0.0655,
2.8,2.284599183,0.0655,
2.9,2.285304144,0.0655,
3,2.44976,0.15449,
3.1,2.286305961,0.0655,
3.2,2.286657121,0.0655,
3.3,2.286935489,0.0655,
3.4,2.287156153,0.0655,
3.5,2.287331076,0.0655,
3.6,2.287469738,0.0655,
3.7,2.287579657,0.0655,
3.8,2.287666791,0.0655,
3.9,2.287735862,0.0655,
4,2.287790616,0.0655,
4.1,2.287834019,0.0655,
4.2,2.287868426,0.0655,
4.3,2.2878957,0.0655,
4.4,2.287917321,0.0655,
4.5,2.287934459,0.0655,
4.6,2.287948045,0.0655,
4.7,2.287958815,0.0655,
4.8,2.287967352,0.0655
分析：
微分方程中
Lu脚本代码（使用OpenLu+lugslmath计算，可从www.forcal.net下载）：
!!!using["luopt","math"]; //使用命名空间
{
dy=0.25*p1*exp(-p1*t)*x,
0 //必须返回0
};
{
//最后一个参数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[@目标函数] //也可以使用此语句
};
196.3535183967741         2.775743871883714         20.18977085052178         0.9084940395557356
以下Lu代码可绘制图形：
!!!using["luopt","math","win"]; //使用命名空间
{
dy=0.25*p1*exp(-p1*t)*x,
0 //必须返回0
};
{
//最后一个参数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
};
{
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取矩阵的列
tyz=gsl_ode[@f,nil,0.0,tA,ra1(0,0), 1e-6, 1e-6, 2, 1e-6,100],  //获得最终结果
outa[tyz],                  //outa输出结果
luShareX2[tyArray, tyz]     //绘制共享X轴视图
};
ChartWnd[@init];                //显示窗口并初始化
图形：

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

time        CLA        CHMF
-0.292699813        296.95        371.15
0        0        0.104205377
30        0.004714244        0.097714144
60        0.011371834        0.091620741
90        0.017292886        0.083395946
120        0.024604931        0.076880763
180        0.034326394        0.067098618
240        0.043770544        0.057635083
300        0.050232723        0.049835334
420        0.0630261        0.035480465
540        0.074273358        0.024649385
-1.187859399        296.95        371.15
0        0.000998896        0.104359628
7        0.010982214        0.093943616
15        0.021126176        0.082156616
30        0.037770911        0.064592594
45        0.049576229        0.050176205
60        0.062465132        0.038429704
90        0.070549574        0.025632208
120        0.080625062        0.014099057
150        0.086861135        0.008531942
180        0.092640851        0.003482618
-0.306744318        296.95        371.15
0        0        1.064590325
30        0.035820956        0.99951873
60        0.091231251        0.937544756
90        0.134419728        0.863743863
120        0.175771334        0.812009778
180        0.261121043        0.714875158
240        0.322707477        0.612848744
300        0.374201765        0.547687408
360        0.43844506        0.459826754
480        0.542643307        0.350285652
-1.051427816        296.95        371.15
0        0.005972939        1.009605954
15        0.111613907        0.807385517
30        0.239243562        0.653177373
45        0.323494645        0.551727619
60        0.401448109        0.416047561
90        0.515684189        0.306615783
120        0.611196296        0.16394207
150        0.637330955        0.128679898
180        0.657593008        0.075611085
240        0.736999768        0.019444307
-1.134063515        297.9        414.45
0.883333333        0        0.103346707
1.883333333        0.001512522        0.100892698
2.883333333        0.01367608        0.088755571
3.883333333        0.020184345        0.082041377
4.883333333        0.044262628        0.054390229
5.883333333        0.056989435        0.041914107
6.883333333        0.077149736        0.020978875
7.883333333        0.085874714        0.011700555
8.883333333        0.090409481        0.002783166
10.88333333        0.097292737        0
-0.316500027        297.9        414.45
0        0        1.070970111
2        0.003948991        1.056056057
3        0.020186072        1.037126132
5        0.072206172        0.953098035
7        0.193384568        0.788352458
10        0.331724523        0.630124042
13        0.46896552        0.434878839
16        0.553992395        0.340920015
20        0.64828259        0.221513537
25        0.720649283        0.129803699
-1.134691646        297.9        414.45
0.883333333        0.000951425        1.015091849
1.883333333        0.022369218        0.964249856
2.883333333        0.109778689        0.841350168
3.883333333        0.200071993        0.735440236
4.883333333        0.265436739        0.647223865
5.883333333        0.534716911        0.339015497
7.883333333        0.746720695        0.087667964
9.883333333        0.80490024        0.026375542
11.88333333        0.814507102        0.003078599
13.88333333        0.844622073        0
-0.069553776        297.9        454.05
0        0        0.105069261
0.883333333        0        0.1033609
2.883333333        0.001240168        0.100456984
5.883333333        0.031816596        0.059602151
7.883333333        0.041520146        0.046256377
8.883333333        0.055889344        0.029139086
9.883333333        0.057040816        0.028828737
11.88333333        0.071843951        0.01235586
15.88333333        0.077975291        0.001686074
-0.132438303        297.9        454.05
0        0        0.107494764
0.883333333        0        0.107741493
4.883333333        0.034722983        0.064167986
6.883333333        0.069460981        0.023443174
7.883333333        0.073060273        0.009201692
8.883333333        0.07989269        0.007082479
9.883333333        0.081267505        0.002823071
10.88333333        0.087028096        0.001622815
12.88333333        0.086598397        0
-0.069742593        297.9        454.05
0        0        1.045091436
0.883333333        0        1.019821604
2.883333333        0.027121965        0.980370094
5.883333333        0.237319106        0.612502376
7.883333333        0.327882492        0.463633654
8.883333333        0.45866343        0.262292609
9.883333333        0.476993327        0.245222172
11.88333333        0.544647928        0.113480966
15.88333333        0.599978991        0.03028414
-0.132277891        297.9        454.05
0        0        1.073088085
1        0        1.042241346
4.883333333        0.246301959        0.661045394
6.883333333        0.517594712        0.290417203
7.883333333        0.559102723        0.22683689
8.883333333        0.662520515        0.082134805
9.883333333        0.701897526        0.056834349
10.88333333        0.723673938        0.009283515
12.88333333        0.73410445        0.00186463
-0.301540898        297.9        414.45
0.883333333        0        0.101558639
1.883333333        0.00563862        0.096383724
3.883333333        0.017255774        0.084089401
5.883333333        0.026288775        0.072388068
7.883333333        0.039693626        0.058040049
10.88333333        0.055995293        0.04294237
13.88333333        0.06437459        0.032954514
16.88333333        0.076142889        0.021339873
20.88333333        0.083887249        0.010987306
25.88333333        0.09228189        0
-0.529462294        296.95000         371.15000
0        0        0.104396538
30        0.01086087        0.092633826
60        0.025903593        0.077374404
90        0.035707982        0.066417356
120        0.046876147        0.055103546
180        0.059555879        0.03988887
240        0.071565642        0.027640828
300        0.076601728        0.021842596
420        0.088010322        0.010088388
540        0.095707003        0.003371728
-2.07036895        296.95000         371.15000
0        0.005354629        0.101113608
5        0.015596272        0.089273904
10        0.031459676        0.070240021
20        0.048850215        0.050642749
30        0.064007378        0.034673716
45        0.080303313        0.017800588
60        0.083870308        0.014048253
80        0.091867922        0.00226478
100        0.096076813        0.000809425
120        0.094754675        0
-0.531088765        296.95000         371.15000
0        0        1.085663867
30        0.092227585        0.959687156
60        0.196777828        0.815955815
90        0.26771054        0.678294715
120        0.339669208        0.584243991
180        0.44493791        0.388930384
240        0.546276325        0.31232313
300        0.61758637        0.234488071
360        0.656355053        0.157608349
480        0.705466486        0.074621884
-1.96250305        296.95000         371.15000
0        0.009942844        0.943763649
7        0.103997441        0.792109828
15        0.26409945        0.602246938
30        0.422141533        0.384870483
45        0.559939068        0.204035091
60        0.625774343        0.110517741
90        0.704614433        0.03486973
120        0.730085821        0.011077271
180        0.738694061        0
300        0.725429428        0
-0.543478572        297.90000         414.45000
0.883333333        0        0.103371092
1.883333333        0        0.101492209
2.883333333        0.010860825        0.0913963
3.883333333        0.018569505        0.083314629
4.883333333        0.031600628        0.069524342
5.883333333        0.038345432        0.061930036
6.883333333        0.05561085        0.044439081
8.883333333        0.068330924        0.029366472
10.88333333        0.079149507        0.017914315
13.88333333        0.094041084        0.001441428
-2.047061922        297.90000         414.45000
0.883333333        0        0.103068737
1.383333333        0.000594026        0.101344014
1.883333333        0.004078897        0.099061111
2.383333333        0.010968284        0.091356357
2.883333333        0.025168999        0.075723863
3.383333333        0.034499333        0.067241358
3.883333333        0.056895681        0.044080416
4.883333333        0.076428187        0.02052236
5.883333333        0.079637269        0.01927162
6.883333333        0.100845748        0
-0.544306924        297.90000         414.45000
0.883333333        0        1.009874204
1.883333333        0        1.009745063
2.883333333        0        1.010440284
3.883333333        0.011090946        0.994311593
4.883333333        0.047614676        0.926359288
6.883333333        0.048051427        0.929046366
8.883333333        0.127466089        0.817776336
10.88333333        0.23769185        0.685748319
12.88333333        0.300362249        0.620264623
15.88333333        0.474651916        0.397967837
-2.067390913        297.90000         414.45000
0.883333333        0.003733023        0.971632292
1.383333333        0.006034621        0.962671549
1.883333333        0.039335525        0.916364421
2.383333333        0.072289459        0.887838939
2.883333333        0.077965623        0.874929815
3.883333333        0.41942945        0.461344221
4.883333333        0.68586173        0.136323751
5.883333333        0.769533999        0.044770056
6.883333333        0.807748957        0.002235684
7.883333333        0.792673928        0

展开全文
weixin_36242516 2021-04-18 06:23:57
• 5星
1KB weixin_42683394 2021-09-29 03:36:33
• Lu 微分方程参数优化（拟合） 微分方程 优化

Lu 微分方程参数优化（拟合） 例子1数学模型： dy/dt=k*y*z+0.095*b*z dz/dt=-b*z-0.222*z 实验数据（ti; yi）： ti yi0.012.3291015630.683.8512927831.14.500939361.636.7491722472.079.1120628722.679....
欢迎访问 Lu程序设计
Lu 微分方程参数优化（拟合）
例子1 数学模型：
dy/dt=k*y*z+0.095*b*z      dz/dt=-b*z-0.222*z
实验数据（ti; yi）：

ti        yi0.01    2.3291015630.68    3.8512927831.1     4.500939361.63    6.7491722472.07    9.1120628722.67    9.6916573663.09    11.169280133.64    10.914518234.65    16.442792045.1     18.296158855.58    21.639890256.11    25.786114216.63    26.342826337.06    26.505812877.62    27.639518238.66    35.027576269.04    35.55620359.63    36.10393415

已知z在t=0.01时的初值为5000，求k和b的最优值？
Lu代码1：

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

结果（前面的数为最优参数，最后一个数是极小值。下同）：忽略Lu警告信息

1.408312606232919e-004 -1.465812213944895e-004 24.29124283367381

Lu代码2：优化（拟合）并绘图

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

结果：

例子2 数学模型：需要拟合的微分方程组为:
df1/dt=a1*(x-f2)^a2,         df2/dt=a3*df1/dt-a4*f2^a5
其中：x=0.3。拟合参数为a1,a2,a3,a4,a5。试验获取的数据为（t,f1）；t=0时，f2=0.15。

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

Lu代码1：

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

结果（没有多次运行，或许还有更优解）：

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

Lu代码2：优化（拟合）并绘图

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

结果：

例子3 数学模型：
有关于y1,y2,y3,y4的微分方程：k1 = 3.5*10^4; k2 = 1.0*10^3; k3 = a*ka+b*kb+c*kc; r1 = k1*y1*y2; r2 = k2*y1*y3; r3 = k3*y1*y4; dy1/dt = -r1-r2-r3; dy2/dt = -r1; dy3/dt = -r2+r1; dy4/dt = -r3+r2; 其中ka，kb，kc为拟合参数；常数a、b、c为已知数据；k1,k2,k3,r1,r2,r3为中间变量；t=0时，y1=30.0*10^-6; y2=10.2*10^-6; y3=4.1*10^-6; y4=6.7。 有4组数据：
当[a b c]=[100 5 20]时 t=[0, 10, 20, 30, 40] y1=[30.0*10^-6, 17.5*10^-6, 15*10^-6, 14*10^-6, 13.5*10^-6 ] 当[a b c]=[120 15 25]时 t=[0 10 20 30 40] y1=[30.0*10^-6, 12.5*10^-6, 10*10^-6, 8*10^-6, 7*10^-6 ]
当[a b c]=[130 25 30]时 t=[0 10 20 30 40] y1=[30.0*10^-6, 11.5*10^-6, 9*10^-6, 7.5*10^-6, 6*10^-6 ] 当[a b c]=[140 30 35]时 t=[0 10 20 30 40] y1=[30.0*10^-6, 10.5*10^-6, 8.5*10^-6, 6.5*10^-6, 5*10^-6 ]
求ka，kb，kc？
分析：观察式子 k3 = a*ka+b*kb+c*kc ，当[a b c]实验数据少于3组时，将有无穷组最优的ka，kb，kc。现在有4组数据，故只有一组最优解。
Lu代码：

!!!using["luopt","math"]; f(t,y1,y2,y3,y4,z1,z2,z3,z4,pp : k1,k2,k3,r1,r2,r3 : a,b,c,ka,kb,kc)= {     k1 = 3.5e4,     k2 = 1.0e3,     k3 = a*ka+b*kb+c*kc,     r1 = k1*y1*y2,     r2 = k2*y1*y3,     r3 = k3*y1*y4,     z1 = -r1-r2-r3,     z2 = -r1,     z3 = -r2+r1,     z4 = -r3+r2,     0 }; 目标函数(kka,kkb,kkc : i,s,yy : tArray,y11Array,y12Array,y13Array,y14Array,max,a,b,c,ka,kb,kc)= {     ka=kka, kb=kkb, kc=kkc,     s=0,     a=100, b=5, c=20, //[a b c]=[100 5 20]    yy=gsl_ode[@f,nil,0.0,tArray,ra1(30.0e-6, 10.2e-6, 4.1e-6, 6.7), 1e-6, 1e-6, 2, 1e-6,100],     i=0, while{++i<max, s=s+[yy(i,1)-y11Array(i)]^2},     a=120, b=15, c=25, //[a b c]=[120, 15, 25]    yy=gsl_ode[@f,nil,0.0,tArray,ra1(30.0e-6, 10.2e-6, 4.1e-6, 6.7), 1e-6, 1e-6, 2, 1e-6,100],     i=0, while{++i<max, s=s+[yy(i,1)-y12Array(i)]^2},     a=130, b=25, c=30, //[a b c]=[130, 25, 30]    yy=gsl_ode[@f,nil,0.0,tArray,ra1(30.0e-6, 10.2e-6, 4.1e-6, 6.7), 1e-6, 1e-6, 2, 1e-6,100],     i=0, while{++i<max, s=s+[yy(i,1)-y13Array(i)]^2},     a=140, b=30, c=35, //[a b c]=[140, 30, 35]    yy=gsl_ode[@f,nil,0.0,tArray,ra1(30.0e-6, 10.2e-6, 4.1e-6, 6.7), 1e-6, 1e-6, 2, 1e-6,100],     i=0, while{++i<max, s=s+[yy(i,1)-y14Array(i)]^2},     sqrt[s/(max*4-4)] }; main(::tArray,y11Array,y12Array,y13Array,y14Array,max)= {     max=5,     tArray=ra1[0, 10, 20, 30, 40],     y11Array=ra1[30.0e-6, 17.5e-6, 15e-6, 14e-6, 13.5e-6],     y12Array=ra1[30.0e-6, 12.5e-6, 10e-6, 8e-6, 7e-6],     y13Array=ra1[30.0e-6, 11.5e-6, 9e-6, 7.5e-6, 6e-6],     y14Array=ra1[30.0e-6, 10.5e-6, 8.5e-6, 6.5e-6, 5e-6],     Opt1[@目标函数, optrange,-1e5,1e5,-1e5,1e5,-1e5,1e5] //Opt1函数优化    //Opt[@目标函数, optrange,-1e5,1e5,-1e5,1e5,-1e5,1e5] //也可以使用此语句};

结果（ka,kb,kc，目标函数终值）：

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

验证（t，实验值，拟合值）：
当[a b c]=[100 5 20]时
10. 1.75e-005 1.73918e-005 20. 1.5e-005  1.54895e-005 30. 1.4e-005  1.40236e-005 40. 1.35e-005 1.28548e-005
当[a b c]=[120 15 25]时
10. 1.25e-005 1.45484e-005 20. 1.e-005   1.09127e-005 30. 8.e-006   8.29305e-006 40. 7.e-006   6.35421e-006
当[a b c]=[130 25 30]时
10. 1.15e-005 1.29864e-005 20. 9.e-006   8.73647e-006 30. 7.5e-006  5.94617e-006 40. 6.e-006   4.07251e-006
当[a b c]=[140 30 35]时
10. 1.05e-005 1.30756e-005 20. 8.5e-006  8.85422e-006 30. 6.5e-006  6.06631e-006 40. 5.e-006   4.18281e-006
例子4 数学模型：
微分方程组为：dn1/dt = a*n1*n2-b*n1; dn2/dt = c*n1*n2; 其中a，b，c为拟合参数；t=0时，n1=n0; n2=n0。
数据：
t    1    2     3     4     5     6
n1   1    3.5   2.8   2.4   1.5   0.9
分析：由于n0未知，故共有4个拟合参数：n0,a,b,c。
Lu代码：

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

结果（n0,a,b,c，误差）：

4.842109347483411e-002 74.58161099198458 0.3525934220396638 -0.7028769944300635 0.1446429325392884

验证（t，实验值，拟合值）：

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

例子5 动力学方程参数估计：
微分方程组为：dC/dt=-k0*exp{-Ea/[R*(273.15+T)]}*(A^m)*(B^n)*(C^o) 单次实验T,A,B都是固定的，而C随时间变化并测量出来，所以T,A,B可以说是可变常数，R=8.314，拟合参数为k0,Ea,m,n,o，实验数据如下： T(℃)  A  B(MPa)  t/min  C 260 0.05 4 0   3759.7 260 0.05 4 10  2772.45 260 0.05 4 20  2100 260 0.05 4 30  1679.03 260 0.05 4 40  1356.12 260 0.05 4 50  1122.72 260 0.05 4 60  818.6 260 0.05 4 70  643.28 260 0.05 4 80  485.26 260 0.05 4 90  361.16 260 0.05 4 100 290.17 260 0.1  4 0   3759.7 260 0.1  4 10  2450 260 0.1  4 20  1803.37 260 0.1  4 30  1246.87 260 0.1  4 40  895 260 0.1  4 50  496.14 260 0.1  4 60  410.13 260 0.1  4 70  179.72 260 0.1  4 80  142.38 260 0.1  4 90  94.73 260 0.1  4 100 93.51 260 0.15 4 0   3759.7 260 0.15 4 10  2300 260 0.15 4 20  1424.68 260 0.15 4 30  884.85 260 0.15 4 40  546.36 260 0.15 4 50  336.83 260 0.15 4 60  203.33 260 0.15 4 70  117.97 260 0.15 4 80  66.95 260 0.15 4 90  54.49 260 0.15 4 100 49.41 260 0.2  4 0   3759.7 260 0.2  4 10  2150 260 0.2  4 20  1200 260 0.2  4 30  650 260 0.2  4 40  380 260 0.2  4 50  230 260 0.2  4 60  130 260 0.2  4 70  80 260 0.2  4 80  65 260 0.2  4 90  50 260 0.2  4 100 45 200 0.1  4 0   3759.7 200 0.1  4 10  2893.56 200 0.1  4 20  2037.91 200 0.1  4 30  1319.34 200 0.1  4 40  1005.13 200 0.1  4 50  773.12 200 0.1  4 60  562.77 200 0.1  4 70  384.27 200 0.1  4 80  267.1 200 0.1  4 90  171.38 200 0.1  4 100 157.51 220 0.1  4 0   3759.7 220 0.1  4 10  2600 220 0.1  4 20  1691.77 220 0.1  4 30  1126.16 220 0.1  4 40  902.03 220 0.1  4 50  628.67 220 0.1  4 60  371.19 220 0.1  4 70  240.16 220 0.1  4 80  139.84 220 0.1  4 90  116.87 220 0.1  4 100 83.8 240 0.1  4 0   3759.7 240 0.1  4 10  2350.51 240 0.1  4 20  1450 240 0.1  4 30  1000 240 0.1  4 40  708.77 240 0.1  4 50  375.68 240 0.1  4 60  271.72 240 0.1  4 70  175.88 240 0.1  4 80  123.3 240 0.1  4 90  90 240 0.1  4 100 66.69 260 0.1  4 0   3759.7 260 0.1  4 10  2080 260 0.1  4 20  1200 260 0.1  4 30  884.85 260 0.1  4 40  546.36 260 0.1  4 50  336.83 260 0.1  4 60  203.33 260 0.1  4 70  117.97 260 0.1  4 80  66.95 260 0.1  4 90  54.49 260 0.1  4 100 49.41 280 0.1  4 0   3759.7 280 0.1  4 10  1935.79 280 0.1  4 20  1007.44 280 0.1  4 30  758.08 280 0.1  4 40  500 280 0.1  4 50  315.13 280 0.1  4 60  190.44 280 0.1  4 70  100 280 0.1  4 80  50 280 0.1  4 90  40 280 0.1  4 100 30 260 0.1  4 0   3759.7 260 0.1  4 10  2080 260 0.1  4 20  1200 260 0.1  4 30  749.04 260 0.1  4 40  500 260 0.1  4 50  280 260 0.1  4 60  203.33 260 0.1  4 70  117.97 260 0.1  4 80  66.95 260 0.1  4 90  54.49 260 0.1  4 100 49.41 260 0.1  3.5 0   3759.7 260 0.1  3.5 10  2244.58 260 0.1  3.5 20  1484.05 260 0.1  3.5 30  749.04 260 0.1  3.5 40  547.23 260 0.1  3.5 50  331.18 260 0.1  3.5 60  267.75 260 0.1  3.5 70  198.36 260 0.1  3.5 80  130.51 260 0.1  3.5 90  100.47 260 0.1  3.5 100 93.91 260 0.1  3 0   3759.7 260 0.1  3 10  2483.29 260 0.1  3 20  1669.27 260 0.1  3 30  1162.42 260 0.1  3 40  827.86 260 0.1  3 50  517.85 260 0.1  3 60  357.71 260 0.1  3 70  274.21 260 0.1  3 80  185.53 260 0.1  3 90  129.69 260 0.1  3 100 88.35 260 0.1  2.5 0   3759.7 260 0.1  2.5 10  2527.4 260 0.1  2.5 20  1857.92 260 0.1  2.5 30  1219.77 260 0.1  2.5 40  850.2 260 0.1  2.5 50  628.67 260 0.1  2.5 60  515.01 260 0.1  2.5 70  359.75 260 0.1  2.5 80  265.34 260 0.1  2.5 90  172.36 260 0.1  2.5 100 149.93 260 0.1  2 0   3759.7 260 0.1  2 10  2722.63 260 0.1  2 20  2089.13 260 0.1  2 30  1485.29 260 0.1  2 40  1062.05 260 0.1  2 50  831.48 260 0.1  2 60  766.41 260 0.1  2 70  547.1 260 0.1  2 80  461.41 260 0.1  2 90  347.28 260 0.1  2 100 244.06
分析：一共有14组实验数据，另外注意到时间序列都一样[0,10,20,30,40,50,60,70,80,90,100]，这样可以简化代码。
Lu代码：

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

结果（k0,Ea,m,n,o,误差）：

0.5911010179120556 11993.11965328726 0.5943806176877708 0.5842721259729485 1.094210900092406 2087691.254199586

例子6 求解微分方程组的参数：
求解K1，K2，a, m 方程组为： CB'=-K1*CB^a*(k1/k2*CB^a)^m CH'=K2*(K1/K2*CB^a)^m t   CB     CH 60  10.9237 0 90  10.7462 0 120 10.4357 0.03778 135 10.1695 0.0432 150 9.7481 0.1203 165 9.2346 0.21242 180 8.6613 0.34579 195 8.0058 0.56225 210 7.2423 0.83487 225 6.4188 1.12793 240 5.5353 1.38079 255 4.5768 1.869 270 4.0146 2.5 285 3.5703 3.01 300 3.1158 3.54452 330 2.4438 4.312 360 1.9878 4.70402 390 1.6668 4.8548 代码：

!!!using["luopt","math"]; //使用命名空间f(t,CB,CH,dCB,dCH,pp::K1,K2,a,m)= {     dCB=-K1*CB^a*(K1/K2*CB^a)^m,     dCH=K2*(K1/K2*CB^a)^m,     0 }; 目标函数(_K1,_K2,_a,_m : i,s,tyz : tyArray,tA,max,K1,K2,a,m)= {     K1=_K1, K2=_K2, a=_a, m=_m, //传递优化变量，函数f中要用到K1,K2,a,m    tyz=gsl_ode[@f,nil,0.0,tA,ra1(10.9237,0), 1e-6, 1e-6, 2, 1e-6,100],     i=0, s=0, while{++i<max, s=s+[tyz(i,1)-tyArray(i,1)]^2+[tyz(i,2)-tyArray(i,2)]^2},     s }; main(::tyArray,tA,max)= {     tyArray=matrix{ //存放实验数据"60 10.9237 0 90 10.7462 0 120 10.4357 0.03778 135 10.1695 0.0432 150 9.7481 0.1203 165 9.2346 0.21242 180 8.6613 0.34579 195 8.0058 0.56225 210 7.2423 0.83487 225 6.4188 1.12793 240 5.5353 1.38079 255 4.5768 1.869 270 4.0146 2.5 285 3.5703 3.01 300 3.1158 3.54452 330 2.4438 4.312 360 1.9878 4.70402 390 1.6668 4.8548"     },     len[tyArray,0,&max], tA=tyArray(all:0), //用len函数取矩阵的行数，tA取矩阵的列    Opt1[@目标函数,optwaysimdeep, optwayconfra] //Opt1函数全局优化    //Opt[@目标函数, optwaysimdeep, optwayconfra] //也可以使用此语句};

结果（K1,K2,a,m，目标函数值）： 2.007221403771439e-002 3.477931047132456e-002 0.7598653691173753 -1.205432420176749 18.90800296978626
绘图：

!!!using("win","math"); f(t,CB,CH,dCB,dCH,pp :: K1,K2,a,m)= {     dCB=-K1*CB^a*(K1/K2*CB^a)^m,     dCH=K2*(K1/K2*CB^a)^m,     0 }; init(x,tx::K1,K2,a,m)=     x=matrix[ " 60 10.9237 0 90 10.7462 0 120 10.4357 0.03778 135 10.1695 0.0432 150 9.7481 0.1203 165 9.2346 0.21242 180 8.6613 0.34579 195 8.0058 0.56225 210 7.2423 0.83487 225 6.4188 1.12793 240 5.5353 1.38079 255 4.5768 1.869 270 4.0146 2.5 285 3.5703 3.01 300 3.1158 3.54452 330 2.4438 4.312 360 1.9878 4.70402 390 1.6668 4.8548 "     ],     tx=x(all:0),     K1=2.007221403771439e-002, K2=3.477931047132456e-002, a=0.7598653691173753, m=-1.205432420176749,     luShareX2(x, gsl_ode[@f,nil,0.0,tx,ra1(10.9237,0), 1e-6, 1e-6, 2, 1e-6,100]); ChartWnd[@init];   //显示窗口并初始化

分析：本例所求解的应是最优解，但图形显示数据与曲线不一致，故怀疑数据与模型不匹配。
例子7 常微分方程组的参数拟合：
根据实验数据拟合出细胞对水的渗透参数lp和对保护剂的渗透系数ps。 微分方程组是这样的： da/dtt=lp*r*temp*(0.289*（v0-vb)/a-0.289+b/a-4.809)*4.83579*c^0.6667 db/dt=ps*(4.809-b/a)*4.83579*c^0.6667 dc/dt=da/dt+vcpa*db/dt 其中，a,b,c是因变量，t是自变量，其他的r,temp,v0,vb,vcpa均为常量,如下：r=0.08206,temp=298,v0=0.9402*10^(-12),vb=0.2642*10^(-12),vcpa=63.75*10^(-6) 实验只能测出a和c值，b值无法测出。
实验数据：
t  a      c 0  0.719  1 21 0.2503 0.5314 29 0.1412 0.4222 38 0.1893 0.4703 46 0.2268 0.5078 54 0.2552 0.5362 71 0.2818 0.5628 87 0.3122 0.5932 104 0.3434 0.6244 120 0.3724 0.6534 137 0.3957 0.6767 153 0.4068 0.6878 170 0.428 0.709 195 0.4514 0.7324 228 0.4744 0.7554 261 0.4977 0.7787 294 0.5154 0.7964 335 0.5344 0.8154 377 0.5424 0.8234 426 0.5709 0.8519 476 0.587 0.868 526 0.6022 0.8832
分析：因b值无法测出，故追加b的初值_b为优化参数。下面代码中本来没有[tyz(i,0)-tac(i,0)]^2，但发现所求解存在异常，故添加此代码。
代码：

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

结果（lp,ps,_b,目标函数值）：
2.473962346617309e-003 -3.626422822558594e-003 1.308637154884405 0.1199412216587657
绘图：

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

例子8 常微分方程组的参数拟合：
dz/dt=1/tv*(((lp*t)^av)*(z^(-av-1))-((lp*t)^(-2*av))*(z^(2*av-1))) ds/dt=(ue*(ae-1)*(lp*t)^(ae-2)+ue*(2*ae+1)*(lp*t)^(-2*ae-2)+uv*(av-1)*(lp*t)^(av-2)*z^(-av)+uv*(2*av+1)*(lp*t)^(-2*av-2)*z^(2*av))*lp-(uv*av*(lp*t)^(av-1)*z^(-av-1)+2*uv*av*(lp*t)^(-2*av- 1)*z^(2*av-1))*dz/dt ue,ae,uv,av,tv为待拟合参数 z为中间变量 lp=0.00354 t-s为实验数据： t=282.5,361.0,439.4,517.9,596.4,674.8,753.3,831.8,910.2,988.7; s=0,29099,36228,40626,44290,47686,50964,54181,57363,60521;
分析：因z值无法测出，故追加z的初值_z为优化参数。
代码：

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

结果（ue,ae,uv,av,tv,_z,目标函数值）：
-9114964.001609007 -0.4453565389059405 7335590868.341645 -2.754919559106283e-004 -10485633.9903882 2.654977892066841 244.5944885704022
-249.2203422970038 -2.000969463745669 1822318656671169 -5.117661642056397e-013 -30625.54305602086 2.426847712155262e-010 151.3043733914658（此组解不能确定是否有效）
绘图：

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

版权所有© Lu程序设计 2002-2021，保留所有权利 E-mail: forcal@sina.com  QQ：630715621最近更新: 2021年07月29日
展开全文
forcal 2021-07-30 19:51:44
• 1stOpt计算：Elapsed Time (Hr:Min:Sec:Msec): 00:00:19:284Root of Mean Square Error (RMSE): 1.18378166169196Sum of Square Error: 7.00669511279087Correlation Coef. (R): 0.999185213441573R-Square: 0....

1stOpt计算：
Elapsed Time (Hr:Min:Sec:Msec): 00:00:19:284
Root of Mean Square Error (RMSE): 1.18378166169196
Sum of Square Error: 7.00669511279087
Correlation Coef. (R): 0.999185213441573
R-Square: 0.998371090760282
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
• 美赛整理之带参数的常微分方程拟合问题研究 matlab 算法

shengzimao 2020-10-03 21:34:44
• 非线性微分方程组多参数拟合 测试用例

chimchim12138 2020-06-24 10:15:50
• weixin_30616309 2021-04-18 14:37:41
• weixin_39622217 2021-04-18 14:38:15
• weixin_31889919 2021-04-18 14:38:23
• 6.22MB qq_41934573 2021-08-26 22:42:44
• weixin_31032795 2021-04-18 08:58:08
• 含有未知中间变量同时需要传递其他中间变量的微分方程参数拟合 经验分享

wlfc 2021-08-11 15:15:10
• weixin_36337403 2021-04-21 04:28:54
• weixin_42309722 2021-04-20 06:53:28
• 初值未知且需要传递中间变量的微分方程参数拟合 经验分享

wlfc 2021-08-20 20:06:16
• weixin_33056525 2021-04-21 04:29:52
• weixin_35319894 2021-04-18 08:58:04
• weixin_35695434 2020-12-23 12:25:09
• weixin_35141129 2021-04-22 11:37:30
• weixin_28837357 2021-05-08 15:47:22
• 带参数传递的微分方程拟合 经验分享

wlfc 2021-08-10 10:52:42
• weixin_33277564 2021-04-23 20:30:15
• weixin_34791959 2021-04-21 15:37:00
• weixin_33013787 2021-04-21 15:36:58
• 隐函数式微分方程参数拟合 matlab 动态规划 经验分享

wlfc 2021-08-09 15:25:10
• weixin_39851809 2021-04-20 08:45:41
• 二阶微分方程拟合数据得到方程参数 c++ c语言

weixin_38425472 2020-06-19 12:28:29

...