1stopt高版本代跑 （微分方程与代数方程参数拟合）
7参数拟合问题，含有微分方程与代数方程，请求高版本代跑一下程序。
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
返回小木虫查看更多
是用5组数据拟合出1组参数。不知如何编辑帖子，只能在回复中说明了。
计算很费时间，也不知道楼主设定的参数范围是否合理，下面结果参考下，有时间自己计算或许有更好的结果：
计算用时(时:分:秒:微秒): 00:33:45:772
均方差(RMSE): 0.0495089955802209
残差平方和(SSR): 0.453461019022031
相关系数(R): 0.909905962406545
相关系数之平方(R^2): 0.827928860422981
修正R平方(Adj. R^2): 0.816833445348956
确定系数(DC): 0.494630751059566
F统计(F-Statistic): -0.0452805144567521
参数 最佳估算
-------------------- -------------
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 初值 (数据文件-1) 0.451633585283195
lambv 初值 (数据文件-2) 0.92664566345588
lambv 初值 (数据文件-3) 0.923788301682772
lambv 初值 (数据文件-4) 0.955286830686244
lambv 初值 (数据文件-5) 0.972523739080962
====== 结果输出 ======
文件: 数据文件-1
No t 目标 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
文件: 数据文件-2
No t 目标 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
文件: 数据文件-3
No t 目标 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
文件: 数据文件-4
No t 目标 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
文件: 数据文件-5
No t 目标 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 有什么区别？
前面是一阶微分
用OpenLu求解：
分析：需在微分方程求解中传递中间变量lamb，gsl支持的Lu扩展数学库中求解微分方程的函数gsl_ode提供了传递中间变量的功能。
进一步讨论：微分方程求解传递中间变量时，按时间t的改变量进行插值，可提高求解精度。例如，进行线性插值：lamb=lamb0+(lamb1-lamb0)*[(t-t0)/(t1-t0)]。
另外，微分方程初值lambv未知，将其追加为拟合变量，因有5组数据，故追加5个初值lambv1,lambv2,lambv3,lambv4,lambv5。
Lu代码：
结果：
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
讨论：模型似乎与数据不匹配，拟合效果不好，且lambdif<=0与lambdif>0时公式不一样，图形上有明显转折。另外，数据稍有变动，微分方程求解结果就相差很大，例如，如果中间变量lamb未作线性插值，直接lamb=t_lamb，拟合得最优解为：
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的用法互换使用时，微分方程竟然无法求解，似乎本例微分方程对数据相当敏感。
3楼求解的数据，因不知道1stopt是如何在微分方程求解中使用变量lamb的，故不知如何验证。如果将数据带入以上Lu代码求目标函数值，是非常大的。估计反过来也一样。
故微分方程求解方法不一样，获得的最优值是不一样的。对本例来说，即便变量lamb进行线性插值或未作插值，也有着天壤之别。
关于变量lamb进行线性插值或未作插值，微分方程计算差别大的问题，有了新发现。
问题在于gsl_ode中的算法选择，以上用的是gsl_rkf45（Runge-Kutta-Fehlberg 法），如果改为gsl_rk4 （四阶经典 Runge-Kutta 法），便没有问题。
另外，以上Lu代码的两个结果，经matlab相同的ode45算法验证，结果一致。