1stopt高版本代跑 (微分方程與代數(shù)方程參數(shù)擬合)
7參數(shù)擬合問題,含有微分方程與代數(shù)方程,請求高版本代跑一下程序。
Parameter fac1=[1e-5,],fac2=[1e-5,],fac3=[1e-5,],fac4=[1e-5,],fac5=[,-1e-5],fac6=[,-1e-5],fac7=[1e-5,];
Variable t,lamb,lambdif,pb;
ODEFunction lambv'=if(lambdif<=0,(lamb/3/fac7)*(fac1*((lamb/lambv)^(fac2-1)-(lambv/lamb)^(0.5*fac2+1))+fac3*((lamb/lambv)^(fac4-1)-(lambv/lamb)^(0.5*fac4+1))+fac5*((lamb/lambv)^(fac6-1)-(lambv/lamb)^(0.5*fac6+1))),(lamb/(0.001*3*fac7))*(fac1*((lamb/lambv)^(fac2-1)-(lambv/lamb)^(0.5*fac2+1))+fac3*((lamb/lambv)^(fac4-1)-(lambv/lamb)^(0.5*fac4+1))+fac5*((lamb/lambv)^(fac6-1)-(lambv/lamb)^(0.5*fac6+1))));
Function pb=fac1*(lamb^(fac2-1)*lambv^(-fac2)-lamb^(-0.5*fac2-1)*lambv^(0.5*fac2))+fac3*(lamb^(fac4-1)*lambv^(-fac4)-lamb^(-0.5*fac4-1)*lambv^(0.5*fac4))+fac5*(lamb^(fac6-1)*lambv^(-fac6)-lamb^(-0.5*fac6-1)*lambv^(0.5*fac6));
Data;
//0.1Hz
0.00 0.960000 0.012566 0.057751
0.20 0.962433 0.012467 0.063943
0.40 0.964893 0.012172 0.063200
0.60 0.967303 0.011684 0.063611
0.80 0.969577 0.011012 0.065567
1.00 0.971710 0.010166 0.065701
1.20 0.973627 0.009160 0.065639
1.40 0.975377 0.008010 0.064222
1.60 0.976853 0.006733 0.062565
1.80 0.978063 0.005351 0.061493
2.00 0.979000 0.003883 0.058928
2.20 0.979637 0.002355 0.058219
2.40 0.979960 0.000789 0.055017
2.60 0.979960 -0.000789 0.052617
2.80 0.979660 -0.002355 0.047199
3.00 0.979037 -0.003883 0.042338
3.20 0.978137 -0.005351 0.032877
3.40 0.976940 -0.006733 0.025563
3.60 0.975463 -0.008010 0.018672
3.80 0.973743 -0.009160 0.010299
4.00 0.971830 -0.010166 0.004338
4.20 0.969690 -0.011012 -0.004050
4.40 0.967423 -0.011684 -0.010569
4.60 0.965040 -0.012172 -0.019444
4.80 0.962587 -0.012467 -0.026517
5.00 0.960073 -0.012566 -0.031025
5.20 0.957553 -0.012467 -0.039910
5.40 0.955093 -0.012172 -0.044413
5.60 0.952710 -0.011684 -0.050038
5.80 0.950440 -0.011012 -0.053758
6.00 0.948293 -0.010166 -0.059845
6.20 0.946357 -0.009160 -0.062465
6.40 0.944640 -0.008010 -0.063932
6.60 0.943153 -0.006733 -0.064208
6.80 0.941933 -0.005351 -0.066548
7.00 0.940993 -0.003883 -0.063021
7.20 0.940363 -0.002355 -0.059452
7.40 0.940043 -0.000789 -0.055989
7.60 0.940033 0.000789 -0.051224
7.80 0.940323 0.002355 -0.039363
8.00 0.940947 0.003883 -0.025688
8.20 0.941860 0.005351 -0.009293
8.40 0.943073 0.006733 0.005272
8.60 0.944540 0.008010 0.016966
8.80 0.946260 0.009160 0.028023
9.00 0.948193 0.010166 0.034484
9.20 0.950300 0.011012 0.042393
9.40 0.952567 0.011684 0.048518
9.60 0.954957 0.012172 0.050146
9.80 0.957423 0.012467 0.055289
10.00 0.960000 0.012566 0.057751
Data;
//0.5Hz
0.00 0.960000 0.062831853 0.057592
0.04 0.962137 0.062336405 0.062721
0.08 0.964590 0.060857875 0.067115
0.12 0.967043 0.05841958 0.069002
0.16 0.969327 0.055059973 0.070735
0.20 0.971487 0.050832037 0.072267
0.24 0.973477 0.04580245 0.067785
0.28 0.975240 0.04005053 0.064896
0.32 0.976757 0.03366699 0.065203
0.36 0.978013 0.026752502 0.061654
0.40 0.978987 0.01941611 0.060950
0.44 0.979673 0.011773515 0.055396
0.48 0.980030 0.003945245 0.052443
0.52 0.980090 -0.003945245 0.048458
0.56 0.979813 -0.011773515 0.042921
0.60 0.979240 -0.01941611 0.029484
0.64 0.978350 -0.026752502 0.019637
0.68 0.977177 -0.03366699 0.005958
0.72 0.975753 -0.04005053 -0.008437
0.76 0.974063 -0.04580245 -0.021190
0.80 0.972140 -0.050832037 -0.033505
0.84 0.970023 -0.055059973 -0.046834
0.88 0.967757 -0.05841958 -0.059828
0.92 0.965407 -0.060857875 -0.073131
0.96 0.962903 -0.062336405 -0.084756
1.00 0.960407 -0.062831853 -0.097015
1.04 0.957857 -0.062336405 -0.108150
1.08 0.955407 -0.060857875 -0.116848
1.12 0.952977 -0.05841958 -0.126935
1.16 0.950650 -0.055059973 -0.132818
1.20 0.948500 -0.050832037 -0.140228
1.24 0.946553 -0.04580245 -0.143198
1.28 0.944793 -0.04005053 -0.147031
1.32 0.943243 -0.03366699 -0.148441
1.36 0.941980 -0.026752502 -0.146859
1.40 0.941007 -0.01941611 -0.143694
1.44 0.940333 -0.011773515 -0.137398
1.48 0.939983 -0.003945245 -0.130185
1.52 0.939920 0.003945245 -0.118369
1.56 0.940163 0.011773515 -0.100830
1.60 0.940763 0.01941611 -0.078713
1.64 0.941643 0.026752502 -0.054448
1.68 0.942843 0.03366699 -0.031340
1.72 0.944277 0.04005053 -0.009725
1.76 0.945943 0.04580245 0.004953
1.80 0.947883 0.050832037 0.017859
1.84 0.949980 0.055059973 0.027375
1.88 0.952207 0.05841958 0.037944
1.92 0.954630 0.060857875 0.044558
1.96 0.957083 0.062336405 0.050863
2.00 0.960000 0.062831853 0.057592
Data;
//1Hz
0.00 0.959110 0.125664 0.101944
0.02 0.961607 0.124673 0.101015
0.04 0.964177 0.121716 0.092652
0.06 0.966597 0.116839 0.090688
0.08 0.968997 0.110120 0.087186
0.10 0.971203 0.101664 0.081292
0.12 0.973220 0.091605 0.074067
0.14 0.975063 0.080101 0.072352
0.16 0.976643 0.067334 0.064963
0.18 0.977983 0.053505 0.062168
0.20 0.979003 0.038832 0.059446
0.22 0.979740 0.023547 0.056947
0.24 0.980190 0.007890 0.053218
0.26 0.980303 -0.007890 0.051380
0.28 0.980107 -0.023547 0.049209
0.30 0.979583 -0.038832 0.050486
0.32 0.978767 -0.053505 0.045529
0.34 0.977637 -0.067334 0.042285
0.36 0.976253 -0.080101 0.040140
0.38 0.974577 -0.091605 0.040742
0.40 0.972660 -0.101664 0.028174
0.42 0.970597 -0.110120 0.015944
0.44 0.968337 -0.116839 -0.000255
0.46 0.965933 -0.121716 -0.015326
0.48 0.963463 -0.124673 -0.032873
0.50 0.960933 -0.125664 -0.051744
0.52 0.958353 -0.124673 -0.068495
0.54 0.955860 -0.121716 -0.084269
0.56 0.953367 -0.116839 -0.097367
0.58 0.951017 -0.110120 -0.107719
0.60 0.948813 -0.101664 -0.117079
0.62 0.946750 -0.091605 -0.120800
0.64 0.944930 -0.080101 -0.128407
0.66 0.943330 -0.067334 -0.127666
0.68 0.942017 -0.053505 -0.125147
0.70 0.940943 -0.038832 -0.120878
0.72 0.940243 -0.023547 -0.112602
0.74 0.939793 -0.007890 -0.102641
0.76 0.939683 0.007890 -0.090411
0.78 0.939893 0.023547 -0.066376
0.80 0.940393 0.038832 -0.039621
0.82 0.941243 0.053505 -0.011216
0.84 0.942383 0.067334 0.013463
0.86 0.943800 0.080101 0.032020
0.88 0.945480 0.091605 0.047094
0.90 0.947350 0.101664 0.060936
0.92 0.949430 0.110120 0.076479
0.94 0.951703 0.116839 0.086050
0.96 0.954087 0.121716 0.093122
0.98 0.956563 0.124673 0.102631
1.00 0.959110 0.125664 0.101944
Data;
//2Hz
0.00 0.957683 0.251327 0.103554
0.02 0.962873 0.243431 0.099815
0.04 0.967843 0.220240 0.089030
0.06 0.972333 0.183210 0.079861
0.08 0.976073 0.134668 0.066694
0.10 0.978863 0.077664 0.058112
0.12 0.980390 0.015781 0.053472
0.14 0.980647 -0.047094 0.048643
0.16 0.979643 -0.107010 0.043628
0.18 0.977403 -0.160202 0.036629
0.20 0.974053 -0.203328 0.016363
0.22 0.969813 -0.233678 -0.030502
0.24 0.964933 -0.249346 -0.082703
0.26 0.959757 -0.249346 -0.129067
0.28 0.954600 -0.233678 -0.165265
0.30 0.949830 -0.203328 -0.192737
0.32 0.945690 -0.160202 -0.204648
0.34 0.942393 -0.107010 -0.208398
0.36 0.940157 -0.047094 -0.197350
0.38 0.939237 0.015781 -0.169471
0.40 0.939640 0.077664 -0.110091
0.42 0.941337 0.134668 -0.036682
0.44 0.944197 0.183210 0.015458
0.46 0.948013 0.220240 0.052942
0.48 0.952623 0.243431 0.084528
0.50 0.957683 0.251327 0.103554
Data;
//5Hz
0.00 0.953477 0.628319 0.135643
0.02 0.966487 0.508320 0.098394
0.04 0.976950 0.194161 0.063313
0.06 0.981010 -0.194161 0.047463
0.08 0.977333 -0.508320 0.039043
0.10 0.966537 -0.628319 -0.039222
0.12 0.953520 -0.508320 -0.220273
0.14 0.943087 -0.194161 -0.271013
0.16 0.938930 0.194161 -0.210556
0.18 0.942740 0.508320 0.020174
0.20 0.953477 0.628319 0.135643
返回小木蟲查看更多
京公網(wǎng)安備 11010802022153號
是用5組數(shù)據(jù)擬合出1組參數(shù)。不知如何編輯帖子,只能在回復(fù)中說明了。
計算很費(fèi)時間,也不知道樓主設(shè)定的參數(shù)范圍是否合理,下面結(jié)果參考下,有時間自己計算或許有更好的結(jié)果:
計算用時(時:分:秒:微秒): 00:33:45:772
均方差(RMSE): 0.0495089955802209
殘差平方和(SSR): 0.453461019022031
相關(guān)系數(shù)(R): 0.909905962406545
相關(guān)系數(shù)之平方(R^2): 0.827928860422981
修正R平方(Adj. R^2): 0.816833445348956
確定系數(shù)(DC): 0.494630751059566
F統(tǒng)計(F-Statistic): -0.0452805144567521
參數(shù) 最佳估算
-------------------- -------------
fac1 1.01254846648036E-5
fac2 1.07296406274351E-5
fac3 0.286048998561719
fac4 8.89086313048348
fac5 -2.08243476928144E-5
fac6 -2.45497878879962E-5
fac7 107.530424216471
lambv 初值 (數(shù)據(jù)文件-1) 0.451633585283195
lambv 初值 (數(shù)據(jù)文件-2) 0.92664566345588
lambv 初值 (數(shù)據(jù)文件-3) 0.923788301682772
lambv 初值 (數(shù)據(jù)文件-4) 0.955286830686244
lambv 初值 (數(shù)據(jù)文件-5) 0.972523739080962
====== 結(jié)果輸出 ======
文件: 數(shù)據(jù)文件-1
No t 目標(biāo) pb 計算 pb
1 0.2 0.063943 0.0627235820485261
2 0.4 0.0632 0.00646253000540754
3 0.6 0.063611 0.00156611706497645
4 0.8 0.065567 0.00105921527521241
5 1 0.065701 0.00095368444589652
6 1.2 0.065639 0.00085592351193353
7 1.4 0.064222 0.000778839835616696
8 1.6 0.062565 0.000658251865379009
9 1.8 0.061493 0.000541653484650114
10 2 0.058928 0.000423191745924583
11 2.2 0.058219 0.000294168942682817
12 2.4 0.055017 0.000159422910405084
13 2.6 0.052617 0.000136773366862757
14 2.8 0.047199 -0.00105268803396297
15 3 0.042338 -0.00351793879116287
16 3.2 0.032877 -0.00707040873424754
17 3.4 0.025563 -0.0117828674508465
18 3.6 0.018672 -0.0175802986889611
19 3.8 0.010299 -0.0243093426264643
20 4 0.004338 -0.0317669097854501
21 4.2 -0.00405 -0.0400812355343877
22 4.4 -0.010569 -0.0488560952075601
23 4.6 -0.019444 -0.0580469902999115
24 4.8 -0.026517 -0.0674742283975919
25 5 -0.031025 -0.0771040259784608
26 5.2 -0.03991 -0.0867247216012876
27 5.4 -0.044413 -0.0960835875455815
28 5.6 -0.050038 -0.105118939887549
29 5.8 -0.053758 -0.113695200019873
30 6 -0.059845 -0.121777294447135
31 6.2 -0.062465 -0.1290271491769
32 6.4 -0.063932 -0.135416438972139
33 6.6 -0.064208 -0.140904395383305
34 6.8 -0.066548 -0.145348785980609
35 7 -0.063021 -0.148700162618943
36 7.2 -0.059452 -0.150842944454711
37 7.4 -0.055989 -0.151778124582847
38 7.6 -0.051224 -0.01652565610898
39 7.8 -0.039363 -0.00144911289491686
40 8 -0.025688 0.000135247971510671
41 8.2 -0.009293 0.000404872736131531
42 8.4 0.005272 0.000555403032197519
43 8.6 0.016966 0.000674209041159786
44 8.8 0.028023 0.000788419021982872
45 9 0.034484 0.000889201844346569
46 9.2 0.042393 0.000966512277130337
47 9.4 0.048518 0.00103485709553578
48 9.6 0.050146 0.00108639436658552
49 9.8 0.055289 0.00111672335905731
50 10 0.057751 0.00115869277957147
文件: 數(shù)據(jù)文件-2
No t 目標(biāo) pb 計算 pb
1 0.04 0.062721 0.0946147992848782
2 0.08 0.067115 0.0630037940241945
3 0.12 0.069002 0.0448355069469828
4 0.16 0.070735 0.0335492340759896
5 0.2 0.072267 0.0263892663797244
6 0.24 0.067785 0.0215614871861531
7 0.28 0.064896 0.0179897110362559
8 0.32 0.065203 0.0151338228900227
9 0.36 0.061654 0.0126788857000415
10 0.4 0.06095 0.0104175136701096
11 0.44 0.055396 0.0082617212577555
12 0.48 0.052443 0.0060648242576651
13 0.52 0.048458 0.00599444700667593
14 0.56 0.042921 0.00488635992644049
15 0.6 0.029484 0.0025991789799334
16 0.64 0.019637 -0.000947071849688621
17 0.68 0.005958 -0.00561185498777706
18 0.72 -0.008437 -0.0112623864009205
19 0.76 -0.02119 -0.0179522832922644
20 0.8 -0.033505 -0.0255445468958885
21 0.84 -0.046834 -0.033880447627889
22 0.88 -0.059828 -0.0427787065916436
23 0.92 -0.073131 -0.0519828167083645
24 0.96 -0.084756 -0.0617663608696495
25 1 -0.097015 -0.0714961481636833
26 1.04 -0.10815 -0.0814167422624327
27 1.08 -0.116848 -0.090930808088205
28 1.12 -0.126935 -0.100353676097057
29 1.16 -0.132818 -0.109365876637923
30 1.2 -0.140228 -0.117682418699296
31 1.24 -0.143198 -0.125205345829746
32 1.28 -0.147031 -0.131998385368345
33 1.32 -0.148441 -0.137973459362134
34 1.36 -0.146859 -0.142831078281293
35 1.4 -0.143694 -0.146559670383817
36 1.44 -0.137398 -0.149123747180655
37 1.48 -0.130185 -0.150425655833094
38 1.52 -0.118369 -0.0993213281073989
39 1.56 -0.10083 -0.0635879705228149
40 1.6 -0.078713 -0.0389312062381749
41 1.64 -0.054448 -0.0221528861378749
42 1.68 -0.03134 -0.0105377068227763
43 1.72 -0.009725 -0.00254866233203608
44 1.76 0.004953 0.00309949161446674
45 1.8 0.017859 0.00737215909938563
46 1.84 0.027375 0.0104411466004864
47 1.88 0.037944 0.0125000759856942
48 1.92 0.044558 0.0144615073672508
49 1.96 0.050863 0.0155256109286916
50 2 0.057592 0.0176108786266836
文件: 數(shù)據(jù)文件-3
No t 目標(biāo) pb 計算 pb
1 0.02 0.101015 0.132654995221037
2 0.04 0.092652 0.110298802223814
3 0.06 0.090688 0.0929646295885082
4 0.08 0.087186 0.0799575512153793
5 0.1 0.081292 0.0693909086463457
6 0.12 0.074067 0.0606559664328367
7 0.14 0.072352 0.05334763588412
8 0.16 0.064963 0.0468230409453347
9 0.18 0.062168 0.040973924273115
10 0.2 0.059446 0.0353791710462901
11 0.22 0.056947 0.0301026095015497
12 0.24 0.053218 0.0250575552329378
13 0.26 0.05138 0.0251674018911858
14 0.28 0.049209 0.0243655627273676
15 0.3 0.050486 0.0222336664685919
16 0.32 0.045529 0.0189214195432357
17 0.34 0.042285 0.014345001696673
18 0.36 0.04014 0.00875354948417075
19 0.38 0.040742 0.00200075355986696
20 0.4 0.028174 -0.00570035751252772
21 0.42 0.015944 -0.0139628075703394
22 0.44 -0.000255 -0.0229866025636025
23 0.46 -0.015326 -0.0325561823056165
24 0.48 -0.032873 -0.0423598255848552
25 0.5 -0.051744 -0.0523745541888915
26 0.52 -0.068495 -0.0625622258501379
27 0.54 -0.084269 -0.0723847749226578
28 0.56 -0.097367 -0.0821895541178059
29 0.58 -0.107719 -0.091417272068946
30 0.6 -0.117079 -0.100060643452097
31 0.62 -0.1208 -0.108142872330871
32 0.64 -0.128407 -0.115266212914614
33 0.66 -0.127666 -0.121523166344025
34 0.68 -0.125147 -0.126652005522404
35 0.7 -0.120878 -0.130841852830688
36 0.72 -0.112602 -0.13356291013483
37 0.74 -0.102641 -0.135302082607774
38 0.76 -0.090411 -0.110340729026219
39 0.78 -0.066376 -0.088396952462021
40 0.8 -0.039621 -0.0694027754728905
41 0.82 -0.011216 -0.0527193662819628
42 0.84 0.013463 -0.0382031240076955
43 0.86 0.03202 -0.0255542156233595
44 0.88 0.047094 -0.0145134464002296
45 0.9 0.060936 -0.00505930729303474
46 0.92 0.076479 0.003160128499167
47 0.94 0.08605 0.0103087461625105
48 0.96 0.093122 0.0163018339469903
49 0.98 0.102631 0.021298446969059
50 1 0.101944 0.0254243992669727
文件: 數(shù)據(jù)文件-4
No t 目標(biāo) pb 計算 pb
1 0.02 0.099815 0.0249938602394822
2 0.04 0.08903 0.0362477257843847
3 0.06 0.079861 0.0434393972397149
4 0.08 0.066694 0.0462246283651297
5 0.1 0.058112 0.0453273028378734
6 0.12 0.053472 0.0403392185226596
7 0.14 0.048643 0.0410587046413378
8 0.16 0.043628 0.0369173345581378
9 0.18 0.036629 0.0277190459885104
10 0.2 0.016363 0.0140372566049685
11 0.22 -0.030502 -0.00316751277523413
12 0.24 -0.082703 -0.0228350590851213
13 0.26 -0.129067 -0.0435623757494601
14 0.28 -0.165265 -0.0641034130327071
15 0.3 -0.192737 -0.0830264770473762
16 0.32 -0.204648 -0.0994059691320063
17 0.34 -0.208398 -0.112427112575224
18 0.36 -0.19735 -0.121244564464212
19 0.38 -0.169471 -0.10135367343319
20 0.4 -0.110091 -0.0803523636093328
21 0.42 -0.036682 -0.0587651925747703
22 0.44 0.015458 -0.0372386917923555
23 0.46 0.052942 -0.0166020545489229
24 0.48 0.084528 0.00238922277897762
25 0.5 0.103554 0.0190173377906893
文件: 數(shù)據(jù)文件-5
No t 目標(biāo) pb 計算 pb
1 0.02 0.098394 -0.0192474838395384
2 0.04 0.063313 0.0183935545163442
3 0.06 0.047463 0.0346054449370822
4 0.08 0.039043 0.0196053654096042
5 0.1 -0.039222 -0.0239116624539082
6 0.12 -0.220273 -0.0756364740518034
7 0.14 -0.271013 -0.116767151536529
8 0.16 -0.210556 -0.108189529369281
9 0.18 0.020174 -0.0745205350857752
10 0.2 0.135643 -0.0242572692209939,
lambv'=if(lambdif<=0,(lamb/3/fac7)*(fac1*((lamb/lambv)^(fac2-1)-(lambv/lamb)^(0.5*fac2+1))+fac3*((lamb/lambv)^(fac4-1)-(lambv/lamb)^(0.5*fac4+1))+fac5*((lamb/lambv)^(fac6-1)-(lambv/lamb)^(0.5*fac6+1))),(lamb/(0.001*3*fac7))*(fac1*((lamb/lambv)^(fac2-1)-(lambv/lamb)^(0.5*fac2+1))+fac3*((lamb/lambv)^(fac4-1)-(lambv/lamb)^(0.5*fac4+1))+fac5*((lamb/lambv)^(fac6-1)-(lambv/lamb)^(0.5*fac6+1))));
這個式子中,等號前的 lambv' 和等號后的 lambv 有什么區(qū)別?
前面是一階微分
用OpenLu求解:
分析:需在微分方程求解中傳遞中間變量lamb,gsl支持的Lu擴(kuò)展數(shù)學(xué)庫中求解微分方程的函數(shù)gsl_ode提供了傳遞中間變量的功能。
進(jìn)一步討論:微分方程求解傳遞中間變量時,按時間t的改變量進(jìn)行插值,可提高求解精度。例如,進(jìn)行線性插值:lamb=lamb0+(lamb1-lamb0)*[(t-t0)/(t1-t0)]。
另外,微分方程初值lambv未知,將其追加為擬合變量,因有5組數(shù)據(jù),故追加5個初值lambv1,lambv2,lambv3,lambv4,lambv5。
Lu代碼:
結(jié)果:
4.559193187270802e-004 0.1385329162525733 2.105891636951903 1.069500926752445 -0.1172195254220325 -2.229532073322357e-004 61.97780044873748 0.943941790941042 0.9369235375692753 0.9210261290423717 0.9321935108690476 0.9479390743267424 0.4826192722906011
討論:模型似乎與數(shù)據(jù)不匹配,擬合效果不好,且lambdif<=0與lambdif>0時公式不一樣,圖形上有明顯轉(zhuǎn)折。另外,數(shù)據(jù)稍有變動,微分方程求解結(jié)果就相差很大,例如,如果中間變量lamb未作線性插值,直接lamb=t_lamb,擬合得最優(yōu)解為:
19.08839920909583 3.165001986998071e-005 0.4606063592598733 5.102787301134413 -2.111635805852794e-005 -3.981718621914052 79.73597242036105 0.9448641243550566 0.9381929974925054 0.9260059193870247 0.9399688967878371 0.9552532765459085 0.4726755802116919
兩個解與變量lamb的用法互換使用時,微分方程竟然無法求解,似乎本例微分方程對數(shù)據(jù)相當(dāng)敏感。
3樓求解的數(shù)據(jù),因不知道1stopt是如何在微分方程求解中使用變量lamb的,故不知如何驗證。如果將數(shù)據(jù)帶入以上Lu代碼求目標(biāo)函數(shù)值,是非常大的。估計反過來也一樣。
故微分方程求解方法不一樣,獲得的最優(yōu)值是不一樣的。對本例來說,即便變量lamb進(jìn)行線性插值或未作插值,也有著天壤之別。
關(guān)于變量lamb進(jìn)行線性插值或未作插值,微分方程計算差別大的問題,有了新發(fā)現(xiàn)。
問題在于gsl_ode中的算法選擇,以上用的是gsl_rkf45(Runge-Kutta-Fehlberg 法),如果改為gsl_rk4 (四階經(jīng)典 Runge-Kutta 法),便沒有問題。
另外,以上Lu代碼的兩個結(jié)果,經(jīng)matlab相同的ode45算法驗證,結(jié)果一致。