Julia語言 範例1-13 已知不明分行物體, 經光測儀器之量測,印出Hermite 差除表

Julia語言 範例1-13 已知不明分行物體, 經光測儀器之量測,印出Hermite 差除表

範例1-13 已知不明分行物體, 經光測儀器之量測其結果如下
=========================================
時間                位置               飛行速度
 x                       f(x)                  f ' (x)
 0                       50.00               50.00
 2                       216.67             94.44
 4                       410.00             98.00
=========================================
印出Hermite 差除表
當xa=1.5時 求 Hn(xa)之值

程式

using Printf

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

i z(i)   f(i)   f(i-1,i)   f(i-2,i-1,i)...

0  0.0 50.0000

1  0.0 50.0000  50.0000

2  2.0 216.6700 83.3350 16.6675

3  2.0 216.6700 94.4400  5.5525 -5.5575

4  4.0 410.0000 96.6650  1.1125 -1.1100  1.1119

5  4.0 410.0000 98.0000  0.6675 -0.2225  0.2219 -0.2225

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 , 2.0 ,4.0]

f=[50.00 ,  216.67 ,  410.00]

ff=[50.00 ,  94.44  ,  98.00]

n=length(x)

z=[0.0 for i=1:2*n]

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

xa=1.0
p=1
H=q[1][1]

for i=2:2*n
    p *= (xa-z[i-1])
    H = H+p*q[i][i]
   
    #println(q[i][i])
    #println(p)
    #println(H)
    #println()
   
end
println()

println("Xa=",xa," Hn(xa)=" ,H, "---Hermite方法")
H1=(100*xa + (50.0/(xa+1)))
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 50.0000
2 0.0000 50.0000  50.0000
3 2.0000 216.6700 83.3350
3 2.0000 216.6700 83.3350 16.6675
4 2.0000 216.6700 94.4400 5.5525 -5.5575
5 4.0000 410.0000 96.6650 1.1125 -1.1100 1.1119
6 4.0000 410.0000 98.0000 0.6675 -0.2225 0.2219 -0.2225

Xa=1.0 Hn(xa)=124.004375---Hermite方法
Xa=1.0 f(xa)=125.0---真實值
Hn(x)與fn(x)的誤差0.995625000000004