一、龙格-库塔方法的公式推导

摘自《数值计算方法（MATLAB版）》

二、龙格-库塔方法在水库调洪演算的应用

一般库水面坡降很小，忽略动库容影响，近似看成静水面，水库蓄水量V只随坝前水位Z而变。若假定水库水位为水平起落，则水库调洪演算的实质，乃是对下列微分方程求解，即：
$$\frac{dV(Z)}{dt}=Q(t)-R(t)$$
若已知第n时段内的预报人库平均流量$Q_n$，n时段初的水$Z_{n-1}$与库容$V_{n-1}$，时段初的泄流量$R(Z_{n-1})$，泄流设备的开启状态，则应用定步长四阶龙格一库塔法求解，可得n时段末的库容$V_n$，即：
$$V_n=V_{n-1}+\frac{1}{6}[k_1+2(k_2+k_3)+k_4]$$
其中，
{ k 1 = Δ t { Q n − R [ Z ( V n − 1 ) ] } k 2 = Δ t { Q n − R [ Z ( V n − 1 + k 1 / 2 ) ] } k 3 = Δ t { Q n − R [ Z ( V n − 1 + k 2 / 2 ) ] } k 4 = Δ t { Q n − R [ Z ( V n − 1 + k 3 ) ] } \left\{\begin{align*} &k_1 = \Delta t\{ Q_n-R[Z(V_{n-1})]\} \\ &k_2 = \Delta t\{ Q_n-R[Z(V_{n-1}+k_1/2)]\} \\ &k_3 = \Delta t\{ Q_n-R[Z(V_{n-1}+k_2/2)]\}\\ &k_4 = \Delta t\{ Q_n-R[Z(V_{n-1}+k_3)]\} \end{align*}\right. ⎩ ⎨ ⎧​​k1​=Δt{ Qn​−R[Z(Vn−1​)]}k2​=Δt{ Qn​−R[Z(Vn−1​+k1​/2)]}k3​=Δt{ Qn​−R[Z(Vn−1​+k2​/2)]}k4​=Δt{ Qn​−R[Z(Vn−1​+k3​)]}​
式中
Z，由库容在水库水位一库容关系曲线上用三次自然样条插值法求得；
R，由水库水位及泄流设备的开启状态确定；
$\Delta t$，第n时段的时段长，$\Delta t$越小则计算精度越高，但计算时间相应加长。

求得$V_n$后，即可求得$V_n$，从而可得水库水位，库容，泄量随时间的变化过程。
调洪演算时段步长的选取可依据计算中要求达到的水位变幅精度确定，特别对复式泄洪建筑物或带有胸墙的泄水堰，需通过水位确定泄洪建筑物的开启及出流情况。

三、C-2 水库调洪演算的数值解程序

说明

这是水利水电工程设计计算程序集中调洪演算的数值计算方法。
已知水库的水位–水面面积关系，洪水过程线，对于每一种调洪方案（包括泄流条件、调洪方式、泄水建筑物参数）由调洪起始水位依次计算，直至洪水过程结束，计算机输出各时段末之水位、泄洪洞流量、溢洪道流量、水库出库总流量等。并用彩色曲线绘制洪水过程线、泄洪过程线和水库水位变化线。
这里已知水位面积曲线，递推采用水位~

C-2使用中应注意的问题

1、关于‘时段间隔T’
调洪演算的数值解从数学上说，是解一个微分方程，其自变量微分dt，具体到工程上，就是水文数据洪水过程线的时段间隔T（单位为秒），它的长短直接影响着调洪演算的成败。时段间隔选取的过大，往往成为计算失败的主要原因。有时计算失败，有时虽然能够计算通过，但从计算成果的图形上，可以看到跳跃式的锯齿形，其数据当然是不能用的了。从各地反馈来的数据中，有1800秒、3600秒，还有5400秒的。有的工程5400秒也能满足要求，可有的工程3600秒就不能满足要求。究其原因，还要看洪水过程线的起伏程度。变化剧烈的，就要求时段间隔小一些，变化不剧烈的，就可以允许时段间隔大一些。
如果时段间隔大了，可以将洪水过程线内插，同时时段间隔T减半。这并不影响洪水过程线的精度，是一个有效的解决办法。这个工作在Excel表上做起来，非常容易。

2、充分利用C-2程序设置的泄流条件。
以前在用图解法、半图解法做调洪演算时，都要制作工程的泄流曲线。所谓泄流曲线，就是对水库中特定的泄流建筑物（溢洪道、泄洪洞等），特定的尺寸，特定的底高程，特定的运行方式，计算水位与泄流量的对应关系，这就是泄流曲线，以便在具体调洪时使用。只要泄流条件有任何变动，泄流曲线的数值就要改变。这也是过去调洪演算比较麻烦，工作量特大的原因之一。C-2程序将这部分工作都用具体的参数（尺寸、流量系数等）所取代。你只要改动其中的任何参数，都会形成不同的调洪方案。这对于优化设计方案非常有用，改变几个数，再计算一遍，不过是几秒的时间，应该充分利用。C-2程序也设置了泄流曲线的输入方法，那是为校核过去手算已经有数据的工程而用的，新的工程设计没有必要自找这个麻烦。

附 龙格-库塔python代码

    def runge_kutta(self, f, x0, y0, h, n):
        """
        使用龙格-库塔方法求解一阶常微分方程的数值解
        :param f: 函数 f(x, y) 的句柄，表示 y' = f(x, y)
        :param x0: 自变量的初始值
        :param y0: 因变量的初始值
        :param h: 步长
        :param n: 迭代次数
        :return: 返回一个二元组，包含自变量和因变量的数组
        """
        x = np.zeros(n + 1)
        y = np.zeros(n + 1)
        x[0] = x0
        y[0] = y0
        for i in range(n):
            k1 = f(x[i], y[i])
            k2 = f(x[i] + h / 2, y[i] + h / 2 * k1)
            k3 = f(x[i] + h / 2, y[i] + h / 2 * k2)
            k4 = f(x[i] + h, y[i] + h * k3)
            y[i + 1] = y[i] + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
            x[i + 1] = x[i] + h
        return x, y

参考文献

• 数值计算方法（MATLAB版)，马昌凤、柯艺芬编著
• 水库信息化工程理论与实践，高英利，韩建忠编著
• 水库调洪演算的数值解程序，张校正