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]