2010年10月16日土曜日

藻ミジンコver7.2

・池に、藻類、ミジンコ、魚がいる。また、底の泥に、タニシ、ユスリカがいる。
・リンは１ppm。
・生物が死んだら、底の泥となる。
・外部からの流入は無し。
青　：リン
緑　：藻類
紫　：泥
水色：ミジンコ

※リンは1/40（mg/L）
　リン以外は、全て炭素（mg/L）
　タニシ、ユスリカの初期値はミジンコの１０分の１、魚の初期値はミジンコの１００分の１と小さく設定したため、グラフ上にはほとんど現われていない。

プログラムはコメント欄に記載
投稿者 時刻:

2 件のコメント:

  1. kf332010年10月16日 22:53

    //藻ミジンコver7.2
    clear;
    xbasc;

    //--------------------------------------------------------------[ 初期値 ]
    // [ 変数 ]
    P=2; A=1; Z=0.1; F=0.001; D=(1/10000); T=0.01; Y=0.01; POM=(1/10000); // P：栄養塩（リン）※リンはCOD換算（リン0.2mgのとき、4mgと表記）、A:藻類、Z：ミジンコ、F：魚、D：底泥のCOD量 (mg/L)、T：タニシ、Y：ユスリカ、POM：懸濁有機物（魚の糞等）

    // [ 定数 ]
    V=1; //V:湖沼の単位体積(1L)
    DA=100; DZ=60; DF=1000; DT=900; DY=900; //DA：藻類の寿命(水深m/沈降深m・日), DZ：ミジンコの寿命, DF：魚の寿命，DT：タニシの寿命、DY：ユスリカの寿命（日）
    qZ=0.01; qF=0.1; qDF=0.1; //qZ：ミジンコ１mg/Lあたりのろ過速度、qF：魚１gm/Lあたりのろ過速度(ろ過L/L・日)、qDF：魚１gm/Lあたりの泥ろ過速度(ろ過L/L・日)
    qT=0.01; qDT=0.01; qDY=0.01;//qT：タニシの1mg/Lあたりのろ過速度(ろ過L/L・日)、qDT：タニシの1mg/Lあたりの泥摂食速度(泥g/日)、QDY：ユスリカの1mg/Lあたりの泥摂食速度(泥g/日)(ろ過L/L・日)
    K=0.4;//K：藻類の栄養塩（リン）吸収飽和定数
    Mmax=0.2;//Mmax：藻類の最大比増殖速度
    EfZ=0.05; EfF=0.05; EfT=0.05;//EfZ：ミジンコが食べた藻類の変換効率、EfF：魚が食べたミジンコ、ユスリカの変換効率、EfT：タニシが食べた藻類の変換効率
    EfDT=0.01; EfDY=0.01;//EfDT：タニシが食べた泥の変換効率、EfDY：ユスリカが食べた泥の変換効率
    Elution=0.01;//Elution：泥からの栄養塩(リン)溶出率 (※リンはCOD換算)
    Year=0.5;//シミュレーション期間（年）

    x=[P
    A
    Z
    F
    D
    T
    Y
    POM];



    //--------------------------------------------------------------[　関数　]

    function xdot=func(t,x);
    P=x(1,1); A=x(2,1); Z=x(3,1); F=x(4,1); D=x(5,1); T=x(6,1); Y=x(7,1); POM=x(8,1);
    Temp=10
    //-----------------------------------------------------------------------------------------------------------------------
    // P A Z F D T Y POM
    //-----------------------------------------------------------------------------------------------------------------------
    xmat=[ 0, A*Mmax*P/(K+P)*(Temp/10*exp((10-Temp)/10))^3, 0, 0, 0, 0, 0, 0;
    0, 0, A*Z*qZ/V*EfF, 0, A/DA, A*T*qT/V*EfT, 0, A*Z*qZ/V*(1-EfZ)+A*T*qT/V*(1-EfT);
    0, 0, 0, Z*F*qF/V*EfF, Z/DZ, 0, 0, Z*F*qF/V*(1-EfF);
    0, 0, 0, 0, F/DF, 0, 0, 0;
    Elution*D, 0, 0, 0, 0, T*D*qDT/V*EfDT, Y*D*qDY/V*EfDY, T*D*qDT/V*(1-EfDT)+Y*D*qDY/V*(1-EfDY);
    0, 0, 0, 0, T/DT, 0, 0, 0;
    0, 0, 0, Y*F*qDF/V*EfF, Y/DY, 0, 0, F*Y*qDF/V*(1-EfF);
    0, 0, 0, 0, POM, 0, 0, 0; ]
    //-----------------------------------------------------------------------------------------------------------------------

    xdot(1)=-sum(xmat(1,:)) + sum(xmat(:,1));
    xdot(2)=-sum(xmat(2,:)) + sum(xmat(:,2));
    xdot(3)=-sum(xmat(3,:)) + sum(xmat(:,3));
    xdot(4)=-sum(xmat(4,:)) + sum(xmat(:,4));
    xdot(5)=-sum(xmat(5,:)) + sum(xmat(:,5));
    xdot(6)=-sum(xmat(6,:)) + sum(xmat(:,6));
    xdot(7)=-sum(xmat(7,:)) + sum(xmat(:,7));
    xdot(8)=-sum(xmat(8,:)) + sum(xmat(:,8));

    endfunction

    //--------------------------------------------------------------[　計算　]
    ta=0;tstep=0.01;tb=365*Year;
    t=ta:tstep:tb;
    xx=ode([P;A;Z;F;D;T;Y;POM], 0, t, func);
    //--------------------------------------------------------------[プロット]
    clf();
    xset('color',1);//色
    plot2d(t,xx(1,:),2); //P
    plot2d(t,xx(2,:),3); //A
    plot2d(t,xx(3,:),4); //Z
    plot2d(t,xx(4,:),5); //F
    plot2d(t,xx(5,:),6); //D
    plot2d(t,xx(6,:),7); //T
    plot2d(t,xx(7,:),8); //Y

    返信削除
  2. 匿名2012年11月28日 9:27

    Thank you for thе good writeup. It in fаct was
    а аmusement account it. Look аdvanceԁ tο far
    аdded agreeable from you! Bу the way, hoω can
    we communicate?
    My page :: pikavippii

    返信削除

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