例題3-3 不等距的函數f(x)的微分近似值

例題3-3 不等距的函數f(x)的微分近似值
先使用牛頓向前的內插多項式

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝====
x        f(x)          第一項P'n(x)      前二項P'n(x)　　前三項P'n(x)
0.5     0.4794 
0.6     0.5646 
0.8     0.7174 
1.05   0.8674     
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝====

# 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 = [ 0.5 , 0.6 , 0.8 , 1.05 ];
# y[][] is used for divided difference
# table where y[][0] is used for input
y[0][0] = 0.4794
y[1][0] = 0.5646
y[2][0] = 0.7174
y[3][0] = 0.8674;

# 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.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/Ex3-3-1.py ===========

0.5 0.4794 0.852 -0.2933 -0.1293

0.6 0.5646 0.764 -0.3644

0.8 0.7174 0.6

1.05 0.8674

*********** Show By Backward **********

0.5 0.4794

0.6 0.5646 0.852

0.8 0.7174 0.764   -0.2933

1.05 0.8674 0.6 　-0.3644 -0.1293

>>> 

P'n(x)= 0.852  +  -0.2933 * [(x-0.6) + (x- 0.5) ]  +  -0.1293  *[ (x-0.8) + (x-0.6) + (x- 0.5)


x=0.5

第一項P'n(x)      　　
P'n(0.5)=0.852

前二項P'n(x)
P'n(0.5)=0.852 + -0.2933 * [ (0.5-0.6) + (0.5-0.5) = 0.8813

前三項P'n(x)
P'n(0.5)= 0.852  +  -0.2933 * [(0.5-0.6) + (0.5-  0.5) ]  +  -0.1293  *[ (0.5-0.8) + ( 0.5 - 0.6) + ( 0.5 - 0.5) = 0.87745

x=0.6
.第一項P'n(x)      　　
P'n(0.6)=0.852

前二項P'n(x)
P'n(0.6)=0.852 + -0.2933 * [ (0.6-0.6) + (0.6-0.5) = 0.82267

前三項P'n(x)
P'n(0.6)= 0.852  +  -0.2933 * [(0.6-0.6) + (0.6 - 0.5) ]  +  -0.1293  *[ (0.6-0.8) + ( 0.6 - 0.6) + ( 0.6 - 0.5) = 0.82526

x=0.8
.第一項P'n(x)      　　
P'n(0.8)=0.852

前二項P'n(x)
P'n(0.8)=0.852 + -0.2933 * [ (0.8-0.6) + (0.8-0.5) = 0.70534

前三項P'n(x)
P'n(0.8)= 0.852  +  -0.2933 * [(0.8-0.6) + (0.8 - 0.5) ]  +  -0.1293  *[ (0.8-0.8) + ( 0.8 - 0.6) + ( 0.8 - 0.5) = 0.69758

................................................................
................................................................