歡迎訪問 Forcal程序設計

Forcal 微分方程參數優化(擬合)

推薦:Lu 微分方程參數優化(擬合)

例子1 數學模型:

      dy/dt=k*y*z+0.095*b*z
      dz/dt=
-b*z-0.222*z

    實驗數據(ti; yi):

 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

    已知z在t=0.01時的初值為5000,求k和b的最優值?

    Forcal代碼1:

!using["IMSL","fcopt","math"]; //使用命名空間
實驗數據(::tyArray,max)=
{
    max=18,                   
//實驗數據組數
    tyArray=arrayinitns{max,2,
//存放實驗數據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"
    }.free()
};
f(t,y,z,dy,dz::k,b)=
{
    dy=k*y*z+0.095*b*z,
    dz=-b*z-0.222*z
};
目標函數(_k,_b : i,s,IDO,y,z,yy,t,tt : tyArray,max,k,b)=
{
    k=_k,b=_b,                
//傳遞優化變量,函數f中要用到k,b
    i=0, s=0, IDO=1, t=tyArray[0,0], y=tyArray[0,1], z=5000,
    (++i<max).while{
        tt=tyArray[i,0], yy=tyArray[i,1],
        IVPRK[&IDO,HFor("f"),&t,tt,1e-6,0,&y,&z],
        s=s+[y-yy]^2
    },
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,&y,&z],
    s
};
SimOpt[HFor("目標函數")];     
//求n維極值的單形調優法

    結果(前面的數為最優參數,最后一個數是極小值。下同):

1.408334652276297e-004 -1.467403743595575e-004 24.2912485903651

    Forcal代碼2:優化(擬合)并繪圖

!using["IMSL","fcopt","math","fc2d"];
實驗數據(::tyArray,max,tA,tyA)=
{
    max=18,
    tyArray=arrayinitns{max,2,
        "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"
    }.free(),
    tA=tyArray
(neg:0).free(), tyA=tyArray(neg:1).free()  //獲得列向量
};
f(t,y,z,dy,dz::k,b)=
{
    dy=k*y*z+0.095*b*z,
    dz=-b*z-0.222*z
};
目標函數(_k,_b : i,s,IDO,y,z,yy,t,tt : tyArray,max,k,b)=
{
    k=_k,b=_b,
    i=0, s=0, IDO=1, t=tyArray[0,0], y=tyArray[0,1], z=5000,
    (++i<max).while{
        tt=tyArray[i,0], yy=tyArray[i,1],
        IVPRK[&IDO,HFor("f"),&t,tt,1e-6,0,&y,&z],
        s=s+[y-yy]^2
    },
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,&y,&z],
    s
};
yi
(x,y,n::tA,tyA,max)= x=tA, y=tyA, n=max;  //繪制實驗數據點(ti,yi)的函數
yyi(tt : y,z,t : tyArray)=                  //繪制擬合曲線的函數
{
    t=tyArray[0,0], y=tyArray[0,1], z=5000,
    IVPRK[1,HFor("f"),&t,tt,1e-6,0,&y,&z],
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,y,z],
    y
};
main(::k,b)=
{
    SimOpt[HFor("目標函數"), optstarter,&k,&b,0]
,
   
Plot{Ix : 0,10,
     Igarray : HFor("yi"), Acolor,Vred,
Adot, Athick,8, //繪制函數yi的圖形,將繪制散點圖
     Iufun : HFor("yyi"), Acolor,Vblue, Adots,100
      //繪制函數yyi的圖形
    }

};

    結果:

例子2 數學模型:三房室模型給藥后,藥物在房室1和2,2和3間轉運,并經房室3排出。假定轉運均符合一級速率模型,則各房室藥物量(X1,X2,X3)與時間t的關系如下:
        dX1/dt=-K12*X1+K21*X2
        dX2/dt=(-K21-K23)*X2+K12*X1+K32*X3
        dX3/dt=(-K32-K30)*X3+K23*X2
    已知初值t=0時:X1=100, X2=X3=0 。
    實驗值:
        t=0: X1=100, X2=0, X3=0
        t=1: X1=90,  X2=8, X3=2
        t=3: X1=73,  X2=19,X3=7
        t=5: X1=60,  X2=14,X3=23
    求擬合參數K12,K21,K23,K32,和K30。
    限制條件:K12,K21,K23,K32,和K30均大于零;X1+X2+X3小于等于100(總量不能超過給藥量)。
    擬合可以用非線形最小二乘法,即求取minΣ(y實驗值-y真實值)2時的K值。
    當然這個問題的結構相對簡單,可以求出各方程積分式,然后進行擬合。但是對于更復雜的系統,如十房室以上,積分式將非常復雜,手工求方程積分式幾近不可能。

    Forcal代碼1:

!using["IMSL","fcopt","math","sys"];
實驗數據(::tArray,max)=
{
    max=4,
    tArray=arrayinit{2,max,4 :
        0, 100, 0,  0,
        1, 90,  8,  2,
        3, 73,  19, 7,
        5, 60,  14, 23
    }.free()
};
f(t,X1,X2,X3,dX1,dX2,dX3::k12,k21,k23,k32,k30)=
{
    dX1=-k12*X1+k21*X2,
    dX2=(-k21-k23)*X2+k12*X1+k32*X3,
    dX3=(-k32-k30)*X3+k23*X2
};
J(_k12,_k21,_k23,_k32,_k30 : IDO,t1,t2,X1,X2,X3,X11,X22,X33,i,s : k12,k21,k23,k32,k30,tArray,max)=
    //目標函數
{
    k12=_k12, k21=_k21, k23=_k23, k32=_k32, k30=_k30,
    tArray.GA(0,&t1,&X1,&X2,&X3),
    IDO=1, s=0, i=0, (++i<max).while{
        tArray.GA(i*4,&t2,&X11,&X22,&X33),
        IVPRK[&IDO,HFor("f"),&t1,t2,1e-6,0,&X1,&X2,&X3],
        s=s+(X1-X11)^2+(X2-X22)^2+(X3-X33)^2
    },
    IVPRK[3,HFor("f"),&t1,t2+0.1,1e-6,0,&X1,&X2,&X3],
    s
};
S(_k12,_k21,_k23,_k32,_k30 : IDO,t1,t2,X1,X2,X3,i,s : k12,k21,k23,k32,k30,tArray,max)=   
//約束函數
{
    k12=_k12, k21=_k21, k23=_k23, k32=_k32, k30=_k30,
    tArray.GA(0,&t1,&X1,&X2,&X3),
    IDO=1, s=1, i=0, (++i<max).while{
        IVPRK[&IDO,HFor("f"),&t1,t2,1e-6,0,&X1,&X2,&X3],
        s=s & (X1+X2+X3<=100),
     //計算約束條件
        i++
    },
    IVPRK[3,HFor("f"),&t1,t2+0.1,1e-6,0,&X1,&X2,&X3],
    s
};
main(:k12,k21,k23,k32,k30)=
{
    k12=0.1, k21=0.1, k23=0.1, k32=0.1, k30=0.1,   
//參數初值
    //獲得一組符合約束的初值
    ROptIni[HFor("J"),HFor("S"),optmax,1000, optstarter : &k12, &k21, &k23, &k32, &k30, 0, optrange : 0,1e10, 0,1e10, 0,1e10, 0,1e10, 0,1e10],
    //優化
    RGOpt[HFor("J"),HFor("S"),optmax,2000, optmode,1, optstarter : &k12, &k21, &k23, &k32, &k30, 0, optrange : 0,1e10, 0,1e10, 0,1e10, 0,1e10, 0,1e10]
};

    結果:

0.1154404925322101 8.73875042635457e-002 0.3435968132778263 0. 8.11974164395054e-003 33.56570615692139

    Forcal代碼2:優化(擬合)并繪圖

!using["IMSL","fcopt","math","sys","fc2d"];
實驗數據(::tArray,max,tX,tX1,tX2,tX3)=
{
    max=4,
    tArray=arrayinit{2,max,4 :
        0, 100, 0,  0,
        1, 90,  8,  2,
        3, 73,  19, 7,
        5, 60,  14, 23
    }.free(),
    tX=tArray
(neg:0).free(), tX1=tArray(neg:1).free(), tX2=tArray(neg:2).free(), tX3=tArray(neg:3).free()
};
f(t,X1,X2,X3,dX1,dX2,dX3::k12,k21,k23,k32,k30)=
{
    dX1=-k12*X1+k21*X2,
    dX2=(-k21-k23)*X2+k12*X1+k32*X3,
    dX3=(-k32-k30)*X3+k23*X2
};
J(_k12,_k21,_k23,_k32,_k30 : IDO,t1,t2,X1,X2,X3,X11,X22,X33,i,s : k12,k21,k23,k32,k30,tArray,max)=

{
    k12=_k12, k21=_k21, k23=_k23, k32=_k32, k30=_k30,
    tArray.GA(0,&t1,&X1,&X2,&X3),
    IDO=1, s=0, i=0, (++i<max).while{
        tArray.GA(i*4,&t2,&X11,&X22,&X33),
        IVPRK[&IDO,HFor("f"),&t1,t2,1e-6,0,&X1,&X2,&X3],
        s=s+(X1-X11)^2+(X2-X22)^2+(X3-X33)^2
    },
    IVPRK[3,HFor("f"),&t1,t2+0.1,1e-6,0,&X1,&X2,&X3],
    s
};
S(_k12,_k21,_k23,_k32,_k30 : IDO,t1,t2,X1,X2,X3,i,s : k12,k21,k23,k32,k30,tArray,max)=
{
    k12=_k12, k21=_k21, k23=_k23, k32=_k32, k30=_k30,
    tArray.GA(0,&t1,&X1,&X2,&X3),
    IDO=1, s=1, i=0, (++i<max).while{
        IVPRK[&IDO,HFor("f"),&t1,t2,1e-6,0,&X1,&X2,&X3],
        s=s & (X1+X2+X3<=100),
        i++
    },
    IVPRK[3,HFor("f"),&t1,t2+0.1,1e-6,0,&X1,&X2,&X3],
    s
};
X1
(x,y,n::tX,tX1,max)= x=tX, y=tX1, n=max;  //繪制實驗數據點(t,X1)的函數
X2(x,y,n::tX,tX2,max)= x=tX, y=tX2, n=max;  //繪制實驗數據點(t,X2)的函數
X3(x,y,n::tX,tX3,max)= x=tX, y=tX3, n=max;  //繪制實驗數據點(t,X3)的函數
XX1(tt : X1,X2,X3,t : tArray)=              //繪制擬合曲線X1的函數
{
    tArray.GA(0,&t,&X1,&X2,&X3),
    IVPRK[1,HFor("f"),&t,tt,1e-6,0,&X1,&X2,&X3],
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,X1,X2,X3],
    X1
};
XX2(tt : X1,X2,X3,t : tArray)= 
            //繪制擬合曲線X2的函數
{
    tArray.GA(0,&t,&X1,&X2,&X3),
    IVPRK[1,HFor("f"),&t,tt,1e-6,0,&X1,&X2,&X3],
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,X1,X2,X3],
    X2
};
XX3(tt : X1,X2,X3,t : tArray)= 
            //繪制擬合曲線X3的函數
{
    tArray.GA(0,&t,&X1,&X2,&X3),
    IVPRK[1,HFor("f"),&t,tt,1e-6,0,&X1,&X2,&X3],
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,X1,X2,X3],
    X3
};
XXX(tt : X1,X2,X3,t : tArray)= 
            //繪制擬合曲線X1+X2+X3的函數
{
    tArray.GA(0,&t,&X1,&X2,&X3),
    IVPRK[1,HFor("f"),&t,tt,1e-6,0,&X1,&X2,&X3],
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,X1,X2,X3],
    X1+X2+X3
};
main(::k12,k21,k23,k32,k30)=
{
    k12=0.1, k21=0.1, k23=0.1, k32=0.1, k30=0.1,
    ROptIni[HFor("J"),HFor("S"),optmax,1000, optstarter : &k12, &k21, &k23, &k32, &k30, 0, optrange : 0,1e10, 0,1e10, 0,1e10, 0,1e10, 0,1e10],
    RGOpt[HFor("J"),HFor("S"),optmax,2000, optmode,1, optstarter : &k12, &k21, &k23, &k32, &k30, 0, optrange : 0,1e10, 0,1e10, 0,1e10, 0,1e10, 0,1e10],
   
Plot{Ix : 0,5, Iy : 0,110, Iydynamic : 0,
     Igarray : HFor("X1"), Acolor,Vred,
Adot, Athick,8, Alegend,"X1實驗數據",   //繪制函數X1的圖形,將繪制散點圖
     Igarray : HFor("X2"), Acolor,Vblue,
Adot, Athick,8, Alegend,"X2實驗數據",  //繪制函數X2的圖形,將繪制散點圖
     Igarray : HFor("X3"), Acolor,Vgreen,
Adot, Athick,8, Alegend,"X3實驗數據", //繪制函數X3的圖形,將繪制散點圖
     Iufun : HFor("XX1"), Acolor,Vred, Adots,100, Alegend,"X1模擬數據"
,         //繪制函數XX1的圖形
     Iufun : HFor("XX2"), Acolor,Vblue, Adots,100, Alegend,"X2模擬數據"
,        //繪制函數XX2的圖形
     Iufun : HFor("XX3"), Acolor,Vgreen, Adots,100, Alegend,"X3模擬數據"
,       //繪制函數XX3的圖形
     Iufun : HFor("XXX"), Acolor,RGB(255,128,0), Adots,100
, Athick,8, Alegend,"X1+X2+X3模擬數據"  //繪制函數XXX的圖形
    }

};

    結果:

例子3 數學模型:需要擬合的微分方程組為:

        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

    Forcal代碼1:

!using["IMSL","fcopt","math"];
init(::tArray,max)=
{
    max=101,
    tArray=arrayinit{2,max,2 :
        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
    }.free()
};
f(t,f1,f2,df1,df2 : x : a1,a2,a3,a4,a5)=
{
    x=0.3,
    df1=a1*(x-f2)^a2,
    df2=a3*df1-a4*f2^a5
};
J(_a1,_a2,_a3,_a4,_a5 : IDO,i,s,f1,f2,f11,t1,t2 : tArray,max,a1,a2,a3,a4,a5)=
{
    a1=_a1, a2=_a2, a3=_a3, a4=_a4, a5=_a5,
    t1=tArray(0,0), f1=tArray(0,1), f2=0.15,
    IDO=1, i=0, s=0,
    (++i<max).while{
        t2=tArray(i,0), f11=tArray(i,1),
        IVPRK[&IDO,HFor("f"),&t1,t2,1e-6,0,&f1,&f2],
        s=s+[f1-f11]^2
    },
    IVPRK[3,HFor("f"),&t1,t2+0.1,1e-6,0,&f1,&f2],
    s
};
ClearImslErr(),            
//清空IMSL錯誤輸出
ERSET(0,0,0),              
//關閉IMSL所有警告
Opt[HFor("J")],     
       //Opt函數全局優化
ERSET(0,2,2), ERSET(0,1,0); //恢復IMSL警告

    結果(沒有多次運行,或許還有更優解):

27594.47852777774 8.968224513288728 18.32857783239671 -4.477342854126431e-006 4.851885736868509e-015 1.515524149642018e-009

    Forcal代碼2:優化(擬合)并繪圖

!using["IMSL","fcopt","math","fc2d"];
init(::tArray,max,tA,f1A)=
{
    max=101,
    tArray=arrayinit{2,max,2 :
        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
    }.free(),
    tA=tArray
(neg:0).free(), f1A=tArray(neg:1).free()
};
f(t,f1,f2,df1,df2 : x : a1,a2,a3,a4,a5)=
{
    x=0.3,
    df1=a1*(x-f2)^a2,
    df2=a3*df1-a4*f2^a5
};
J(_a1,_a2,_a3,_a4,_a5 : IDO,i,s,f1,f2,f11,t1,t2 : tArray,max,a1,a2,a3,a4,a5)=
{
    a1=_a1, a2=_a2, a3=_a3, a4=_a4, a5=_a5,
    t1=tArray(0,0), f1=tArray(0,1), f2=0.15,
    IDO=1, i=0, s=0,
    (++i<max).while{
        t2=tArray(i,0), f11=tArray(i,1),
        IVPRK[&IDO,HFor("f"),&t1,t2,1e-6,0,&f1,&f2],
        s=s+[f1-f11]^2
    },
    IVPRK[3,HFor("f"),&t1,t2+0.1,1e-6,0,&f1,&f2],
    s
};
f1i
(x,y,n::tA,f1A,max)= x=tA, y=f1A, n=max;  //繪制實驗數據點(t,f1)的函數
ff1i
(tt : f1,f2,t : tArray)=                 //繪制擬合曲線的函數
{
    t=tArray[0,0], f1=tArray[0,1], f2=0.15,
    IVPRK[1,HFor("f"),&t,tt,1e-6,0,&f1,&f2],
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,f1,f2],
    f1
};
main(::a1,a2,a3,a4,a5)=
{
    ClearImslErr(),            
//清空IMSL錯誤輸出
    ERSET(0,0,0),              
//關閉IMSL所有警告
    Opt[HFor("J"), optstarter,&a1,&a2,&a3,&a4,&a5,0],
  //Opt函數全局優化
    Plot{Iclear, Ix : 0,1000,
     Igarray : HFor("f1i"), Acolor,Vred,
Adot,         //繪制函數f1i的圖形,將繪制散點圖
     Iufun : HFor("ff1i"), Acolor,Vblue, Adots,100
    //繪制函數ff1i的圖形
    }
,
    ERSET(0,2,2), ERSET(0,1,0)
 //恢復IMSL警告
};

    結果:

例子4 數學模型:

有關于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]
y(1)=[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]
y(1)=[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組數據,故只有一組最優解。

Forcal代碼:

!using["IMSL","fcopt","math"];
實驗數據(::tArray,y11Array,y12Array,y13Array,y14Array,max)=
{
    max=5,
    tArray=arrayinit[1,max : 0, 10, 20, 30, 40].free(),
    y11Array=arrayinit[1,max : 30.0e-6, 17.5e-6, 15e-6, 14e-6, 13.5e-6].free(),
    y12Array=arrayinit[1,max : 30.0e-6, 12.5e-6, 10e-6, 8e-6, 7e-6].free(),
    y13Array=arrayinit[1,max : 30.0e-6, 11.5e-6, 9e-6, 7.5e-6, 6e-6].free(),
    y14Array=arrayinit[1,max : 30.0e-6, 10.5e-6, 8.5e-6, 6.5e-6, 5e-6].free()
};
f(t,y1,y2,y3,y4,z1,z2,z3,z4 : 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
};
目標函數(kka,kkb,kkc : i,s,IDO,t,tt,y1,y2,y3,y4,z1 : 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]
    i=0, IDO=1, t=tArray[0], y1=30.0e-6, y2=10.2e-6, y3=4.1e-6, y4=6.7,
    (++i<max).while{
        tt=tArray[i], z1=y11Array[i],
        IVPRK[&IDO,HFor("f"),&t,tt,1e-6,0,&y1,&y2,&y3,&y4],
        s=s+(y1-z1)^2
    },
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,&y1,&y2,&y3,&y4],
    a=120, b=15, c=25,
//[a b c]=[120, 15, 25]
    i=0, IDO=1, t=tArray[0], y1=30.0e-6, y2=10.2e-6, y3=4.1e-6, y4=6.7,
    (++i<max).while{
        tt=tArray[i], z1=y12Array[i],
        IVPRK[&IDO,HFor("f"),&t,tt,1e-6,0,&y1,&y2,&y3,&y4],
        s=s+(y1-z1)^2
    },
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,&y1,&y2,&y3,&y4],
    a=130, b=25, c=30,
//[a b c]=[130, 25, 30]
    i=0, IDO=1, t=tArray[0], y1=30.0e-6, y2=10.2e-6, y3=4.1e-6, y4=6.7,
    (++i<max).while{
        tt=tArray[i], z1=y13Array[i],
        IVPRK[&IDO,HFor("f"),&t,tt,1e-6,0,&y1,&y2,&y3,&y4],
        s=s+(y1-z1)^2
    },
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,&y1,&y2,&y3,&y4],
    a=140, b=30, c=35,
//[a b c]=[140, 30, 35]
    i=0, IDO=1, t=tArray[0], y1=30.0e-6, y2=10.2e-6, y3=4.1e-6, y4=6.7,
    (++i<max).while{
        tt=tArray[i], z1=y14Array[i],
        IVPRK[&IDO,HFor("f"),&t,tt,1e-6,0,&y1,&y2,&y3,&y4],
        s=s+(y1-z1)^2
    },
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,&y1,&y2,&y3,&y4],
    sqrt[s/(max*4-4)]
};
ClearImslErr(),            
//清空IMSL錯誤輸出
ERSET(0,0,0),              
//關閉IMSL所有警告
Opt[HFor("目標函數"), optrange,-1e5,1e5,-1e5,1e5,-1e5,1e5],
//Opt函數優化
ERSET(0,2,2), ERSET(0,1,0);
//恢復IMSL警告

結果(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

例子5 數學模型:

微分方程組為:

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。

Forcal代碼:

!using["IMSL","fcopt","math"];
實驗數據(::tArray,n1Array,max)=
{
    max=6,
    tArray=arrayinitns{max : "1 2 3 4 5 6"}.free(),
    n1Array=arrayinitns{max : "1 3.5 2.8 2.4 1.5 0.9"}.free()
};
f(t,n1,n2,dn1,dn2::a,b,c)=
{
    dn1=a*n1*n2-b*n1,
    dn2=c*n1*n2
};
目標函數(n0,aa,bb,cc : i,s,IDO,n1,n2,nn1,t,tt : tArray,n1Array,max,a,b,c)=
{
    a=aa, b=bb, c=cc,
    i=0, s=0, IDO=1, t=0, n1=n0, n2=n0,
    (i<max).while{
        tt=tArray[i], nn1=n1Array[i],
        IVPRK[&IDO,HFor("f"),&t,tt,1e-6,0,&n1,&n2],
        s=s+[n1-nn1]^2,
        i++
    },
    IVPRK[3,HFor("f"),&t,tt+0.1,1e-6,0,&n1,&n2],
    sqrt[s/max]
};
ClearImslErr(),            
//清空IMSL錯誤輸出
ERSET(0,0,0),              
//關閉IMSL所有警告
Opt[HFor("目標函數"), optwaysimdeep, optwayconfra, optrange,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10],
//Opt函數優化
ERSET(0,2,2), ERSET(0,1,0);
//恢復IMSL警告

結果(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


版權所有© Forcal程序設計 2002-2011,保留所有權利
E-mail: [email protected]
  QQ:630715621
最近更新: 2012年11月02日

欧冠杯