多项式插值（数值计算方法）Matlab实现

一. 原理介绍

1. 关于插值的定义及基本原理可以参照如下索引
   插值原理（数值计算方法）
2. 前面已经介绍过插值原理的唯一性表述，对于分立的数据点，方程组：

$$\begin{array}{rl}& P(x_0)=y_0\Rightarrow a_0+a_1{x}_0+a_2{x}_0^2+\cdots +a_n{x}_0^n=y_0,\\ & P(x_1)=y_1\Rightarrow a_0+a_1{x}_1+a_2{x}_1^2+\cdots +a_n{x}_1^n=y_1,\\ & ⋮\\ & P(x_n)=y_n\Rightarrow a_0+a_1{x}_n+a_2{x}_n^2+\cdots +a_n{x}_n^n=y_n.\end{array}$$

恒有解，多项式插值的目标即为在这一过程中求解系数$$a_0、a_1、...、a_n⟺\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_n\end{array}\right]$$

3. 即解方程组：
   $$\left[\begin{array}{ccccc}1& x_0& x_0^2& \cdots & x_0^n\\ 1& x_1& x_1^2& \cdots & x_1^n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_n& x_n^2& \cdots & x_n^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_n\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ ⋮\\ y_n\end{array}\right].$$

关于该方程组的解法在线性代数中有多种，这里主要提及两种：
①高斯消元法
②克莱姆法则
程序设计过程中一般有封装好的库函数，如果为了考虑减少库依赖和提高程序运行效率及占用可能会用到上述方法（这里就不详细展开了）

二. 程序设计

1. 构建矩阵

% 构造Vandermonde矩阵A
A = zeros(n, n);
for i = 1:n
    for j = 1:n
        A(i, j) = x_data(i)^(j-1);  % Vandermonde矩阵
    end
end

Ⅰ 构建一个$(n\times n)$的矩阵 A 来描述多项式矩阵：
$$\left[\begin{array}{ccccc}1& x_0& x_0^2& \cdots & x_0^n\\ 1& x_1& x_1^2& \cdots & x_1^n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_n& x_n^2& \cdots & x_n^n\end{array}\right]$$

其中：
$$A[i][j]=x_i^{j-1}$$

式中第一列的1是通过$x_i^0$得到的。

y_data = data(:, 2);

Ⅱ 构建系数矩阵 B，即原始数据对应的 y 值：

2. 求解矩阵方程

% 解线性方程组 A * coefficients = y_data
coefficients = A \ y_data;

注释：该部分通过反斜杠运算符 \ 计算线性方程组 $A\cdot coefficients=y_{d}ata$的解。
①当方程组超定时（方程数大于未知数个数），返回最小二乘解，即最小化残差平方和 $\Vert A\cdot coefficients-y_{data}{\Vert }^2$ ；
②当方程组适定时，返回精确解；
③当方程组欠定时返回最小范数解。

求解出矩阵形式形如
6.0000
-7.8333
4.5000
-0.6667
从上至下为最低次项到最高次项系数

3. 作出多项式函数

% 生成插值多项式的x和y值
x_vals = linspace(min(x_data) - 1, max(x_data) + 1, 500);
y_vals = polyval(flip(coefficients), x_vals);  % 计算插值多项式的y值

注释: 前面注释提到coefficients数组中的系数对应从左到右为最低到最高次项系数，而函数polyval()要求输入具有逆序的项系数：flip函数将系数的顺序反转，将变为从最高次到最低次项系数
y_vals = polyval(flip(coefficients), x_vals) 将计算每一个 x_val 对应的多项式值，并返回一个 y_vals 数组，包含每个 x_val 对应的 y 值。

4. 绘制插值曲线

% 绘制插值曲线
figure;
plot(x_vals, y_vals, 'b-', 'DisplayName', '插值曲线');
hold on;
scatter(x_data, y_data, 'ro', 'DisplayName', '数据点');
title('插值多项式');
xlabel('X轴');
ylabel('Y轴');
legend;
grid on;

5. 完整代码

% 输入数据 (x, y)
data = [
    1,2
    2,3
    3,5
    4,4
];

% 提取x和y值
x_data = data(:, 1);
y_data = data(:, 2);
n = length(data);

% 构造Vandermonde矩阵A
A = zeros(n, n);
for i = 1:n
    for j = 1:n
        A(i, j) = x_data(i)^(j-1);  % Vandermonde矩阵
    end
end

% 解线性方程组 A * coefficients = y_data
coefficients = A \ y_data;

% 输出插值多项式的系数
disp('插值多项式的系数:');
disp(coefficients);

% 生成插值多项式的x和y值
x_vals = linspace(min(x_data) - 1, max(x_data) + 1, 500);
y_vals = polyval(flip(coefficients), x_vals);  % 计算插值多项式的y值

% 绘制插值曲线
figure;
plot(x_vals, y_vals, 'b-', 'DisplayName', '插值曲线');
hold on;
scatter(x_data, y_data, 'ro', 'DisplayName', '数据点');
title('插值多项式');
xlabel('X轴');
ylabel('Y轴');
legend;
grid on;

三. 图例

这要求我们的输入数据都具有上述形式：

data = [
    x_1, y_1
    x_2, y_2
    x_3, y_3
    x_4, y_4
    ...]

最后我们插值一组随机生成的测试数据

data = [
    7.264384, 3.931292
    1.943873, 6.218903
    8.384019, 2.584103
    5.672210, 9.032674
    0.294315, 4.726018
    6.129531, 7.912846
    9.516347, 1.478264
    3.824679, 5.596042
]

实际应用时应避免数据点过多导致的多项式次数过高

希望能够帮到迷途之中的你，知识有限，如有学术错误请及时指正，感谢大家的阅读

(^^)／▽ ▽＼(^^)