Julia語言例題1-12 已知3點座標及導數 ,建立 Hermite 向後差除表

Julia語言例題1-12 已知3點座標及導數 ,建立 Hermite 向後差除表

例題1-12 已知3點座標及導數如下

=================================
i       xi        f(xi)          f'(xi)
0      0.0      1.0000       1.0000
1      1.0      2.7183        2.7183
2      2.0      2.7183        2.7183
==================================
建立 Hermite 向後差除表
Hn(1.5) 之值 ,  H'n(1.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=[0.0 , 1.0 ,2.0]
f=[1.00 ,  2.7183 ,  7.3891]
ff=[1.00 ,  2.7183 ,  7.3891]

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 ]]

           
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)

xa=1.5
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=(exp(xa))
println("真實解 f(x)=exp(x)")
println("Xa=",xa," f(xa)=" ,H1 , "---真實值")
println("Hn(x)與fn(x)的誤差" ,abs(H1-H))


輸出畫面

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

1 0.0000 1.0000
2 0.0000 1.0000  1.0000
3 1.0000 2.7183  1.7183
3 1.0000 2.7183  1.7183 0.7183
4 1.0000 2.7183  2.7183 1.0000 0.2817
5 2.0000 7.3891  4.6708 1.9525 0.4762 0.0973
6 2.0000 7.3891  7.3891 2.7183 0.7658 0.1448 0.0238

Xa=1.5 Hn(xa)=4.481125---Hermite方法
真實解 f(x)=exp(x)
Xa=1.5 f(xa)=4.4816890703380645---真實值
Hn(x)與fn(x)的誤差0.000564070338064937