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

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

已知下面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]之值

using Printf

#======================================================

function [H]=hermite(X,x,f,fd);
m=length(x);
for i=1:m
    z(2*i-1)=x(i);
    Q(2*i-1,1)=f(i);
    z(2*i)=x(i);
    Q(2*i,1)=f(i);
    Q(2*i,2)=fd(i);
    if i~=1
        Q(2*i-1,2)=(Q(2*i-1,1)-Q(2*i-2,1))/(z(2*i-1)-z(2*i-2));
    end;
end;
for i=3:2*m
    for j=3:i
        Q(i,j)=(Q(i,j-1)-Q(i-1,j-1))/(z(i)-z(i-j+1));
    end
end

p=1;
H=Q(1,1);
for i=2:2*m
    p=p.*(X-z(i-1));
    H=H+p.*Q(i,i);
end;

======================================================#
println("Hermite 差除表")
x=[ 1.0 ,2.0 ,3.0  ,4.0 ]
f=[0.0 ,  0.693 , 1.099 , 1.386]
ff=[1.00 , 0.5 , 0.3333 , 0.25 ]

xx=[1.5 , 2.5 , 3.5 , 4.5]
n=length(x)
z=[0.0 for i=1:2*n]

# q 需為 f 的 2n * 2n 倍
q= [[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
    [ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
    [ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
    [ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
    [ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
    [ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
    [ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
    [ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ]]

           
println("i\tz(i)\tf(i)\t\tf(i-1,i)\tf(i-2,i-1,i)...\n");

for i=1:n
    z[2*i-1]=x[i]
    z[2*i]=x[i]
   
    q[2*i-1][1]=f[i]
    q[2*i][1]=f[i]
    q[2*i][2]=ff[i]
   
    if (i!=1)
        q[2*i-1][2]= (q[2*i-1][1]- q[2*i-2][1]) / (z[2*i-1]-z[2*i-2])
    end
    #println(i,z,q)
end

for i=1:n
    if(i==1)
        s=@sprintf("%d\t%0.4f\t%0.4f\n",i,z[i],q[i][1])
        print(s)
    else (i==2)
        s=@sprintf("%d\t%0.4f\t%0.4f\t\t%0.4f\n",i,z[i],q[i][1],q[i][2])
        print(s)
    end
end


for i=3:2*n
    s=@sprintf("%d\t%0.4f\t%0.4f\t\t%0.4f",i,z[i],q[i][1],q[i][2])
    print(s)
    for j=3:i
        q[i][j]=(q[i][j-1]-q[i-1][j-1])/(z[i]-z[i-j+1])
        s=@sprintf("\t%0.4f",q[i][j])
        print(s)
    end
    println()
end

#println(q)


for m=1:length(xx)
    xa=xx[m]
    println()
    println("============================================")
    println("i=",m,"  xa=",xa)
 
    p=1
    H=q[1][1]
    for i=2:2*n
        p *= (xa-z[i-1])
        H = H+p*q[i][i]
    end
    println()

    println("Xa=",xa," Hn(xa)=" ,H, "---Hermite方法")
    H1=(log(xa))
    println("真實解 f(x)= ln(x)")
    println("Xa=",xa," f(xa)=" ,H1 , "---真實值")
    println("Hn(x)與fn(x)的誤差" ,abs(H1-H))
end

輸出畫面

Hermite 差除表
i z(i) f(i)  f(i-1,i) f(i-2,i-1,i)...

1 1.0000 0.0000
2 1.0000 0.0000  1.0000
3 2.0000 0.6930  0.6930
4 2.0000 0.6930  0.5000
3 2.0000 0.6930  0.6930 -0.3070
4 2.0000 0.6930  0.5000 -0.1930 0.1140
5 3.0000 1.0990  0.4060 -0.0940 0.0495 -0.0323
6 3.0000 1.0990  0.3333 -0.0727 0.0213 -0.0141 0.0091
7 4.0000 1.3860  0.2870 -0.0463 0.0132 -0.0040 0.0034 -0.0019
8 4.0000 1.3860  0.2500 -0.0370 0.0093 -0.0019 0.0011 -0.0008 0.0004

============================================
i=1  xa=1.5

Xa=1.5 Hn(xa)=0.4057314453125---Hermite方法
真實解 f(x)= ln(x)
Xa=1.5 f(xa)=0.4054651081081644---真實值
Hn(x)與fn(x)的誤差0.00026633720433560937

============================================
i=2  xa=2.5

Xa=2.5 Hn(xa)=0.9164583984375---Hermite方法
真實解 f(x)= ln(x)
Xa=2.5 f(xa)=0.9162907318741551---真實值
Hn(x)與fn(x)的誤差0.0001676665633448815

============================================
i=3  xa=3.5

Xa=3.5 Hn(xa)=1.2529150390625001---Hermite方法
真實解 f(x)= ln(x)
Xa=3.5 f(xa)=1.252762968495368---真實值
Hn(x)與fn(x)的誤差0.00015207056713206768

============================================
i=4  xa=4.5

Xa=4.5 Hn(xa)=1.5076044921874998---Hermite方法
真實解 f(x)= ln(x)
Xa=4.5 f(xa)=1.5040773967762742---真實值
Hn(x)與fn(x)的誤差0.0035270954112256447