May 11, 2026

miniyuan

插值法

---

题 5.1

基函数：

$$\begin{aligned} l_0(x) &= \frac{(x-3)(x-4)(x-6)}{(1-3)(1-4)(1-6)} = -\frac{(x-3)(x-4)(x-6)}{30} \\ l_1(x) &= \frac{(x-1)(x-4)(x-6)}{(3-1)(3-4)(3-6)} = \frac{(x-1)(x-4)(x-6)}{6} \\ l_2(x) &= \frac{(x-1)(x-3)(x-6)}{(4-1)(4-3)(4-6)} = -\frac{(x-1)(x-3)(x-6)}{6} \\ l_3(x) &= \frac{(x-1)(x-3)(x-4)}{(6-1)(6-3)(6-4)} = \frac{(x-1)(x-3)(x-4)}{30} \end{aligned}$$

从而：

$$\begin{aligned} \varphi_3(x) &= \sum_{i=0}^3 y_i l_i(x) \\ &= \frac{7}{30}(x-3)(x-4)(x-6) + \frac{5}{6}(x-1)(x-4)(x-6) \\ &\quad- \frac{4}{3}(x-1)(x-3)(x-6) + \frac{7}{15}(x-1)(x-3)(x-4) \\ &= \frac{1}{5}(x^{3}-13x^{2}+69x-92) \end{aligned}$$

从而：

$$f(2) \approx \varphi_3(2) = 0.4$$

题 5.2

$$p(t) = b-\frac{3(b-t)^2}{b-a}+\frac{2(b-t)^3}{(b-a)^2}$$

对 $t$ 求导得：

$$p'(t) = 6\frac{b-t}{b-a}-6(\frac{b-t}{b-a})^2$$

从而

$$p(a) = b-3(b-a)+2(b-a) = a, \quad p(b) = b$$ $$p'(a) = 0, \quad p'(b) = 0$$

令 $u = \dfrac{b-t}{b-a}$，在 $[a,b]$ 上 $u \in [0, 1]$。从而：

$$|p'(t)| = |6u - 6u^2| \le \frac{3}{2}$$

当且仅当顶点 $u=\dfrac{1}{2}$ 处取等，也即 $t=\dfrac{b-a}{2}$ 处取等。证毕。

题 5.3

$$f(x)= \begin{cases} 3+x-9x^{2}, & x\in[0,1] \\ a+b(x-1)+c(x-1)^{2}+d(x-1)^{3}, & x\in[1,2] \end{cases}$$ $$f'(x)= \begin{cases} 1-18x, & x\in[0,1] \\ b-2c(x-1)+3d(x-1)^{2}, & x\in[1,2] \end{cases}$$ $$f''(x)= \begin{cases} -18, & x\in[0,1] \\ 2c+6d(x-1), & x\in[1,2] \end{cases}$$

又因为：

$$\begin{cases} f(1^-)=f(1^+) \\ f'(1^-)=f'(1^+) \\ f''(1^-)=f''(1^+) \end{cases}$$

所以：

$$\begin{cases} -5=a \\ -17=b \\ -18=2c \end{cases}$$

解得：

$$a=-5,\quad b=-17,\quad c=-9$$

从而：

$$f''(x)= \begin{cases} -18, & x\in[0,1] \\ -18+6d(x-1), & x\in[1,2] \end{cases}$$

故：

$$\begin{aligned} \int_0^2[f''(x)]^2\,\mathrm{d}x &=324+\int_1^2\bigl[-18+6d(x-1)\bigr]^2\,\mathrm{d}x \\ &=324+\int_0^1\bigl(324-216du+36d^2u^2\bigr)\,\mathrm{d}u \\ &=12d^2-108d+648 \end{aligned}$$

故当 $\int_0^2[f''(x)]^2\,\mathrm{d}x$ 取极小值时，有：

$$d=\frac{9}{2}$$

题 5.4

取 $s=t=n+1$，并固定 $q_{n+1}=1,\;p_{n+1}=L$。此时 $\displaystyle\lim_{x\to\infty}R(x)=L$ 自动满足。注意这对 $L=0$ 也是成立的。

设：

$$P(x) = LQ(x) + S(x), \quad R(x) = L + \frac{S(x)}{Q(x)}$$

其中 $\deg Q = n+1,\; \deg S \le n$。

尝试构造：

$$\begin{cases} Q(x_i) = 1, & i=1, 2, \cdots, n \\ S(x_i) = y_i - L, & i=1, 2, \cdots, n \\ \deg Q = n+1 \\ \deg S \le n \end{cases}$$

对于 $Q(x)$，可取：

$$Q(x) = 1+x\prod_{j=1}^{n}(x-x_j), \quad \deg Q = n+1$$

对于 $S(x)$，可由 Lagrange 插值得：

$$S(x) = \sum_{i=1}^{n}(y_i-L)\,l_i(x), \quad \deg S = n-1$$

易知此时满足提示要求！从而：

$$\begin{aligned} Q(x) &= 1+x\prod_{j=1}^{n}(x-x_j) \\ S(x) &= \sum_{i=1}^{n}(y_i-L)\,l_i(x) \\ P(x) &= LQ(x) + S(x) \end{aligned}$$

题 5.5

算法思路

在田径运动中，常规百米成绩包含起跑反应时间。为得到运动员纯跑100米的最短时间，需从分段计时数据中分离出起跑反应时 $t_{\text{react}}$。

设 $x_i = 10, 20, \dots, 100\,\text{m}$ 为分段距离，$t_i$ 为到达 $x_i$ 的累计时间。一般思路是插值 $x(t)$ 并求解 $x(t_{\text{react}}) = 0$，但这涉及非线性求根，存在初值敏感和根不唯一的问题。

所以采用反向思路，将 $t$ 视为 $x$ 的函数，即 $t = t(x)$，通过插值已知点 $(x_i, t_i)$ 后直接外推 $x=0$ 处的值 $t(0)$，从而避免求根运算。

算法设计

基本方法

对数据点 $\{(x_i, t_i)\}_{i=1}^{10}$ 构造三次样条插值 $S(x)$，满足：

• $S(x_i) = t_i$
• $S \in C^2[0, 100]$

利用样条函数的自然外推性质，直接计算：

$$t_{\text{react}}^{\text{(est)}} = S(0)$$

纯跑时间估计为：

$$T_{\text{pure}} = t_{100} - t_{\text{react}}^{\text{(est)}}$$

改进：局部样条外推

全局样条使用全部 10 个数据点，但 10–100 m 阶段 $t(x)$ 接近线性，而 0–10 m 起跑阶段呈强非线性。全局样条被大尺度线性段”拉直”，导致外推至 $x=0$ 时严重高估。

因此仅使用前 $m$ 个数据点（即 $x = 10, 20, \dots, 10m\,\text{m}$）构造局部样条，再外推至 $x=0$。

具体步骤：

1. 对 $m = 3, 4, \dots, 9$，分别用前 $m$ 个点构造三次样条 $S_m(x)$；
2. 计算外推值 $t_{\text{react}}^{(m)} = S_m(0)$；
3. 选取使 $|t_{\text{react}}^{(m)} - t_{\text{react}}^{\text{(lit)}}|$ 最小的 $m$ 作为最优参数。

实验结果

采用三组百米分段数据（单位：s），距离点统一为 $x = [10, 20, \dots, 100]\,\text{m}$：

Dataset	累计时间 $t_i$	文献反应时 $t_{\text{react}}^{\text{(lit)}}$
1	1.89, 2.88, 3.78, 4.64, 5.47, 6.29, 7.10, 7.92, 8.75, 9.58	0.146 s
2	1.85, 2.87, 3.78, 4.65, 5.50, 6.32, 7.14, 7.96, 8.79, 9.69	0.165 s
3	1.727, 2.825, 3.849, 4.852, 5.839, 6.841, 7.837, 8.863, 9.903, 11.015	0.186 s

Dataset	使用点数 $m$	估计反应时	文献值	误差 (ms)	纯跑时间
1	9	0.1780 s	0.146 s	32.0	9.402 s
2	3	0.1614 s	0.165 s	3.6	9.529 s
3	3	0.1792 s	0.186 s	6.8	10.836 s

• 局部最优样条平均绝对误差：14.1 ms
• 局部最优样条最大绝对误差：32.0 ms（Dataset 1）

代码

import numpy as np
from scipy.interpolate import CubicSpline
import matplotlib.pyplot as plt

# 三组数据
datasets = [
    {
        'split_times': np.array([1.89, 2.88, 3.78, 4.64, 5.47, 6.29, 7.10, 7.92, 8.75, 9.58]),
        't_react_lit': 0.146,
        'label': 'Dataset 1'
    },
    {
        'split_times': np.array([1.85, 2.87, 3.78, 4.65, 5.50, 6.32, 7.14, 7.96, 8.79, 9.69]),
        't_react_lit': 0.165,
        'label': 'Dataset 2'
    },
    {
        'split_times': np.array([1.541, 2.639, 3.663, 4.666, 5.653, 6.655, 7.651, 8.677, 9.717, 10.829]) + 0.186,
        't_react_lit': 0.186,
        'label': 'Dataset 3'
    }
]

x = np.array([10, 20, 30, 40, 50, 60, 70, 80, 90, 100])
u = np.sqrt(x)

print("=" * 60)
print("方案：t(u) 三次样条外推 —— 前m个点最优值")
print("=" * 60)

fig, axes = plt.subplots(1, 3, figsize=(18, 5))
fig.suptitle('Optimal m for Cubic Spline Extrapolation t(u) → t(0)',
             fontsize=14, fontweight='bold')

summary_best = []

for idx, data in enumerate(datasets):
    split_times = data['split_times']
    t_react_lit = data['t_react_lit']
    label = data['label']

    # 遍历 m = 2 到 10，找最优的
    best_m = None
    best_t_react = None
    best_error = np.inf

    print(f"\n{label}:")
    print(f"  {'m':>3} {'t_react (s)':>12} {'误差 (ms)':>10}")
    print(f"  " + "-" * 27)

    for m in range(2, 11):
        u_m = u[:m]
        t_m = split_times[:m]

        cs = CubicSpline(u_m, t_m)
        t_react_m = cs(0)
        error_m = abs(t_react_m - t_react_lit) * 1000

        print(f"  {m:>3} {t_react_m:>12.4f} {error_m:>10.1f}")

        if error_m < best_error:
            best_error = error_m
            best_m = m
            best_t_react = t_react_m

    T_pure = split_times[-1] - best_t_react
    summary_best.append((label, best_m, best_t_react, best_error, T_pure))

    print(
        f"  → 最优 m = {best_m}, t_react = {best_t_react:.4f}s, 误差 = {best_error:.1f}ms")

    # 用最优 m 绘图
    u_m = u[:best_m]
    t_m = split_times[:best_m]
    cs = CubicSpline(u_m, t_m)

    ax = axes[idx]
    ax.scatter(u_m, t_m, color='red', s=60, zorder=5, label=f'First {best_m} points')
    ax.scatter(u[best_m:], split_times[best_m:], color='gray', s=40, zorder=4,
               alpha=0.4, label='Unused points')
    ax.scatter([0], [0], color='gray', marker='x', s=80, label='Start line')

    u_fine = np.linspace(0, 10, 500)
    ax.plot(u_fine, cs(u_fine), 'b-', lw=2, label='Cubic spline t(u)')
    ax.axvline(x=0, color='gray', ls='-', lw=0.5)
    ax.scatter([0], [best_t_react], color='blue', marker='*', s=150, zorder=6,
               label=f'Est. $t_{{react}}$={best_t_react:.3f}s')
    ax.scatter([0], [t_react_lit], color='orange', marker='D', s=60, zorder=6,
               label=f'True $t_{{react}}$={t_react_lit:.3f}s')

    ax.text(5, max(split_times)*0.15,
            f'm={best_m}\nError: {best_error:.1f} ms',
            fontsize=10, ha='center',
            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.7))

    ax.set_xlabel('u = √x (m^{1/2})', fontsize=11)
    ax.set_ylabel('Cumulative Time t (s)', fontsize=11)
    ax.set_title(f'{label}\nOptimal m = {best_m}', fontsize=12)
    ax.legend(loc='lower right', fontsize=7)
    ax.set_xlim(-0.5, 10.5)
    ax.grid(True, alpha=0.3)

# 汇总
print("\n" + "=" * 60)
print("汇总对比")
print("=" * 60)
print(f"{'Dataset':<12} {'m':>3} {'Est. t_react':>12} {'True t_react':>12} {'Error (ms)':>10} {'T_pure':>8}")
print("-" * 60)

for label, best_m, best_t_react, best_error, T_pure in summary_best:
    data = datasets[[d['label'] for d in datasets].index(label)]
    print(
        f"{label:<12} {best_m:>3} {best_t_react:>12.4f} {data['t_react_lit']:>12.4f} {best_error:>10.1f} {T_pure:>8.4f}")

errors = [s[3] for s in summary_best]
print(f"\n平均误差: {np.mean(errors):.1f} ms")
print(f"最大误差: {np.max(errors):.1f} ms")

plt.tight_layout(rect=[0, 0.03, 1, 0.92])
plt.savefig('./sprint_optimal_m.png', dpi=200, bbox_inches='tight')
plt.show()

print("\n图表已保存")