python---龙格库塔法

原创

已于 2026-01-21 17:04:42 修改 · 1.1k 阅读

标签

#算法

于 2026-01-20 16:34:32 首次发布

一、简介

龙格库塔方法（Runge-Kutta methods）是一类迭代方法，用于求解常微分方程的数值解。它是数值分析中最重要的方法之一。

常见类型

欧拉法（一阶龙格库塔）
中点法（二阶龙格库塔）
经典龙格库塔法（RK4，四阶）

RK4（经典四阶龙格库塔）

最常用的是四阶龙格库塔方法，其局部截断误差为O(h⁵)，全局误差为O(h⁴)。

公式如下：

k1 = f(tn, yn)
k2 = f(tn + h/2, yn + h*k1/2)
k3 = f(tn + h/2, yn + h*k2/2)
k4 = f(tn + h, yn + h*k3)

yn+1 = yn + h/6 * (k1 + 2k2 + 2k3 + k4)

基本实现：

示例1: 求解 dy/dt = y, y(0)=1 (解析解: y=e^t)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp  # 用于比较

def runge_kutta_4(f, y0, t_span, h):
    """
    四阶龙格库塔方法
    
    参数:
        f: 微分方程 dy/dt = f(t, y)
        y0: 初始值
        t_span: 时间范围 [t0, tf]
        h: 步长
    
    返回:
        t: 时间点数组
        y: 解数组
    """
    t0, tf = t_span
    n_steps = int((tf - t0) / h)
    
    t = np.linspace(t0, tf, n_steps + 1)
    y = np.zeros(n_steps + 1)
    y[0] = y0
    
    for i in range(n_steps):
        k1 = f(t[i], y[i])
        k2 = f(t[i] + h/2, y[i] + h*k1/2)
        k3 = f(t[i] + h/2, y[i] + h*k2/2)
        k4 = f(t[i] + h, y[i] + h*k3)
        
        y[i+1] = y[i] + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
    
    return t, y

# 示例1: 求解 dy/dt = y, y(0)=1 (解析解: y=e^t)
def f1(t, y):
    return y

# 使用RK4求解
t_rk4, y_rk4 = runge_kutta_4(f1, 1, [0, 2], 0.1)

# 与解析解比较
y_exact = np.exp(t_rk4)

# 使用SciPy的solve_ivp作为参考
sol = solve_ivp(f1, [0, 2], [1], method='RK45', t_eval=t_rk4)

# 绘图
plt.figure(figsize=(10, 6))
plt.plot(t_rk4, y_rk4, 'b-', linewidth=2, label='RK4_result')
plt.plot(t_rk4, y_exact, 'r--', linewidth=2, label='Analytical_result $e^t$')
plt.plot(t_rk4, sol.y[0], 'g:', linewidth=2, label='SciPy RK45')
plt.xlabel('time t')
plt.ylabel('y(t)')
plt.title('runge_kutta dy/dt = y')
plt.legend()
plt.grid(True, alpha=0.3)

# 保存图像
plt.savefig('runge_kutta_example1.png', dpi=300, bbox_inches='tight')
plt.savefig('runge_kutta_example1.pdf', bbox_inches='tight')  # 保存为矢量图
print("图像已保存为 'runge_kutta_example1.png' 和 'runge_kutta_example1.pdf'")

plt.show()

# 计算误差
error_rk4 = np.abs(y_rk4 - y_exact)
error_scipy = np.abs(sol.y[0] - y_exact)

print(f"\n误差分析:")
print(f"RK4最大误差: {np.max(error_rk4):.6f}")
print(f"SciPy最大误差: {np.max(error_scipy):.6f}")
print(f"RK4平均误差: {np.mean(error_rk4):.6f}")
print(f"SciPy平均误差: {np.mean(error_scipy):.6f}")

运行结果：

二、自适应龙格库塔法

自适应龙格库塔法通过动态调整步长来控制误差：

核心思想：使用两个不同阶数的龙格库塔方法计算同一个步长的结果，通过比较两者的差异来估计误差
步长控制：根据误差估计自动调整步长（误差大时减小步长，误差小时增大步长）
常用方法：Runge-Kutta-Fehlberg (RKF45)、Dormand-Prince (DOPRI5)

RKF45 方法（Runge-Kutta-Fehlberg）

RKF45使用一个4阶方法和一个5阶方法来估计误差，计算公式如下：

k1 = h * f(t, y)
k2 = h * f(t + h/4, y + k1/4)
k3 = h * f(t + 3h/8, y + (3k1 + 9k2)/32)
k4 = h * f(t + 12h/13, y + (1932k1 - 7200k2 + 7296k3)/2197)
k5 = h * f(t + h, y + (439k1/216 - 8k2 + 3680k3/513 - 845k4/4104))
k6 = h * f(t + h/2, y - 8k1/27 + 2k2 - 3544k3/2565 + 1859k4/4104 - 11k5/40)

# 四阶估计
y4 = y + 25k1/216 + 1408k3/2565 + 2197k4/4104 - k5/5

# 五阶估计
y5 = y + 16k1/135 + 6656k3/12825 + 28561k4/56430 - 9k5/50 + 2k6/55

# 误差估计
error = |y5 - y4|

实现代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import time

class AdaptiveRungeKutta:
    """
    Adaptive Runge-Kutta-Fehlberg (RKF45) method implementation
    """
    
    def __init__(self, rtol=1e-6, atol=1e-8, hmin=1e-10, hmax=1.0, safety=0.9):
        """
        Initialize adaptive RK solver
        
        Parameters:
        rtol: relative tolerance
        atol: absolute tolerance
        hmin: minimum step size
        hmax: maximum step size
        safety: safety factor for step size adjustment
        """
        self.rtol = rtol
        self.atol = atol
        self.hmin = hmin
        self.hmax = hmax
        self.safety = safety
        
    def rkf45_step(self, f, t, y, h):
        """
        Perform one RKF45 step
        
        Returns:
        y_new: new solution (4th order)
        y_high: higher order solution (5th order)
        error: estimated error
        """
        # RKF45 coefficients
        a2, a3, a4, a5, a6 = 1/4, 3/8, 12/13, 1, 1/2
        b21 = 1/4
        b31, b32 = 3/32, 9/32
        b41, b42, b43 = 1932/2197, -7200/2197, 7296/2197
        b51, b52, b53, b54 = 439/216, -8, 3680/513, -845/4104
        b61, b62, b63, b64, b65 = -8/27, 2, -3544/2565, 1859/4104, -11/40
        
        # 4t