Julia語言 範例1-11 牛頓向後與向前差除法 求P(1.5)之值

Julia語言 範例1-11 牛頓向後與向前差除法

xa=1.5

x       f(x)
===========
1.0    0.0
2.0    0.693
3.0    1.099
4.0    1.386
============
求P(1.5)之值

using Printf

function newtonform(x::Array{Float64,1}, d::Array{Float64,1}, xa::Float64)
    #
    # Evaluates the Newton form of the
    # interpolating polynomial, with abscisses
    # in x and divided differences in d at xa.
    #
    n = length(d)
    result = d[n]
    for i=n-1:-1:1
        result = result*(xa - x[i]) + d[i]
    end
    return result
end



function newton_err(x::Array{Float64,1}, da::Float64 , xa::Float64)
    result = da
    for i=n-1:-1:1
        result = result*(xa - x[i])
    end
    return result
end


x=[ 1.0 , 2.0 , 3.0 ,4.0]
f= [  [0.00     ,0.0 , 0.0   ,0.0  ,0.0 ,0.0 ],
      [0.693    ,0.0 , 0.0   ,0.0  ,0.0 ,0.0 ],
      [1.099    ,0.0 , 0.0   ,0.0  ,0.0 ,0.0 ],
      [1.386    ,0.0 , 0.0   ,0.0  ,0.0 ,0.0 ] ]
y=f

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]
for i=1:n
    d[i]=f[1][i]
end 

s=@sprintf("牛頓前向差除表")
println(s,d)

println()
println()
xa=1.5
p = newtonform(x,d,xa)
s=@sprintf("%0.5f",p)
print("\nNewton內插法的差除表")
println("P(",xa,")=",s)

da=d[n]
p = newton_err(x,da,xa)
s=@sprintf("%0.8f",abs(p))
print("\nNewton內插法的誤差")
println("P(",xa,")=", s)


#================================================
 for(i=1;i<=n;i++)
   {
      printf("%4.2lf   %8.5lf",x[i],f[i][0]);
      for(j=1;j<=i;j++)
      {
f[i][j]=(f[i][j-1]-f[i-1][j-1])/(x[i]-x[i-j]);
printf("  %8.5lf",f[i][j]);
      }
      printf("\n");
   }

===============================================#
println()
println()

for i=2:n
    for j=2:i
        y[i][j]=(y[i][j-1]-y[i-1][j-1]) / (x[i]-x[i-j+1])
        # println(y[i][j])
    end
end 

println("  Divided Difference Table: ")
println("=============================")
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:i
        s=@sprintf("\t%8.5f",y[i][j])
        print(s)
    end
    println()
 
end


d=[0.0 for i=1:n]
for i=1:n
    d[i]=y[n][i]
end 

s=@sprintf("牛頓後向差除表")
println(s,d)




輸出畫面

  Divided Difference Table: 
=============================
i x(i)  f(i)  f(i,i+1) f(i,i+1.i+2),  ......................
1  1.00000  0.00000  0.69300 -0.14350  0.02800
2  2.00000  0.69300  0.40600 -0.05950
3  3.00000  1.09900  0.28700
4  4.00000  1.38600
牛頓前向差除表[0.0, 0.693, -0.1435, 0.028]



Newton內插法的差除表P(1.5)=0.39287

Newton內插法的誤差P(1.5)=0.01050000


  Divided Difference Table: 
=============================
i x(i)  f(i)  f(i,i+1) f(i,i+1.i+2),  ......................
1  1.00000  0.00000
2  2.00000  0.69300  0.69300
3  3.00000  1.09900  0.40600 -0.14350
4  4.00000  1.38600  0.28700 -0.05950  0.02800
牛頓後向差除表[1.386, 0.287, -0.0595, 0.028]