插值（四）——Hermite插值（C++）

龙行泽雨

于 2024-02-14 14:26:35 发布

阅读量860

点赞数 13

分类专栏：计算方法文章标签： c++ 开发语言线性代数

本文链接：https://blog.csdn.net/m0_50793313/article/details/136113262

版权

计算方法专栏收录该内容

12 篇文章 1 订阅

订阅专栏

hermite插值

实际问题中，对所构造的插值多项式，不仅要求函数值重合，而且若干阶导数也重合，即： $p(x_i)=f(x_i),p'(x_i)=f'(x_i),\cdots,p^{(n)}(x_i)=f^{(n)}(x_i)$ 。
一般情况都是使用二阶hermite插值，即要求： $H(x_i)=f(x_i),H'(x_i)=f'(x_i)$ ，构造其多项式为：
$H(x)=\sum_{i=0}^n g_i(x)f(x_i)+\sum_{i=0}^n h_i(x)f'(x_i)$
其中
$g_i(x)=[1-2(x-x_i)\sum_{i \neq j}\frac{1}{x_i-x_j}]l_i^2(x) h_i(x)=(x-x_i)l^2_i(x)$
这里的 $l_i(x)$ 是Lagrange插值的基函数。

代码

这里以二阶hermite插值为例，给出代码。

//二阶hermite插值
#include <iostream>
#include <cmath>
#define DerivativeIntervalAccuracy 0.0001//求导区间精度

//求拉格朗日插值函数的基函数li(x)
double Lagrange_li(double *xi,double *yi,double x,int i,int n)
{

    double li = 1.0;
    for (int j = 0; j < n;j++)
    {
        if(i!=j)
        {
            li*=((x-xi[j])/(xi[i]-xi[j]));    
        }
    }
    return li;
}
double hermite(double *xi,double *yi,double *f_dx,double x,int n)
{
    double li;
    double hermite = 0.0f;
    for (int i = 0; i < n;i++)
    {
        li=Lagrange_li(xi, yi, x, i, n);
        double gi = 0.0f;
        for (int j = 0; j < n;j++)
        {
            if(i!=j)
            {
                gi += 1 / (xi[i] - xi[j]);
            }            
        }
        gi = (1 - 2 * (x - xi[i]) * gi) * li * li;
        double hi = 2 * (x - xi[i]) * li * li;
        hermite += gi * yi[i] + hi * f_dx[i];
    }
    return hermite;
}
//近似求导，传递函数进来
double f_derivative(double (*fx)(double),double x)
{
    return (fx(x + DerivativeIntervalAccuracy) - fx(x))/DerivativeIntervalAccuracy;
}

double fx(double x)
{
    //这里填入被插值函数
    //如果插值区间包括分母为0的情况需要手动调整
    return sqrt(x + 3);
}

int main()
{
    int n;//插值点数
    std::cout<<"请输入插值点数：";
    std::cin>>n;
    double error[3]={0.0f,0.0f,0.0f};//误差评价
    double *x = new double [n];
    double *y = new double [n];
    double *d_fx = new double[n];//导数
    double a=3, b=10;//定义插值区间
    //生成插值数据
    for (int i = 0;i<n;i++)
    {
        x[i] = a + (b-a)*i/(n-1);
        d_fx[i]=f_derivative(fx,x[i]);
        y[i] = fx(x[i]);
        //std::cout<<"["<<x[i]<<","<<y[i]<<"]:"<<d_fx[i]<<std::endl;
    }

    //std::cout << std::endl;
    // 验证误差，抽取区间内的100个点进行误差检测
    for (int i = 0; i < 100; i++)
    {
        double x_temp = a + (b - a) * (i+1) / 100;
        double y_temp = fx(x_temp);
        double y_temp2 =hermite(x,y,d_fx,x_temp,n);
        if (fabs(y_temp - y_temp2) > error[0])
        {
            error[0] = fabs(y_temp - y_temp2);
            // std::cout << y_temp2 << std::endl;
        }
        error[1] += fabs((y_temp - y_temp2) / y_temp);
        error[2] += ((y_temp - y_temp2)) * ((y_temp - y_temp2));
    }
    //err[0]得到的是在区间内最大的误差
    //err[1]得到的是平均最大相对误差
    //err[2]是均方根误差
    error[1] = error[1]/100;
    error[2] = sqrt(error[2]/100);
    for(int i = 0;i < 3;i++)
    {
        std::cout<<"误差"<<i+1<<"："<<error[i]<<std::endl;
    }
    delete[] x;
    delete[] y;
    delete[] d_fx;
    return 0;
}

结果分析

在这里插入图片描述

其精度比牛顿插值法略差，但是hermite最重要的性质在于它的n阶导数也可以重合，此处是二阶Hermite插值，所以它的导数也会重合，所以这里插入一段代码查看导数的重合情况：

for (int j = 0; j < n;j++)
{
    if(x[j]==x_temp)
    {
        std::cout<<"x="<<x[j]<<":"<<fabs(y_temp - y_temp2)<<" d_f(x)_error:"<<
        fabs((hermite(x,y,d_fx,x_temp+DerivativeIntervalAccuracy,n)-hermite(x,y,d_fx,x_temp,n))/DerivativeIntervalAccuracy)-d_fx[j]
        <<std::endl;
        break;
    }
}