在Mathematica中，使用鸟枪法求解在无穷远处的边值常微分方程
边界值问题最简单的例子是定义在区间a ≤ x ≤ b上的二阶常微分方程y′′ = f (x, y, y′)，其边界条件为y(a) = α和y(b) = β，其中α和β是给定的常数。需要注意的是，由于自变量通常表示空间坐标，因此常用x表示（此处符号'≡d/dx）。
我们将探讨几种计算上述二阶边值问题解的方法。第一种方法利用时间步进算法（注意此处"时间"指代x变量），具体步骤如下：
步骤分解
- 方程转换将y′′ = f (x, y, y′)改写为一阶常微分方程组：u′ = f (x, u)，其中u = (y, y′)ᵀ（符号ᵀ表示列向量转置）
- 初值问题求解
- 选取初始斜率S作为y′(a)的猜测值
- 使用时间步进算法求解初值问题：u′ = f (x, u)，初始条件u(a) = (α, S)ᵀ从x = a积分至x = b
- 检查x = b处的解值：
- 若计算值大于目标值β → 初始斜率S导致过冲
- 若计算值小于目标值β → 初始斜率S导致欠冲
- 斜率范围确定
- 找到另一个初始斜率值，使其在x = b处产生相反的过冲/欠冲现象
- 此时获得两个初始斜率值S₁和S₂，分别对应过冲和欠冲
- 根查找算法应用设̃y(x; S)表示使用初始斜率S计算的近似解，本质上是寻找使表达式̃y(b; S) - β为零的根S*。由于：
- ̃y(b; S) - β的显式公式通常未知
- 但已知根S*位于S₁和S₂之间因此适合使用以下方法：
- 二分法（利用区间套原理）
- 割线法（基于线性插值）来迭代求解满足边界条件的初始斜率S*
(*Define a finite version of "infinity"*)
inf = 5;(*Define the differential equation and its initial conditions, parameterized by the initial gradient y'[0] == dy0. For simplicity, set y[0] == 1*)
deqn = {y''[x] - x y[x] == 0, y[0] == 1, y'[0] == dy0};(*Compute the numerical solution parameterised by dy0*)
ydysol = ParametricNDSolve[deqn, y, {x, 0, inf}, dy0][[1]](*Find the value of the initial gradient dy0 that makes the solution go to zero at "infinity". choose to minimise y[inf]^2 w.r.t. dy0*)
dysol = FindMinimum[((y[dy0] /. ydysol)[inf])^2 // Evaluate, {dy0, -1}]
(*{7.51024*10^-27, {dy0 -> -0.729012}}*)Plot[(y[dy0] /. ydysol /. dysol[[2]])[x] // Evaluate, {x, 0, inf}]