先使用牛頓向前的內插多項式
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝====
x f(x) 第一項P'n(x) 前二項P'n(x) 前三項P'n(x)
0.5 0.4794
0.6 0.5646
0.8 0.7174
1.05 0.8674
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝====
P'(n) 第1項 f0,1
P'(n) 第2項 f0,1+f0,1,2 *( (x-x[1]) + (x-x[2]) )
P'(n) 第3項 f0,1+f0,1,2 *( (x-x[1]) + (x-x[2]) )
+ f0,1,2,3 * ( (x-x[1])*(x-x[2]) ) + ( (x-x[2])*(x-x[3]) ) + ( (x-x[3])*(x-x[1]) )
x=[ 0.5 , 0.6 , 0.8 , 1.05]
f= [ [0.4794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[0.5646 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[0.7174 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[0.8674 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ] ]
n=length(x)
println(" Divided Difference Table: ")
println("=============================")
for j=2:n
for i=1:n-j+1
f[i][j]=(f[i+1][j-1]-f[i][j-1])/(x[i+j-1]-x[i])
end
end
print("i\tx(i)\t\tf(i)\t\tf(i,i+1)\tf(i,i+1.i+2), ......................\n")
for i=1:n
s=@sprintf("%d\t%8.5f",i,x[i])
print(s)
for j=1:n-i+1
s=@sprintf("\t%8.5f",f[i][j])
print(s)
end
println()
end
d=[0.0 for i=1:n-1]
for i=1:n-1
d[i]=f[1][i+1]
end
s=@sprintf("\n牛頓一階微分前向差除表")
println(s,d)
xa=[ 0.5 , 0.6 , 0.8 , 1.05]
pn=[0.0 , 0.0 , 0.0 ,0.0]
s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
pn[i]=d[1]
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end
s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
xb=xa[i]
dx=0.0
for j=1:2
#println(xb," ", x[j])
dx=dx+(xb-x[j])
end
#println(dx)
dx=d[2]*dx
pn[i]=pn[i]+dx
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end
s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
xb=xa[i]
dx=0.0
dx=(xb-x[1])*(xb-x[2])+ (xb-x[2])*(xb-x[3])+ (xb-x[3])*(xb-x[1])
dx=d[3]*dx
pn[i]=pn[i]+dx
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end
println("\n\n")
for i=1:length(xa)
xb=xa[i]
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , 誤差 f'(x)-P'(n) = %0.5f",i,xa[i],abs(cos(xb)-pn[i]) )
println(s)
end
輸出畫面
i f(i) f(i,i+1) f(i,i+1.i+2), ...................... 1 0.50000 0.47940 0.85200 -0.29333 -0.12929 2 0.60000 0.56460 0.76400 -0.36444 3 0.80000 0.71740 0.60000 4 1.05000 0.86740 牛頓一階微分前向差除表[0.852, -0.293333, -0.129293] THE RESULTS OF INTERPOLATION: ======================================================== xa[1]= 0.500 , P'n( 0.500 ) = 0.85200 ======================================================== xa[2]= 0.600 , P'n( 0.600 ) = 0.85200 ======================================================== xa[3]= 0.800 , P'n( 0.800 ) = 0.85200 ======================================================== xa[4]= 1.050 , P'n( 1.050 ) = 0.85200 THE RESULTS OF INTERPOLATION: ======================================================== xa[1]= 0.500 , P'n( 0.500 ) = 0.88133 ======================================================== xa[2]= 0.600 , P'n( 0.600 ) = 0.82267 ======================================================== xa[3]= 0.800 , P'n( 0.800 ) = 0.70533 ======================================================== xa[4]= 1.050 , P'n( 1.050 ) = 0.55867 THE RESULTS OF INTERPOLATION: ======================================================== xa[1]= 0.500 , P'n( 0.500 ) = 0.87745 ======================================================== xa[2]= 0.600 , P'n( 0.600 ) = 0.82525 ======================================================== xa[3]= 0.800 , P'n( 0.800 ) = 0.69758 ======================================================== xa[4]= 1.050 , P'n( 1.050 ) = 0.49434
======================================================== xa[1]= 0.500 , 誤差 f'(x)-P'(n) = 0.00013 ======================================================== xa[2]= 0.600 , 誤差 f'(x)-P'(n) = 0.00008 ======================================================== xa[3]= 0.800 , 誤差 f'(x)-P'(n) = 0.00087 ======================================================== xa[4]= 1.050 , 誤差 f'(x)-P'(n) = 0.00323
沒有留言:
張貼留言