数值分析实验：Hermite插值

1 题目

题目::人从甲地开车去乙地，每隔一段时间对行车距离和速率进行一次采样，得到在 $n+1$ 个采样时刻点 $t_i$ 的里程 $s_i$ 和速率 $v_i（i=0,1,...,n）$。

要求::编程构造埃尔米特插值多项式 $H_{2n+1}(t)$

满足::
对所有i=0, 1, …, n

$\{\begin{array}{ccc}{H}_{2n+1}(t_i)=s_i& & \\ H_{2n+1}^{\prime }(t_i)=v_i& & \end{array}$

并据此计算$m$个给定时刻$ht$及其里程$hs$和速率$hv$。

函数接口定义::
void Hermite_Interpolation( int N, double t[], double s[], double v[], int m, double ht[], double hs[], double hv[] );

其中$N$为采样点个数（注意这个$N$不是公式中的最大下标$n$，而是等于$n+1$）

采样时刻点$t_i$、里程$s_i$、速率$v_i$分别通过$t、s、v$传入；

$m$是需要估算的给定时刻的个数，$ht$传入给定的时刻点，相应计算出的里程和速率应分别存储在$hs$和$hv$中。

---

2 c++代码实现

2.0 程序目录

伪代码

main.cpp

Hermite_Interpolation.h

input.txt

output.txt

2.1 伪代码

===========================================================================================================

伪代码：Hermite_Interpolation

输入：$N$ 个采样时刻点 $t_i$ 的里程 $s_i$ 和速率 $v_i（i=0,1,...,n）$, $m$个给定时刻

输出：这$m$个给定时刻$ht$及其里程$hs$和速率$hv$。

===========================================================================================================
00：init $hs[]$, $hv[]$, $l[]$, $p[]$, $pp[]$, $k$

01：for 　$k\leftarrow 0$　 to 　 $m$

02：　　$xx\leftarrow ht[k]$

03：　　for 　$i\leftarrow 0$　 to 　$N$

04：　　　　$l_i\leftarrow {\prod }_{j=0,j̸=i}^N\frac{x-t_j}{t_i-t_j}$

05：　　　　$p_i\leftarrow {\sum }_{j=0,j̸=i}^N\frac{l_i}{x-t_j};$

06：　　　　$p{p}_i\leftarrow {\sum }_{j=0,j̸=i}^N\frac{1}{t_i-t_j};$

07：　　　　$h{1}_i\leftarrow [1-2\cdot p{p}_i\cdot (x-t_i)]\cdot l_i^2;$

08：　　　　$h{2}_i\leftarrow (x-t_i)\cdot l_i^2;$

09：　　end

10：　　$h{s}_k\leftarrow {\sum }_{i=0}^N(s_i\cdot h{1}_i+v_i\cdot h{2}_i);$

11：　　if 　find　$a$　s.t.　$h{t}_k==t_a$ 　then

12：　　　　$h{v}_k\leftarrow v_a$

13：　　else

14：　　　　$h{v}_k\leftarrow {\sum }_{i=0}^N{s}_i\cdot [2\cdot l_i\cdot p_i-2\cdot p{p}_i\cdot l_i-4\cdot p{p}_i\cdot (x-t_i)\cdot l_i\cdot p_i]+v_i\cdot [l_i^2+2\cdot l_i\cdot p_i\cdot (x-t_i)]$

15：　　end if

16：end

17：return　 $hs[],hv[]$

===========================================================================================================

2.2 main.cpp

#include<stdio.h>
#include<iostream>
using namespace std;

#define MAXN 5 /* 最大采样点个数 */
#define MAXM 10 /* 最大估算点个数 */

void Hermite_Interpolation( int N, double t[], double s[], double v[], int m, double ht[], double hs[], double hv[] );

int main()
{
    int N, m;
    double t[MAXN], s[MAXN], v[MAXN]; /* 用于构造的数据 */
    double ht[MAXM], hs[MAXM], hv[MAXM]; /* 用于估算的数据 */
    int i;

    freopen("input.txt", "r", stdin);

    while ( scanf("%d", &N) != EOF ) {
            //可以测试多组数据
            //输入N:n+1个采样点

        for ( i=0; i<N; i++ )
            scanf("%lf", &t[i]);
        for ( i=0; i<N; i++ )
            scanf("%lf", &s[i]);
        for ( i=0; i<N; i++ )
            scanf("%lf", &v[i]);

        scanf("%d", &m);
            //m个需要估算的给定时刻

        for ( i=0; i<m; i++ )
            scanf("%lf", &ht[i]);

        //重点::利用Hermite插值多项式进行估算
        Hermite_Interpolation(N, t, s, v, m, ht, hs, hv );

        //输出对应时刻里程
        for ( i=0; i<m; i++ )
            printf("%.4lf ", hs[i]);

        printf("\n");

        //输出对应时刻速率
        for ( i=0; i<m; i++ )
            printf("%.4lf ", hv[i]);

        printf("\n\n");
    }
    fclose(stdin);

  return 0;
}

#include "Hermite_Interpolation.h"

2.3 Hermite_Interpolation.h

#ifndef HERMITE_INTERPOLATION_H_INCLUDED
#define HERMITE_INTERPOLATION_H_INCLUDED
#include<cstdio>

void Hermite_Interpolation( int N, double t[], double s[], double v[], int m, double ht[], double hs[], double hv[] )
{
    double l[MAXM], f[MAXM], h1[MAXM], h2[MAXM], t_x, a[10];

    //遍历每一组要求时刻
    for(int k=0;k<m;k++)
    {
        //初始化
        t_x=ht[k];
        hs[k]=0;
        hv[k]=0;

        //遍历每个采样点
        for(int i=0;i<N;i++)
        {
            //初始化
            l[i] = 1;
            a[i] = 0;

            //遍历每个采样点,求l[i],用于求hs
            for(int j=0;j<N;j++)
            {
                if(i != j)
                {
                    l[i] *= (t_x-t[j])/(t[i]-t[j]);
                    a[i] += 1/(t[i] - t[j]);
                }
            }

            //初始化
            f[i]=0;

            //遍历每个采样点,求f[i],用于求hv
            for(int j=0;j<N;j++)
            {
                if(i!=j)
                    f[i] += l[i]/(t_x-t[j]);
            }

            h1[i] = (1-2*a[i]*(t_x-t[i]))*l[i]*l[i];
            h2[i] = (t_x-t[i])*l[i]*l[i];
            hs[k] = hs[k] + s[i]*h1[i] + v[i]*h2[i];

            int e=0;
            for(;e<N;e++)
            {
                if(t_x==t[e])
                break;
            }

            if(t_x==t[e])
                hv[k]=v[e];
            else
                hv[k]=hv[k] + (2*l[i]*f[i] - 2*a[i]*l[i]*l[i] - 4*a[i]*(t_x-t[i])*l[i]*f[i]) * s[i] + (l[i]*l[i] + 2*l[i]*f[i]*(t_x-t[i])) * v[i];
        }
    }
}


#endif // HERMITE_INTERPOLATION_H_INCLUDED

2.4 input.txt

2
0.0 1.0
0.0 1.0
0.0 0.0
5
0.0 0.2 0.5 0.8 1.0
 
3
0.0 0.5 1.0
100.0 170.0 200.0
30.0 150.0 0.0
5
0.0 0.25 0.5 0.75 1.0
 
5
0.0 1.0 2.0 3.0 4.0
0.0 60.0 160.0 260.0 300.0
5.0 70.0 100.0 120.0 20.0  
10
0.5 1.0 1.5 2.0 2.5 3.0 3.5 3.8 3.95 4.0

2.5 output.txt

0.0000 0.1040 0.5000 0.8960 1.0000 
0.0000 0.9600 1.5000 0.9600 0.0000 
 
100.0000 127.9297 170.0000 195.9766 200.0000 
30.0000 165.4688 150.0000 52.9688 0.0000 
 
30.2222 60.0000 105.9303 160.0000 206.3438 260.0000 307.9764 305.7687 299.9796 300.0000 
62.6024 70.0000 109.0488 100.0000 92.9745 120.0000 41.2374 -44.8421 -16.2783 20.0000 

---

3 难点解析 – Hermite_Interpolation.h

Hermite插值满足 在节点上等于给定函数值 ，而且 在节点上的导数值也等于给定的导数值 。对于高阶导数的情况，Hermite插值多项式比较复杂，在实际情况中，常常遇到的是 函数值与一阶导数给定的情况 。

在此情况下，对给定的n个节点x1,x2,…,xn的Hermite插值多项式充分理解 二重hermite插值多项式的一般形式 , 并在此基础上,编程实现

尽可能 节约空间 ,并 降低算法复杂 的条件下实现hermite插值多项式( 公式见伪代码 ),最终利用求出的插值多项式求出在$m$个$ht[]$时刻下的$hs[]$与$hv[]$

---

---

Author

lance
2018.10.16