MATLAB之多项式插值

一、算法原理

1、插值问题定义

当精确函数 y = f(x) 非常复杂或未知时，在区间【a，b】上一系列节点x0,x1,x2,......xn处测得函数值f(x0)、f(x1)、f(x2)......f(xn)；

由此构造一个简单的简单易算的近似函数g(x)≈f(x)，满足条件g(xi)=f(xi) (i=0,1,2...,n)。这个问题就被称为插值问题。

其中，g(x)被称为插值函数；

x0,x1,x2,......xn位插值节点；

g(x)≈f(x)为插值条件；

区间【a,b】为插值区间；

2、插值函数分类

插值函数有多项式插值、拉格朗日插值、牛顿插值、埃米特插值，最小二乘法插值等，这里我们来介绍多项式插值。其他的插值方法请关注博主的其他文章。

3、多项式插值原理

函数解析式未知，但已知一些列点的函数值。如下表所示，对于n+1个点，我们可以找到一个次数不超过n的插值多项式。

$f=an*X^n+a(n-1)*X^(n-1)+....+a1*X+a0$ ，称f为x的n次插值多项式。

x0	x1	x2	x3	......	xn
f(x0)	f(x1)	f(x2)	f(x3)	......	f(xn)

将表中的n+1个点带入，即可得到n+1个插值多项式，

an*X0^n+a(n-1)*X0^(n-1)+....+a1*X0+a0=f(x1);

an*X1^n+a(n-1)*X1^(n-1)+....+a1*X1+a0=f(x1);

......

an*Xn^n+a(n-1)*Xn^(n-1)+....+a1*Xn+a0=f(x1);

在该插值多项式中x，f(x)已知,未知系数为a0 a1 a2...an;

将其写成矩阵形式：

$\begin{bmatrix}x_{0}^n & x_{0}^{n-1} & ... & x_{0}^{1}& 1 \\ x_{1}^n & x_{1}^{n-1}& ... &x_{1}^{1}& 1 \\ ... & ... & ... &... & ...\\ x_{(n-1)}^n & x_{n-1}^{n-1}& ...&x_{n-1}^{1}& 1 \\ x_{n}^n & x_{n}^{n-1} & ... & x_{n}^1 & 1 \end{bmatrix}$ $*\begin{bmatrix}a_{n} \\ a_{n-1} \\ ... \\ a_{1} \\ a_{0} \end{bmatrix}$ $=\begin{bmatrix}f(x_{0}) \\ f(x_{1}) \\ ... \\ f(x_{n-1}) \\ f(x_{n}) \end{bmatrix}$

x组成的行列式为范德蒙德行列式 =可写成vander(x0 x1 ... x(n-1) xn) ,系数an...a1 a0 为列矩阵，f(x) 为列矩阵

注：因为x0....xn互不相同，多以该方程组的系数行列式不为零，方程有唯一解。

二、matlab代码

clc
clear
close
x=[1 2 3 4 5 6];   %一系列插值点
y=[2 5 6 9 14 26];
plot(x,y,'o');  %绘制出散点图
hold on
A=vander(x); %利用范德蒙矩阵生成插值多项式系数矩阵
%[x1^n x1^(n-1)......x1 1     [an
%  x2^n x2^(n-1)......x2 1     a(n-1)
%   ......                  *  ......    =y
%   ......                     ......
%  xk^n xk^(n-1)......xk 1]    a0]
y=y';
B=[A,y]; %生成增广矩阵
C=rref(B);%  化最简型，解方程组
D=C(:,end);
fprintf('\n  这是一个%d阶的多项式，从高阶到低阶排列！\n',length(y)-1);
x=x(1):0.01:x(end); %绘制拟合曲线
y=polyval(D,x)% 写出拟合曲线方程
plot(x,y) %绘图

三、matlab自带插值函数

多项式插值函数 polyfit

P=polyfit(x,y,n) %n为拟合阶数，P为多项式的系数,x和y为原始数据点。

[p,S,mu] = polyfit(x,y,n)

S是规模为1×1的结构数组，包括系数矩阵的QR分解的上三角阵，自由度，拟合误差平方和的算术平方根。

mu为居中和缩放参数，mu (1)是平均值，mu (2)是标准差。这种定心和缩放变换改善了多项式和拟合算法的数值特性。

matlab代码

%% 多项式插值 polyfit   polyval
close
clc
clear
x=linspace(0,2*pi,10); %取y=sinx上的十个数据点
y=sin(x);
P=polyfit(x,y,7) %7为拟合阶数，p为多项式的系数
poly2str(P,'x'); %组装成插值多项式
x1=linspace(0,2.5*pi,100);
y1=polyval(P,x1);%用x1代替插值多项式中的x
plot(x,y,'o',x1,y1)x=linspace(1000,2000,10);  %初始10个数据点
y=1e5*[0.02 0.04 0.05 0.055 0.5 0.6 0.9 1.52 2.3 3.6];
plot(x,y,'o')%数据点绘图
hold on
x1=linspace(1000,2000,100);
[P,~,mu]=polyfit(x,y,6); %mu为居中和缩放参数 具体查看函数定义
y1=polyval(P,x1,[],mu);
plot(x1,y1)