Julia語言 例題1-1 若已知下面二點 請問 x=1.5時 則f(1.5)= ??

Julia語言  例題1-1 若已知下面二點

   x              y=f(x)
=================
  1.0            0.000
  2.0            0.693 
請問 x=1.5時 則f(1.5)= ??


程式
using Printf

function lagrange(x::Array{Float64,1},f::Array{Float64,1},xa::Float64)
    #
    # implements the interpolation algorithm of Newton
    #
    # ON ENTRY :
    # x abscisses, given as a column vector;
    # f ordinates, given as a column vector;
    # xa point where to evaluate the interpolating
    # polynomial through (x[i],f[i]).
    #
    # ON RETURN :
    # d divided differences, computed from and f;
    # p value of the interpolating polynomial at xa.
    #
    # EXAMPLE :
    n = length(x)
    tmp2=0.0
    for k=1:n
        tmp1=1.0
        for i=1:n
            if (i != k)
                tmp1 *= (xa-x[i]) / (x[k]-x[i])
            end 
        end
        tmp2=tmp2+tmp1*f[k] 
    end 
     
    return tmp2
end


x = [1.0 , 2.0 ]
f = [0.0 , 0.693 ]
xa = 1.5
result1 = lagrange(x,f,xa)
println("Lagrange 內插法理則 ")
println("x=    " , x)
println("f(x)= " , f)
 
s = @sprintf("Pn(x)=%0.5f" , result1 )
println(s)
s = @sprintf("f(xa)=%0.5f" , log(xa) )
println(s)
s = @sprintf("誤差 =%0.5f" , abs(log(xa)-result1) )
println(s)


輸出畫面

Lagrange 內插法理則 
x=    [1.0, 2.0]
f(x)= [0.0, 0.693]
Pn(x)=0.34650
f(xa)=0.40547
誤差 =0.05897

 

Julia語言例題1-5 已知函數f(x)=ln(x)請使用Lagrange 內插法 求f(1.5) ,f(2.5) f(3.5)的P(x)與f(x) 分別用x=1.5 ,2.5 ,3.5 計算f(x)的e(x)誤差值

Julia語言例題1-5 已知函數f(x)=ln(x) 

xi     f(xi)

==============

1.0    0.000

2.0    0.693

3.0    1.099

4.0    1.386

請使用Lagrange 內插法 求f(1.5) ,f(2.5) f(3.5)的P(x)與f(x) 

分別用x=1.5 ,2.5 ,3.5 計算f(x)的e(x)誤差值

using Printf

function lagrange(x::Array{Float64,1},f::Array{Float64,1},xa::Float64)
    #
    # implements the interpolation algorithm of Newton
    #
    # ON ENTRY :
    # x abscisses, given as a column vector;
    # f ordinates, given as a column vector;
    # xa point where to evaluate the interpolating
    # polynomial through (x[i],f[i]).
    #
    # ON RETURN :
    # d divided differences, computed from and f;
    # p value of the interpolating polynomial at xa.
    #
    # EXAMPLE :
    n = length(x)
    tmp2=0.0
    for k=1:n
        tmp1=1.0
        for i=1:n
            if (i != k)
                tmp1 *= (xa-x[i]) / (x[k]-x[i])
            end    
        end
        tmp2=tmp2+tmp1*f[k]    
    end    
        
    return tmp2
end


function er(xa::Float64,x0::Float64,x1::Float64,x2::Float64,x3::Float64)
    tmp1=0.0
    tmp2=0.0
    sum=1.0
    xm=(x[1]+x[length(x)])/2
    #println(xm)
    for i=1:length(x)
        sum=sum*i
    end 
    #println(sum)
    tmp1=(xa-x0)*(xa-x1)*(xa-x2)*(xa-x3)/sum
    #println(tmp1)
    tmp2=tmp1*(-6/xm^4)
    return tmp2
end

x= [1.0 ,2.0 , 3.0 , 4.0] 
f= [0.0,0.693,1.099,1.386]
xb = [ 1.5 , 2.5 ,3.5]

for i = 1:length(xb)
    xa=xb[i]
    result1 = lagrange(x,f,xa)
    println("Lagrange 內插法理則 ")
    println("x=    " , x)
    println("f(x)= " , f)
    println("xa=   " , xa)
    
   
    s = @sprintf("Pn(x)=%0.5f" , result1 )
    println(s)
    s = @sprintf("f(xa)=%0.5f" , log(xa) )
    println(s)
    s = @sprintf("誤差 =%0.5f" , abs(log(xa)-result1) )
    println(s)
    println()

end

for i=1:length(xb)
    result2=0.0
    xa=xb[i]
    x0=x[1]
    x1=x[2]
    x2=x[3]
    x3=x[4]
    result2= er(xa, x0 ,x1 , x2 , x3)
    s=@sprintf("x= %0.2f 誤差e(x)=%0.5f",xa ,result2)
    println(s)
end    

輸出畫面

Lagrange 內插法理則 
x=    [1.0, 2.0, 3.0, 4.0]
f(x)= [0.0, 0.693, 1.099, 1.386]
xa=   1.5
Pn(x)=0.39287
f(xa)=0.40547
誤差 =0.01259

Lagrange 內插法理則 
x=    [1.0, 2.0, 3.0, 4.0]
f(x)= [0.0, 0.693, 1.099, 1.386]
xa=   2.5
Pn(x)=0.92138
f(xa)=0.91629
誤差 =0.00508

Lagrange 內插法理則 
x=    [1.0, 2.0, 3.0, 4.0]
f(x)= [0.0, 0.693, 1.099, 1.386]
xa=   3.5
Pn(x)=1.24687
f(xa)=1.25276
誤差 =0.00589

x= 1.50 誤差e(x)=0.00600
x= 2.50 誤差e(x)=-0.00360
x= 3.50 誤差e(x)=0.00600

Julia語言例題EX4-7-B利用梯形法求雙重積分 f (x,y) = x* exp(y) a=0.0 , b=1.0 , c(x) =0.0 d(x)= x 取 n=2 n-4 並比較誤差值

Julia語言例題EX4-7-B利用梯形法求雙重積分 f (x,y) = x* exp(y) a=0.0 , b=1.0 , c(x) =0.0 d(x)= x 取 n=2  n=4 並比較誤差值


#==================================================
/* ex4-7-B.jl based on Trapezoidal Rule to
 * compute the double integral.
 */
===================================================#
using Printf

function F(x::Float64, y::Float64) #//  f(x,y)= x exp(y)
    return  (x*exp(y))
end

function C(x::Float64) #//  c(x)= 0.0
    return  (0.0)
end

function D(x::Float64) #//  d(x)= x
    return  (x)
end


function gy(n::Int64)
    sum1=0.0
    sum2=0.0
    for i=0:n
        for j=1:n-1
            if(i%2==0)
                sum2=sum2+F(x[i+1],y[i+1][j+1])
            else
                sum1=sum1+F(x[i+1],y[i+1][j+1])
            end 
        end
        g1[i+1]=(hy[i+1]/3)*(F(x[i+1],y[i+1][1])+F(x[i+1],y[i+1][n+1])+2*sum2+4*sum1)
        println(i+1,"----",g1[i+1])
        sum1=0.0
        sum2=0.0
    end
    return g1
end


s=@sprintf("辛普森積分計算雙重積分 n=2 ")
println(s)

x= [0.0 for i=1:5 ]
g1=[0.0 for i=1:5 ]
g2=[0.0 for i=1:5 ]
hy=[0.0 for i=1:5 ]

f(i) = [0.0 for i=1:5]
y= f.([0.0 for i=1:5])

sum1=0.0
sum2=0.0

n=2
a=0.0
b=1.0

hx=(b-a)/n
ts=0.0
for i=0:n
    x[i+1]=a+i*hx
    hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
    for j=0:n
        y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
        #println(y)
    end
end

#println(x)
#println(hy)
#println(y)

g2=gy(n)
sum1=0.0
sum2=0.0

println("\n\n")
for i=1:n-1
    if(i%2==0)
        sum2=sum2+g2[i+1]
    else
        sum1=sum1+g2[i+1]
    end
 
end 

ts= (hx/3) * (g1[1] + g1[n+1] + 4*sum1 + 2*sum2)
s=@sprintf("辛普森積分計算雙重積分結果 T%d=%0.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%0.6lf\n",(1/2))
println(s)
tn=abs( 1/2 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)


s=@sprintf("梯形積分計算雙重積分 n=4 ")
println(s)

x= [0.0 for i=1:10 ]
g1=[0.0 for i=1:10 ]
g2=[0.0 for i=1:10 ]
hy=[0.0 for i=1:10 ]

f(i) = [0.0 for i=1:10]
y= f.([0.0 for i=1:10])


sum1=0.0
sum2=0.0

n=4
a=0.0
b=1.0

hx=(b-a)/n
ts=0.0
for i=0:n
    x[i+1]=a+i*hx
    hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
    for j=0:n
        y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
        #println(y)
    end
end

#println(x)
#println(hy)
#println(y)

g2=gy(n)
sum1=0.0
sum2=0.0

println("\n\n")
for i=1:n-1
    if(i%2==0)
        sum2=sum2+g2[i+1]
    else
        sum1=sum1+g2[i+1]
    end
 
end 

ts= (hx/3) * (g1[1] + g1[n+1] + 4*sum1 + 2*sum2)
s=@sprintf("辛普森積分計算雙重積分結果 T%d=%0.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%0.6lf\n",(1/2))
println(s)
tn=abs( 1/2 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)

輸出畫面

辛普森積分計算雙重積分 n=2 
1----0.0
2----0.3243676223937956
3----1.16928739497655



辛普森積分計算雙重積分結果 T2=0.411126

實際值=0.500000

誤差值=0.088874

梯形積分計算雙重積分 n=4 
1----0.0
2----0.0828099899078667
3----0.21652191332178472
4----0.9741611859661217
5----1.151481269705011



辛普森積分計算雙重積分結果 T4=0.484367

實際值=0.500000

誤差值=0.015633

Julia語言例題EX4-7-A利用梯形法求雙重積分 f (x,y) = x* exp(y) a=0.0 , b=1.0 , c(x) =0.0 d(x)= x 取 n=2 與 n=4 二者並比較誤差值

Julia語言例題EX4-7-A利用梯形法求雙重積分
f (x,y) = x* exp(y)   a=0.0 , b=1.0 , c(x) =0.0   d(x)= x
取 n=2 與 n=4 二者並比較誤差值

#==================================================
/* ex4-7-A.jl based on Trapezoidal Rule to
 * compute the double integral.
 */
===================================================#
using Printf

function F(x::Float64, y::Float64) #//  f(x,y)= x exp(y)
    return  (x*exp(y))
end

function C(x::Float64) #//  c(x)= 0.0
    return  (0.0)
end

function D(x::Float64) #//  d(x)= x
    return  (x)
end


function gy(n::Int64)
    sum=0.0;
    for i=0:n
        for j=1:n-1
          sum=sum+F(x[i+1],y[i+1][j+1])
          #println(y[i][j])
        end
        g1[i+1]=(0.5*hy[i+1])*(F(x[i+1],y[i+1][1])+F(x[i+1],y[i+1][n+1])+2*sum)
        println(i+1,"----",g1[i+1])
        sum=0.0;
    end
    return g1
end


s=@sprintf("梯形積分計算雙重積分 n=2 ")
println(s)

x= [0.0 for i=1:20 ]
g1=[0.0 for i=1:20 ]
g2=[0.0 for i=1:20 ]
hy=[0.0 for i=1:20 ]

f(i) = [0.0 for i=1:20]
y= f.([0.0 for i=1:20])

sum=0.0
n=2
a=0.0
b=1.0

hx=(b-a)/n
ts=0.0
for i=0:n
    x[i+1]=a+i*hx
    hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
    for j=0:n
        y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
        #println(y)
    end
end


g2=gy(n)
sum1=0.0
println("\n\n")
for i=1:n-1
    sum1=sum1+g2[i+1];
    #println(g2[i+1])
end 

ts= (hx/2) * (g1[1] + g1[n+1] + 2*sum1)
s=@sprintf("梯形積分計算雙重積分結果 T%d=%0.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%0.6lf\n",(1/2))
println(s)
tn=abs( 1/2 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)


s=@sprintf("梯形積分計算雙重積分 n=4 ")
println(s)

x= [0.0 for i=1:20 ]
g1=[0.0 for i=1:20 ]
g2=[0.0 for i=1:20 ]
hy=[0.0 for i=1:20 ]

f(i) = [0.0 for i=1:20]
y= f.([0.0 for i=1:20])

sum=0.0
n=4
a=0.0
b=1.0

hx=(b-a)/n
ts=0.0
for i=0:n
    x[i+1]=a+i*hx
    hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
    for j=0:n
        y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
        #println(y)
    end
end


g2=gy(n)
sum1=0.0
println("\n\n")
for i=1:n-1
    sum1=sum1+g2[i+1];
    #println(g2[i+1])
end 

ts= (hx/2) * (g1[1] + g1[n+1] + 2*sum1)
s=@sprintf("梯形積分計算雙重積分結果 T%d=%0.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%0.6lf\n",(1/2))
println(s)
tn=abs( 1/2 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)


輸出畫面

梯形積分計算雙重積分 n=2 
1----0.0
2----0.3260482565047257
3----1.7539310924648253



梯形積分計算雙重積分結果 T2=0.601507

實際值=0.500000

誤差值=0.101507

梯形積分計算雙重積分 n=4 
1----0.0
2----0.07102946671483652
3----0.3247828699826771
4----0.8402029213086307
5----1.7272219045575166



梯形積分計算雙重積分結果 T4=0.524907

實際值=0.500000

誤差值=0.024907

Julia語言 例題1-2 利用Lagrange 已知下列數據 求 xa=1.5 的值並計算誤差值=?

Julia語言 例題1-2 利用Lagrange 已知下列數據 求 xa=1.5 的值並計算誤差值=?
 x      f(x)
=======
1.0    0.0
2.0    0.693
3.0    1.099

程式
using Printf

function lagrange(x::Array{Float64,1},f::Array{Float64,1},xa::Float64)
    #
    # implements the interpolation algorithm of Newton
    #
    # ON ENTRY :
    # x abscisses, given as a column vector;
    # f ordinates, given as a column vector;
    # xa point where to evaluate the interpolating
    # polynomial through (x[i],f[i]).
    #
    # ON RETURN :
    # d divided differences, computed from and f;
    # p value of the interpolating polynomial at xa.
    #
    # EXAMPLE :
    n = length(x)
    tmp2=0.0
    for k=1:n
        tmp1=1.0
        for i=1:n
            if (i != k)
                tmp1 *= (xa-x[i]) / (x[k]-x[i])
            end 
        end
        tmp2=tmp2+tmp1*f[k] 
    end 
     
    return tmp2
end


x = [1.0 , 2.0 , 3.0 ]
f = [0.0 , 0.693 , 1.099]
xa = 1.5
result1 = lagrange(x,f,xa)
println("Lagrange 內插法理則 ")
println("x=    " , x)
println("f(x)= " , f)
 
s = @sprintf("Pn(x)=%0.5f" , result1 )
println(s)
s = @sprintf("f(xa)=%0.5f" , log(xa) )
println(s)
s = @sprintf("誤差 =%0.5f" , abs(log(xa)-result1) )
println(s)
 
輸出畫面

Lagrange 內插法理則 
x=    [1.0, 2.0, 3.0]
f(x)= [0.0, 0.693, 1.099]
Pn(x)=0.38237
f(xa)=0.40547
誤差 =0.02309

Julia語言 例題1-1 利用Lagrange 已知下列數據 求 xa=1.5 的值並計算誤差值=?

Julia語言 例題1-1 利用Lagrange 已知下列數據 求 xa=1.5 的值並計算誤差值=?
 x       f(x)
==========
1.0     0.0
2.0     0.693

程式
using Printf

function lagrange(x::Array{Float64,1},f::Array{Float64,1},xa::Float64)
    #
    # implements the interpolation algorithm of Newton
    #
    # ON ENTRY :
    # x abscisses, given as a column vector;
    # f ordinates, given as a column vector;
    # xa point where to evaluate the interpolating
    # polynomial through (x[i],f[i]).
    #
    # ON RETURN :
    # d divided differences, computed from and f;
    # p value of the interpolating polynomial at xa.
    #
    # EXAMPLE :
    n = length(x)
    tmp2=0.0
    for k=1:n
        tmp1=1.0
        for i=1:n
            if (i != k)
                tmp1 *= (xa-x[i]) / (x[k]-x[i])
            end 
        end
        tmp2=tmp2+tmp1*f[k] 
    end 
     
    return tmp2
end


x = [1.0 , 2.0 ]
f = [0.0 , 0.693 ]
xa = 1.5
result1 = lagrange(x,f,xa)
println("Lagrange 內插法理則 ")
println("x=    " , x)
println("f(x)= " , f)
 
s = @sprintf("Pn(x)=%0.5f" , result1 )
println(s)
s = @sprintf("f(xa)=%0.5f" , log(xa) )
println(s)
s = @sprintf("誤差 =%0.5f" , abs(log(xa)-result1) )
println(s)
 

輸出畫面

Lagrange 內插法理則 
x=    [1.0, 2.0]
f(x)= [0.0, 0.693]
Pn(x)=0.34650
f(xa)=0.40547
誤差 =0.05897

Julia語言 例題1-4 利用Lagrange 已知下列數據 求 xa=1.5 的值並計算誤差值=?

Julia語言 例題1-4 利用Lagrange 已知下列數據 求 xa=1.5
的值並計算誤差值=?  真實值f(x)=ln(x) , f(1.5)=ln(1.5=0.405
===============
(1)
xa=1.5
1.0  0.0
2.0  0.693
3.0  1.099
4.0  1.386
===============
(2)
xa=1.5
1.0  0.000
1.2  0.182
1.4  0.336
2.0  0.693
===============
(3)
xa=1.5
1.0  0.000
1.2  0.182
1.7  0.531
2.0  0.693
2.2  0.788
2.7  0.993
3.0  1.099
3.2  1.163
3.7  1.308
4.0  1.386
===============
程式
using Printf

function lagrange(x1::Array{Float64,1},f1::Array{Float64,1},xa::Float64)
    #
    # implements the interpolation algorithm of Newton
    #
    # ON ENTRY :
    # x abscisses, given as a column vector;
    # f ordinates, given as a column vector;
    # xa point where to evaluate the interpolating
    # polynomial through (x[i],f[i]).
    #
    # ON RETURN :
    # d divided differences, computed from and f;
    # p value of the interpolating polynomial at xa.
    #
    # EXAMPLE :
    n = length(x1)
    tmp2=0.0
    for k=1:n
        tmp1=1.0
        for i=1:n
            if (i != k)
                tmp1 *= (xa-x1[i]) / (x1[k]-x1[i])
            end 
        end
        tmp2=tmp2+tmp1*f1[k] 
    end 
     
    return tmp2
end


x = [ [1.0 , 2.0 ,3.0 , 4.0],
      [1.0 , 1.2 ,1.4 , 2.0],
      [1.0 , 1.2 , 1.7 , 2.0  , 2.2 ,  2.7  , 3.0 , 3.2  , 3.7  , 4.0]]
   
f = [ [0.0 , 0.693 , 1.099 , 1.386],
      [0.0 , 0.182 , 0.336 , 0.693],
      [0.0 , 0.182 , 0.531 , 0.693  , 0.788  ,  0.993 , 1.099 , 1.163 , 1.308 , 1.386] ]
     
xa = 1.5
for m=1:3
    x1=x[m]
    f1=f[m]
    result1 = lagrange(x1,f1,xa)
    println("Lagrange 內插法理則 ")
    println("x=    " , x1)
    println("f(x)= " , f1)
 
    s = @sprintf("Pn(x)=%0.5f" , result1 )
    println(s)
    s = @sprintf("f(xa)=%0.5f" , log(xa) )
    println(s)
    s = @sprintf("誤差 =%0.5f" , abs(log(xa)-result1) )
    println(s)
 
    println()
end

輸出畫面

Lagrange 內插法理則 
x=    [1.0, 2.0, 3.0, 4.0]
f(x)= [0.0, 0.693, 1.099, 1.386]
Pn(x)=0.39287
f(xa)=0.40547
誤差 =0.01259

Lagrange 內插法理則 
x=    [1.0, 1.2, 1.4, 2.0]
f(x)= [0.0, 0.182, 0.336, 0.693]
Pn(x)=0.40447
f(xa)=0.40547
誤差 =0.00100

Lagrange 內插法理則 
x=    [1.0, 1.2, 1.7, 2.0, 2.2, 2.7, 3.0, 3.2, 3.7, 4.0]
f(x)= [0.0, 0.182, 0.531, 0.693, 0.788, 0.993, 1.099, 1.163, 1.308, 1.386]
Pn(x)=0.40614
f(xa)=0.40547
誤差 =0.00068

Julia語言 範例1-7 牛頓內插法 差除表

Julia語言 範例1-7 牛頓內插法 差除表
已知 下面3點座標值利用牛頓內插法求差除表
x = [0.0, 1.0 , 2.0 ]
f = [0.0 , -3.0 , 0.0]

==========程式==============
using Printf

function divdif(x::Array{Float64,1},f::Array{Float64,1})
    #
    # Returns the vector of divided differences for the
    # Newton form of the interpolating polynomial.
    #
    # ON ENTRY:
    # x abscisses, given as a column vector;
    # f ordinates, given as a column vector.
    #
    # ON RETURN:
    # d divided differences
    #
    #
    n = length(x)
    d = deepcopy(f)
    for i=2:n
        for j=1:i-1
            d[i] = (d[j] - d[i])/(x[j] - x[i])
        end
     
    end
    return d
end

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

function newton(x::Array{Float64,1},f::Array{Float64,1},xx::Float64)
    #
    # implements the interpolation algorithm of Newton
    #
    # ON ENTRY :
    # x abscisses, given as a column vector;
    # f ordinates, given as a column vector;
    # xx point where to evaluate the interpolating
    # polynomial through (x[i],f[i]).
    #
    # ON RETURN :
    # d divided differences, computed from and f;
    # p value of the interpolating polynomial at xx.
    #
    # EXAMPLE :
    divided = divdif(x,f)
    result = newtonform(x,divided,xx)
    return divided, result
end


x = [0.0, 1.0 , 2.0 ]
f = [0. , -3.0 , 0.0]
xx = 1.5
d, p = newton(x,f,xx)
print("Newton內插法的差除表")
println(d)



輸出畫面

Newton內插法的差除表[0.0, -3.0, 3.0]

Julia語言 例題1-8 已知5點的座標值 使用牛頓內插法求差除表？

Julia語言 例題1-8 已知5點的座標值 使用牛頓內插法求差除表？
x = [0.0, 1.0  , 2.0  , 3.0  , 4.0 ]
f = [0.0  , -3.0 , 0.0 , 15.0, 48.0]
並計算P(1.5)與P(2.5)之值=?

程式

using Printf

function divdif(x::Array{Float64,1},f::Array{Float64,1})
    #
    # Returns the vector of divided differences for the
    # Newton form of the interpolating polynomial.
    #
    # ON ENTRY:
    # x abscisses, given as a column vector;
    # f ordinates, given as a column vector.
    #
    # ON RETURN:
    # d divided differences
    #
    #
    n = length(x)
    d = deepcopy(f)
    for i=2:n
        for j=1:i-1
            d[i] = (d[j] - d[i])/(x[j] - x[i])
        end
    end
    return d
end

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

function newton(x::Array{Float64,1},f::Array{Float64,1},xx::Float64)
    #
    # implements the interpolation algorithm of Newton
    #
    # ON ENTRY :
    # x abscisses, given as a column vector;
    # f ordinates, given as a column vector;
    # xx point where to evaluate the interpolating
    # polynomial through (x[i],f[i]).
    #
    # ON RETURN :
    # d divided differences, computed from and f;
    # p value of the interpolating polynomial at xx.
    #
    # EXAMPLE :
    divided = divdif(x,f)
    result = newtonform(x,divided,xx)
    return divided, result
end


x = [0.0, 1. , 2. , 3. , 4. ]
f = [0. , -3.0 , 0. , 15.0, 48.0]
xx = 2.5
d, p = newton(x,f,xx)
print("Newton內插法的差除表")
println(d)
s = @sprintf("Newton內插法的Pn(%0.2f)值= %0.5f " , float(xx) , float(p) )
println(s)

xx = 1.5
d, p = newton(x,f,xx)
print("Newton內插法的差除表")
println(d)
s = @sprintf("Newton內插法的Pn(%0.2f)值= %0.5f " , float(xx) , float(p) )
println(s)

輸出結果

Newton內插法的差除表[0.0, -3.0, 3.0, 1.0, -0.0]
Newton內插法的Pn(2.50)值= 5.62500 
Newton內插法的差除表[0.0, -3.0, 3.0, 1.0, -0.0]
Newton內插法的Pn(1.50)值= -2.62500 

Julia語言 求Lagrange 的值 [ 1.5 , 2.5 ,3.5] 與誤差

Julia語言 求Lagrange 的值 [ 1.5 , 2.5 ,3.5] 與誤差 

x= [1.0 ,2.0 , 3.0 , 4.0] 
y= [0.0,0.693,1.099,1.386]
xa = [ 1.5 , 2.5 ,3.5]

程式

using Printf


#This function returns another function, which is the Lagrange Interpolant of the values xvals and yvals.
function LagrangeInterpolantGenerator(xvals,yvals)
    function LagrangeInterpolant(x)
        numvalstoevaluate = length(x)
        numvalstoevaluate == 1 ? output = 0 : output = zeros(numvalstoevaluate)
        for k = 1:numvalstoevaluate
            N = length(xvals)
            LagrangePolynomials = ones(N)
            for i in 1:N  
                for j in [1:i-1;i+1:N]     #Surprisingly, this works even in the i=1 and i=N cases.
                    LagrangePolynomials[i] = LagrangePolynomials[i].*(x[k]-xvals[j])./(xvals[i]-xvals[j])
                end
            end
            numvalstoevaluate == 1 ? output = sum(LagrangePolynomials.*yvals) : output[k] = sum(LagrangePolynomials.*yvals)
        end
        return output
    end
    return LagrangeInterpolant
end

#Examples
x= [1.0 ,2.0 , 3.0 , 4.0] 
y= [0.0,0.693,1.099,1.386]
xa = [ 1.5 , 2.5 ,3.5]

for i = 1:3

    interpolantfunc = LagrangeInterpolantGenerator(x,y)
    xb=xa[i]
    a=interpolantfunc(xb)   #returns  0.34650 
    s = @sprintf("Lagrange插值法 P(%0.2f) = %0.5f " , float(xb) , float(a) )
    println(s)

end

輸出結果

Lagrange插值法 P(1.50) = 0.39287 
Lagrange插值法 P(2.50) = 0.92138 
Lagrange插值法 P(3.50) = 1.24687 

using Printf


#This function returns another function, which is the Lagrange Interpolant of the values xvals and yvals.
function LagrangeInterpolantGenerator(xvals,yvals)
    function LagrangeInterpolant(x)
        numvalstoevaluate = length(x)
        numvalstoevaluate == 1 ? output = 0 : output = zeros(numvalstoevaluate)
        for k = 1:numvalstoevaluate
            N = length(xvals)
            LagrangePolynomials = ones(N)
            for i in 1:N  
                for j in [1:i-1;i+1:N]     #Surprisingly, this works even in the i=1 and i=N cases.
                    LagrangePolynomials[i] = LagrangePolynomials[i].*(x[k]-xvals[j])./(xvals[i]-xvals[j])
                end
            end
            numvalstoevaluate == 1 ? output = sum(LagrangePolynomials.*yvals) : output[k] = sum(LagrangePolynomials.*yvals)
        end
        return output
    end
    return LagrangeInterpolant
end

#Examples
x= [1.0 ,2.0 , 3.0 , 4.0] 
y= [0.0,0.693,1.099,1.386]
xa = [ 1.5 , 2.5 ,3.5]

for i = 1:3

    interpolantfunc = LagrangeInterpolantGenerator(x,y)
    xb=xa[i]
    a=interpolantfunc(xb)   #returns  0.34650 
    s = @sprintf("Lagrange插值法 P(%0.2f) = %0.5f  " , float(xb) , float(a) )
    println(s)
    xc = abs( float(a) - float(log(xb)) )
    s = @sprintf("Lagrange插值法 P(%0.2f) 與 真實值的誤差 = %0.5f " , float(xb) ,  float(xc)   )
    println(s)
   

end

Lagrange插值法 P(1.50) = 0.39287 
Lagrange插值法 P(1.50) 與 真實值的誤差 = 0.01259 
Lagrange插值法 P(2.50) = 0.92138 
Lagrange插值法 P(2.50) 與 真實值的誤差 = 0.00508 
Lagrange插值法 P(3.50) = 1.24687 
Lagrange插值法 P(3.50) 與 真實值的誤差 = 0.00589