Edison.X. Blog

雖然向前差除和向後差除在計算順序上似乎有些不同，但對程式設計而言，事實上前差法和後差法是一樣的東西，只是表格輸出不同而已。說明圖如下所示。

/*
        Newton Forward and Backward Interpolation Algorithm
        Assume function is -4x^4 + 2x^3-3x^2+4x-1

        In fact, forward and backward is the same for Interpolation.
        If we use the 2 dim array to store the divided-difference value,
        we can find that the arrays of forward and backward are the same,
        just different display for the two arrays.

        事實上, 前除與後除是一樣的東西。若使用二維陣列儲存差分值，
        觀查其結果會發現，前除與後除之凍列是一樣的，只是顯示方式不同而已。

                                                                   -by EdisonX
*/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define N 4

// 前差方式顯示
void ShowForward(double **f, unsigned n)
{
        unsigned i, j;
        printf("\n **** Show By Forward **** \n");
        for(i=0; i<n; i++){
                for(j=0; j<n-i; j++) printf("%5.0lf ", f[i][j]);
                printf("\n");
        }
}

// 後差方式顯示
void ShowBackward(double **f, unsigned n)
{
        unsigned i, j;
        printf("\n **** Show By Backward **** \n");
        for(i=0; i<n; i++){
                for(j=0; j<=i; j++) printf("%5.0lf ", f[i-j][j]);
                printf("\n");
        }
}

double func(double x)
{
        double sum=0.0;
        double x2=x*x, x3=x2*x,x4=x3*x;
        sum = -4*x4+2*x3-3*x2+4*x-1;
        return sum;
}

double NewtonInterpolation(double x,double *xn, double *yn, unsigned n)
{
        double mult=0.0;
        double **f; // 差分矩陣
        double P; // result
        unsigned i, j;

        // 配置記憶體
        f = (double**)malloc(sizeof(double*)*n);
        for(i=0; i<n; i++) f[i]=(double*)malloc(sizeof(double)*n);
        for(j=0; j<n; j++) f[j][0] = yn[j];

        // 建立差分表
        for(j=1; j<n; j++)
                for(i=0; i<n-j; i++)
                        f[i][j]= (f[i+1][j-1]-f[i][j-1])/(xn[i+j]-xn[i]);

        // 代入公式
        for(i=0, P=0.0; i<n; i++){
                for(j=0, mult=f[0][i]; j<i && mult!=0.0; j++){
                        mult*=(x-xn[j]);
                }
                P+=mult;
        }

        // 顯示前差與後差表格
        ShowForward(f, n), ShowBackward(f, n), printf("\n");

        // 釋放記憶體
        for(i=0; i<n; i++) free(f[i]);
        free(f);
        return P;
}

int main()
{
        unsigned i;
        double x=1.5;

        double xn[N+1]={1.234, 3.782, 4.563, 5.987, -1.157};
        double yn[N+1];

        for(i=0; i<N+1; i++) yn[i]=func(xn[i]); // reference point
        printf("correct: %lf\n", func(x));
        printf("Newton: %lf\n",NewtonInterpolation(x, xn, yn, N+1));
        return 0;
}

執行結果

correct: -15.250000

 **** Show By Forward ****
   -6  -288  -241   -60    -4
 -739 -1089  -527   -51
-1589 -2251  -277
-4795  -668
  -20

 **** Show By Backward ****
   -6
 -739  -288
-1589 -1089  -241
-4795 -2251  -527   -60
  -20  -668  -277   -51    -4

Newton: -15.250000