2011年3月6日日曜日

生物７種の計算例



（2011/03/09）
　図ではA1（緑藻）、A２（大型珪藻）としていますが、プログラム等ではA1（大型珪藻）、A２（緑藻）でとしています。
投稿者 時刻:

2 件のコメント:

  1. kf332011年3月6日 23:22

    //初期値等

    //---------------------------------------------------------------------------------------------[ Constants ]

    Mmax1 = 0.2 , Mmax2 = 0.3
    K1= 0.5 , K2=1
    Ts1 = 10 , Ts2 = 20

    dZ1=1/60 , dZ2=1/60 , dF1=1/900 , dF2=1/900 , dF3=1/900
    fZ1=0.1 , fZ2=0.1 , fF1=0.1 , fF2=0.1 , fF3_1=0.01 , fF3_2=0.1
    eZ1=0.2 , eZ2=0.2 , eF1=0.1 , eF2=0.1 , eF3=0.1

    //---------------------------------------------------------------------------------------------[ Initial Condition ]

    Po = 2 // Phoshorous
    A1o = 0.1 // oogatakeisou
    A2o = 0.1 // ryokusou
    Z1o = 0.01 // harinagamijinko
    Z2o = 0.01 // zoumijinko
    F1o = 0.001 // himemasu
    F2o = 0.002 // wakasagi
    F3o = 0.0005 // sakuramasu
    POMo = 0.001 // POM
    Do = 0.0001 // Detritas


    InitialCondition = [ Po ; A1o ; A2o ; Z1o ; Z2o ; F1o ; F2o ; F3o ; POMo ; Do ]

    //---------------------------------------------------------------------------------------------[ Matrix Number ]

    NN=10

    // ----------------------------------------------[Graphics]

    Time = 3000
    Ymax = 2.2 / 1
    Color = [ 1 2 3 4 5 6 7 2 ]

    返信削除
  2. kf332011年3月6日 23:23

    P=u1(1,1)
    A1=u1(2,1)
    A2=u1(3,1)
    Z1=u1(4,1)
    Z2=u1(5,1)
    F1=u1(6,1)
    F2=u1(7,1)
    F3=u1(8,1)
    POM=u1(9,1)
    D=u1(10,1)

    Temp=( 22.9/2 ) * sin( 2*3.14/365 * u2(1,1) ) + 16.6 ;y2(1,1)=Temp

    Kt1= ( ( Temp / Ts1 ) * exp( ( Ts1 -Temp ) / Ts1 ) ) ^3
    Kt2= ( ( Temp / Ts2 ) * exp( ( Ts2 -Temp ) / Ts2 ) ) ^3
    //------------------------------------------------------------------ [ Matrix ]

    //-------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- --------------
    // P A1 A2 Z1 Z2 F1 F2 F3 POM D
    //-------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- --------------
    Mat= [ 0 A1*Mmax1*P/(K1+P)*Kt1 A2*Mmax2*P/(K2+P)*Kt2 0 0 0 0 0 0 0 ;
    0 0 0 A1*Z1*fZ1*eZ1 0 0 0 0 (A1*Z1*fZ1*(1-eZ1))+A1*(1/100) 0 ;
    0 0 0 A2*Z1*fZ2*eZ2 A2*Z2*fZ2*eZ2 0 0 0 (A2*Z1*fZ1*(1-eZ1))+(A2*Z2*fZ2*(1-eZ2))+A2*(1/100) 0 ;
    0 0 0 0 0 Z1*F1*fF1*eF1 Z1*F2*fF2*eF2 0 (Z1*F1*fF1*(1-eF1))+(Z1*F2*fF2*(1-eF2))+Z1*(dZ1) 0 ;
    0 0 0 0 0 0 Z2*F2*fF2*eF2 0 (Z2*F2*fF2*(1-eF2))+Z2*(dZ2) 0 ;
    0 0 0 0 0 0 0 F1*F3*fF3_1*eF3 F1*dF1 0 ;
    0 0 0 0 0 0 0 F2*F3*fF3_2*eF3 F2*dF2 0 ;
    0 0 0 0 0 0 0 0 F3*dF3 0 ;
    0 0 0 0 0 0 0 0 0 POM*0.5 ;
    D*(1/100) 0 0 0 0 0 0 0 0 0 ;
    ]
    //-------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- --------------


    //----------------------------------------------------------------- [ Matrix ---> xdot ]
    for i= 1 : NN ; y1(i,1)= - sum(Mat(i , :) ) + sum(Mat(: , i) ) ; end

    返信削除

登録: コメントの投稿 (Atom)