clear;xbasc;//--------------------------------------------------------------[ 初期値 ]// [ 変数 ]P=2; A1=0.1; A2=0.1; Z1=0.01; Z2=0.01; F1=0.0001; F2=0.0002; D=0.001; POM=0.001; // P：栄養塩（リン）※リンはCOD換算（リン0.2mgのとき、4mgと表記）、// A：藻類、Z：ミジンコ、F：魚、D：底泥のCOD量 (mg/L)、POM：懸濁有機物（魚の糞等）x=[ P A1 A2 Z1 Z2 F1 F2 POM D]; //--------------------------------------------------------------[ 関数 ]function xdot=func(t,x); P=x(1,1); A1=x(2,1); A2=x(3,1); Z1=x(4,1); Z2=x(5,1); F1=x(6,1); F2=x(7,1); POM=x(8,1); D=x(9,1); Temp=12*sin(2*3.1415*t/365-5.2)+19//-----------------------------------------------------------------------------------------------------------------------// P A1 A2 Z1 Z2 F1 F2 POM D //-----------------------------------------------------------------------------------------------------------------------xmat=[ 0, A1*0.2*P/(0.1+P), A2*0.3*P/(0.3+P), 0, 0, 0, 0, 0, 0; 0, 0, 0, A1*Z1*0.1*0.1, 0, 0, 0, (A1*Z1*0.1*0.9) + A1*(1/300), 0; 0, 0, 0, A2*Z1*0.1*0.1, A2*Z2*0.1*0.1, 0, 0, (A2*Z1*0.1*0.9) + A2*(1/300), 0; 0, 0, 0, 0, 0, Z1*F1*0.1*0.1, Z1*F2*0.1*0.1,(Z1*F1*0.1*0.9)+(Z1*F2*0.1*0.9)+Z1*(1/90), 0; 0, 0, 0, 0, 0, 0, Z2*F1*0.1*0.1, (Z2*F1*0.1*0.9)+Z2*(1/90), 0; 0, 0, 0, 0, 0, 0, 0, F1*(1/900), 0; 0, 0, 0, 0, 0, F2*F1*0.1*0.1, 0, (F2*F1*0.1*0.9)+F2*(1/900), 0; 0, 0, 0, 0, 0, 0, 0, 0, POM*0.5; D*(1/100), 0, 0, 0, 0, 0, 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));xdot(9)=-sum(xmat(9,:)) + sum(xmat(:,9)); endfunction //--------------------------------------------------------------[ 計算 ]ta=0;tstep=0.1;tb=365*2;t=ta:tstep:tb;xx=ode([P;A1;A2;Z1;Z2;F1;F2;POM;D], 0, t, func);//--------------------------------------------------------------[プロット]clf();subplot(1,3,1); plot2d(t,xx(1,:),4); //Pplot2d(t,xx(2,:),5); //A1plot2d(t,xx(3,:),3); //A2subplot(1,3,2); plot2d(t,xx(4,:),2); //Z1plot2d(t,xx(5,:),6); //Z2legend('Zooplankton(Large)','Zooplankton(Small)',1);subplot(1,3,3); plot2d(t,xx(6,:),5); //F1plot2d(t,xx(7,:),7); //F2legend('Fish(Himemasu)','Fish(Wakasagi)',1);subplot(1,3,1); plot2d(t,xx(8,:),1); //POMplot2d(t,xx(9,:),7,); //Dlegend('P','Algae(Diatom)','Algae(Green Algae)','POM','Sediment',1);
// Ver8_2clear;xbasc;//--------------------------------------------------------------[ 初期値 ]// [ 変数 ]P=2; A1=0.1; A2=0.1; Z1=0.1; Z2=0.1; F1=0.005; F2=0.01; D=0.001; POM=0.001; // P：栄養塩（リン）※リンはCOD換算（リン0.2mgのとき、4mgと表記）、// A：藻類、Z：ミジンコ、F：魚、D：底泥のCOD量 (mg/L)、POM：懸濁有機物（魚の糞等）x=[ P A1 A2 Z1 Z2 F1 F2 POM D]; //--------------------------------------------------------------[ 関数 ]function xdot=func(t,x); P=x(1,1); A1=x(2,1); A2=x(3,1); Z1=x(4,1); Z2=x(5,1); F1=x(6,1); F2=x(7,1); POM=x(8,1); D=x(9,1); Temp=12*sin(2*3.1415*t/365-5.2)+19//-----------------------------------------------------------------------------------------------------------------------// P A1 A2 Z1 Z2 F1 F2 POM D //-----------------------------------------------------------------------------------------------------------------------xmat=[ 0, A1*0.2*P/(0.2+P), A2*0.3*P/(0.3+P), 0, 0, 0, 0, 0, 0; 0, 0, 0, A1*Z1*0.1*0.1, 0, 0, 0, (A1*Z1*0.1*0.9) + A1*(1/300), 0; 0, 0, 0, A2*Z1*0.1*0.1, A2*Z2*0.1*0.1, 0, 0, (A2*Z1*0.1*0.9)+(A2*Z2*0.1*0.9) + A2*(1/300), 0; 0, 0, 0, 0, 0, Z1*F1*0.1*0.1, Z1*F2*0.1*0.1,(Z1*F1*0.1*0.9)+(Z1*F2*0.1*0.9)+Z1*(1/120), 0; 0, 0, 0, 0, 0, 0, Z2*F1*0.1*0.1, (Z2*F1*0.1*0.9)+Z2*(1/120), 0; 0, 0, 0, 0, 0, 0, 0, F1*(1/900), 0; 0, 0, 0, 0, 0, F2*F1*0.05*0.1, 0, (F2*F1*0.05*0.9)+F2*(1/900), 0; 0, 0, 0, 0, 0, 0, 0, 0, POM*0.5; D*(1/100), 0, 0, 0, 0, 0, 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));xdot(9)=-sum(xmat(9,:)) + sum(xmat(:,9)); endfunction //--------------------------------------------------------------[ 計算 ]ta=0;tstep=0.1;tb=365*5;t=ta:tstep:tb;xx=ode([P;A1;A2;Z1;Z2;F1;F2;POM;D], 0, t, func);//--------------------------------------------------------------[プロット]clf();subplot(1,3,1); plot2d(t,xx(1,:),4); //Pplot2d(t,xx(2,:),5); //A1plot2d(t,xx(3,:),3); //A2subplot(1,3,2); plot2d(t,xx(4,:),2); //Z1plot2d(t,xx(5,:),6); //Z2legend('Zooplankton(Large)','Zooplankton(Small)',1);subplot(1,3,3); plot2d(t,xx(6,:),5); //F1plot2d(t,xx(7,:),7); //F2legend('Fish(Himemasu)','Fish(Wakasagi)',1);subplot(1,3,1); plot2d(t,xx(8,:),1); //POMplot2d(t,xx(9,:),7,); //Dlegend('P','Algae(Diatom)','Algae(Green Algae)','POM','Sediment',1);
clear;
返信削除xbasc;
//--------------------------------------------------------------[ 初期値 ]
// [ 変数 ]
P=2; A1=0.1; A2=0.1; Z1=0.01; Z2=0.01; F1=0.0001; F2=0.0002; D=0.001; POM=0.001;
// P：栄養塩（リン）※リンはCOD換算（リン0.2mgのとき、4mgと表記）、
// A：藻類、Z：ミジンコ、F：魚、D：底泥のCOD量 (mg/L)、POM：懸濁有機物（魚の糞等）
x=[ P
A1
A2
Z1
Z2
F1
F2
POM
D];
//--------------------------------------------------------------[ 関数 ]
function xdot=func(t,x);
P=x(1,1); A1=x(2,1); A2=x(3,1); Z1=x(4,1); Z2=x(5,1); F1=x(6,1); F2=x(7,1); POM=x(8,1); D=x(9,1);
Temp=12*sin(2*3.1415*t/365-5.2)+19
//-----------------------------------------------------------------------------------------------------------------------
// P A1 A2 Z1 Z2 F1 F2 POM D
//-----------------------------------------------------------------------------------------------------------------------
xmat=[ 0, A1*0.2*P/(0.1+P), A2*0.3*P/(0.3+P), 0, 0, 0, 0, 0, 0;
0, 0, 0, A1*Z1*0.1*0.1, 0, 0, 0, (A1*Z1*0.1*0.9) + A1*(1/300), 0;
0, 0, 0, A2*Z1*0.1*0.1, A2*Z2*0.1*0.1, 0, 0, (A2*Z1*0.1*0.9) + A2*(1/300), 0;
0, 0, 0, 0, 0, Z1*F1*0.1*0.1, Z1*F2*0.1*0.1,(Z1*F1*0.1*0.9)+(Z1*F2*0.1*0.9)+Z1*(1/90), 0;
0, 0, 0, 0, 0, 0, Z2*F1*0.1*0.1, (Z2*F1*0.1*0.9)+Z2*(1/90), 0;
0, 0, 0, 0, 0, 0, 0, F1*(1/900), 0;
0, 0, 0, 0, 0, F2*F1*0.1*0.1, 0, (F2*F1*0.1*0.9)+F2*(1/900), 0;
0, 0, 0, 0, 0, 0, 0, 0, POM*0.5;
D*(1/100), 0, 0, 0, 0, 0, 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));
xdot(9)=-sum(xmat(9,:)) + sum(xmat(:,9));
endfunction
//--------------------------------------------------------------[ 計算 ]
ta=0;tstep=0.1;tb=365*2;
t=ta:tstep:tb;
xx=ode([P;A1;A2;Z1;Z2;F1;F2;POM;D], 0, t, func);
//--------------------------------------------------------------[プロット]
clf();
subplot(1,3,1);
plot2d(t,xx(1,:),4); //P
plot2d(t,xx(2,:),5); //A1
plot2d(t,xx(3,:),3); //A2
subplot(1,3,2);
plot2d(t,xx(4,:),2); //Z1
plot2d(t,xx(5,:),6); //Z2
legend('Zooplankton(Large)','Zooplankton(Small)',1);
subplot(1,3,3);
plot2d(t,xx(6,:),5); //F1
plot2d(t,xx(7,:),7); //F2
legend('Fish(Himemasu)','Fish(Wakasagi)',1);
subplot(1,3,1);
plot2d(t,xx(8,:),1); //POM
plot2d(t,xx(9,:),7,); //D
legend('P','Algae(Diatom)','Algae(Green Algae)','POM','Sediment',1);
// Ver8_2
返信削除clear;
xbasc;
//--------------------------------------------------------------[ 初期値 ]
// [ 変数 ]
P=2; A1=0.1; A2=0.1; Z1=0.1; Z2=0.1; F1=0.005; F2=0.01; D=0.001; POM=0.001;
// P：栄養塩（リン）※リンはCOD換算（リン0.2mgのとき、4mgと表記）、
// A：藻類、Z：ミジンコ、F：魚、D：底泥のCOD量 (mg/L)、POM：懸濁有機物（魚の糞等）
x=[ P
A1
A2
Z1
Z2
F1
F2
POM
D];
//--------------------------------------------------------------[ 関数 ]
function xdot=func(t,x);
P=x(1,1); A1=x(2,1); A2=x(3,1); Z1=x(4,1); Z2=x(5,1); F1=x(6,1); F2=x(7,1); POM=x(8,1); D=x(9,1);
Temp=12*sin(2*3.1415*t/365-5.2)+19
//-----------------------------------------------------------------------------------------------------------------------
// P A1 A2 Z1 Z2 F1 F2 POM D
//-----------------------------------------------------------------------------------------------------------------------
xmat=[ 0, A1*0.2*P/(0.2+P), A2*0.3*P/(0.3+P), 0, 0, 0, 0, 0, 0;
0, 0, 0, A1*Z1*0.1*0.1, 0, 0, 0, (A1*Z1*0.1*0.9) + A1*(1/300), 0;
0, 0, 0, A2*Z1*0.1*0.1, A2*Z2*0.1*0.1, 0, 0, (A2*Z1*0.1*0.9)+(A2*Z2*0.1*0.9) + A2*(1/300), 0;
0, 0, 0, 0, 0, Z1*F1*0.1*0.1, Z1*F2*0.1*0.1,(Z1*F1*0.1*0.9)+(Z1*F2*0.1*0.9)+Z1*(1/120), 0;
0, 0, 0, 0, 0, 0, Z2*F1*0.1*0.1, (Z2*F1*0.1*0.9)+Z2*(1/120), 0;
0, 0, 0, 0, 0, 0, 0, F1*(1/900), 0;
0, 0, 0, 0, 0, F2*F1*0.05*0.1, 0, (F2*F1*0.05*0.9)+F2*(1/900), 0;
0, 0, 0, 0, 0, 0, 0, 0, POM*0.5;
D*(1/100), 0, 0, 0, 0, 0, 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));
xdot(9)=-sum(xmat(9,:)) + sum(xmat(:,9));
endfunction
//--------------------------------------------------------------[ 計算 ]
ta=0;tstep=0.1;tb=365*5;
t=ta:tstep:tb;
xx=ode([P;A1;A2;Z1;Z2;F1;F2;POM;D], 0, t, func);
//--------------------------------------------------------------[プロット]
clf();
subplot(1,3,1);
plot2d(t,xx(1,:),4); //P
plot2d(t,xx(2,:),5); //A1
plot2d(t,xx(3,:),3); //A2
subplot(1,3,2);
plot2d(t,xx(4,:),2); //Z1
plot2d(t,xx(5,:),6); //Z2
legend('Zooplankton(Large)','Zooplankton(Small)',1);
subplot(1,3,3);
plot2d(t,xx(6,:),5); //F1
plot2d(t,xx(7,:),7); //F2
legend('Fish(Himemasu)','Fish(Wakasagi)',1);
subplot(1,3,1);
plot2d(t,xx(8,:),1); //POM
plot2d(t,xx(9,:),7,); //D
legend('P','Algae(Diatom)','Algae(Green Algae)','POM','Sediment',1);