2010年11月28日日曜日

計算例（琵琶湖）



--- 珪藻
--- Urogrena
--- 緑藻








プログラム（コメント欄）
投稿者 時刻:

1 件のコメント:

  1. kf332010年11月28日 0:29

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

    umax1=0.21 ; umax2=0.223; umax3=0.2;
    umax4=umax3+0.0322;

    Ts1=10, Ts2=20, Ts3=20; Ts4=24;

    D=1/20;
    Year=8;

    Kp=0.1;

    Xs=[2; 0.1; 0.1; 0.05; 0.05];
    Xin=[2; 0; 0; 0; 0];


    function Xdot=f(t,X) ;

    P=X(1)
    A1=X(2)
    A2=X(3)
    A3=X(4)
    A4=X(5)

    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
    Kt3=(Temp/Ts3*exp((Ts3-Temp)/Ts3))^3
    Kt4=(Temp/Ts4*exp((Ts4-Temp)/Ts4))^3

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

    Xmat=[
    0, A1*umax1*P/(Kp+P)*Kt1,A2*umax2*P/(Kp+P)*Kt2,A3*umax3*P/(Kp+P)*Kt3,A4*umax4*P/(Kp+P)*Kt4;
    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) );
    ]

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

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

    // Xdot=[dX * (1-D) + Xin * D]

    endfunction

    //　初期値


    //時間

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

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

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

    scf(0)
    // plot2d(t,XX(1,:),4);
    plot2d(t,XX(2,:),2);
    plot2d(t,XX(3,:),5);
    //plot2d(t,XX(4,:),3);
    //plot2d(t,XX(5,:),3);
    plot2d(t,XX(4,:)+XX(5,:),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,:),1,1,'Diatom@Uroglena@gleen-algae'); //x y z

    返信削除

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