插值方法

合集 - 数学建模算法(27)

1.整数规划2024-07-072.非线性规划2024-10-113.非线性规划之飞行管理问题2024-10-114.遗传算法与直接搜索2024-10-115.图论基础：最小路径问题2024-07-076.网络最大流量问题2024-07-077.最小费用最大流问题2024-07-058.旅行商（TSP）问题2024-07-059.微分方程模型2024-07-0910.数据预处理2024-10-11

11.插值方法2024-07-06

12.拟合算法2024-07-0813.数据初始化2024-10-1114.层次分析法2024-07-1315.模糊综合评价2024-07-1416.CRITIC指标客观赋权方法2024-10-1117.Topsis评价法2024-10-1118.熵权法2024-10-1119.灰色关联度分析2024-10-1120.聚类分析2024-10-1121.主成分分析2024-10-1122.因子分析2024-10-1123.判别分析2024-10-1124.典型相关分析2024-10-1125.对应分析2024-10-1126.多维标度法2024-10-1127.亮度变化与空间滤波2024-10-11

收起

插值是什么

在工程中，我们经常要用一条曲线将一些点依次连接起来，称为插值。

插值的可行性

插值法定理:对n+1个不同的节点有唯一多项式$\varphi (x)=a_0+a_{1}x+\cdots +a_n{x}^n$,使得${\varphi }_n(x_j)=y_j(j=0,1,2,\cdots ,n)$

证明：将$x_0$到$x_n$带入多项式能得到一个线性方程组，AX=Y.其中A中元素为$x_i$，X中为$a_i$，Y中为$y_i$.

A的形式是一个范德蒙行列式，一定可逆，所以线性方程有解，多项式一定存在。

插值函数

已知n+1个点$(x_i,y_i)(i=0,1,\cdots ,n)$下面讲各种插值函数。

分段线性插值函数

简单说就是将相邻的点两两用直线连接起来，如此形成一条折线，折线的方程为$I_n(x)=\sum _{i=0}^n{y}_i{l}_i(x)$，满足$I_n(x_i)=y_i$其中

$$l_i(x)=\{\begin{array}{c}\frac{x-x_{i-1}}{x_i-x_{i-1}},x\in [x_{i-1},x_i],i\ne 0\\ \frac{x-x_{i+1}}{x_i-x_{i+1}},x\in [x_i,x_{i+1}],i\ne n\\ 0,\mathrm{其}\mathrm{他}\end{array}$$

相当于是在求每个点的y值时将前面的点累加起来。

拉格朗日插值多项式

$L_n(x)=\sum _{i=0}^n{y}_i{l}_i(x)=\sum _{i=0}^n{y}_i(\prod _{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j})$
对于每一个$x_i$,$l_i(x_i)\ne 0,l_j(x_i)=0(j\ne i)$.
这样求$I_n(x_i)=y_i$.

function y=lagrange(x0,y0,x)
%x0,y0为节点
%x是插值点
n = length(x0);m = length(x);
 for i = 1:m
   z = x(i);
   s = 0.0;
   for k = 1:n
       p = 1.0;
        for j = 1:n
            if j ~= k
               p = p * (z - x0(j)) / (x0(k) - x0(j));
            end
        end
        s = p*y0(k)+s;
   end
   y(i)=s;
  end

牛顿插值

拉格朗日插值法中存在的“增加一个节点时整个计算工作重新开始”的问题。
牛顿插值法是拉格朗日插值法的改进方法，运用迭代的思路求函数，且可以节省乘除法运算次数, 同时在Newton插值多项式中用到差分与差商等概念，又与数值计算的其他方面有密切的关系。

function [A,y]= newtonzi(X,Y,x)
%   Newton插值函数
%   X为已知数据点的x坐标
%   Y为已知数据点的y坐标
%   x为插值点的x坐标
%   函数返回A差商表
%   y为各插值点函数值
n=length(X); m=length(x);
for t=1:m
    z=x(t); A=zeros(n,n);A(:,1)=Y';
    s=0.0; y=0.0; c1=1.0;
    for  j=2:n
       for i=j:n
           A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));
       end
    end
    C=A(n,n);
    for k=1:n
        p=1.0;
        for j=1:k-1
            p=p*(z-X(j));
        end
        s=s+A(k,k)*p;        
    end
    ss(t)=s;
end
    y=ss;
    A=[X',A];    
end

Hermite插值法

牛顿和拉格朗日插值法都有一些缺陷，它们只考虑了给出的点的拟合特性，没有考虑原函数的其他性质如导数与原函数是否适配。
有时候我们需要插值函数与原函数不仅需要一阶导数相同，还需要更高阶。
设原函数的导数为m(i),插值函数为H_i. 则Hermite函数满足

$$\{\begin{array}{c}H(x_i)=y_i\\ H^{\prime }(x_i)=m_i\end{array}$$

具体算法见：https://blog.csdn.net/SanyHo/article/details/106849323
下面给出代码.

function f = Hermite( x,y,y_1,x0 )
%Hermite插值函数
%   x为已知数据点的x坐标
%   y为已知数据点的y坐标
%   y_1为数据点y值导数
%   x0为插值点的x坐标
syms t;
f = 0.0;
if(length(x) == length(y))
    if(length(x) == length(y_1))
        n = length(x);
    else
        disp('y和y的导数维数不相等');
        renturn;
    end
else
    disp('x和y的维数不相等');
    return;
end
%以上为输入判断和确定“n”的值
for i = 1:n
    h =  1.0;
    a =  0.0;
    for j = 1:n
        if(j ~= i)
            h = h*(t-x(j))^2/((x(i)-x(j))^2);%求得值为(li(x))^2
            a = a + 1/(x(i)-x(j));   %求得ai（x）表达式之中的累加部分
        end
    end
    f = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));
    if(i == n)
        if(nargin == 4)
            f = subs(f, 't' , x0);  %输出结果
        else
            f = vpa(f,6);  %输出精度为有效数字为6位的函数表达式
        end
    end
end

样条插值

给定区间[a,b]的一个划分
$\mathrm{\Delta }:a=x_0<x_1<\cdots <x_{n-1}<x_n=b$
如果S(x)满足：1. 在每一个小区间上$[x_i,x_{i+1}](i=0,1,\cdots ,n-1)$，S(x)为m次多项式
2. 它有m-1阶连续导数。
则S(x)为m次样条函数。
计算比较繁琐，就不详细说了，感兴趣看

% 定义三次样条函数的系数矩阵、插值宽度、差商表、g值和M值
function [D,h,A,g,M]=threesimple(X,Y)
% X: 已知的插值节点数组
% Y: 插值节点对应的函数值数组

    n = length(X); % 计算插值节点的数量
    A = zeros(n,n); % 初始化差商表为n*n的零矩阵
    A(:,1) = Y'; % 第一列填充插值点的函数值
    D = zeros(n-2,n-2); % 初始化系数矩阵D为(n-2)*(n-2)的零矩阵
    g = zeros(n-2,1); % 初始化g值为(n-2)*1的零向量

    % 构建差商表A
    for j = 2:n
        for i = j:n
            A(i,j) = (A(i,j-1) - A(i-1,j-1)) / (X(i) - X(i-j+1));
        end
    end

    % 计算插值宽度h
    for i = 1:n-1
        h(i) = X(i+1) - X(i);
    end

    % 构建系数矩阵D和初始化g值
    for i = 1:n-2
        D(i,i) = 2;
        g(i) = (6 / (h(i+1) + h(i))) * (A(i+2,2) - A(i+1,2));
    end

    % 填充D矩阵的非对角线元素
    for i = 2:n-2
        u(i) = h(i) / (h(i) + h(i+1));
        n(i-1) = h(i) / (h(i-1) + h(i));
        D(i-1,i) = n(i-1);
        D(i,i-1) = u(i);
    end

    % 解线性方程组求解M值
    M = D \ g;
    M = [0; M; 0]; % 边界条件，M的首尾元素设为0

end

% 定义三次样条插值函数
function s = threesimple1(X,Y,x)
% X: 已知的插值节点数组
% Y: 插值节点对应的函数值数组
% x: 需要计算插值值的点

    [D,h,A,g,M] = threesimple(X,Y); % 调用threesimple函数计算相关参数
    n = length(X); % 插值节点的数量
    m = length(x); % 需要计算插值值的点的数量

    % 计算每个插值点的函数值
    for t = 1:m
        for i = 1:n-1
            if (x(t) <= X(i+1)) && (x(t) >= X(i))
                % 计算三次样条插值公式的各个部分
                p1 = M(i,1) * (X(i+1) - x(t)) ** 3 / (6 * h(i));
                p2 = M(i+1,1) * (x(t) - X(i)) ** 3 / (6 * h(i));
                p3 = (A(i,1) - M(i,1) / 6 * (h(i) ** 2)) * (X(i+1) - x(t)) / h(i);
                p4 = (A(i+1,1) - M(i+1,1) / 6 * (h(i) ** 2)) * (x(t) - X(i)) / h(i);
                s(t) = p1 + p2 + p3 + p4; % 计算插值点的函数值
                break; % 找到对应的区间后跳出循环
            else
                s(t) = 0; % 如果x不在X的范围内，插值结果为0
            end
        end
    end
end

matlab 工具包（香）

三次Hermite插值法

p = pchip(x,y,new_x)

三次样条插值

p = spline(x,y,new_x)

n维插值
2维为例，函数为二元函数

p = interp2(x,y,z,new_x,new_y,method)
%method具体有 'linear','cublic','spline','nearest'