範例1-9 已知6點座標,依照牛頓差除法值計算P(0.2)之值 (向前與向後差除法)

範例1-9 已知6點座標,依照牛頓差除法值計算P(0.2)之值

已知6點座標

n=5
0.0  -6.0
0.1  -5.89483
0.3  -5.65014
0.6  -5.17788
1.0  -4.28172
1.1  -3.99583
依照牛頓差除法值計算P(0.2)之值

# Python3 program for implementing
# Newton divided difference formula

# Function to find the product term
def proterm(i, value, x):
    pro = 1;
    for j in range(i):
        pro = pro *(value - x[j]);
    return pro;


# Function for calculating
# divided difference table
def dividedDiffTable(x, y, n):
    for i in range(1, n):
        for j in range(n -  i):
            y[j][i] = ((y[j][i - 1] -  y[j + 1][i -  1]) / (x[j] -  x[i + j]));
    return y;

# Function for applying Newton’s
# divided difference formula
def applyFormula(value, x, y, n):
    sum = y[0][0];
    for i in range(1, n):
        sum = sum + (proterm(i, value, x) * y[0][i]);
    return sum;

# Function for displaying divided
# difference table
def printDiffTable(y, n):
    for i in range(0,n):
        print(round(x[i], 4),"\t", end = ""); 
        for j in range(0,n -  i):
            print( round(y[i][j], 4),"\t", end = "");
        print("");


# 後差方式顯示
def ShowBackward(y , n):
        print("\n*********** Show By Backward **********");
        for i in range(0,n):
                print(round(x[i], 4),"\t", end = ""); 
                for j in range(0,i+1):
                    print(round(y[i-j][j],4) ,"\t", end = "");
                print("");
       
# Driver Code
# number of inputs given
#=======main subroutine==================
n = 5+1
# create zero array
#from numpy import zeros
y = [[0 for i in range(10)]  for j in range(10)];

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
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
y=dividedDiffTable(x, y, n);

# displaying divided difference table
printDiffTable(y, n);
ShowBackward(y , n)
# value to be interpolated
value = 0.2;

# printing the value
print("\nValue at", value,  "is",round(applyFormula(value, x, y, n), 6))

# This code is contributed by mits

輸出畫面
======= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/Ex1-9-6.py ===========
0.0 -6.0 1.0517 0.5725 0.215 0.063 0.0142
0.1 -5.8948 1.2234 0.7015 0.278 0.0786
0.3 -5.6501 1.5742 0.9517 0.3566
0.6 -5.1779 2.2404 1.237
1.0 -4.2817 2.8589
1.1 -3.9958

*********** Show By Backward **********
0.0 -6.0
0.1 -5.8948 1.0517
0.3 -5.6501 1.2234 0.5725
0.6 -5.1779 1.5742 0.7015 0.215
1.0 -4.2817 2.2404 0.9517 0.278 0.063
1.1 -3.9958 2.8589 1.237 0.3566 0.0786 0.0142

Value at 0.2 is -5.778599
>>>