範例1-10 , 1-11 牛頓向後與向前差除法

範例1-10 , 1-11 牛頓向後與向前差除法
n=3
xa=1.5
1.0  0.0
2.0  0.693
3.0  1.099
4.0  1.386

求P(1.5)之值

# 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 = 3+1
# create zero array
#from numpy import zeros
y = [[0 for i in range(10)]  for j in range(10)];

x = [ 1.0 , 2.0 , 3.0 , 4.0 ];
# y[][] is used for divided difference
# table where y[][0] is used for input
y[0][0] = 0.000;
y[1][0] = 0.693;
y[2][0] = 1.099;
y[3][0] = 1.386;

# calculating divided difference table
y=dividedDiffTable(x, y, n);

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

# 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-10.py ==========
1.0 0.0         0.693 -0.1435 0.028
2.0 0.693 0.406 -0.0595
3.0 1.099 0.287
4.0 1.386

*********** Show By Backward **********
1.0 0.0
2.0 0.693 0.693
3.0 1.099 0.406 -0.1435
4.0 1.386 0.287 -0.0595 0.028

Value at 1.5 is 0.392875
>>> 

向前差除法:

P(x)=0.0 + 0.693(x-1.0) -0.1453 (x-1.0)(x-2.0) + 0.028(x-1.0)(x-2.0)(x-3.0)

P(1.5)=0.392875


向後差除法:

P(x)=1.386+ 0.287(x-4.0) -0.0595 (x-4.0)(x-3.0) + 0.028(x-4.0)(x-3.0) (x-2.0)

P(1.5)=0.392875