四阶龙格-库塔方法求解微分方程

简介

提供一个使用四阶龙格-库塔方法 [1] 求解二阶微分方程的计算器，并提供一个基于解析方法的计算器作为对比和教育用途。
首先介绍如何将四阶龙格-库塔方法应用于一阶微分方程，然后将相同的方法应用于二阶线性微分方程。最后，提供一个互动教程，以研究在龙格-库塔方法中使用的步长 \( h \) 对近似解的影响。

应用于一阶微分方程的四阶龙格-库塔方法

四阶龙格-库塔方法是一种常用于求解形如以下形式的常微分方程（ODE）的数值方法： \[ \dfrac{dy}{dt} = f(y, t) \] 其中初始条件为 \( y(0) = y_0 \)。四阶龙格-库塔方法基于在未来时刻 \( t + h \) （其中 \( h \) 是步长）时的信息来近似解 \( y(t) \)。该方法使用四个中间步骤来计算在区间 \([t, t+h]\) 内不同点的斜率，并将这些斜率组合成加权平均斜率，随后用该斜率估计 \( y(t+h) \)。以下是四阶龙格-库塔方法的步骤：
步骤1： \[ k_1 = h \cdot f(y_n, t_n) \] 此处，\( y_n \) 和 \( t_n \) 是当前 \( y \) 和 \( t \) 的值。
步骤2： \[ k_2 = h \cdot f\left(y_n + \dfrac{k_1}{2}, t_n + \dfrac{h}{2}\right) \] 步骤3： \[ k_3 = h \cdot f\left(y_n + \dfrac{k_2}{2}, t_n + \dfrac{h}{2}\right) \] 步骤4： \[ k_4 = h \cdot f(y_n + k_3, t_n + h) \] 步骤5：
更新 \( y \) 和 \( t \) 四阶龙格-库塔方法结合上述量（斜率）来近似区间 \([t, t+h]\) 末端的 \( y \) 值。使用以下公式计算这些斜率的加权平均值： \[ \dfrac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \] 然后将该平均斜率乘以 \( h \)，并加到 \( y_n \) 以更新 \( t_{n+1} \) 处的 \( y \) 值。 \[ y_{n+1} = y_n + h (\dfrac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)) \] \[ t_{n+1} = t_n + h \]

应用于二阶微分方程的龙格-库塔方法

对于形如以下形式的二阶微分方程： \[ \dfrac{d^2y}{dt^2} = f(y, \dfrac{dy}{dt}, t) \] 初始条件为： \[ y(0) = y_0 \] \[ \dfrac{dy}{dt}(0) = \dfrac{dy}{dt}_0 \] 我们可以通过引入新变量将其转换为两个一阶微分方程的系统：
设 \[ \dfrac{dy}{dt} = v \] 因此，二阶微分方程可以写成如下系统： \[ \dfrac{dy}{dt} = v \] \[ \dfrac{dv}{dt} = \dfrac{d^2y}{dt^2} = f(y, v, t) \] 然后我们可以使用四阶龙格-库塔方法求解这个一阶微分方程组，如下所示：
步骤1：
计算 \( k_{1y} \) 和 \( k_{1v} \) \[ k_{1y} = h \cdot v \] \[ k_{1v} = h \cdot f(y_n, v_n, t_n) \]
步骤2：
计算 \( k_{2y} \) 和 \( k_{2v} \) \[ k_{2y} = h \cdot (v_n + \dfrac{k_{1v}}{2}) \] \[ k_{2v} = h \cdot f\left(y_n + \dfrac{k_{1y}}{2}, v_n + \dfrac{k_{1v}}{2}, t_n + \dfrac{h}{2}\right) \]
步骤3：
计算 \( k_{3y} \) 和 \( k_{3v} \) \[ k_{3y} = h \cdot (v_n + \dfrac{k_{2v}}{2}) \] \[ k_{3v} = h \cdot f\left(y_n + \dfrac{k_{2y}}{2}, v_n + \dfrac{k_{2v}}{2}, t_n + \dfrac{h}{2}\right) \]
步骤4：
计算 \( k_{4y} \) 和 \( k_{4v} \) \[ k_{4y} = h \cdot (v_n + k_{3v}) \] \[ k_{4v} = h \cdot f(y_n + k_{3y}, v_n + k_{3v}, t_n + h) \]
步骤5：
更新 \( y \) 和 \( v \) \[ y_{n+1} = y_n + \dfrac{1}{6}(k_{1y} + 2k_{2y} + 2k_{3y} + k_{4y}) \] \[ v_{n+1} = v_n + \dfrac{1}{6}(k_{1v} + 2k_{2v} + 2k_{3v} + k_{4v}) \] \[ t_{n+1} = t_n + h \]

龙格-库塔方法计算器

我们现在提供一个计算器，使用上述描述的四阶龙格-库塔方法来近似求解如下形式的二阶微分方程： \[ a \dfrac{d^2y}{dt^2} + b \dfrac{dy}{dt} + c y = 0 \] 同时也用解析方法求解该微分方程，并进行比较。
上面的微分方程也会通过解析方法求解，以便进行比较。
设 \[ \dfrac{dy}{dt} = v \] 因此我们需要将四阶龙格-库塔方法应用于以下系统。 \[ \dfrac{dy}{dt} = v \] \[ \dfrac{dv}{dt} = \dfrac{d^2y}{dt^2} = \dfrac{1}{a} (- b \; v - c \; y ) \] 该计算器使用四阶龙格-库塔方法，并同时解析求解相同的微分方程，并在图形上进行比较。

输入系数和初始条件

请输入系数 \(a, b \) 和 \( c \) 以及初始条件 \( y(0) \) 和 \( dy/dt(0) \) 以及步长 \( h \)，然后求解并比较所得到的数值解和解析解的图形。
注意
1) 您可以使用鼠标在图形上悬停，以获得图上点的读数，从而进行比较。
读取图形上的点
2) 使用鼠标可以一次显示一张图或两张图，通过点击右侧图例中的图形来切换。

3) 使用龙格-库塔方法的误差依赖于步长 \( h \)。本计算器允许更改 \( h \) 并比较数值方法和解析方法，并研究 \( h \) 对四阶龙格-库塔方法中使用的近似解的误差的影响。
4) 您还可以更改计算中使用的步数 \( n \)。图形的范围从 \( x = 0 \) 到 \( x = n \times h \)。

a:
b:
c:
y(0):
dy/dt(0):
步长: \( h \)
步数: \( n \)

解:

互动教程: 步长 \( h \) 对龙格-库塔方法近似解的影响

输入以下常数 \( b , c \) 和 \( d \) 的值：
\( a = 1 \) , \( b = 0.5 \) 和 \( c = 5 \)，这应当产生一个振幅递减的振荡函数。
现在输入从 \( h = 0.01 \) 开始的 \( h \) 值，然后逐步增加到1。哪些 \( h \) 值会在与解析方法进行比较时产生更好的近似解？