写在前面：特征线法是基本，必须要懂。本身这个问题不涉及迭代，只是节点数据的传递，所以比较简单。本程序的思路为：设置两个二维矩阵H（x,t)、Q(x,t),先根据初始条件赋值，再计算内部节点和边界点的数据值。通过两个矩阵计算H、Q的值会更为直观。
clear clc format long; disp('管道瞬变流动分析.Page59-1.'); D=0.209; g=9.81; L=5000; H0=100; v0=0.5; i=0.002;%基本数据赋值。 a=1100;%赋值基本数据条件。 T=input('请输入总时间T的值（单位：s）:');%输入总时间值，可根据需要调整。 n=input('请输入步长的个数:');%输入步长个数，可根据需要调整。 dx=L/n; dt=dx/a; Ad=pi*D*D/4; Q0=v0*Ad; B=a/g/Ad; f=2*g*i*D/v0/v0; R=f*dx/2/g/D/Ad/Ad;%推导数据赋值。 for x=1:1:n+1 H(x,1)=H0-(x-1)*i*dx; Q(x,1)=Q0; end%设置初始条件 t=0; k=1; while (t t=t+dt; k=k+1; for x=2:1:n%计算内部节点的值 Cp=H(x-1,k-1)+B*Q(x-1,k-1)-R*Q(x-1,k-1)*abs(Q(x-1,k-1)); Cm=H(x+1,k-1)-B*Q(x+1,k-1)+R*Q(x+1,k-1)*abs(Q(x+1,k-1)); H(x,k)=(Cp+Cm)/2; Q(x,k)=(H(x,k)-Cm)/B; end H(1,k)=H0;%计算前边界H Q(1,k)=(H(1,k)-(H(2,k-1)-B*Q(2,k-1)+R*Q(2,k-1)*abs(Q(2,k-1))))/B;%计算前边界Q H(n+1,k)=H(n,k-1)+B*Q(n,k-1)-R*Q(n,k-1)*abs(Q(n,k-1));%计算后边界H Q(n+1,k)=0;%计算后边界Q end disp('计算已完成。'); disp('H、Q矩阵的行数，即空间节点数为：'); disp(size(H,1)) disp('H、Q矩阵的列数，即时间节点数为：'); disp(size(H,2)) m=input('请输入在允许范围内的空间节点序号：'); c=input('请输入在允许范围内的时间节点序号：'); disp('该状态下的H值和Q值分别为：') disp(H(m,c)) disp(Q(m,c)) disp('该状态点与上游端点的距离(m)为：') disp((m-1)*dx) disp('该状态点距阀门关闭的时间(s)为：') disp((c-1)*dt) [X,Y]=meshgrid([1:size(H,1)],[1:size(H,2)]); surf(X',Y',H) rotate3d%画H的三维图像