習題1-7 已知下面4點座標求Hermite內插法的插除表並 計算Hn(x) = [1.5 , 2.5 , 3.5 , 4.5]之值

習題1-7 已知下面4點座標

x      f(x)       f'(x)=ff(x)
=================
1.0  0.000     1.0000
2.0  0.693     0.5000
3.0  1.099     0.3333

4.0  1.386     0.2500
==================
若已知 f(x) = ln(x)
求Hermite內插法的插除表並
計算Hn(x) = [1.5 , 2.5 , 3.5 , 4.5]之值


'''
/* pr1-7.py 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 H(xa) for different values of xa.
 */
 input  data
n=3  m=4
1.5  2.5  3.5  4.5

1.0  0.000  1.0000
2.0  0.693  0.5000
3.0  1.099  0.3333
4.0  1.386  0.2500

i  z(i)   f(i)   f(i-1,i)   f(i-2,i-1,i)...
0  1.0  0.0000
1  1.0  0.0000  1.0000
2  2.0  0.6930  0.6930 -0.3070
3  2.0  0.6930  0.5000 -0.1930  0.1140
4  3.0  1.0990  0.4060 -0.0940  0.0495 -0.0323
5  3.0  1.0990  0.3333 -0.0727  0.0213 -0.0141  0.0091
6  4.0  1.3860  0.2870 -0.0463  0.0132 -0.0040  0.0034 -0.0019
7  4.0  1.3860  0.2500 -0.0370  0.0093 -0.0019  0.0011 -0.0008  0.0004
Hn(1.50)= 0.4057
Hn(2.50)= 0.9165
Hn(3.50)= 1.2529
Hn(4.50)= 1.5076

 '''

xa=[1.5 , 2.5 , 3.5 ,  4.5 ]
n=4    #4  point n=4
m=len(xa)
x=[1.0  , 2.0 , 3.0 ,4.0]
f=[0.0  , 0.693 , 1.099 , 1.386]
ff=[1.0  , 0.5 , 0.3333 , 0.25 ]

z= [0.0 for i in range(2*len(x)+1)]     # [n] 矩陣
q=  [ [0.0 for i in range(2*len(x)+1)]  for j in range(2*len(x)+1) ]       #[n][n],


print('xa=',end='')
for i  in range (len(xa)) :   
    print( round( xa[i],4), ' ,', end='')

print('\n') 
print('x', '\t' , 'f' , '\t' , 'ff' )
print('===================')
for i  in range (len(f)) :   
    print( round( x[i],4), "\t",round( f[i],4),"\t",round( ff[i],4) )

print('\n\n')
print("i\tz(i)\tf(i)\tf(i-1,i)\tf(i-2,i-1,i)...\n");
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])
    elif(i==0):
        print("{%d}\t{%3.1f}\t{%7.4f}\n" %(i,z[i],q[i][0]))

i=1
print("{%d}\t{%3.1f}\t{%7.4f}\t{%7.4f}\n" %(i,z[i],q[i][0],q[i][1]))
       
 
for i in range (2 , 2*n ):
    print("{%d}\t{%3.1f}\t{%7.4f}\t{%7.4f}" %(i,z[i],q[i][0],q[i][1]),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{%7.4f}" %(q[i][j]) ,end='')
   
    print("\n")
 
for i in range (len(xa)):
    k=0
    dx=1
    h=q[0][0]
    for j in range (1 , 2*n+1):
        if(j%2==0):
            k=k-1

        dx=dx*(xa[i]-x[k])+0.0
        h=h+q[j][j]*dx
        k=k+1

    print("Hn(%4.2lf)=%7.4f\n" %(xa[i],h))


======== RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/Pr1-7.py ============
xa=1.5  ,2.5  ,3.5  ,4.5  ,

x f ff
===================
1.0 0.0 1.0
2.0 0.693 0.5
3.0 1.099 0.3333
4.0 1.386 0.25



i z(i) f(i) f(i-1,i) f(i-2,i-1,i)...

{0} {1.0} { 0.0000}

{1} {1.0} { 0.0000} { 1.0000}

{2} {2.0} { 0.6930} { 0.6930} {-0.3070}

{3} {2.0} { 0.6930} { 0.5000} {-0.1930} { 0.1140}

{4} {3.0} { 1.0990} { 0.4060} {-0.0940} { 0.0495} {-0.0323}

{5} {3.0} { 1.0990} { 0.3333} {-0.0727} { 0.0213} {-0.0141} { 0.0091}

{6} {4.0} { 1.3860} { 0.2870} {-0.0463} { 0.0132} {-0.0040} { 0.0034} {-0.0019}

{7} {4.0} { 1.3860} { 0.2500} {-0.0370} { 0.0093} {-0.0019} { 0.0011} {-0.0008} { 0.0004}

Hn(1.50)= 0.4057

Hn(2.50)= 0.9165

Hn(3.50)= 1.2529

Hn(4.50)= 1.5076

>>>