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

Julia語言例題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   
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝====
P'(n) 第1項 f0,1

P'(n) 第2項 f0,1+f0,1,2 *( (x-x[1]) + (x-x[2]) )

P'(n) 第3項 f0,1+f0,1,2 *( (x-x[1]) + (x-x[2]) )
                    + f0,1,2,3 *  ( (x-x[1])*(x-x[2]) ) + ( (x-x[2])*(x-x[3]) ) + ( (x-x[3])*(x-x[1]) )

using Printf


x=[ 0.5 , 0.6 , 0.8 , 1.05]
f= [  [0.4794    ,0.0 , 0.0   ,0.0  ,0.0 ,0.0 ],
      [0.5646    ,0.0 , 0.0   ,0.0  ,0.0 ,0.0 ],
      [0.7174    ,0.0 , 0.0   ,0.0  ,0.0 ,0.0 ],
      [0.8674    ,0.0 , 0.0   ,0.0  ,0.0 ,0.0 ] ]

n=length(x)
println("  Divided Difference Table: ")
println("=============================")
for j=2:n
    for i=1:n-j+1
        f[i][j]=(f[i+1][j-1]-f[i][j-1])/(x[i+j-1]-x[i])
    end
end

print("i\tx(i)\t\tf(i)\t\tf(i,i+1)\tf(i,i+1.i+2),  ......................\n")
for i=1:n
    s=@sprintf("%d\t%8.5f",i,x[i])
    print(s)
    for j=1:n-i+1
        s=@sprintf("\t%8.5f",f[i][j])
        print(s)
    end
    println()

end


d=[0.0 for i=1:n-1]
for i=1:n-1
    d[i]=f[1][i+1]
end

s=@sprintf("\n牛頓一階微分前向差除表")
println(s,d)

xa=[ 0.5 , 0.6 , 0.8 , 1.05]
pn=[0.0 , 0.0 , 0.0 ,0.0]
s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
    pn[i]=d[1]
    println("========================================================")
    s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
    println(s)
end


s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
    xb=xa[i]
    dx=0.0
    for j=1:2
        #println(xb,"  ", x[j])
        dx=dx+(xb-x[j]) 
    end
    #println(dx)
    dx=d[2]*dx
    pn[i]=pn[i]+dx
    println("========================================================")
    s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
    println(s)
end


s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
    xb=xa[i]
    dx=0.0
    dx=(xb-x[1])*(xb-x[2])+ (xb-x[2])*(xb-x[3])+ (xb-x[3])*(xb-x[1])
    dx=d[3]*dx
    pn[i]=pn[i]+dx
    println("========================================================")
    s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
    println(s)
end


println("\n\n")
for i=1:length(xa)
    xb=xa[i]
    println("========================================================")
    s=@sprintf("xa[%1d]= %0.3f , 誤差 f'(x)-P'(n) = %0.5f",i,xa[i],abs(cos(xb)-pn[i]) )
    println(s)
end

輸出畫面

i f(i)  f(i,i+1) f(i,i+1.i+2),  ......................
1  0.50000  0.47940  0.85200 -0.29333 -0.12929
2  0.60000  0.56460  0.76400 -0.36444
3  0.80000  0.71740  0.60000
4  1.05000  0.86740

牛頓一階微分前向差除表[0.852, -0.293333, -0.129293]

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.500 , P'n( 0.500 ) = 0.85200
========================================================
xa[2]= 0.600 , P'n( 0.600 ) = 0.85200
========================================================
xa[3]= 0.800 , P'n( 0.800 ) = 0.85200
========================================================
xa[4]= 1.050 , P'n( 1.050 ) = 0.85200

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.500 , P'n( 0.500 ) = 0.88133
========================================================
xa[2]= 0.600 , P'n( 0.600 ) = 0.82267
========================================================
xa[3]= 0.800 , P'n( 0.800 ) = 0.70533
========================================================
xa[4]= 1.050 , P'n( 1.050 ) = 0.55867

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.500 , P'n( 0.500 ) = 0.87745
========================================================
xa[2]= 0.600 , P'n( 0.600 ) = 0.82525
========================================================
xa[3]= 0.800 , P'n( 0.800 ) = 0.69758
========================================================
xa[4]= 1.050 , P'n( 1.050 ) = 0.49434

========================================================
xa[1]= 0.500 , 誤差 f'(x)-P'(n) = 0.00013
========================================================
xa[2]= 0.600 , 誤差 f'(x)-P'(n) = 0.00008
========================================================
xa[3]= 0.800 , 誤差 f'(x)-P'(n) = 0.00087
========================================================
xa[4]= 1.050 , 誤差 f'(x)-P'(n) = 0.00323