Python详细实现龙格-库塔算法

闲人编程

已于 2025-03-29 21:40:47 修改

阅读量1.9k

点赞数 35

分类专栏： python 文章标签： python 算法开发语言 RC4 RC5 龙格-库塔法常微分

于 2024-11-06 07:15:00 首次发布

本文链接：https://blog.csdn.net/qq_42568323/article/details/143527399

版权

python 专栏收录该内容

123 篇文章

订阅专栏

龙格-库塔算法详解及案例分析

1. 引言

龙格-库塔(Runge-Kutta)方法是求解常微分方程(ODE)初值问题的一类重要数值方法，由德国数学家Carl Runge和Wilhelm Kutta在20世纪初提出。这类方法通过构造多个中间点的斜率估计，加权平均得到更高精度的数值解，在科学计算和工程应用中具有广泛用途。

本文将系统介绍龙格-库塔算法的数学原理、实现细节，并通过三个典型案例展示其实际应用。每个案例包含完整的Python实现、mermaid流程图以及结果可视化。

2. 算法原理

2.1 常微分方程初值问题

考虑一阶常微分方程初值问题：
$\frac{dy}{dt} = f(t,y), \quad y(t_0) = y_0$

2.2 经典四阶龙格-库塔方法(RK4)

最常用的四阶龙格-库塔方法计算步骤如下：

计算斜率k1：
$k_1 = f(t_n, y_n)$
计算斜率k2：
$k_2 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_1)$
计算斜率k3：
$k_3 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_2)$
计算斜率k4：
$k_4 = f(t_n + h, y_n + hk_3)$
更新解：
$y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$

其中 $h$ 为步长，该方法具有四阶精度，局部截断误差为 $O(h^5)$ 。

3. 算法实现

3.1 基础实现类

class RungeKutta4:
    def __init__(self, f, y0, t0, t_end, h):
        """
        初始化RK4求解器
        
        参数:
        f: 微分方程右端函数，形式为f(t, y)
        y0: 初始条件
        t0: 初始时间
        t_end: 终止时间
        h: 步长
        """
        self.f = f
        self.y = y0
        self.t = t0
        self.t_end = t_end
        self.h = h
        self.t_values = [t0]
        self.y_values = [y0]
        
    def step(self):
        """执行单步RK4计算"""
        k1 = self.f(self.t, self.y)
        k2 = self.f(self.t + self.h/2, self.y + self.h/2 * k1)
        k3 = self.f(self.t + self.h/2, self.y + self.h/2 * k2)
        k4 = self.f(self.t + self.h, self.y + self.h * k3)
        
        self.y += (self.h / 6) * (k1 + 2*k2 + 2*k3 + k4)
        self.t += self.h
        
        self.y_values.append(self.y)
        self.t_values.append(self.t)
        
    def solve(self):
        """求解整个时间区间的ODE"""
        while self.t < self.t_end:
            self.step()
        return np.array(self.t_values), np.array(self.y_values)

3.2 流程图

4. 案例分析

4.1 案例1：简单指数衰减

4.1.1 问题描述

求解指数衰减方程：
$\frac{dy}{dt} = -ky, \quad y(0) = y_0$
解析解为 $y(t) = y_0 e^{-kt}$ 。

4.1.2 代码实现

import numpy as np
import matplotlib.pyplot as plt

# 定义微分方程
def exponential_decay(t, y, k=0.3):
    return -k * y

# 参数设置
y0 = 5.0
t0 = 0.0
t_end = 20.0
h = 0.5

# 创建求解器实例
rk4 = RungeKutta4(lambda t, y: exponential_decay(t, y), y0, t0, t_end, h)

# 数值求解
t_vals, y_vals = rk4.solve()

# 解析解
t_exact = np.linspace(t0, t_end, 100)
y_exact = y0 * np.exp(-0.3 * t_exact)

# 绘图
plt.figure(figsize=(10,6))
plt.plot(t_vals, y_vals, 'o-', label='RK4 Numerical')
plt.plot(t_exact, y_exact, 'r-', label='Exact Solution')
plt.xlabel('Time')
plt.ylabel('y(t)')
plt.title('Exponential Decay: Numerical vs Exact Solution')
plt.legend()
plt.grid(True)
plt.show()

4.1.3 结果分析

该案例展示了RK4方法对简单ODE的高精度求解能力，数值解与解析解几乎完全重合，验证了算法的正确性。

4.2 案例2：Lotka-Volterra捕食者-猎物模型

4.2.1 问题描述

Lotka-Volterra模型描述捕食者与猎物种群动态：
$\begin{cases} \frac{dx}{dt} = \alpha x - \beta xy \\ \frac{dy}{dt} = \delta xy - \gamma y \end{cases}$
其中 $x$ 为猎物数量， $y$ 为捕食者数量。

4.2.2 代码实现

def lotka_volterra(t, z, alpha=1.1, beta=0.4, delta=0.1, gamma=0.4):
    x, y = z
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return np.array([dxdt, dydt])

# 初始条件：10个猎物，5个捕食者
z0 = np.array([10, 5])
t0 = 0
t_end = 50
h = 0.05

# 向量化RK4实现
class VectorRK4(RungeKutta4):
    def step(self):
        k1 = self.f(self.t, self.y)
        k2 = self.f(self.t + self.h/2, self.y + self.h/2 * k1)
        k3 = self.f(self.t + self.h/2, self.y + self.h/2 * k2)
        k4 = self.f(self.t + self.h, self.y + self.h * k3)
        
        self.y += (self.h / 6) * (k1 + 2*k2 + 2*k3 + k4)
        self.t += self.h
        
        self.y_values.append(self.y.copy())
        self.t_values.append(self.t)

# 求解
rk4_lv = VectorRK4(lotka_volterra, z0, t0, t_end, h)
t_vals, z_vals = rk4_lv.solve()
x_vals = z_vals[:,0]
y_vals = z_vals[:,1]

# 绘图
plt.figure(figsize=(12,5))

plt.subplot(1,2,1)
plt.plot(t_vals, x_vals, label='Prey (x)')
plt.plot(t_vals, y_vals, label='Predator (y)')
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('Population Dynamics Over Time')
plt.legend()
plt.grid(True)

plt.subplot(1,2,2)
plt.plot(x_vals, y_vals)
plt.xlabel('Prey Population')
plt.ylabel('Predator Population')
plt.title('Phase Portrait')
plt.grid(True)

plt.tight_layout()
plt.show()

4.2.3 结果分析

结果展示了捕食者与猎物种群数量的周期性振荡，以及典型的相图极限环，符合生态学预期行为。

4.3 案例3：刚性问题(Stiff Equation)

4.3.1 问题描述

求解刚性方程：
$\frac{dy}{dt} = -15y, \quad y(0) = 1$
该问题对步长选择敏感，适合测试算法的稳定性。

4.3.2 代码实现

def stiff_equation(t, y):
    return -15 * y

# 不同步长比较
h_values = [0.1, 0.2, 0.25]
t_end = 1.0

plt.figure(figsize=(10,6))

for h in h_values:
    rk4 = RungeKutta4(stiff_equation, 1.0, 0.0, t_end, h)
    t_vals, y_vals = rk4.solve()
    plt.plot(t_vals, y_vals, 'o-', label=f'h={h}')

# 解析解
t_exact = np.linspace(0, t_end, 100)
y_exact = np.exp(-15 * t_exact)
plt.plot(t_exact, y_exact, 'k-', label='Exact')

plt.yscale('log')
plt.xlabel('Time')
plt.ylabel('y(t) (log scale)')
plt.title('Stiff Equation Solution with Different Step Sizes')
plt.legend()
plt.grid(True)
plt.show()

4.3.3 结果分析

该案例展示了RK4方法在解决刚性问题时的表现。当步长过大时(h=0.25)，数值解出现不稳定振荡，说明对刚性方程需要谨慎选择步长或使用专门方法。

5. 算法变体与扩展

5.1 自适应步长RK方法

通过比较不同阶数方法的解来估计误差，动态调整步长：

class AdaptiveRK:
    def __init__(self, f, y0, t0, t_end, h0=0.1, tol=1e-6):
        self.f = f
        self.y = y0
        self.t = t0
        self.t_end = t_end
        self.h = h0
        self.tol = tol
        self.t_values = [t0]
        self.y_values = [y0]
        
    def rk4_step(self, t, y, h):
        k1 = self.f(t, y)
        k2 = self.f(t + h/2, y + h/2 * k1)
        k3 = self.f(t + h/2, y + h/2 * k2)
        k4 = self.f(t + h, y + h * k3)
        return y + (h / 6) * (k1 + 2*k2 + 2*k3 + k4)
    
    def rk5_step(self, t, y, h):
        # 嵌入的RK5方法
        k1 = self.f(t, y)
        k2 = self.f(t + h/4, y + h/4 * k1)
        k3 = self.f(t + 3*h/8, y + 3*h/8 * k2)
        k4 = self.f(t + 12*h/13, y + 12*h/13 * k3)
        k5 = self.f(t + h, y + h * k4)
        k6 = self.f(t + h/2, y + h/2 * k5)
        return y + h * (16*k1/135 + 6656*k2/12825 + 28561*k3/56430 - 9*k4/50 + 2*k5/55)
    
    def step(self):
        while True:
            # 计算RK4和RK5
            y4 = self.rk4_step(self.t, self.y, self.h)
            y5 = self.rk5_step(self.t, self.y, self.h)
            
            # 估计误差
            error = np.linalg.norm(y5 - y4)
            
            # 调整步长
            if error < self.tol:
                break
            self.h *= 0.9 * (self.tol / error)**0.2
            
        # 接受步长
        self.y = y4
        self.t += self.h
        self.t_values.append(self.t)
        self.y_values.append(self.y)
        
        # 调整下一步步长
        self.h *= 0.9 * (self.tol / error)**0.25
        
    def solve(self):
        while self.t < self.t_end:
            if self.t + self.h > self.t_end:
                self.h = self.t_end - self.t
            self.step()
        return np.array(self.t_values), np.array(self.y_values)