本题考虑对于定义在[−1,1]上的一个光滑函数𝑓(𝑥)的三次样条插值的使用。下面所说的误差 都是指绝对误差。 (a) （10 分）仿照课本推导出关于额外给定边界点处（即−1和1）三次样条插值多项式的一 次导数值时其在各插值点上的二次导数值应该满足的线性方程组。请给出推导过程。 (b) （10 分）令三次样条插值多项式在−1和1处的导数为0，基于上一问中的结果使用𝑛 = 2 4个子区间插值一个定义在[−1,1]上的函数𝑓(𝑥) = 𝑠𝑖𝑛(4𝑥 2 ) + 𝑠𝑖𝑛2 (4𝑥)并使用 semilogy 图通过在2000个等距点上取真实值画出你构造的三次样条插值的逐点误差。 (c) （15 分）使用不同的𝑛，令𝑛 = 2 4 , 2 5 , ⋯ , 2 10重复上一问，取关于不同𝑛的2000个等距点 上的误差的最大值，用 loglog 图描述插值区间上最大误差值随𝑛变化的情况（横轴是𝑛）。 用比较新的matlab语言来写，我只要matlab程序

(a) 求解三次样条插值多项式的一次导数值时其在各插值点上的二次导数值应该满足的线性方程组，可以仿照课本中的推导方法，得到以下线性方程组：

$$
\begin{cases}
h_1M_0 + 2(h_1+h_2)M_1 + h_2M_2 = 3\frac{f(x_1)-f(x_0)}{h_1^2} - 3\frac{f(x_0)-f_{-1}}{h_0^2} \\
h_2M_1 + 2(h_2+h_3)M_2 + h_3M_3 = 3\frac{f(x_3)-f(x_2)}{h_3^2} - 3\frac{f(x_2)-f(x_1)}{h_2^2}
\end{cases}
$$

其中，$M_0$ 和 $M_3$ 分别为边界点 $x_{-1}$ 和 $x_4$ 处的二次导数值，$h_i = x_{i+1} - x_i$，$f_{-1}$ 和 $f_4$ 分别为边界点处的函数值。

(b) 根据题目中的要求，在边界点 $x_{-1}$ 和 $x_4$ 处，三次样条插值多项式的一阶导数值为0。因此，我们可以先求出满足这个条件的 $M_0$ 和 $M_3$ 的值，然后在每个子区间上分别插值，得到整个区间上的三次样条插值多项式。

具体地，我们可以先计算出 $M_0$ 和 $M_3$ 的值：

$$
\begin{cases}
2M_0 + M_1 = \frac{3}{h_1^2}(f(x_1)-f_{-1}) \\
h_2M_1 + 2(h_2+h_3)M_2 + h_3M_3 = 3\frac{f(x_3)-f(x_2)}{h_3^2} - 3\frac{f(x_2)-f(x_1)}{h_2^2} \\
M_2 + 2M_3 = \frac{3}{h_3^2}(f_4-f(x_3))
\end{cases}
$$

然后，在每个子区间 $[x_i, x_{i+1}]$ 上，我们可以使用课本中的方法，得到三次样条插值多项式的表达式：

$$
S_i(x) = \frac{M_{i-1}}{6h_i}(x_{i+1}-x)^3 + \frac{M_i}{6h_i}(x-x_i)^3 + \frac{f(x_i)}{h_i}(x_{i+1}-x) + \frac{f(x_{i+1})}{h_i}(x-x_i)
$$

最后，我们可以在2000个等距点上计算出真实值和插值值的差，得到逐点误差，并使用 semilogy 图来表示。具体的 MATLAB 代码如下：

% (b) 插值
n = 2^4; % 子区间数
x = linspace(-1, 1, n+1); % 等距插值点
h = diff(x);
f = sin(4*x.^2) + sin(2*4*x);
f1 = [0 (3./h(2:end)).*(f(3:end)-f(2:end-1))-(3./h(1:end-1)).*(f(2:end-1)-f(1:end-2)) 0]; % 一次导数值
M = zeros(1, n+1);
A = zeros(n-1, n-1);
B = zeros(n-1, 1);
for i = 1:n-1
    A(i, i) = 2*(h(i)+h(i+1));
    if i > 1
        A(i, i-1) = h(i);
    end
    if i < n-1
        A(i, i+1) = h(i+1);
    end
    B(i) = f1(i+1);
end
M(2:end-1) = A\B;
S = zeros(1, 2000);
x_interp = linspace(-1, 1, 2000);
for i = 1:n
    idx = x_interp >= x(i) & x_interp <= x(i+1);
    S(idx) = (M(i)/(6*h(i))*(x_interp(idx)-x(i+1)).^3 + M(i+1)/(6*h(i))*(x(i+1)-x_interp(idx)).^3 + f(i)/h(i)*(x_interp(idx)-x(i+1)) + f(i+1)/h(i)*(x(i+1)-x_interp(idx)));
end
error = abs(S - (sin(4*x_interp.^2) + sin(2*4*x_interp)));
semilogy(x_interp, error);

(c) 接下来，我们可以对不同的 $n$ 值（$n=2^4,2^5,\cdots,2^{10}$）重复上述插值过程，并在2000个等距点上计算出逐点误差的最大值。然后，我们可以使用 loglog 图来描述插值区间上最大误差值随 $n$ 变化的情况。

具体的 MATLAB 代码如下：

% (c) 改变 n 的值，计算最大误差
n_values = 2.^(4:10);
max_errors = zeros(size(n_values));
for k = 1:length(n_values)
    n = n_values(k);
    x = linspace(-1, 1, n+1); % 等距插值点
    h = diff(x);
    f = sin(4*x.^2) + sin(2*4*x);
    f1 = [0 (3./h(2:end)).*(f(3:end)-f(2:end-1))-(3./h(1:end-1)).*(f(2:end-1)-f(1:end-2)) 0]; % 一次导数值
    M = zeros(1, n+1);
    A = zeros(n-1, n-1);
    B = zeros(n-1, 1);
    for i = 1:n-1
        A(i, i) = 2*(h(i)+h(i+1));
        if i > 1
            A(i, i-1) = h(i);
        end
        if i < n-1
            A(i, i+1) = h(i+1);
        end
        B(i) = f1(i+1);
    end
    M(2:end-1) = A\B;
    S = zeros(1, 2000);
    x_interp = linspace(-1, 1, 2000);
    for i = 1:n
        idx = x_interp >= x(i) & x_interp <= x(i+1);
        S(idx) = (M(i)/(6*h(i))*(x_interp(idx)-x(i+1)).^3 + M(i+1)/(6*h(i))*(x(i+1)-x_interp(idx)).^3 + f(i)/h(i)*(x_interp(idx)-x(i+1)) + f(i+1)/h(i)*(x(i+1)-x_interp(idx)));
    end
    error = abs(S - (sin(4*x_interp.^2) + sin(2*4*x_interp)));
    max_errors(k) = max(error);
end
loglog(n_values, max_errors);

运行上述程序可以得到逐点误差和最大误差随 $n$ 变化的图像。