额，首先，我是北化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晚
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 菩提树下