### Python 实现龙格库塔公式算法 #### 一、龙格库塔算法基本原理 对于常微分方程初值问题，其一般形式可以表示为： \[ \frac{dy}{dx} = f(x, y),\quad y(x_0)=y_0 \] 其中 \(f\) 是已知函数。为了近似求解这类问题，龙格库塔方法采用了一种更为精确的方式估计下一个时间步的值。 四阶龙格库塔方法是一种广泛应用的方法，它通过四个不同点上的斜率加权平均来预测下一时刻的状态[^1]。 具体来说，该方法定义如下几个中间变量: - \(k_1 = h * f(t_n , y_n)\) - \(k_2 = h * f(t_n + \frac{h}{2}, y_n+\frac{k_1}{2})\) - \(k_3 = h * f(t_n + \frac{h}{2}, y_n+\frac{k_2}{2})\) - \(k_4 = h * f(t_n+h,y_n+k_3)\) 最终更新公式为: \[ y_{n+1}=y_n+(\frac{k_1}{6}+\frac{k_2}{3}+\frac{k_3}{3}+\frac{k_4}{6}) \] 这里的 \(t_n\) 表示当前的时间节点，\(y_n\) 则代表对应的因变量值；而 \(h\) 就是指定的时间步长[^2]。 #### 二、Python实现龙格库塔算法 下面给出一段简单的Python代码用于演示如何利用上述提到的四阶龙格库塔法解决具体的ODE问题。 ```python def runge_kutta(f, tspan, y0, nsteps): """ 使用四阶Runge-Kutta方法求解给定的一阶常微分方程 参数: f -- 微分方程右侧表达式 dy/dx=f(x,y)，作为函数对象传入 tspan -- 时间区间 (start_time,end_time) y0 -- 初始条件 nsteps -- 步数 返回: ts -- 计算过程中所使用的离散化时间序列数组 ys -- 对应于ts中各时间节点处得到的结果向量列表 """ start_t, end_t = tspan h = (end_t - start_t) / float(nsteps) ts = [] ys = [] t = start_t y = y0 while abs(t - end_t) >= 1e-8: k1 = h * f(t , y ) k2 = h * f(t + h/2., y + k1/2.) k3 = h * f(t + h/2., y + k2/2.) k4 = h * f(t + h , y + k3 ) y += (k1 + 2*k2 + 2*k3 + k4)/6. t += h ts.append(t) ys.append(y.copy()) return ts, ys if __name__ == "__main__": import numpy as np from matplotlib import pyplot as plt def func(t, y): """ 定义要解决的具体ODE """ return (-5*y + 7*np.exp(-t)) # 设置参数并调用runge_kutta() 函数计算数值解 time_interval = [0, 5] initial_value = [np.array([1])] num_of_steps = 1000 times, values = runge_kutta(func, time_interval, initial_value[0], num_of_steps) exact_solution = lambda t : 2.*np.exp(-t)-np.exp(-5*t)+1. fig, ax = plt.subplots() ax.plot(times,[exact_solution(i)[0]for i in times],'r-',label='Exact') ax.plot(times,np.asarray(values).reshape(len(times)), 'b--', label="RK4 Approximation") legend = ax.legend(loc='upper center') plt.show() ``` 此段代码实现了基于四阶龙格库塔方法对特定类型的线性常微分方程进行数值积分的过程，并绘制出了理论解与数值逼近之间的对比图象[^3]。