Julia語言例題EX4-7-B利用梯形法求雙重積分 f (x,y) = x* exp(y) a=0.0 , b=1.0 , c(x) =0.0 d(x)= x 取 n=2 n-4 並比較誤差值

Julia語言例題EX4-7-B利用梯形法求雙重積分 f (x,y) = x* exp(y) a=0.0 , b=1.0 , c(x) =0.0 d(x)= x 取 n=2  n=4 並比較誤差值


#==================================================
/* ex4-7-B.jl based on Trapezoidal Rule to
 * compute the double integral.
 */
===================================================#
using Printf

function F(x::Float64, y::Float64) #//  f(x,y)= x exp(y)
    return  (x*exp(y))
end

function C(x::Float64) #//  c(x)= 0.0
    return  (0.0)
end

function D(x::Float64) #//  d(x)= x
    return  (x)
end


function gy(n::Int64)
    sum1=0.0
    sum2=0.0
    for i=0:n
        for j=1:n-1
            if(i%2==0)
                sum2=sum2+F(x[i+1],y[i+1][j+1])
            else
                sum1=sum1+F(x[i+1],y[i+1][j+1])
            end 
        end
        g1[i+1]=(hy[i+1]/3)*(F(x[i+1],y[i+1][1])+F(x[i+1],y[i+1][n+1])+2*sum2+4*sum1)
        println(i+1,"----",g1[i+1])
        sum1=0.0
        sum2=0.0
    end
    return g1
end


s=@sprintf("辛普森積分計算雙重積分 n=2 ")
println(s)

x= [0.0 for i=1:5 ]
g1=[0.0 for i=1:5 ]
g2=[0.0 for i=1:5 ]
hy=[0.0 for i=1:5 ]

f(i) = [0.0 for i=1:5]
y= f.([0.0 for i=1:5])

sum1=0.0
sum2=0.0

n=2
a=0.0
b=1.0

hx=(b-a)/n
ts=0.0
for i=0:n
    x[i+1]=a+i*hx
    hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
    for j=0:n
        y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
        #println(y)
    end
end

#println(x)
#println(hy)
#println(y)

g2=gy(n)
sum1=0.0
sum2=0.0

println("\n\n")
for i=1:n-1
    if(i%2==0)
        sum2=sum2+g2[i+1]
    else
        sum1=sum1+g2[i+1]
    end
 
end 

ts= (hx/3) * (g1[1] + g1[n+1] + 4*sum1 + 2*sum2)
s=@sprintf("辛普森積分計算雙重積分結果 T%d=%0.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%0.6lf\n",(1/2))
println(s)
tn=abs( 1/2 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)


s=@sprintf("梯形積分計算雙重積分 n=4 ")
println(s)

x= [0.0 for i=1:10 ]
g1=[0.0 for i=1:10 ]
g2=[0.0 for i=1:10 ]
hy=[0.0 for i=1:10 ]

f(i) = [0.0 for i=1:10]
y= f.([0.0 for i=1:10])


sum1=0.0
sum2=0.0

n=4
a=0.0
b=1.0

hx=(b-a)/n
ts=0.0
for i=0:n
    x[i+1]=a+i*hx
    hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
    for j=0:n
        y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
        #println(y)
    end
end

#println(x)
#println(hy)
#println(y)

g2=gy(n)
sum1=0.0
sum2=0.0

println("\n\n")
for i=1:n-1
    if(i%2==0)
        sum2=sum2+g2[i+1]
    else
        sum1=sum1+g2[i+1]
    end
 
end 

ts= (hx/3) * (g1[1] + g1[n+1] + 4*sum1 + 2*sum2)
s=@sprintf("辛普森積分計算雙重積分結果 T%d=%0.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%0.6lf\n",(1/2))
println(s)
tn=abs( 1/2 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)

輸出畫面

辛普森積分計算雙重積分 n=2 
1----0.0
2----0.3243676223937956
3----1.16928739497655



辛普森積分計算雙重積分結果 T2=0.411126

實際值=0.500000

誤差值=0.088874

梯形積分計算雙重積分 n=4 
1----0.0
2----0.0828099899078667
3----0.21652191332178472
4----0.9741611859661217
5----1.151481269705011



辛普森積分計算雙重積分結果 T4=0.484367

實際值=0.500000

誤差值=0.015633