2010年12月5日日曜日

琵琶湖２（緑藻５種類）





--- 珪藻

--- Urogrena

--- 緑藻
　・下の線
　　　個々の緑藻
　・上の線
　　　５種の合計



投稿者 時刻:

1 件のコメント:

  1. kf332010年12月5日 21:13

    //　競合 (２０年分 珪藻・Uroglena・緑藻)

    umax1=0.25 ; umax2=0.20; umax3=0.22;

    Ts1=10; Ts2=20;
    //---------
    Tsg1=20; Tsg2=21; Tsg3=22; Tsg4=23; Tsg5=24;

    D=1/20;
    Du=1/22.2;
    Year=20;

    Kp=0.1;

    Xs= [2; 0.5; 0.3; 0.1; 0.1; 0.1; 0.1; 0.1; 0];
    Xin=[2; 0; 0; 0; 0; 0; 0; 0; 0];

    X=Xs;

    function Xdot=f(t,X) ;

    //---------
    P=X(1)
    A1=X(2)
    A2=X(3)
    //---------
    Ag1=X(4)
    Ag2=X(5)
    Ag3=X(6)
    Ag4=X(7)
    Ag5=X(8)

    //---------

    Temp=15*sin(2*3.1415*t/365-158)+(19)
    Kt1=(Temp/Ts1*exp((Ts1-Temp)/Ts1))^3
    Kt2=(Temp/Ts2*exp((Ts2-Temp)/Ts2))^3
    //---------
    Ktg1=(Temp/Tsg1*exp((Tsg1-Temp)/Tsg1))^3
    Ktg2=(Temp/Tsg2*exp((Tsg2-Temp)/Tsg2))^3
    Ktg3=(Temp/Tsg3*exp((Tsg3-Temp)/Tsg3))^3
    Ktg4=(Temp/Tsg4*exp((Tsg4-Temp)/Tsg4))^3
    Ktg5=(Temp/Tsg5*exp((Tsg5-Temp)/Tsg5))^3

    V1=A1*umax1*P/(Kp+P)*Kt1
    V2=A2*umax2*P/(Kp+P)*Kt2
    //---------
    Vg1=Ag1*umax3*P/(Kp+P)*Ktg1
    Vg2=Ag2*umax3*P/(Kp+P)*Ktg2
    Vg3=Ag3*umax3*P/(Kp+P)*Ktg3
    Vg4=Ag4*umax3*P/(Kp+P)*Ktg4
    Vg5=Ag5*umax3*P/(Kp+P)*Ktg5


    //------------------------------------------------------

    Xmat=[
    0, V1, V2, Vg1, Vg2, Vg3, Vg4, Vg5, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0
    ]

    dX=[
    -sum( Xmat(1,:) ) + sum( Xmat(:,1) );
    -sum( Xmat(2,:) ) + sum( Xmat(:,2) );
    -sum( Xmat(3,:) ) + sum( Xmat(:,3) );
    -sum( Xmat(4,:) ) + sum( Xmat(:,4) );
    -sum( Xmat(5,:) ) + sum( Xmat(:,5) );
    -sum( Xmat(6,:) ) + sum( Xmat(:,6) );
    -sum( Xmat(7,:) ) + sum( Xmat(:,7) );
    -sum( Xmat(8,:) ) + sum( Xmat(:,8) );
    -sum( Xmat(9,:) ) + sum( Xmat(:,9) );
    ]

    //------------------------------------------------------

    Xdot=dX -X.*[D;D;Du;D;D;D;D;D;D] + Xin.*[D;D;Du;D;D;D;D;D;D]



    endfunction

    //　初期値


    //時間

    t0=365*0; t1=365 *Year; dt=0.1;
    t= t0 : dt : t1;

    //微分方程式solver
    XX=ode(Xs,t0,t,f);

    // グラフの表示
    clf(0);

    scf(0)
    // plot2d(t,XX(1,:),4);
    plot2d(t,XX(2,:),2);
    plot2d(t,XX(3,:),5);
    //-------------------------
    plot2d(t,XX(4,:),3,'g--');
    plot2d(t,XX(5,:),3,'g--');
    plot2d(t,XX(6,:),3,'g--');
    plot2d(t,XX(7,:),3,'g--');
    plot2d(t,XX(8,:),3,'g--');

    plot2d(t,XX(4,:)+XX(5,:)+XX(6,:)+XX(7,:)+XX(8,:),3);



    // legend(['P';'Asterionella';'Fragilaria';'Peridinium';'Cyanobacteria'],1);

    legend(['Keisou';'Uroglena';'Ryokusou'],1);

    scf(1)
    xset('color',2);//色
    param3d(XX(2,:),XX(3,:),XX(4,:)+XX(5,:)+XX(6,:)+XX(7,:)+XX(8,:),1,1,'Diatom@Uroglena@gleen-algae'); //x y z

    返信削除

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