Julia語言 例題4-4使用
求 f(x)= (1/x) a=1.0 , b=2.0 求n=??
辛普森積分法需要 n 次 n=8
梯形積分法需要 n 次 n=41
x2
為了方便，符號簡化: ∫ dx = ∫ , h = x2-x0
x0
Thm:
Given [x0,x2] , x1 = (x0+x2)/2, assume f in C^4[x0,x2]
∫f = (h/3) * ( f(x0) + 4f(x1) + f(x2) ) - (h^5/90)*f^(4)(ξ)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~
這項是標準的辛普森法 這項是誤差
#========================================================
/* ex4-4.jl based on Simpson's Rule to compute
* definite integral with domain [a,b] and
* n even-grid. n must be even.
*/
========================================================#
using Printf
function F4(x::Float64) #// 四次微分後函數
return (24/x^5)
end
function F1(x::Float64) #// 欲微分函數
return (1/x)
end
function F2(x::Float64) #// 二次微分後函數
return (2/x^3)
end
a=1.0
b=2.0
TOL=0.0001
xa=F4(a)
xb=F4(b)
if xa>xb
k=xa
else
k=xb
end
n1= k* (b-a)^5
n2= n1 / (180*TOL)
n3= sqrt(n2)
n4=sqrt(n3)
n=Int32(round(n4+0.5))
if mod(n,2)!=0
n=n+1
end
println("---------------------------------")
println("辛普森積分法需要",n,"次")
#驗證
println("---------------------------------")
sn=0.0
a=1.0
b=2.0
m=n/2
h=(b-a)/n
sum1=0.0
sum2=0.0
for i=1:2*m-1
x=a+i*h;
if(i%2==0)
sum2=sum2+F1(x)
sn= (h/3.0)*(F1(a)+F1(b)+2.0*sum2+4.0*sum1)
#s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F1(x))
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f err=%0.6f ",i,x,F1(abs(x)), abs(sn-(log(2)-log(1))))
println(s)
else
sum1=sum1+F1(x)
sn= (h/3.0)*(F1(a)+F1(b)+2.0*sum2+4.0*sum1)
#s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F1(x))
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f err=%0.6f ",i,x,F1(abs(x)), abs(sn-(log(2)-log(1))))
println(s)
end
end
print("辛普森積分法 ")
s=@sprintf("S%d=%lf 誤差=%0.6f\n",n,sn ,abs(sn-(log(2)-log(1))))
println(s)
#梯形法
a=1.0
b=2.0
TOL=0.0001
xa=F2(a)
xb=F2(b)
if xa>xb
k=xa
else
k=xb
end
n1= k* (b-a)^5
n2= n1 / (12*TOL)
n4= sqrt(n2)
n=Int32(round(n4+0.5))
println("---------------------------------")
println("梯形積分法需要",n,"次")
println("---------------------------------")
sn=0.0
a=1.0
b=2.0
h=(b-a)/n
x=a
result=0.0
for i=1:n-1
x=x+h
result=result+F1(abs(x))
tn=(h/2.0)*(F1(abs(a))+F1(abs(b))+2.0*result)
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f err=%0.6f ",i,x,F1(abs(x)), abs(tn-(log(2)-log(1))))
println(s)
end
#tn=(h/2.0)*(F1(abs(a))+F1(abs(b))+2.0*result)
print("梯形積分法 ")
s=@sprintf("S%d=%lf 誤差=%0.6f\n",n,tn ,abs(tn-(log(2)-log(1))))
println(s)
輸出畫面
---------------------------------
辛普森積分法需要8次
---------------------------------
i= 1 , x=1.13 ---- f(x)=0.888889 err=0.482499
i= 2 , x=1.25 ---- f(x)=0.800000 err=0.415832
i= 3 , x=1.38 ---- f(x)=0.727273 err=0.294620
i= 4 , x=1.50 ---- f(x)=0.666667 err=0.239065
i= 5 , x=1.63 ---- f(x)=0.615385 err=0.136501
i= 6 , x=1.75 ---- f(x)=0.571429 err=0.088882
i= 7 , x=1.88 ---- f(x)=0.533333 err=0.000007
辛普森積分法 S8=0.693155 誤差=0.000007
---------------------------------
梯形積分法需要41次
---------------------------------
i= 1 , x=1.02 ---- f(x)=0.976190 err=0.651045
i= 2 , x=1.05 ---- f(x)=0.953488 err=0.627789
i= 3 , x=1.07 ---- f(x)=0.931818 err=0.605062
i= 4 , x=1.10 ---- f(x)=0.911111 err=0.582840
i= 5 , x=1.12 ---- f(x)=0.891304 err=0.561101
i= 6 , x=1.15 ---- f(x)=0.872340 err=0.539824
i= 7 , x=1.17 ---- f(x)=0.854167 err=0.518991
i= 8 , x=1.20 ---- f(x)=0.836735 err=0.498582
i= 9 , x=1.22 ---- f(x)=0.820000 err=0.478582
i=10 , x=1.24 ---- f(x)=0.803922 err=0.458975
i=11 , x=1.27 ---- f(x)=0.788462 err=0.439744
i=12 , x=1.29 ---- f(x)=0.773585 err=0.420876
i=13 , x=1.32 ---- f(x)=0.759259 err=0.402357
i=14 , x=1.34 ---- f(x)=0.745455 err=0.384176
i=15 , x=1.37 ---- f(x)=0.732143 err=0.366318
i=16 , x=1.39 ---- f(x)=0.719298 err=0.348775
i=17 , x=1.41 ---- f(x)=0.706897 err=0.331533
i=18 , x=1.44 ---- f(x)=0.694915 err=0.314584
i=19 , x=1.46 ---- f(x)=0.683333 err=0.297917
i=20 , x=1.49 ---- f(x)=0.672131 err=0.281524
i=21 , x=1.51 ---- f(x)=0.661290 err=0.265395
i=22 , x=1.54 ---- f(x)=0.650794 err=0.249522
i=23 , x=1.56 ---- f(x)=0.640625 err=0.233897
i=24 , x=1.59 ---- f(x)=0.630769 err=0.218512
i=25 , x=1.61 ---- f(x)=0.621212 err=0.203361
i=26 , x=1.63 ---- f(x)=0.611940 err=0.188435
i=27 , x=1.66 ---- f(x)=0.602941 err=0.173729
i=28 , x=1.68 ---- f(x)=0.594203 err=0.159237
i=29 , x=1.71 ---- f(x)=0.585714 err=0.144951
i=30 , x=1.73 ---- f(x)=0.577465 err=0.130867
i=31 , x=1.76 ---- f(x)=0.569444 err=0.116978
i=32 , x=1.78 ---- f(x)=0.561644 err=0.103279
i=33 , x=1.80 ---- f(x)=0.554054 err=0.089765
i=34 , x=1.83 ---- f(x)=0.546667 err=0.076432
i=35 , x=1.85 ---- f(x)=0.539474 err=0.063274
i=36 , x=1.88 ---- f(x)=0.532468 err=0.050287
i=37 , x=1.90 ---- f(x)=0.525641 err=0.037467
i=38 , x=1.93 ---- f(x)=0.518987 err=0.024809
i=39 , x=1.95 ---- f(x)=0.512500 err=0.012309
i=40 , x=1.98 ---- f(x)=0.506173 err=0.000037
梯形積分法 S41=0.693184 誤差=0.000037
沒有留言:
張貼留言