範例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
沒有留言:
張貼留言