Matlab 实现几种常用插值法

插值法 在数值算例中有广泛的应用。

插值的作法,很直觀地來講,就是,(1) 先從表列值來獲得函數 f(x),再 (2) 用函數 f(x) 求出我們所要的任何特定 x 之 f(x) 函數值。然而,比較精密且系統化的數值方法卻不是用這兩個步驟來進行插值,原因是前述兩階段方法對於插值的精密度並沒有控制,效率較差,也比較會有進位誤差。一般在做插值法,是從欲插值點 x 附近的幾個表列點 xi 開始,建立插值函數 f(x),並且也試著網羅更多表列點來插值,看隨著項數變多誤差會不會變小,如此找出最適合的函數 f(x)。

Reference from Interpolation and Extrapolation,这是一个对插值法探讨的比较好的网站,也提供了相对的程式供使用,非常推荐参考。

函数逼近 在很多方面与插值法比较类似,也具有非常广泛的用处,尤其在近些年的前沿问题深度学习中有着重要的作用:什么是深度学习?(从函数逼近论的角度来理解) - 中国科学技术大学。函数逼近我也会在深度学习中进行更多的讨论。

本文结构

  • 插值法
    • Lagrange 插值法
    • Newton 插值法
    • Hermite 插值法
    • 三次样条插值

在这篇 post 中,将不会对各种插值法的原理进行太多的讨论,只是将相关资料与 Matlab Code 附于此处,供读者参考。能力有限,可能会有一些错误在里面,希望可以通过留言告知,也请多多理解。

若公式无法正常加载,请刷新重试。

正文

Lagrange 插值法

Lagrange 插值法 - Wiki

Lagrange 插值多项式(n = 2 时): $$ p_2(x)=y_0 \times \frac{(x-x_1) \times (x - x_2)}{(x_0 - x_1) \times (x_0-x_2)}+y_1 \times \frac{(x-x_0) \times (x - x_2)}{(x_0 - x_0) \times (x_0-x_2)}+y_2 \times \frac{(x-x_0) \times (x - x_1)}{(x_0 - x_0) \times (x_0-x_1)} $$ 余项(次数为 n 时):

$$ f(x)=p_n(x)+\frac{f^{(n+1)}(\xi)}{(n+1)!}\times w_{n+1}(x), a<\xi <b. $$ 其中:

$$ w_{n+1}(x)=(x-x_0)(x-x_1)\cdots (x-x_n) $$ Matlab 实现:

function [y] = Lagrange (x0, y0)
% Lagrange 插值法
% input: (x0, y0): 插值点
% output: y: Lagrange 插值函数
n = length (x0);
syms x;
s = 0.0;
for k = 1:n
    p = 1.0;
    for j = 1:n
        if j ~= k
            p = p * (x - x0(j)) / (x0(k) - x0(j));
        end % if
    end % for
    s = p * y0(k) + s;
end % for
y = s;
end % function

本函数直接返回 Lagrange 插值多项式,亦有其他输出方法,也可将误差考虑在内,参考 拉格朗日(lagrange)插值及其MATLAB程序

Newton 插值法

牛顿多项式 - Wiki

均差:

  • 一阶均差:$f[x_i,x_j]=\frac{f(x_j)-f(x_i)}{x_j-x_i}$
  • 二阶均差:$f[x_i,x_j,x_k]=\frac{f[x_j,x_k]-f[x_i,x_j]}{x_k-x_i}$
  • n 阶均差:$f[x_0,x_1,…,x_n]=\frac{f[x_1,x_2,…,x_n]-f[x_0,x_1,…,x_{n-1}]}{x_n-x_0}$

Newton 插值多项式:

$$ f(x)=N_n(x)+f[x,x_0,x_1,…,x_n]w_{n+1}(x) $$ 其中:

$$ N_n(x)=f(x_0)+f[x_0,x_1]w_1(x)+\dots+f[x_0,x_1,…,x_n]w_n(x) $$

  • 输入:数 $x_0,x_1,…,x_n,x$;函数值 $f(x_0),f(x_1),…,f(x_n)$ 作为 Q 的第一列元素 $Q_{0,0},Q_{1,0},…,Q_{n,0}$

  • 输出:$f(x)$ 的近似值 $b_0$

Matlab 实现:

function [ y ] = Newton_Insert( x0, y0 )
% Newton 均差插值法
% input: (x0, y0): 插值点
% output: y: f(x) 近似值
n = length(x0);
for i = 0 : n
    Q(i, 0) = y0(i);
end % for
for i = 1 : n
    for j = 1 : i
        Q(i, j) = (Q(i, j - 1) - Q(i - 1, j - 1)) / (x(i) - x(i - j));
    end % for
end % for
b(n) = Q(n, n);
for k = n : -1 : 1
    b(k - 1) = Q(k - 1, k - 1) + b(k) * (x - x(k - 1));
end % for
y = b(0);
end % function

Hermit 插值法

Hermit 插值法 - Wiki

前面介绍的插值公式,都只要求插值多项式在插值基点处取给定的函数值,在许多实际问题中,不仅要求插值多项式与函数在各插值基点值相等,还需要在相应插值基点处有若干阶导数相等,这类问题称为 Hermite 插值问题。

在这里只考虑在各基点均有一阶导数的 Hermite 插值问题。

function H = Hermite( x0, y0, z0 )
% Hermit 插值法
% input: (x0, y0): 插值点; (x0, z0): 对应点一阶导
% output: H: f(x) 近似值
n=length(x0);
syms x;
s=0.0;
    for k=1:n
        p=1;
        for j=1:n
            if j~=k
                p=(p*(x-x0(j)))/(x0(k)-x0(j));
            end % if
        end % for
        b(k)=p*p*(x-x0(k));
        q=0;
        for j=1:n
            if j~=k
                q=q+1/(x0(k)-x0(j));
            end %if
        end % for
        a(k)=p*p*(1-(2*(x-x0(k))*q));
        s=a(k)*y0(k)+b(k)*z0(k)+s;
    end % for
    H=s;
end % function

三次样条插值法

样条插值 - Wiki

样条插值是分段多项式插值的优化,其中有分段线性插值、分段抛物线插值等多种插值方式。为使分段插值函数在插值区间的子区间的端点连续而使用条样插值的方法。

函数:

MathJax BUG

> \begin{equation} 
> S(x)=
> \left\{
> \begin{array}{lr}
> S_1(x), x\in [x_1,x_2],&  \\
> ... & \\
> S_i(x), x\in [x_i,x_{i+1}],&  \\
> ... & \\
> S_n(x), x\in [x_n,x_{n+1}]& \\
> \end{array}
> \right. 
> \end{equation}

$$ a=x_1<x_2<…<x_{n+1}=b $$

是二次连续可微的,右方函数都是不高于三次的多项式或者零多项式,且满足条件:

$$ S(x_j)=f(x_j),j=1,…,n+1. $$ 则称 S(x) 为函数 f(x) 的三次条样插值函数

  • 输入:$x_1,…,x_{n+1};f_1=f(x_1),…,f_{n+1}(x_{n+1});f_{p_1}=f’(x_1),f_{p_2}=f’(x_{n+1}).$
  • 输出:$S_i(x)=\sum\nolimits_{i=1}^{4} A_{k,i}(x-x_i)^{k-1} = 0$ 的系数 $A_{1,i},A_{2,i},A_{3,i},A_{4,i},i=1,…,n$

Matlab 实现:

function [ A ] = Cubic_Spline_Interpolation( x0, y0, fp1, fp2 )
% 完备的三次条样插值法
% input: (x0, y0): 插值点; fp1 = f'(x1); fp2 = f'(x_n+1);
% output: 详见上述描述
n = length(x0) - 1;
for i = 1 : n
    h(i) = x0(i + 1) - x(i);
    b(i) = (y0(i + 1) - y0(i)) / h(i);
end % for
d(1) = b(1) - fp1;
d(n + 1) = fp2 - b(n);
for i = 2 : n
    d(i) = b(i) - b(i - 1);
end % for

% ————————————————————
%用三对角算法求上述方程组的解
% ————————————————————

for i = 1 : n
    A(1, i) = y0(i);
    A(2, i) = b(i) - (h(i) / 6) * (m(i + 1) + 2 * m(i));
    A(3, i) = m(i) / 2;
    A(4, i) = (m(i + 1) - m(i)) / (6 * h(i));
end % for
end % function

结尾

插值法是函数逼近问题的一种解决方案,函数逼近更多的在于理论方面的说明,如最佳一致逼近、最佳平方逼近等等,这里我也没有尝试写一些比较通用的算法与 Code 供使用(Matlab 工具包中提供了非常多的函数逼近工具可供使用)。在后面了解深度学习的过程中,函数逼近问题肯定是无法被避免的。