常微分方程

一阶微分方程

变量分离法

一阶微分方程有如下形式： [ \frac{dy}{dx} = g(x)h(y) ] 操作步骤:

  1. 分离变量：将方程改写为所有 (y) 项与 (dy) 在一起，所有 (x) 项与 (dx) 在一起： [ \frac{1}{h(y)} \, dy = g(x) \, dx ] 注意：若 (h(y) = 0) 有解 (y = c)，则 (y = c) 也是原方程的解。

  2. 两边积分： [ \int \frac{1}{h(y)} \, dy = \int g(x) \, dx + C ] 其中 (C) 为任意常数。

  3. 求解 (y)：若积分可求出原函数，则得到隐式或显式通解。

变量替换法

通过变量代换将方程化为可分离或线性方程。常见类型包括：

齐次方程 形如 $\frac{dy}{dx} = f\left(\frac{y}{x}\right)$的微分方程。
令 (v = \frac{y}{x})，则 (y = vx)，(\frac{dy}{dx} = v + x \frac{dv}{dx})。
代入得： [ v + x \frac{dv}{dx} = f(v) ] 分离变量： [ \frac{dv}{f(v) - v} = \frac{dx}{x} ] 积分即可。

线性组合 形如:(\frac{dy}{dx} = f(ax+by+c))
令 (u = ax+by+c)，则 (\frac{du}{dx} = a + b \frac{dy}{dx} =a+bf(u))，代入后可分离变量。

常数变异法

对于如下形式微分方程: [ \frac{dy}{dx} + P(x)y = Q(x) ] 操作步骤:

  1. 解齐次方程： [ \frac{dy}{dx} + P(x)y = 0 ] 通解为： [ y_h = C e^{-\int P(x) dx} ]

  2. 常数变易：将常数 (C) 替换为函数 (C(x))，设特解： [ y_p = C(x) e^{-\int P(x) dx} ] 则 \(\begin{equation*} \frac{dy_p}{dx} = C'(x) e^{-\int P(x) dx} - P(x) C(x) e^{-\int P(x) dx} \end{equation*}\)

  3. 代入原方程： \(\begin{align*} C'(x) e^{-\int P(x) dx} - P(x) C(x) e^{-\int P(x) dx} + P(x) C(x) e^{-\int P(x) dx} &= Q(x) \\ C'(x) e^{-\int P(x) dx} &= Q(x) \\ C'(x) &= Q(x) e^{\int P(x) dx} \\ C(x) &= \int Q(x) e^{\int P(x) dx} \, dx + C \end{align*}\)
  4. 通解： [ y = e^{-\int P(x) dx} \left( \int Q(x) e^{\int P(x) dx} \, dx + C \right) ]

积分因子法

适用于一阶线性方程： [ \frac{dy}{dx} + P(x)y = Q(x) ] 操作步骤:

  1. 寻找积分因子 (\mu(x))，使得乘以 (\mu(x)) 后方程左边成为某个函数的导数： \(\begin{align*} \mu(x) \frac{dy}{dx} + \mu(x) P(x) y &= \frac{d}{dx} \left[ \mu(x) y \right] \\ \mu(x) \frac{dy}{dx} + \mu(x) P(x) y &= \mu'(x) y + \mu(x) \frac{dy}{dx} \qquad //\text{右边导数展开} \\ \mu(x) P(x) &= \mu'(x) \\ \mu'(x) &= \mu(x) P(x) \end{align*}\)

  2. 解出 (\mu(x))： [ \frac{d\mu}{\mu} = P(x) \, dx \quad \Rightarrow \quad \ln|\mu| = \int P(x) dx ] 取： [ \mu(x) = e^{\int P(x) dx} ]

  3. 乘以积分因子： \(\begin{equation*} \frac{d}{dx} \left[ \mu(x) y \right] = \mu(x) Q(x) \end{equation*}\)

  4. 积分求解： [ \mu(x) y = \int \mu(x) Q(x) \, dx + C ] [ y = \frac{1}{\mu(x)} \left( \int \mu(x) Q(x) \, dx + C \right) ]

总结

  • 变量分离法：适用于方程右侧可分离为 (g(x)h(y)) 的情况。
  • 变量替换法：通过代换（如齐次替换、线性组合、伯努利替换等）转化为可解形式。
  • 常数变异法：用于线性非齐次方程，将齐次解中的常数变为函数求特解。
  • 积分因子法：针对一阶线性方程，构造积分因子化为全微分形式。

ODE数值方法

牛顿插值公式

牛顿插值公式是通过差商概念构造插值多项式的方法。下面我将详细推导这一公式，展示其从基本思想到最终形式的完整过程。

一、问题背景与基本思想

  • 1.1 插值问题 给定 ( n+1 ) 个互异节点 ( (x_0, y_0), (x_1, y_1), \dots, (x_n, y_n) )，寻找一个次数不超过 ( n ) 的多项式 ( P_n(x) ) 满足： [ P_n(x_i) = y_i, \quad i = 0, 1, \dots, n ]

  • 1.2 牛顿插值的基本思路 希望将多项式写成如下形式（牛顿形式）： [ P_n(x) = a_0 + a_1(x-x_0) + a_2(x-x_0)(x-x_1) + \cdots + a_n(x-x_0)(x-x_1)\cdots(x-x_{n-1}) ] 这种形式具有递推性：增加新节点时，只需添加一项并确定新的系数。

二、差商的定义与性质

  • 2.1 差商的定义 差商（Divided Difference）是递归定义的：

  1. 零阶差商（函数值）： \(f[x_i] = f(x_i) = y_i\)

  2. 一阶差商： \(f[x_i, x_j] = \frac{f[x_j] - f[x_i]}{x_j - x_i}\)

  3. 二阶差商： \(\begin{equation*} f[x_i, x_j, x_k] = \frac{f[x_j, x_k] - f[x_i, x_j]}{x_k - x_i} \end{equation*}\)

  4. k阶差商（递归定义）： \(f[x_0, x_1, \dots, x_k] = \frac{f[x_1, \dots, x_k] - f[x_0, \dots, x_{k-1}]}{x_k - x_0}\)

  • 2.2 差商的性质
    1. 对称性：差商值与节点的排列顺序无关
    2. 与导数的关系：若 ( f ) 在包含节点的区间上 ( k ) 次可导，则存在 ( \xi ) 使得： \(f[x_0, x_1, \dots, x_k] = \frac{f^{(k)}(\xi)}{k!}\)

三、牛顿插值公式的推导

  • 3.1 多项式形式设定 设多项式为： [ P_n(x) = a_0 + a_1(x-x_0) + a_2(x-x_0)(x-x_1) + \cdots + a_n(x-x_0)(x-x_1)\cdots(x-x_{n-1}) \tag{1} ]

  • 3.2 逐步确定系数

步骤1：利用 ( P_n(x_0) = y_0 ) 代入 ( x = x_0 )： \(P_n(x_0) = a_0 = y_0 \quad \Rightarrow \quad a_0 = f[x_0]\)

步骤2：利用 ( P_n(x_1) = y_1 ) 代入 ( x = x_1 )： \(P_n(x_1) = a_0 + a_1(x_1 - x_0) = y_1\) 解得： \(a_1 = \frac{y_1 - y_0}{x_1 - x_0} = f[x_0, x_1]\)

步骤3：利用 ( P_n(x_2) = y_2 ) 代入 ( x = x_2 )： [ P_n(x_2) = a_0 + a_1(x_2 - x_0) + a_2(x_2 - x_0)(x_2 - x_1) = y_2 ] 将 ( a_0, a_1 ) 代入整理： [ a_2 = \frac{y_2 - a_0 - a_1(x_2 - x_0)}{(x_2 - x_0)(x_2 - x_1)} ] 将 ( a_0 = f[x_0], a_1 = f[x_0, x_1] ) 代入，通过代数运算可得： \(a_2 = \frac{\frac{y_2 - y_1}{x_2 - x_1} - \frac{y_1 - y_0}{x_1 - x_0}}{x_2 - x_0} = f[x_0, x_1, x_2]\)

步骤4：归纳假设与一般情况 假设前 ( k ) 个系数满足： \(a_j = f[x_0, x_1, \dots, x_j], \quad j = 0, 1, \dots, k\) 现在利用 ( P_n(x_{k+1}) = y_{k+1} ) 来确定 ( a_{k+1} )。

由(1)式： [ P_n(x_{k+1}) = \sum_{j=0}^{k} a_j \prod_{i=0}^{j-1} (x_{k+1} - x_i) + a_{k+1} \prod_{i=0}^{k} (x_{k+1} - x_i) = y_{k+1} ] 其中约定空乘积为1。

解得： [ a_{k+1} = \frac{y_{k+1} - \sum_{j=0}^{k} a_j \prod_{i=0}^{j-1} (x_{k+1} - x_i)}{\prod_{i=0}^{k} (x_{k+1} - x_i)} \tag{2} ]

我们需要证明 ( a_{k+1} = f[x_0, x_1, \dots, x_{k+1}] )。

  • 3.3 证明 ( a_{k+1} = f[x_0, \dots, x_{k+1}] )

考虑两个插值多项式：

  1. ( Q(x) )：基于节点 ( x_0, \dots, x_k ) 的牛顿插值多项式（次数 ≤ k）
  2. ( R(x) )：基于节点 ( x_1, \dots, x_{k+1} ) 的牛顿插值多项式（次数 ≤ k）

由归纳假设： \(Q(x) = f[x_0] + \sum_{j=1}^{k} f[x_0, \dots, x_j] \prod_{i=0}^{j-1} (x - x_i)\) \(R(x) = f[x_1] + \sum_{j=2}^{k+1} f[x_1, \dots, x_j] \prod_{i=1}^{j-1} (x - x_i)\)

关键构造：定义多项式 [ S(x) = \frac{(x - x_0)R(x) - (x - x_{k+1})Q(x)}{x_{k+1} - x_0} \tag{3} ]

验证 ( S(x) ) 的性质：

  1. 次数：分子中最高次项为 ( x^{k+1} )，除以线性因子后，( S(x) ) 次数 ≤ k+1
  2. 插值性：对于 ( i = 0, \dots, k+1 )：
    • 当 ( x = x_0:S(x_0) = \frac{(x_0 - x_0)R(x_0) - (x_0 - x_{k+1})Q(x_0)}{x_{k+1} - x_0} = \frac{-(x_0 - x_{k+1})y_0}{x_{k+1} - x_0} = y_0 )
    • 当 ( x = x_{k+1} )：( S(x_{k+1}) = \frac{(x_{k+1} - x_0)R(x_{k+1}) - (x_{k+1} - x_{k+1})Q(x_{k+1})}{x_{k+1} - x_0} = R(x_{k+1}) = y_{k+1} )
    • 当 ( x = x_j )（( 1 ≤ j ≤ k )）：直接验证得 ( S(x_j) = y_j )

因此 ( S(x) ) 是满足所有 ( k+2 ) 个插值条件的次数 ≤ k+1 的多项式，由唯一性，( S(x) = P_{k+1}(x) )。

比较系数： 展开(3)式，考虑 ( x^{k+1} ) 项的系数：

  • ( R(x) ) 中 ( x^k ) 项系数为 ( f[x_1, \dots, x_{k+1}] )，所以 ( (x-x_0)R(x) ) 中 ( x^{k+1} ) 项系数为 ( f[x_1, \dots, x_{k+1}] )
  • ( Q(x) ) 中 ( x^k ) 项系数为 ( f[x_0, \dots, x_k] )，所以 ( (x-x_{k+1})Q(x) ) 中 ( x^{k+1} ) 项系数为 ( f[x_0, \dots, x_k] )

因此 ( S(x) ) 中 ( x^{k+1} ) 项系数为： \(\frac{f[x_1, \dots, x_{k+1}] - f[x_0, \dots, x_k]}{x_{k+1} - x_0} = f[x_0, x_1, \dots, x_{k+1}]\)

而 ( P_{k+1}(x) ) 中 ( x^{k+1} ) 项系数正是 ( a_{k+1} )，所以： [ a_{k+1} = f[x_0, x_1, \dots, x_{k+1}] ]

这就完成了归纳证明。

四、牛顿插值公式的最终形式

综合以上结果，我们得到牛顿插值公式： \(P_n(x) = f[x_0] + \sum_{k=1}^{n} f[x_0, x_1, \dots, x_k] \cdot \prod_{i=0}^{k-1} (x - x_i)\)

其中：

  • ( f[x_0] = y_0 )
  • ( f[x_0, x_1, \dots, x_k] ) 是 k 阶差商
  • $\prod_{i=0}^{k-1} (x - x_i) = (x-x_0)(x-x_1)\cdots(x-x_{k-1})$

五、余项公式

  • 5.1 差商形式的余项 若增加一个节点 ( x )，则插值多项式变为： [ P_{n+1}(t) = P_n(t) + f[x_0, \dots, x_n, x] \prod_{i=0}^{n} (t - x_i) ] 由于 ( P_{n+1}(x) = f(x) )，令 ( t = x ) 得： [ f(x) = P_n(x) + f[x_0, \dots, x_n, x] \prod_{i=0}^{n} (x - x_i) ] 所以余项为： [ R_n(x) = f(x) - P_n(x) = f[x_0, x_1, \dots, x_n, x] \prod_{i=0}^{n} (x - x_i) \tag{4} ]

  • 5.2 导数形式的余项 若 ( f ) 在包含节点 ( x_0, \dots, x_n, x ) 的区间上 ( n+1 ) 次可导，由差商与导数的关系，存在 ( \xi ) 使得： [ f[x_0, \dots, x_n, x] = \frac{f^{(n+1)}(\xi)}{(n+1)!} ] 代入(4)式得： [ R_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^{n} (x - x_i) \tag{5} ] 这就是牛顿插值余项的导数形式。

六、算法实现与差商表

  • 6.1 差商表的计算 差商可以通过表格递归计算：
( x_i ) ( f[x_i] ) 一阶差商 二阶差商 三阶差商 …
( x_0 ) ( f[x_0] )        
( x_1 ) ( f[x_1] ) ( f[x_0,x_1] )      
( x_2 ) ( f[x_2] ) ( f[x_1,x_2] ) ( f[x_0,x_1,x_2] )    
( x_3 ) ( f[x_3] ) ( f[x_2,x_3] ) ( f[x_1,x_2,x_3] ) ( f[x_0,x_1,x_2,x_3] )  
… … … … … …

计算规则：( f[x_i, \dots, x_{i+k}] = \frac{f[x_{i+1}, \dots, x_{i+k}] - f[x_i, \dots, x_{i+k-1}]}{x_{i+k} - x_i} )

七、与拉格朗日插值的比较

特性 拉格朗日插值 牛顿插值
构造方式 基于基函数线性组合 基于差商的递推构造
计算复杂度 每次计算 O(n²) 预计算 O(n²)，求值 O(n)
添加新节点 需重新计算所有基函数 只需计算新差商，添加一项
数值稳定性 高次插值可能不稳定 相对更稳定（差商有导数意义）
形式特点 对称美观，理论分析方便 递推形式，实用性强

八、总结

牛顿插值公式的推导过程体现了以下关键点：

  1. 递推思想：多项式采用嵌套乘积形式，便于逐步添加新节点
  2. 差商的核心作用：系数恰好是各阶差商，建立了离散差商与多项式系数的联系
  3. 构造性证明：通过巧妙构造辅助多项式，证明了系数与差商的等价性
  4. 余项的自然表达：差商形式余项简洁且便于计算

牛顿插值法不仅提供了高效的计算方法，还揭示了插值问题的深刻数学结构：差商作为离散导数的推广，连接了离散插值与连续微分。

Euler法

Heun法

Runge-Kutta法

Adams-Bashforth法

Adams-Bashforth公式的推导基于拉格朗日插值多项式和积分近似思想，其核心是通过前k步的函数值构造插值多项式，近似被积函数，进而推导出递推公式。以下以二阶（AB2）和四阶（AB4）方法为例，详细拆解推导过程：

1. 基础思想：从常微分方程到积分方程

对于初值问题： [ y’(t) = f(t, y(t)), \quad y(t_0) = y_0 ] 可转化为积分形式： [ y(t_{n+1}) = y(t_n) + \int_{t_n}^{t_{n+1}} f(t, y(t)) \, dt ] Adams-Bashforth方法通过数值积分近似右侧积分项，而积分核(f(t, y(t)))由前k步的函数值通过拉格朗日插值多项式近似。

2. 拉格朗日插值多项式的构造

假设已知前(k)步的函数值(f(t_{n-i}, y_{n-i}))（(i=0,1,\dots,k-1)），在区间([t_n, t_{n+1}])上构造(k-1)次拉格朗日插值多项式(P_{k-1}(t))，逼近(f(t, y(t)))： [ P_{k-1}(t) = \sum_{j=0}^{k-1} f(t_{n-j}, y_{n-j}) \cdot l_j(t) ] 其中基函数(l_j(t))满足： [ l_j(t) = \prod_{\substack{m=0 \ m \neq j}}^{k-1} \frac{t - t_{n-m}}{t_{n-j} - t_{n-m}} ] 基函数在对应节点(t_{n-j})处取值为1，其他节点为0，确保插值多项式精确通过已知点。

3. 积分近似与递推公式推导

将积分近似为：
\(\begin{align*} \int_{t_n}^{t_{n+1}} f(t, y(t)) \, dt &\approx \int_{t_n}^{t_{n+1}} P_{k-1}(t) \, dt \\ \Leftrightarrow y(t_{n+1}) & \approx y(t_n) + \int_{t_n}^{t_{n+1}} P_{k-1}(t) \, dt \\ & = y(t_n) + \int_{t_n}^{t_{n+1}}\sum_{j=0}^{k-1} \left[ f(t_{n-j}, y_{n-j}) \cdot l_j(t)\right] \, dt \quad //带入P_{k-1}(t)\\ & = y(t_n) + \sum_{j=0}^{k-1} \left[ f(t_{n-j}, y_{n-j}) \cdot \int_{t_n}^{t_{n+1}} l_j(t) \, dt \right] \quad //求和与积分交换顺序\\ &= y_n + \sum_{j=0}^{k-1} b_j f(t_{n-j}, y_{n-j}) \quad//令 b_j = \int_{t_n}^{t_{n+1}} l_j(t) \, dt \end{align*}\)

积分系数(b_j = \int_{t_n}^{t_{n+1}} l_j(t) \, dt)，则公式为： \(\begin{align*} b_j &= \int_{t_n}^{t_{n+1}} l_j(t) \, dt \\ &= \int_{t_n}^{t_{n+1}} \prod_{\substack{m=0 \\ m \neq j}}^{k-1} \frac{t - t_{n-m}}{t_{n-j} - t_{n-m}} \, dt \\ &= \int_{t_n}^{t_{n+1}} \prod_{\substack{m=0 \\ m \neq j}}^{k-1} \frac{nh+\tau h - (n-m)h}{(n-j)h - (n-m)h} \, dt \quad //令t_i=i h, t=nh+\tau h ,h为均匀长\\ &= h\int_0^1 \prod_{\substack{m=0 \\ m \neq j}}^{k-1} \frac{m+\tau}{m-j} \, d\tau \quad //dt=hd\tau, \tau \in [0,1]\\ \end{align*}\)

4. 具体阶数公式推导

二阶方法（AB2，k=2）

  • 积分系数计算： [ b_0 = h\int_0^1 \prod_{\substack{m=0 \ m \neq j}}^{k-1} \frac{m+\tau}{m-j} \, d\tau = h\int_0^1 \frac{1+\tau}{1-0} \, d\tau = \frac {3h} {2} \quad //k=2,j=0 ] [ b_1 = h\int_0^1 \prod_{\substack{m=0 \ m \neq j}}^{k-1} \frac{m+\tau}{m-j} \, d\tau = h\int_0^1 \frac{0+\tau}{0-1} \, d\tau = -\frac {h} {2} \quad //k=2,j=1 ]
  • 递推公式： \(\begin{equation*} y_{n+1} = y_n + \frac{h}{2} \left[ 3f(t_n, y_n) - f(t_{n-1}, y_{n-1}) \right] \end{equation*}\)

四阶方法（AB4，k=4）

  • 积分系数计算： \(\begin{align*} b_0 &= h\int_0^1 \prod_{\substack{m=0 \\ m \neq j}}^{k-1} \frac{m+\tau}{m-j} \, d\tau \\ &= h\int_0^1 \frac{1+\tau}{1-0} \cdot \frac{2+\tau}{2-0} \cdot \frac{3+\tau}{3-0} \, d\tau \quad //k=4,j=0 \\ &= \frac {55} {24}h \end{align*}\)

    类似的可以算出$b_1=-\frac {59} {24}h,b_2=\frac {37} {24}h,b_3=-\frac {9} {24}h$

  • 递推公式： \(\begin{equation*} y_{n+1} = y_n + \frac{h}{24} \left[ 55f(t_n, y_n) - 59f(t_{n-1}, y_{n-1}) + 37f(t_{n-2}, y_{n-2}) - 9f(t_{n-3}, y_{n-3}) \right] \end{equation*}\)

5. 关键特性

  • 系数确定：积分系数(b_j)由基函数在区间([t_n, t_{n+1}])上的积分确定，需确保多项式在节点处精确匹配，并通过代数运算简化系数。
  • 阶数与步数关系：s步Adams-Bashforth方法的阶数为s（局部截断误差为(O(h^{s+1}))），阶数与步数相等，这是通过插值多项式的次数（s-1次）和积分精度共同决定的。
  • 显式特性：由于插值多项式仅依赖前k步的函数值，递推公式无需解方程，直接由已知值计算下一步解，属于显式方法。

6. 初始值问题与启动方法

Adams-Bashforth方法非自启动，需前k步的初始值。通常使用同阶或更高阶的Runge-Kutta方法（如四阶Runge-Kutta）提供初始值。例如，四阶AB4需前4步的初始值，由RK4方法计算得到。

7. 误差分析与稳定性

  • 局部截断误差：s阶方法的局部截断误差为(O(h^{s+1}))，随阶数升高而降低。
  • 稳定性区域：显式方法稳定性区域较小，步长过大易发散，常与隐式Adams-Moulton方法组合为预测-校正系统以增强稳定性。

总结

Adams-Bashforth公式的推导本质是拉格朗日插值多项式在数值积分中的应用，通过构造前k步函数值的插值多项式，近似被积函数并积分得到递推关系。其阶数与步数相等，高阶方法精度高但稳定性弱，需结合问题特性选择阶数，并搭配初始值方法和稳定性增强策略。这一方法在工程和科学计算中广泛应用于高精度、光滑函数的常微分方程求解。

Adams–Moulton法