【工程应用】两点三次Hermite(埃尔米特)、分段Hermite插值法

mjiansun

已于 2024-09-18 15:00:24 修改

阅读量10

点赞数

分类专栏：常用数学知识文章标签： matlab 开发语言

于 2024-09-18 10:28:42 首次发布

原文链接：https://blog.csdn.net/Math_Boy_W/article/details/137855845

版权

常用数学知识专栏收录该内容

91 篇文章 69 订阅

订阅专栏

一、两点三次Hermite插值

已知: 两个插值节点机器对应的导数值

其中

二、分段三次Hermite插值

三、北太天元 or Matlab 实现

两点三次Hermite插值

function [output]= Hermite_interp(x0,y0,dy0,x)
% 两点三次 Hermite 插值
% x0 : [x1 x2]
% y0 : [y1 y2]
% dy0 : [y1' y2']
% length(x)=1
%
%   Version:            1.0
%   last modified:      05/16/2023
    h = x0(2)-x0(1);
    h1 = x - x0(1);
    h2 = x - x0(2);
    X1 = (h2/h)^2;
    X2 = (h1/h)^2;
    alpha1 = (1+2*h1/h)*X1;
    alpha2 = (1+2*h2/(-h))*X2;
    beta1 = h1*X1;
    beta2 = h2*X2;
    output = alpha1*y0(1)+alpha2*y0(2)+beta1*dy0(1)+beta2*dy0(2);
end

将上述代码保存为 Hermite_interp.m 文件.

分段三次Hermite插值

function [out1] = piecewise_Hermite_interp(x0,y0,dy0,x)
    % 分段Hermite插值
    % x0 : [x1 x2 ...]
    % y0 : [y1 y2...]
    % dy0 : [y1' y2'...]
    %
    %   Version:            1.0
    %   last modified:      09/13/2023
    %   file need: Hermite_interp.m
    m = length(x);
    n = length(x0);
    y = zeros(1,m);
    for k =1:1:m
        for i =1:1:n-1
        	   if x0(i)<=x(k) && x(k)<= x0(i+1)
        	   	 y(k) = Hermite_interp(x0([i,i+1]),y0([i,i+1]),dy0([i,i+1]),x(k));
        	   end
        end
    end
    out1 = y;
end

将上述代码保存为 piecewise_Hermite_interp 文件.

这里命名比较长,不过也没什么关系,我比较喜欢把意思表示清楚,另外也有代码的自动补全,也没有怎么不方便。

四、简单实现一下

例1

% Hermite_v0 例子1
%   last modified:      04/17/2024
%   file need: Hermite_interp.m   piecewise_Hermite_interp.m
clc;clear all;
%样本点
x0 = [0 pi/2];%x0 = [0 pi];
y0 = sin(x0);
dy0 = cos(x0);

x = 0:pi/20:pi
y = zeros(1,length(x));
for i= 1:1:length(x)
    y(i) = Hermite_interp(x0,y0,dy0,x(i));
end

plot(x,y);
hold on 
    plot(x,sin(x),'r');
    legend('H(x)','sin(x)');
    title('H(x)与sin(x)的图像')
hold off

运行后得到上图，

上面这张图中设置的插值节点是在 0 , pi/2，所以在前半部分比较贴合。

例2

%% 分段Hermite_v0 例子2
clc;clear all;
Runge = @(x)1./(x.^2+1);  % 定义Runge函数
dRunge = @(x) -(2*x)./((x.^2+1).^2); % 定义导函数

x0 = linspace(-5,5,11);
y0 = Runge(x0);
dy0 = dRunge(x0);
x = linspace(-5,5,100);
yy = piecewise_Hermite_interp(x0,y0,dy0,x);
plot(x,yy);
hold on
plot(x,Runge(x),'r');
legend('分段Hermite','Runge函数')
hold off

下面这张图为分段三次Hermite插值在Runge函数上的表现，也算是分段上的一种优势。

五、matlab源码

function yq = pchip(x,y,xq)
%PCHIP  Piecewise Cubic Hermite Interpolating Polynomial.
%   YQ = PCHIP(X,Y,XQ) performs shape-preserving cubic Hermite
%   interpolation using the values Y at sample points X to find
%   interpolated values YQ at the query points XQ.
%       - X must be a vector.
%       - If Y is a vector, Y(j) is the value at X(j).
%       - If Y is a matrix or n-D array, Y(:,...,:,j) is the value at X(j).
%
%   PCHIP chooses slopes at X(j) such that YQ respects monotonicity and is
%   shape-preserving. Namely, on intervals where the Y data is monotonic,
%   so is YQ; at points where the Y data has a local extremum, so does YQ.
%
%   PP = PCHIP(X,Y) returns the piecewise polynomial form PP of the
%   interpolant. You can use PP as an input to PPVAL or UNMKPP.
%
%   Comparison of SPLINE, PCHIP, and MAKIMA:
%       - All three are a form of piecewise cubic Hermite interpolation,
%         but each function computes the slopes of YQ at X(j) differently.
%       - SPLINE chooses slopes at X(j) such that the second derivative of
%         YQ is continuous. Therefore, SPLINE is smoother and more accurate
%         if the Y data represents values of a smooth function.
%       - PCHIP has no overshoots and less oscillation than SPLINE.
%       - MAKIMA has less oscillation than SPLINE but may have overshoots.
%       - PCHIP and MAKIMA are less expensive than SPLINE to set up PP.
%       - All three are equally expensive to evaluate.
%       - SPLINE and MAKIMA generalize to n-D grids. See INTERPN.
%
%   Example: Compare SPLINE, PCHIP, and MAKIMA
%
%       x = [1 2 3 4 5 5.5 7 8 9 9.5 10];
%       y = [0 0 0 0.5 0.4 1.2 1.2 0.1 0 0.3 0.6];
%       xq = 0.75:0.05:10.25;
%       yqs = spline(x,y,xq);
%       yqp = pchip(x,y,xq);
%       yqm = makima(x,y,xq);
%
%       plot(x,y,'ko','LineWidth',2,'MarkerSize',10)
%       hold on
%       plot(xq,yqp,'LineWidth',4)
%       plot(xq,yqs,xq,yqm,'LineWidth',2)
%       legend('(x,y) data','pchip','spline','makima')
%
%   See also INTERP1, MAKIMA, SPLINE, PPVAL, MKPP, UNMKPP.

% References:
%   F. N. Fritsch and R. E. Carlson, "Monotone Piecewise Cubic
%   Interpolation", SIAM J. Numerical Analysis 17, 1980, 238-246.
%   David Kahaner, Cleve Moler and Stephen Nash, Numerical Methods
%   and Software, Prentice Hall, 1988.

%   Copyright 1984-2021 The MathWorks, Inc.

% Check and adjust input data
[x,y,sizey] = chckxy(x,y);

% Compute slopes
h = diff(x);
m = prod(sizey);
delta = diff(y,1,2)./repmat(h,m,1);
slopes = zeros(size(y));
for r = 1:m
    if isreal(delta)
        slopes(r,:) = pchipSlopes(h,delta(r,:));
    else
        realSlopes = pchipSlopes(h,real(delta(r,:)));
        imagSlopes = pchipSlopes(h,imag(delta(r,:)));
        slopes(r,:) = complex(realSlopes,imagSlopes);
    end
end

% Compute piecewise cubic Hermite interpolant for those values and slopes
yq = pwch(x,y,slopes,h,delta);
yq.dim = sizey;

if nargin == 3
    % Evaluate the piecewise cubic Hermite interpolant
    yq = ppval(yq,xq);
end

function d = pchipSlopes(h,del)
% Derivative values for PCHIP.
% d = pchipslopes(x,y,del) computes the first derivatives, d(k) = P'(x(k)).

%  Special case n = 2, use linear interpolation.
n = numel(h)+1;
if n == 2
    d = repmat(del(1),[1 n]);
    return
end

%  Slopes at interior points.
%  d(k) = weighted average of del(k-1) and del(k) when they have the same sign.
%  d(k) = 0 when del(k-1) and del(k) have opposites signs or either is zero.

k = find(sign(del(1:n-2)).*sign(del(2:n-1)) > 0);

w1 = 2*h(k+1)+h(k);
w2 = h(k+1)+2*h(k);
d = zeros(1,n);
d(k+1) = (w1+w2)./(w1./del(k) + w2./del(k+1));

%  Slopes at end points.
%  Set d(1) and d(n) via non-centered, shape-preserving three-point formulae.

d(1) = ((2*h(1)+h(2))*del(1) - h(1)*del(2))/(h(1)+h(2));
if sign(d(1)) ~= sign(del(1))
    d(1) = 0;
elseif (sign(del(1)) ~= sign(del(2))) && (abs(d(1)) > abs(3*del(1)))
    d(1) = 3*del(1);
end

d(n) = ((2*h(n-1)+h(n-2))*del(n-1) - h(n-1)*del(n-2))/(h(n-1)+h(n-2));
if sign(d(n)) ~= sign(del(n-1))
    d(n) = 0;
elseif (sign(del(n-1)) ~= sign(del(n-2))) && (abs(d(n)) > abs(3*del(n-1)))
    d(n) = 3*del(n-1);
end