当前位置: 首页 > 计算模拟 >1stopt高版本代跑 (微分方程与代数方程参数拟合)

1stopt高版本代跑 (微分方程与代数方程参数拟合)

作者 ckm0811
来源: 小木虫 300 6 举报帖子
+关注

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 返回小木虫查看更多

今日热帖
  • 精华评论
  • ckm0811

    是用5组数据拟合出1组参数。不知如何编辑帖子,只能在回复中说明了。

  • dingd

    计算很费时间,也不知道楼主设定的参数范围是否合理,下面结果参考下,有时间自己计算或许有更好的结果:

    计算用时(时:分:秒:微秒): 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,

  • wlfc

    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 有什么区别?

  • ckm0811

    引用回帖:
    4楼: Originally posted by wlfc at 2021-08-18 08:07:39
    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* ...

    前面是一阶微分

  • wlfc

    用OpenLu求解:

    分析:需在微分方程求解中传递中间变量lamb,gsl支持的Lu扩展数学库中求解微分方程的函数gsl_ode提供了传递中间变量的功能。
    进一步讨论:微分方程求解传递中间变量时,按时间t的改变量进行插值,可提高求解精度。例如,进行线性插值:lamb=lamb0+(lamb1-lamb0)*[(t-t0)/(t1-t0)]。

    另外,微分方程初值lambv未知,将其追加为拟合变量,因有5组数据,故追加5个初值lambv1,lambv2,lambv3,lambv4,lambv5。

    Lu代码:

    CODE:
    !!!using["luopt","math"]; //使用命名空间
    f(  t,lambv,dlambv, i
        : lamb
        : t_A, t_lamb, t_lambdif, fac1,fac2,fac3,fac4,fac5,fac6,fac7) =
    {
        lamb=t_lamb[i]+(t_lamb[i+1]-t_lamb[i])*[(t-t_A[i])/(t_A[i+1]-t_A[i])], //线性插值
        dlambv=which(t_lambdif(i)<=0.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)))),
        0 //必须返回0
    };
    g(  tArray1, lambv1
        : lamb,lambv,max, tArray, t_pb, i,s,tf
        : t_A, t_lamb, t_lambdif, fac1,fac2,fac3,fac4,fac5,fac6,fac7)=
    {
        tArray=tArray1, len[tArray,0,&max], t_A=tArray(all:0).reshape(), t_lamb=tArray(all:1).reshape(), t_lambdif=tArray(all:2).reshape(), t_pb=tArray(all:3).reshape(), //t_A等取矩阵的列,并转为一维数组
        tf=gsl_ode[@f, nil, nil, t_A, ra1(lambv1), 1e-6, 1e-6, gsl_rkf45, 1e-6,50],
        i=-1, s=0.0, while{++i<max,
            lamb=t_lamb(i), lambv=tf(i,1),
            s=s+[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)) - t_pb(i)]^2.0
        },
        s
    };
    目标函数(_fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7,lambv1,lambv2,lambv3,lambv4,lambv5
        :
        : fac1,fac2,fac3,fac4,fac5,fac6,fac7, tArray1,tArray2,tArray3,tArray4,tArray5)=
    {
        fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7,  //传递优化变量
        g(tArray1, lambv1) +  g(tArray2, lambv2) +  g(tArray3, lambv3) +  g(tArray4, lambv4) +  g(tArray5, lambv5)
    };
    main(:: tArray1,tArray2,tArray3,tArray4,tArray5)=
    {
        tArray1=matrix{ //存放实验数据//t,lamb,lambdif,pb  //0.1Hz
        "0.00         0.960000         0.012566         0.057751
    0.20         0.962433         0.012467         0.063943
    ... ...  //数据省略
    10.00         0.960000         0.012566         0.057751"
        },

        tArray2=matrix{ //存放实验数据//t,lamb,lambdif,pb  //0.5Hz
        "0.00         0.960000         0.062831853        0.057592
    0.04         0.962137         0.062336405        0.062721
    ... ...  //数据省略
    2.00         0.960000         0.062831853        0.057592
    "
        },

        tArray3=matrix{ //存放实验数据//t,lamb,lambdif,pb  //1Hz
        "0.00         0.959110         0.125664         0.101944
    0.02         0.961607         0.124673         0.101015
    ... ...  //数据省略
    1.00         0.959110         0.125664         0.101944"
        },

        tArray4=matrix{ //存放实验数据//t,lamb,lambdif,pb  //2Hz
        "0.00         0.957683         0.251327         0.103554
    0.02         0.962873         0.243431         0.099815
    ... ...  //数据省略
    0.50         0.957683         0.251327         0.103554"
        },

        tArray5=matrix{ //存放实验数据//t,lamb,lambdif,pb  //5Hz
        "0.00         0.953477         0.628319         0.135643
    0.02         0.966487         0.508320         0.098394
    ... ...  //数据省略
    0.20         0.953477         0.628319         0.135643"
        },

        Opt1[@目标函数,optmax,200,optmode,10,optwaycom, optwaysimdeep, optwayconfra, optrange : 1e-5,1e10; 1e-5,1e10; 1e-5,1e10; 1e-5,1e10; -1e10,-1e-5; -1e10,-1e-5; 1e-5,1e10;  -1e10,1e10;  -1e10,1e10;  -1e10,1e10;  -1e10,1e10;  -1e10,1e10 ] //Opt1函数全局优化
    };

    结果:
    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进行线性插值或未作插值,也有着天壤之别。

  • wlfc

    引用回帖:
    6楼: Originally posted by wlfc at 2021-08-20 21:16:50
    用OpenLu求解:

    分析:需在微分方程求解中传递中间变量lamb,gsl支持的Lu扩展数学库中求解微分方程的函数gsl_ode提供了传递中间变量的功能。
    进一步讨论:微分方程求解传递中间变量时,按时间t的改变量进行插值 ...

    关于变量lamb进行线性插值或未作插值,微分方程计算差别大的问题,有了新发现。
    问题在于gsl_ode中的算法选择,以上用的是gsl_rkf45(Runge-Kutta-Fehlberg 法),如果改为gsl_rk4 (四阶经典 Runge-Kutta 法),便没有问题。

    另外,以上Lu代码的两个结果,经matlab相同的ode45算法验证,结果一致。

猜你喜欢
下载小木虫APP
与700万科研达人随时交流
  • 二维码
  • IOS
  • 安卓