额，首先，我是北化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: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));%S-P模型中BOD的变化规律，阿尔法α改为了a
end
a;
for j=1:N
    for i=1:N
     if j==i
        A(i,j)=1;%中间一个穿孔全是1
     else 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;
         end
      end
   end
end
for i=1:N
     for j=1:N
         if i==j
            B(i,j)=Q(i,1)/Q2(i,1);%B矩阵有效值的计算
         else B(i,j)=0;
         end
      end
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: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
for j=1:N
    for i=1:N
        if j==i-1
           D(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;
        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)=(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，新建两个文件，复制进去，然后命名保存到桌面，效果一模一样。

 调整到，交互页面。全选交互页面所有的代码，右击鼠标，选择第一个（evaluate selection）到这步已经完成了60%了，继续加油！（上图中红字请阅读）

MATLAB是无法识别中文的，这就是不能改文件名的缘由。

那可不可以自己起一个英文名？

理论上可以，但是我觉得你没学MATLAB看不懂报错，你改不回来，我不建议你改。

（程序猿的公理：代码能跑，就别轻易去动，否则你会哭的）

 现在，MATLAB已经蓄势待发了，就等待你提供数据，那么我们开始吧！

数据输入如下：

（注意看我输入的多个数据，用英文中括号包裹，数据中间用英文分号间隔开）

请输入河流断面个数n 5

请输入背景河水流量Q0(m3/s) 20

请输入在断面i处注入河流的污水流量Q(m3/s) [0.3;0.5;0.8;0.5;1]

请输入由断面i注入河流的污水的污染物（例如BOD）浓度L(mg/L) [150;150;200;200;150]

请输入由断面i注入河流的污水溶解氧（DO）浓度O(mg/L) [1;1;0;0;1]

请输入BOD的降解速度常数kd(/d) [0.1;0.12;0.15;0.18;0.15]

请输入大气复氧速度ka(/d) [0.16;0.2;0.28;0.25;0.3]

请输入由断面i到断面i+1经过时间t(d) [0.5;0.8;0.6;0.5;0.5]

请输入在断面i处引出的流水流量Q3(m3/s) [0.5;1;0.3;0;0.4]

请输入河水水温/℃ 25

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

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

回车键之后我们看看输出的六个矩阵是怎么样的？

输出如下：

剩下的时间就是老师询问了，你得回答出其中某一部分代码的含义，或者某个矩阵的运算方法。

加油，祝你满绩咯！

                                                                                                                               -----------2021.11.23晚

                                                                                                                                                  菩提树下呀