北化环工多河段BOD—DO耦合矩阵模型

本文主要是介绍北化环工多河段BOD—DO耦合矩阵模型，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

额，首先，我是北化19级的老学长，你如果是环工，那你寻找的就是这份代码了。

鄙人才疏学浅，C语言不是很懂，而且听说440行，我是直接放弃了。

矩阵模型自然是MATLAB最适合了，MATLAB的翻译就是：矩阵实验室（Matrix Laboratory）

所以我下载了MATLAB，学习了3天。（是我的舍友一直不断的教我）

假设没有MATLAB，那就只能是去机房了，机房有MATLAB。（能用，就是很慢，是2014版本）

废话不多说，上代码。（最好带自己的电脑过去）

第一部分的代码是设计运算的，使用《环境系统分析教程》第二版程声通。第82-84页。按照逻辑和运算顺序输入算法（没学过的可以去菜鸟教程学习一下基础语法，看懂每行作用，老师会问）

代码如下：calculate模块

%%%%%%%%%% 以下是“多河段 BOD-DO 耦合矩阵模型”的程序代码%%%%%%%%%%
function F=calculate(N,Q1,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
% 输入题目中变量值Os=468/(31.6+T);%饱和氧算法
for i=1:1:NQ2(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:Na(i,1)=exp(-kd(i,1)*t(i,1));%S-P模型中BOD的变化规律，阿尔法α改为了a
end
a;
for j=1:Nfor i=1:Nif j==iA(i,j)=1;%中间一个穿孔全是1else if j==i-1%穿孔朝下是正常运算的值A(i,j)=-a(i,1)*(Q1(i,1)-Q3(i,1))/Q2(i,1);%A矩阵有效值的计算else A(i,j)=0;endendend
end
for i=1:Nfor j=1:Nif i==jB(i,j)=Q(i,1)/Q2(i,1);%B矩阵有效值的计算else B(i,j)=0;endend
end
A;
B;
g=zeros(N,1);
g(1,1)=a(1,1)*(Q1(1,1)-Q3(1,1))/Q2(1,1)*L20;% 零维矩阵的建立，书上公式
g;
for i=1:N
y(i,1)=exp(-ka(i,1)*t(i,1));
end
y;
for j=1:Nfor i=1:Nif j==iC(i,j)=1;else if j==i-1C(i,j)=(Q1(i,1)-Q3(i,1))/Q2(i,1)*-y(i,1);else C(i,j)=0;endendend
end
for j=1:Nfor i=1:Nif j==i-1D(i,j)=(Q1(i,1)-Q3(i,1))/Q2(i,1)*kd(i,1)/(ka(i,1)-kd(i,1))*(a(i,1)-y(i,1));else D(i,j)=0;endend
end
C;
D;
for i=1:Nf(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)=(Q1(1,1)-Q3(1,1))/Q2(1,1)*y(1,1)*O20-(Q1(1,1)-Q3(1,1))/Q2(1,1)*kd(1,1)/(ka(1,1)-kd(1,1))*(a(1,1)-y(1,1))*L20;
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+h)-C^-1*D*A^-1*g;
n;
L2=U*L+m;
L2;
O2=V*L+n;
O2;
res=cell(1,6);
res{1}=U;
res{2}=V;
res{3}=m;
res{4}=n;
res{5}=L2;
res{6}=O2;
F=res;

第二部分就是交互阶段，MATLAB已经有了一套完整的算法了，现在需要你输入一些实际的工程值进行模拟运算。此时，加入能够输入数值的窗口。

代码如下：interact模块

N=input('请输入河流断面个数n ');
Q1(1,1)=input('请输入背景河水流量Q0(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) '); 
res=calculate(N,Q1,Q,L,O,kd,ka,t,Q3,T,L20,O20);
fprintf('U响应矩阵：\n');
disp(res{1});
fprintf('V相应矩阵：\n');
disp(res{2});
fprintf('m向量：\n');
disp(res{3});
fprintf('n向量：\n');
disp(res{4});
fprintf('各断面BOD浓度：\n');
disp(res{5});
fprintf('各断面DO浓度：\n');
disp(res{6});

那我们看看实际的输入情况吧，毕竟输入的格式也是非常的有要求，这也是MATLAB门槛所在。

上图：

文件最好放置在桌面，这样也非常方便找到。基本上完成这个操作就完成了40%了，继续加油！

没有也简单，打开MATLAB，新建两个文件，复制进去，然后命名保存到桌面，效果一模一样。