C語言例題1-9 已知六點的座標值 使用牛頓內插法求差除表？

C語言例題1-9 已知六點的座標值 使用牛頓內插法求差除表？

已知6點座標
n=6
x        f(x)
=============
0.0      -6.0
0.1      -5.89483
0.3      -5.65014
0.6      -5.17788
1.0      -4.28172
1.1      -3.99583
=============     

源自於 https://www.geeksforgeeks.org/newtons-divided-difference-interpolation-formula/

Newton’s Divided Difference Interpolation Formula

Interpolation is an estimation of a value within two known values in a sequence of values.

Newton’s divided difference interpolation formula is a interpolation technique used when the interval difference is not same for all sequence of values.

Suppose f(x0), f(x1), f(x2)………f(xn) be the (n+1) values of the function y=f(x) corresponding to the arguments x=x0, x1, x2…xn, where interval differences are not same
Then the first divided difference is given by

The second divided difference is given by

and so on…
Divided differences are symmetric with respect to the arguments i.e independent of the order of arguments.
so, 
f[x0, x1]=f[x1, x0]
f[x0, x1, x2]=f[x2, x1, x0]=f[x1, x2, x0]
By using first divided difference, second divided difference as so on .A table is formed which is called the divided difference table.

Divided difference table:

NEWTON’S DIVIDED DIFFERENCE INTERPOLATION FORMULA

Examples:

Input : Value at 7
       
Output :
      
      Value at 7 is 13.47

程式

// CPP program for implementing 
// Newton divided difference formula 
#include <stdio.h>

  
// Function to find the product term 
double proterm(int i, float value, float x[]) 
{ 
    double pro = 1; 
    for (int j = 0; j < i; j++) { 
        pro = pro * (value - x[j]); 
    } 
    return pro; 
} 
  
// Function for calculating 
// divided difference table 
void dividedDiffTable(float x[], float y[][10], int n) 
{ 
    for (int i = 1; i < n; i++) { 
        for (int j = 0; j < n - i; j++) { 
            y[j][i] = (y[j][i - 1] - y[j + 1] 
                         [i - 1]) / (x[j] - x[i + j]); 
        } 
    } 
} 
  
// Function for applying Newton's 
// divided difference formula 
double applyFormula(float value, float x[], 
                   float y[][10], int n) 
{ 
    double sum = y[0][0]; 
  
    for (int i = 1; i < n; i++) { 
      sum = sum + (proterm(i, value, x) * y[0][i]); 
    } 
    return sum; 
} 
  
// Function for displaying  
// divided difference table 
void printDiffTable(float y[][10],int n) 
{ 
    for (int i = 0; i < n; i++) { 
        for (int j = 0; j < n - i; j++) { 
            printf("    %0.4lf\t", y[i][j]); 
        } 
        printf("\n"); 
    } 
} 
  
// Driver Function 
int main() 
{ 
    // number of inputs given 
    int n = 6; 
    float value, sum, y[10][10]; 
    float x[] = { 0.0 , 0.1 , 0.3 ,0.6 , 1.0  , 1.1 }; 
  
    // y[][] is used for divided difference 
    // table where y[][0] is used for input 
    // -6.0 ,  -5.89483   ,-5.65014  ,-5.17788 ,-4.28172  ,-3.99583 
    y[0][0] = -6.0; 
    y[1][0] = -5.89483; 
    y[2][0] = -5.65014; 
    y[3][0] = -5.17788;
    y[4][0] = -4.28172;
    y[5][0] = -3.99583;
    
  
    // calculating divided difference table 
    dividedDiffTable(x, y, n); 
  
    // displaying divided difference table 
    printDiffTable(y,n); 
  
    // value to be interpolated 
    value = 0.2; 
    double result = applyFormula(value, x, y, n);
    // printing the value 
    printf("\nValue at %.2lf  is %.5lf ",value , result); 
    return 0; 
}

輸出結果

   -6.0000     1.0517     0.5725     0.2150     0.0630     0.0141 
    -5.8948     1.2235     0.7015     0.2780     0.0786 
    -5.6501     1.5742     0.9517     0.3566 
    -5.1779     2.2404     1.2370 
    -4.2817     2.8589 
    -3.9958 

Value at 0.20  is -5.77860