歡迎訪問 Lu程序設計

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

請使用OpenLu(32位版本)演示本文的例子

例子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的最優值?

    Lu代碼1:

!!!using["IMSL","luopt","math"]; //使用命名空間
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,tyz : tyArray,tA,max,k,b)=
{
    k=_k,b=_b,                
//傳遞優化變量,函數f中要用到k,b
    tyz=ode[@f,tA,
ra1(2.329101563,5000)],
    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取矩陣的列
    ClearImslErr(),            
//清空IMSL錯誤輸出
    ERSET(0,0,0),              
//關閉IMSL所有警告
    Opt[@目標函數],     
       //Opt函數全局優化
    ERSET(0,2,2), ERSET(0,1,0)  //恢復IMSL警告
};

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

1.408334652276297e-004 -1.467403743595575e-004 24.2912485903651

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

!!!using["IMSL","luopt","math","win"]; //使用命名空間
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,tyz : tyArray,tA,max,k,b)=
{
    k=_k,b=_b,                 //傳遞優化變量,函數f中要用到k,b
    tyz=ode[@f,tA,ra1(2.329101563,5000)],
    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取矩陣的列
    ClearImslErr(),             //清空IMSL錯誤輸出
    ERSET(0,0,0),               //關閉IMSL所有警告
    k=0.0, b=0.0, Opt[@目標函數, optstarter,&k,&b,0],   //Opt函數全局優化
    ERSET(0,2,2), ERSET(0,1,0), //恢復IMSL警告
    tyz=ode[@f,tyArray(all:0),ra1(2.329101563,5000)],  //獲得最終結果
    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["IMSL","luopt","math"];
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 : tf,i,s : tArray,tA,max,a1,a2,a3,a4,a5)=
{
    a1=_a1, a2=_a2, a3=_a3, a4=_a4, a5=_a5,
    tf=ode[@f,tA,
ra1(0,0.15)],
    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),
    ClearImslErr(),            
//清空IMSL錯誤輸出
    ERSET(0,0,0),              
//關閉IMSL所有警告
    Opt[@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

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

!!!using["IMSL","luopt","math","win"];
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 : tf,i,s : tArray,tA,max,a1,a2,a3,a4,a5)=
{
    a1=_a1, a2=_a2, a3=_a3, a4=_a4, a5=_a5,
    tf=ode[@f,tA,
ra1(0,0.15)],
    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),
    ClearImslErr(),            
//清空IMSL錯誤輸出
    ERSET(0,0,0),              
//關閉IMSL所有警告
    a1=0.0, a2=0.0, a3=0.0, a4=0.0, a5=0.0,
    Opt[@J, optstarter,&a1,&a2,&a3,&a4,&a5,0],
  //Opt函數全局優化
    ERSET(0,2,2), ERSET(0,1,0), //恢復IMSL警告
    tf=ode[@f,tA,ra1(0,0.15)],  //獲得最終結果
    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["IMSL","luopt","math"];
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,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=ode[@f,tArray,
ra1(30.0e-6, 10.2e-6, 4.1e-6, 6.7)],
    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=ode[@f,tArray,
ra1(30.0e-6, 10.2e-6, 4.1e-6, 6.7)],
    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=ode[@f,tArray,
ra1(30.0e-6, 10.2e-6, 4.1e-6, 6.7)],
    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=ode[@f,tArray,
ra1(30.0e-6, 10.2e-6, 4.1e-6, 6.7)],
    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],
    ClearImslErr(),            
//清空IMSL錯誤輸出
    ERSET(0,0,0),              
//關閉IMSL所有警告
    Opt[@目標函數, 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

例子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["IMSL","luopt","math"];
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,nn : tArray,n1Array,max,a,b,c)=
{
    a=aa, b=bb, c=cc,
    nn=ode[@f,tArray,
ra1(n0, n0)],
    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},
    ClearImslErr(),            
//清空IMSL錯誤輸出
    ERSET(0,0,0),              
//關閉IMSL所有警告
    Opt[@目標函數, 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

例子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["IMSL","luopt","math"];        //使用命名空間
f(t,C,dC::k0,Ea,m,n,o,T,A,B)= dC=-k0*exp{-Ea/[8.314*(273.15+T)]}*(A^m)*(B^n)*(C^o);
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=ode[@f,tTime,ra1(tC(0,i))],
        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"
    },
    ClearImslErr(),             //清空IMSL錯誤輸出
    ERSET(0,0,0),               //關閉IMSL所有警告
    Opt[@g],                    //Opt函數全局優化
    ERSET(0,2,2), ERSET(0,1,0)  //恢復IMSL警告
};

結果(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["IMSL","luopt","math"]; //使用命名空間
f(t,CB,CH,dCB,dCH::K1,K2,a,m)=
{
    dCB=-K1*CB^a*(K1/K2*CB^a)^m,
    dCH=K2*(K1/K2*CB^a)^m
};
目標函數(_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=ode[@f,tA,ra1(10.9237,0)],
    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取矩陣的列
    ClearImslErr(),
//清空IMSL錯誤輸出
    ERSET(0,0,0),
//關閉IMSL所有警告
    Opt[@目標函數,optwaysimdeep, optwayconfra],
//Opt函數全局優化
    ERSET(0,2,2), ERSET(0,1,0)
//恢復IMSL警告
};

結果(K1,K2,a,m,目標函數值):

2.007221403771439e-002 3.477931047132456e-002 0.7598653691173753 -1.205432420176749 18.90800296978626

繪圖:

!!!using("IMSL","win","math");
f(t,CB,CH,dCB,dCH::K1,K2,a,m)=
{
    dCB=-K1*CB^a*(K1/K2*CB^a)^m,
    dCH=K2*(K1/K2*CB^a)^m
};
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, ode[@f,tx,ra1(10.9237,0)]);
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為優化參數。

代碼:

!!!using["IMSL","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::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
};
目標函數(_lp,_ps,_b : i,s,tyz : tac,t,max,lp,ps)=
{
    lp=_lp,ps=_ps,
//傳遞優化變量,函數f中要用到lp,ps
    tyz=ode[@f,t,ra1(0.719,_b, 1)],
    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取矩陣的列
    ClearImslErr(),
//清空IMSL錯誤輸出
    ERSET(0,0,0),
//關閉IMSL所有警告
    Opt[@目標函數],
//Opt函數全局優化
    ERSET(0,2,2), ERSET(0,1,0)
//恢復IMSL警告
};

結果(lp,ps,_b,目標函數值):

2.473962346617309e-003 -3.626422822558594e-003 1.308637154884405 0.1199412216587657

繪圖:

!!!using("IMSL","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::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
};
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=ode[@f,tx,ra1(0.719,1.308637154884405, 1)],
    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=1.0,1.277,1.555,1.833,2.111,2.388,2.666,2.944,3.222,3.5;
s=0,41004.0,53344.0,58652.0,61883.0,84456.0,66836,69186.0,71566.0,73995.0;

分析:因z值無法測出,故追加z的初值_z為優化參數。

代碼:

!!!using["IMSL","luopt","math"]; //使用命名空間
f(t,z,s,dz,ds: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
};
目標函數(_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=ode[@f,t,ra1(_z,0)],
    i=0, ss=0, while{++i<max, ss=ss+[ts(i,2)-s(i)]^2},
    sqrt[ss/max]
};
main(::t,s,max)=
{
    t=ra1[1.0,1.277,1.555,1.833,2.111,2.388,2.666,2.944,3.222,3.5],
    s=ra1[0,41004.0,53344.0,58652.0,61883.0,84456.0,66836,69186.0,71566.0,73995.0],
    max=len[t],
    ClearImslErr(),
//清空IMSL錯誤輸出
    ERSET(0,0,0),
//關閉IMSL所有警告
    Opt[@目標函數],
//Opt函數全局優化
    ERSET(0,2,2), ERSET(0,1,0)
//恢復IMSL警告
};

結果(ue,ae,uv,av,tv,_z,目標函數值):

-9114964.001609007 -0.4453565389059405 7335590868.341645 -2.754919559106283e-004 -10485633.9903882 2.654977892066841 244.5944885704022

繪圖:

!!!using("IMSL","win","math");
f(t,z,s,dz,ds: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
};
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', ode[@f,t,ra1(2.654977892066841,0)]), outa[ode[@f,t,ra1(2.654977892066841,0)]];
ChartWnd[@init];
 //顯示窗口并初始化


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

欧冠杯