如果我正确理解了这个问题，我认为这可以通过从泄漏点左右遍历管道来实现。在每个点，将当前水位与管道标高进行比较，要么导致水位保持不变的淹没点，要么形成海滩和新的干峰。插值对于计算海滩的位置是必要的。
下面显示了一个实现。算法的大部分在
traverse函数中。希望评论提供充分的描述。
#!/usr/bin/python3
import numpy as np
import matplotlib.pyplot as pp
# Positions and elevations
n = 25
h = np.random.random((n, ))
x = np.linspace(0, 1, h.size)
# Index of leak
leak = np.random.randint(0, h.size)
# Traverse a pipe with positions (x) and elevations (h) from a leak index
# (leak) in a direction (step, +1 or -1). Return the positions of the changes
# in water level (y) the elevations at these changes (g) and the water level
# values (w).
def traverse(x, h, leak, step):
# End of the index range for the traversal
end = h.size if step == 1 else -1
# Initialise 1-element output arrays with values at the leak
y, g, w = [x[leak]], [h[leak]], [h[leak]]
# Loop from the index adjacent to the leak
for i in range(leak + step, end, step):
if w[-1] > h[i]:
# The new height is less than the old water level. Location i is
# submerged. No new points are created and the water level stays
# the same.
y.append(x[i])
g.append(h[i])
w.append(w[-1])
else:
# The new height is greater than the old water level. We have a
# "beach" and a "dry peak".
# ...
# Calculate the location of the beach as the position where the old
# water level intersects the pipe section from [i-step] to [i].
# This is added as a new point. The elevation and water level are
# the same as the old water level.
# ...
# The if statement is not strictly necessary. It just prevents
# duplicate points being generated.
if w[-1] != h[i-step]:
t = (w[-1] - h[i-step])/(h[i] - h[i-step])
b = x[i-step] + (x[i] - x[i-step])*t
y.append(b)
g.append(w[-1])
w.append(w[-1])
# ...
# Add the dry peak.
y.append(x[i])
g.append(h[i])
w.append(h[i])
# Convert from list to numpy array and return
return np.array(y), np.array(g), np.array(w)
# Traverse left and right
yl, gl, wl = traverse(x, h, leak, -1)
yr, gr, wr = traverse(x, h, leak, 1)
# Combine, reversing the left arrays and deleting the repeated start point
y = np.append(yl[:0:-1], yr)
g = np.append(gl[:0:-1], gr)
w = np.append(wl[:0:-1], wr)
# Output the total volume of water by integrating water level minus elevation
print('Total volume =', np.trapz(w - g, y), 'm^3 per unit cross sectional area')
# Display
pp.plot(x, h, '.-', label='elevation')
pp.plot(y, w, '.-', label='water level')
pp.plot([x[leak]], [h[leak]], 'o', label='leak')
pp.legend()
pp.show()
