數值分析實驗報告 Lab2 埃爾米特（Hermite）插值

一、問題引出

1. 掌握埃爾米特插值演算法原理；
2. 使用C語言程式設計實現埃爾米特插值演算法。

二、實驗準備

• 閱讀《數值分析》——李慶陽 2.4節

三、實驗要求

問題： 某人從甲地開車去乙地，每隔一段時間對行車距離和速率進行一次取樣，得到在 $n+1$

函式介面定義：

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$ 中。

裁判程式如下：

#include<stdio.h>

#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; /*確定給定的N=n+1組資料*/
  double t[MAXN], s[MAXN], v[MAXN]; /* 用於構造的資料 */
  double ht[MAXM], hs[MAXM], hv[MAXM]; /* 用估算的資料 */

  int m; //我們要估算的資料組數
  int i; //類似指向

 //可以說是主要部分了
  while ( scanf("%d", &N) != EOF ) {
  
 //輸入資料，可以說用於插值-構造Hermite插值多項式了
    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);
    for ( i=0; i<m; i++ ){scanf("%lf", &ht[i]);}
      
  //呼叫函式
    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");/*空行鴨*/
    
  }
  return 0;
}

裁判輸入資料：

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

標準輸出資料：

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

四、程式碼總結

//
//  main.cpp
//  Hermite
//
//  Created by apple on 2018/10/16.
//  Copyright © 2018年 apple. All rights reserved.
//

#include<stdio.h>

#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[] )
{
    double l[10],p[10],h1[10],h2[10],x,ll[10],pp[10];
    int kk;
    for(kk=0;kk<m;kk++)
    {
        x=ht[kk];hs[kk]=0;hv[kk]=0;
        int i;
        for(i=0;i<N;i++)
        {
            l[i]=1;ll[i]=1;
            int j;
            for(j=0;j<N;j++)
            {
                if(i!=j)
                {
                    l[i]=l[i]*(x-t[j])/(t[i]-t[j]);
                }
            }
            p[i]=0;pp[i]=0;
            int k;
            for(k=0;k<N;k++)
            {
                if(i!=k)
                {
                    p[i]=p[i]+l[i]/(x-t[k]);
                    pp[i]=pp[i]+ll[i]/(t[i]-t[k]);
                }
            }
            h1[i]=(1-2*pp[i]*(x-t[i]))*l[i]*l[i];
            h2[i]=(x-t[i])*l[i]*l[i];
            hs[kk]=hs[kk]+s[i]*h1[i]+v[i]*h2[i];
            int kkk;
            for(kkk=0;kkk<N;kkk++)
            {
                if(x==t[kkk])
                    break;
                
            }
            if(x==t[kkk])
                hv[kk]=v[kkk];
            else
                hv[kk]=hv[kk]+s[i]*(2*p[i]*l[i]-4*l[i]*p[i]*(x-t[i])*pp[i]-2*pp[i]*l[i]*l[i])+v[i]*(l[i]*l[i]+2*l[i]*p[i]*(x-t[i]));
            
        }
    }
}
int main()
{
  int N; /*確定給定的N=n+1組資料*/
  double t[MAXN], s[MAXN], v[MAXN]; /* 用於構造的資料 */
  double ht[MAXM], hs[MAXM], hv[MAXM]; /* 用估算的資料 */

  int m; //我們要估算的資料組數
  int i; //類似指向

 //可以說是主要部分了
  while ( scanf("%d", &N) != EOF ) {
  
 //輸入資料，可以說用於插值-構造Hermite插值多項式了
    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);
    //時間
    for ( i=0; i<m; i++ ){scanf("%lf", &ht[i]);}
      
  //呼叫函式
    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");/*空行鴨*/
    
  }
  return 0;
}

本blog通過蒐集資料，暫時只給出程式答案。
具體過程在本年度10.30會給出。

@author Richard_vim
@E-mail: [email protected]
@date 10.16