範例1-13 已知不明分行物體, 經光測儀器之量測,印出Hermite 差除表

範例1-13 已知不明分行物體, 經光測儀器之量測其結果如下
=========================================
時間                位置               飛行速度
 x                       f(x)                  f ' (x)
 0                       50.00               50.00
 2                       216.67             94.44
 4                       410.00             98.00
=========================================
印出Hermite 差除表
x=1 時  Hn(1) = ?  H'n(1)=?
 

'''
To generate the coeffficients for
 * Hermite Interpolating Polynomial H on the distinct numbers
 * x0,x1,x2,...with f0,f1,f2,... and f'0,f'1,f'2...
 * based on the Newton backward divided-difference Algorithm
 * and output the divided_difference table for Hn(x).
'''
#x[30],z[30],q[30][30],f[30],ff[30];
x = [0 for i in range(30)]
z = [0 for i in range(30)]
f = [0 for i in range(30)]
q = [[0 for i in range(30)]  for j in range(30)];
ff = [[0 for i in range(30)]  for j in range(30)];
n=2+1
x=[0.0 , 2.0 ,4.0]
f=[50.00 ,  216.67 ,  410.00]
ff=[50.00 ,  94.44  ,  98.00]
print("i","\t", "z(i)","\t","f(i)","\t","f(i-1,i)","\t","f(i-2,i-1,i)","\t", "....................................");

for i in range (0,n):
    z[2*i]=x[i];
    z[2*i+1]=x[i];
    q[2*i][0]=f[i];
    q[2*i+1][0]=f[i];
    q[2*i+1][1]=ff[i];
    if(i !=0):
        q[2*i][1]=(q[2*i][0]-q[2*i-1][0])/(z[2*i]-z[2*i-1]);

    if(i==0):
        print(round(i, 4),"\t",round(z[i], 4), "\t" , round(q[i][0],4) );
    elif(i==1):
        print(round(i, 4),"\t",round(z[i], 4), "\t" , round(q[i][0],4) ,"\t",round(q[i][1],4) );

for i in range (2,2*n):
    print(round(i, 4),"\t",round(z[i], 4), "\t" , round(q[i][0],4) ,"\t",round(q[i][1],4),end="" );
    for j in range(2,i+1):
        q[i][j]=(q[i][j-1]-q[i-1][j-1])/(z[i]-z[i-j]);
        print("\t",round(q[i][j] ,4) , end="");
    print("");

輸出結果
========= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX1-13.py ===========
i z(i)    f(i)   f(i-1,i) f(i-2,i-1,i) ....................................
0 0.0    50.0
1 0.0    50.0   50.0
2 2.0    216.67   83.335 16.6675
3 2.0    216.67   94.44 5.5525 -5.5575
4 4.0    410.0   96.665 1.1125 -1.11 1.1119
5 4.0    410.0   98.0 0.6675 -0.2225 0.2219 -0.2225
>>> 

Hn(x)= 50.0 + 50.0 (x-0.0) + 16.6675 (x-0.0)^2 + -5.5575 (x-0.0)^2 (x-2.0) +
            1.1119(x-0.0)^2 (x-2.0)^2 - 0.2225(x-0.0)^2 (x-2.0)^2(x-4.0)

        =  -0.2225x^5 + 2,8919x^4 -14.455x^3 + 35.79x^2 + 50x  + 50

H'n(x)= 5* -0.2225x^4 +  4* 2,8919x^3 - 3*14.455x^2 + 2*35.79x + 50