射击法

射击法（shooting method）将边值问题改写成带未知初始参数的初值问题：先猜初始斜率或其他缺失状态，向前积分，再根据末端偏差调整猜测。边界残差是关于初始参数的非线性方程，初值积分器和残差求根器共同决定结果。原边值问题有良好解，并不保证单次向前射击在数值上容易。

用未知斜率命中右端

考虑 $$y^{″}+{\omega }^{2}y=0,y(0)=0,y(1)=1,\omega =\pi /2.$$ 缺失的初始斜率记为 $s=y^{\prime }(0)$。初值问题的解为 $y(x;s)=s\sin (\omega x)/\omega$，末端残差为 $$R(s)=y(1;s)-1=\frac{2}{\pi }s-1.$$ 初猜 $s_0=0$ 未达到目标，$R(s_0)=-1$；取 $s_1=1$，末端为 $2/\pi$，仍偏低。两点割线的零点为 $s_2=\pi /2$，给出精确解 $y(x)=\sin (\pi x/2)$。一般非线性方程的残差不再是直线，通常需要多次调整。

ω=π/2。初始位置相同，末端高度由初始斜率控制；射击调整的是缺失初值，不是沿区间任意修改轨迹。

一般二阶问题与残差导数

对 $y^{″}=f(x,y,y^{\prime })$、$y(a)=A,y(b)=B$，令 $v=y^{\prime }$，从 $(y(a),v(a))=(A,s)$ 积分一阶系统 $y^{\prime }=v,v^{\prime }=f$。残差为 $R(s)=y(b;s)-B$。若相关初值解存在到 $b$，且右端对状态足够光滑，解关于参数也可微。

令 $z=\partial y/\partial s,w=\partial v/\partial s$，对方程求参数导数得到变分系统 $$z^{\prime }=w,w^{\prime }=f_{y}z+f_{v}w,z(a)=0,w(a)=1.$$ 因此 $R^{\prime }(s)=z(b)$，Newton 更新为 $s_{new}=s-R(s)/z(b)$。这些导数要沿当前射击轨迹计算；把它们当作固定常数，仅对线性常系数等特殊问题成立。没有导数时可用割线更新 $$s_{new}=s_1-R(s_1)\frac{s_1-s_0}{R(s_1)-R(s_0)}.$$ Atkinson、Han 与 Stewart 的常微分方程教材将边值问题的射击与初值积分误差一同讨论。

端点映射的可解性与敏感性

在主例中，$R^{\prime }(s)=\sin \omega /\omega =2/\pi \ne 0$，所以残差的零点唯一，且端点目标的小改动可由有限的斜率改动补偿。一般若 $R(s_{\ast })=0,R^{\prime }(s_{\ast })\ne 0$，隐函数定理保证附近端点条件改变时，斜率可局部唯一地随之变化；Newton 的局部快速收敛还需足够光滑和初猜接近。

有解变式将频率改为 $\omega =\pi$。此时 $y(1;s)=0$ 对所有 $s$ 都成立：右端要求一时无解，要求零时却有无穷多解。两个不同斜率给出相同残差，割线分母为零并非偶然的程序故障，而是端点映射失去可逆性。靠近这个频率时，$|\sin \omega /\omega |$ 很小，残差中的小误差会对应很大的斜率误差。

另一类困难来自初值方向的增长。对 $y^{″}={\mu }^{2}y$、零左端和右端一，精确斜率为 $s=\mu /\sinh \mu$，而端点敏感度为 $\sinh \mu /\mu$。$\mu$ 很大时，需要用极小斜率控制极大增长；一般状态中的舍入和积分误差会被放大。多重射击将区间拆成较短的初值问题，再通过匹配条件联立求解，可减轻这种长距离增长；有限差分或有限元直接处理全局边界也可能更合适。

可运行的 RK4 与割线射击

代码固定区间 $[0,1]$，接受右端函数 $f(x,y,v)$、两端值及两个斜率猜测。每次用经典 RK4 积分，返回斜率和轨迹，不修改输入。$K$ 次射击、每次 $N$ 步的成本为 $O(KN)$ 次右端求值，保存一条轨迹为 $O(N)$ 空间。程序控制的是数值端点残差，正式使用还需加密积分步长，检验残差求根没有被积分误差限制。

Python 3

import math


def shooting_secant(force, left, right, s0, s1, steps=100, tol=1e-10, maxiter=30):
    left, right, s0, s1, tol = map(float, (left, right, s0, s1, tol))
    if not all(math.isfinite(v) for v in (left, right, s0, s1, tol)) or tol <= 0:
        raise ValueError("finite data and positive tolerance required")
    if not isinstance(steps, int) or steps <= 0 or not isinstance(maxiter, int) or maxiter <= 0:
        raise ValueError("positive integer step and iteration counts required")
    h = 1/steps
    def integrate(slope):
        y, v = left, slope; track = [(0.0, y, v)]
        for i in range(steps):
            x = i*h
            k1 = (v, float(force(x, y, v)))
            k2 = (v+h*k1[1]/2, float(force(x+h/2, y+h*k1[0]/2, v+h*k1[1]/2)))
            k3 = (v+h*k2[1]/2, float(force(x+h/2, y+h*k2[0]/2, v+h*k2[1]/2)))
            k4 = (v+h*k3[1], float(force(x+h, y+h*k3[0], v+h*k3[1])))
            y += h*(k1[0]+2*k2[0]+2*k3[0]+k4[0])/6
            v += h*(k1[1]+2*k2[1]+2*k3[1]+k4[1])/6
            if not math.isfinite(y) or not math.isfinite(v):
                raise ArithmeticError("shooting trajectory became nonfinite")
            track.append(((i+1)*h, y, v))
        return y-right, track
    r0, track0 = integrate(s0)
    if abs(r0) <= tol:
        return s0, track0
    for _ in range(maxiter):
        r1, track1 = integrate(s1)
        if abs(r1) <= tol:
            return s1, track1
        if r1 == r0:
            raise ArithmeticError("zero secant denominator")
        s2 = s1-r1*(s1-s0)/(r1-r0)
        if not math.isfinite(s2):
            raise ArithmeticError("nonfinite shooting parameter")
        s0, r0, s1 = s1, r1, s2
    raise RuntimeError("endpoint residual did not reach tolerance")


w = math.pi/2
slope, trajectory = shooting_secant(lambda x, y, v: -w*w*y, 0, 1, 0, 1)
print(slope, math.pi/2)  # close; remaining difference is RK4 integration error
print(trajectory[-1][1])  # near 1

对非线性问题，猜测轨迹可能在到达末端前爆破，或者残差有多个零点。应根据方程选择可行初猜、步长与求根策略，并报告找到了哪一个分支。一个很小的端点残差也不单独保证整个区间上的轨迹误差小。