XMU《计算方法》实验一 三次样条插值算法

实验一 三次样条插值算法

一、Matlab 代码

clear;
x = input('请输入插值结点的 x：');
y = input('请输入插值结点的 y：');

[x, I] = sort(x);
y = y(I);

if length(y) ~= length(x)
    error('x 和 y 的数量不相等！');
end
n = length(x) - 1;
N = n * 4;

% 函数值约束
A = [];
b = [];
for i = 1:n
    a = zeros(1, N);
    a(4*i-3 : 4*i) = [x(i)^3, x(i)^2, x(i), 1];
    A= [A; a]; b = [b; y(i)];
    a(4*i-3 : 4*i) = [x(i+1)^3, x(i+1)^2, x(i+1), 1];
    A = [A; a]; b = [b; y(i+1)];
end

% 一阶导数连续约束
for i = 2:n
    a = zeros(1, N);
    a(4*i-7 : 4*i) = [3*x(i)^2, 2*x(i), 1, 0, -3*x(i)^2, -2*x(i), -1, 0];
    A = [A; a]; b = [b; 0];
end

% 二阶导数连续约束
for i = 2:n
    a = zeros(1, N);
    a(4*i-7 : 4*i) = [6*x(i), 2, 0, 0, -6*x(i), -2, 0, 0];
    A = [A; a]; b = [b; 0];
end

type = input('请输入边界条件类型（0: 自然边界条件；1: 一阶导；2: 二阶导）：');
if type == 0
    % 自然边界条件约束
    a = zeros(1, N);
    a(1:4) = [6*x(1), 2, 0, 0];
    A = [A; a]; b = [b; 0];
    a = zeros(1, N);
    a(4*n-3 : 4*n) = [6*x(n+1), 2, 0, 0];
    A = [A; a]; b = [b; 0];
elseif type == 1
    y0 = input(strcat('请输入 x=', num2str(x(1)), ' 处的一阶导：'));
    yn = input(strcat('请输入 x=', num2str(x(n+1)), ' 处的一阶导：'));
    a = zeros(1, N);
    a(1:4) = [3*x(1)^2, 2*x(1), 1, 0];
    A = [A; a]; b = [b; y0];
    a = zeros(1, N);
    a(4*n-3 : 4*n) = [3*x(n+1)^2, 2*x(n+1), 1, 0];
    A = [A; a]; b = [b; yn];
else
    y0 = input(strcat('请输入 x=', num2str(x(1)), ' 处的二阶导：'));
    yn = input(strcat('请输入 x=', num2str(x(n+1)), ' 处的二阶导：'));
    a = zeros(1, N);
    a(1:4) = [6*x(1), 2, 0, 0];
    A = [A; a]; b = [b; y0];
    a = zeros(1, N);
    a(4*n-3 : 4*n) = [6*x(n+1), 2, 0, 0];
    A = [A; a]; b = [b; yn];
end

% 解方程并将解对应还原到各个样条三次插值的系数
p = A \ b;
a = []; b = []; c = []; d = [];
for i = 1:n
    a(i) = p(4 * i - 3);
    b(i) = p(4 * i - 2);
    c(i) = p(4 * i - 1);
    d(i) = p(4 * i);
    fprintf('f(x) = %gx^3 + %gx^2 + %gx + %g\tx in [%g, %g]\n', a(i), b(i), c(i), d(i), x(i), x(i+1))
end

% 计算各个给定点的插值结果
X = linspace(x(1), x(n+1), 200);
Y = [];
m = length(X);
for i = 1:m
    for j = 1:n
        if X(i) >= x(j) && X(i) <= x(j+1)
            Y(i) = a(j)*X(i)^3 + b(j)*X(i)^2 + c(j)*X(i) + d(j);
        end
    end
end

% 画图
hold on;
plot(X, Y);
plot(x, y, 'o');

二、算法原理

三次样条插值应该满足以下条件：每一小段都是三次多项式，且整体的二阶导数连续。

思路：

给定插值节点 \((x_0, y_0), (x_1, y_1), \dots, (x_n, y_n)\)。

三次样条插值函数：

\[\begin{align*} f(x)&=a_1 x^3+b_1 x^2+c_1 x+d_1 && x_0\leq x\leq x_1\\ &=a_2 x^3+b_2 x^2+c_2 x+d_2&&x_1\leq x\leq x_2\\ &\vdots\\ &=a_n x^3+b_n x^2+c_n x+d_n&&x_{n-1}\leq x\leq x_n \end{align*} \]

一共需要确定 \(4n\) 个未知数。

限制 1——函数值约束

提供了 \(2n\) 个限制。

\[\begin{align*} a_1 x_0^3+b_1 x_0^2+c_1 x_0+d_1&=y_0\\ a_1 x_1^3+b_1 x_1^2+c_1 x_1+d_1&=y_1\\ &\vdots\\ a_i x_{i-1}^3+b_i x_{i-1}^2+c_i x_{i-1}+d_i&=y_{i-1}\\ a_i x_i^3+b_i x_i^2+c_i x_i+d_i&=y_i\\ &\vdots\\ a_n x_{n-1}^3+b_n x_{n-1}^2+c_n x_{n-1}+d_n&=y_{n-1}\\ a_n x_n^3+b_n x_n^2+c_n x_n+d_n&=y_n \end{align*} \]

限制 2——一阶导数连续约束

提供了 \(n-1\) 个限制。

\[\begin{align*} 3a_1 x_1^2+2b_1 x_1+c_1-3a_2 x_1^2-2b_2 x_1-c_2&=0\\ &\vdots\\ 3a_i x_i^2+2b_i x_i+c_i-3a_{i+1} x_i^2-2b_{i+1} x_i-c_{i+1}&=0\\ &\vdots\\ 3a_{n-1} x_{n-1}^2+2b_{n-1} x_{n-1}+c_{n-1}-3a_n x_{n-1}^2-2b_n x_{n-1}-c_n&=0 \end{align*} \]

限制 3——二阶导数连续

\[\begin{align*} 6a_1 x_1+2b_1-6a_2 x_1-2b_2&=0\\ \vdots\\ 6a_i x_i+2b_i-6a_{i+1} x_i-2b_{i+1}&=0\\ \vdots\\ 6a_{n-1} x_{n-1}+2b_{n-1}-6a_n x_{n-1}-2b_n&=0 \end{align*} \]

限制 4 ——边界条件

在本程序中，我们实现了三种边界条件，分别是边界一阶导约束、边界二阶导约束和自然边界条件。

（1）自然边界条件

\[S''(x_0) = S''(x_n) = 0 \]

化作多项式的形式，即

\[6a_1x_0+2b_1 = 0\\ 6a_nx_n+2b_n = 0 \]

（2）一阶导约束

\[S'(x_0) = f'_0, S'(x_n) = f_n' \]

化作多项式的形式，即

\[3a_1x_0^2 + 2b_1x_0 + c_1 = f_0\\ 3a_nx_n^2 + 2b_nx_n + c_n = f_n \]

（3）二阶导约束

\[S''(x_0) = f''_0， S''(x_n) = f''_n \]

化作多项式的形式，即

\[6a_1x_0+2b_1 = f''_0\\ 6a_nx_n+2b_n = f''_n \]

限制的矩阵形式

因为不同的边界条件对应不同的矩阵最后两行，这里以自然边界条件的矩阵为例。

设 \(\boldsymbol x = [a_1, b_1, c_1, d_1, \cdots, a_n, b_n, c_n, d_n]^T\)。

则上述方程组可以写成 \(\boldsymbol{Ax}=\boldsymbol b\) 的形式，即

\[\begin{bmatrix} x_0^3 & x_0^2 & x_0 & 1 & & & & & \cdots \\ x_1^3 & x_1^2 & x_1 & 1 & & & & & \cdots \\ & & & & x_1^3 & x_1^2 & x_1 & 1 & \cdots \\ & & & & x_2^3 & x_2^2 & x_2 & 1 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ & & & & & & & & \cdots & & & & & x_{n-1}^3 & x_{n-1}^2 & x_{n-1} & 1\\ & & & & & & & & \cdots & & & & & x_{n}^3 & x_{n}^2 & x_{n} & 1\\ 3x_1^2 & 2x_1 & 1 & 0 & -3x_1^2 & -2x_1 & -1 & 0 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ & & & & & & & & \cdots & 3x_{n-1}^2 & 2x_{n-1} & 1 & 0 & -3x_{n-1}^2 & -2x_{n-1} & -1 & 0\\ 6x_1 & 2 & 0 & 0 & -6x_1 & -2 & 0 & 0 & \cdots\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ & & & & & & & & \cdots & 6x_{n-1} & 2 & 0 & 0 & -6x_{n-1} & -2 & 0 & 0\\ 6x_0 & 2 & 0 & 0 & & & & & \cdots\\ & & & & & & & & \cdots & & & & & 6x_n & 2 & 0 & 0 \end{bmatrix} \begin{bmatrix} a_1\\ b_1\\ c_1\\ d_1\\ \vdots\\ a_n\\ b_n\\ c_n\\ d_n \end{bmatrix} = \begin{bmatrix} y_0\\ y_1\\ y_1\\ y_1\\ \vdots\\ y_{n-1}\\ y_n\\ 0\\ \vdots\\ 0\\ 0\\ \vdots\\ 0\\ 0\\ 0 \end{bmatrix} \]

三、程序说明

在“二”的最后，我们得到了一个 \(\boldsymbol{Ax}=\boldsymbol b\) 形式的矩阵。我们只需要在 Matlab 中表示出这个矩阵，通过 Matlab 的 A \ b 运算，即可得到 \(\boldsymbol x\) 的解。

在计算出 \(x\) 之后，我们根据我们设计 \(\boldsymbol x\) 的原理，每相邻的四项属于同一个三次多项式，即可得到我们的三次样条插值函数了。

如果需要计算一个 \(x\) 的插值结果，我们只需要判断 \(x\) 属于样条插值的哪一段，然后使用这一段的三次多项式函数计算，即可得到 \(x\) 的解。

在程序中，我们通过画出图像的方式，直观地表现了三次样条插值的结果。其中，插值节点我们使用圆形点出，插值函数我们使用一条连续实线绘制而成。

值得注意的是，为了提高程序的通用性，我们考虑了输入的插值节点并没有按照 \(x\) 坐标排序的情况，因此我们在程序中重新对所有点的节点按照 \(x\) 坐标进行了排序。

此外，我们考虑到了边界条件的多样性，因此提供了三种不同的边界条件约束，在程序中，我们对用户需要的边界条件约束类型进行提问，一次设计不同的系数矩阵，来计算插值的结果。

程序中，除了解线性方程组和使用 plot 绘制图像的过程，全部独立完成，包括线性方程组的建立，矩阵的描述和表示，以及计算插值，计算绘制图像需要的点坐标的过程，全部是我们自己实现的，没有调用任何的 Matlab 工具。

四、运行结果示例

示例 1

示例 2

示例 3