环评3c语言编写多河段bod-do耦合矩阵模型源程序,BOD-DO矩阵耦合模型 matlab mfiles...

用matlab编写的环境评价与系统分析中重要的 BOD-DO矩阵耦合模型。才疏学浅,请多多指教!

%%%%%%%%%% 以下是“多河段 BOD-DO 耦合矩阵模型”的程序代码%%%%%%%%%%

function F=pi(Q,L,O,kd,ka,t,Q3,T,L20,O20)

% Q-在断面i处注入河流的流量,m3/s

% L，O-由断面i注入河流的污水的污染物(例如BOD)浓度与溶解氧(DO)浓度,mg/L

% kd-BOD的降解速度常数,/d

% ka-大气复氧速度,/d

% t-由断面i到断面i+1经过时间，d

% Q3-在断面i处引出的流水流量,m3/s

% L20,O20-河流背景的BOD和DO浓度，mg/L

% 输入题目中变量值

n=input('请输入河流断面个数n ')

Q1(1,1)=input('请输入背景河水流量Q20(m3/s) ')

Q=input('请输入在断面i处注入河流的污水流量Q(m3/s) ')

L=input('请输入由断面i注入河流的污水的污染物(例如BOD)浓度L(mg/L) ')

O=input('请输入由断面i注入河流的污水溶解氧(DO)浓度O(mg/L) ')

kd=input('请输入BOD的降解速度常数kd(/d) ')

ka=input('请输入大气复氧速度ka(/d) ')

t=input('请输入由断面i到断面i+1经过时间t(d) ')

Q3=input('请输入在断面i处引出的流水流量Q3(m3/s) ')

T=input('请输入河水水温/℃ ')

L20=input('请输入河流背景的BOD浓度L20(mg/L) ')

O20=input('请输入河流背景的DO浓度O20(mg/L) ')

Os=468/(31.6+T) % Os-河流的饱和溶解氧的值 mg/L

Os

for i=1:1:n

Q2(i,1)=Q1(i,1)-Q3(i,1)+Q(i,1); % 根据连续性原理，写出每一个断面的流量Q与BOD的平衡关系

Q1(i+1,1)=Q2(i,1);

end

Q1 % Q1i-由上一个河段流到断面i的河水流量，mg/L

Q2 % Q2i-由断面i向下游河段流出的河水流量，mg/L

Q3

for i=1:1:n

a(i,1)=exp(-kd(i,1)*t(i,1));

end

a

% 计算A矩阵

for j=1:n

for i=1:n

if j==i

A(i,j)=1;

else if j==i-1

A(i,j)=-a(i,1)*(Q1(i,1)-Q3(i,1))/Q2(i,1);

else A(i,j)=0;

end

end

end

end

% 计算B矩阵

for i=1:n

for j=1:n

if i==j

B(i,j)=Q(i,1)/Q2(i,1);

else B(i,j)=0;

end

end

end

A

B

g=zeros(n,1);

g(1,1)=a(1,1)*L20;

g

% BOD-DO耦合矩阵模型

for i=1:n

y(i,1)=exp(-ka(i,1)*t(i,1));

end

y;

% 计算C矩阵

for j=1:n

for i=1:n

if j==i

C(i,j)=1;

else if j==i-1

C(i,j)=(Q1(i,1)-Q3(i,1))/Q2(i,1)*-y(i,1);

else C(i,j)=0;

end

end

end

end

%计算D矩阵

for j=1:n

for i=1:n

if j==i-1

D(i,j)=kd(i,1)/(ka(i,1)-kd(i,1))*(a(i,1)-y(i,1));

else D(i,j)=0;

end

end

end

C

D

for i=1:n

f(i,1)=(Q1(i,1)-Q3(i,1))/Q2(i,1)*Os*(1-y(i,1));

end

f

h=zeros(n,1);

h(1,1)=y(1,1)*O20-kd(1,1)/(ka(1,1)-kd(1,1))*(a(1,1)-y(1,1));

h

U=A^-1*B; %U是BOD对BOD响应矩阵

U

V=-C^-1*D*A^-1*B; %V是溶解氧对BOD的响应

矩阵

V

m=A^-1*g;

m

n=C^-1*B*O+C^-1*(f+g)-C^-1*D*A^-1*g;

n

L2=U*L+m;

L2

O2=V*L+n;

O2

%%%%%输入值%%%%%

%Q=[0.3;0.5;0.8;0.5;1]

%L=[150