Julia語言 例題4-4使用 求 f(x)= (1/x) a=1.0 , b=2.0 求n=?? 辛普森積分法 , 梯形積分法需要幾次

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