Julia語言例題EX4-7-A利用梯形法求雙重積分 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-A利用梯形法求雙重積分
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-A.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)
    sum=0.0;
    for i=0:n
        for j=1:n-1
          sum=sum+F(x[i+1],y[i+1][j+1])
          #println(y[i][j])
        end
        g1[i+1]=(0.5*hy[i+1])*(F(x[i+1],y[i+1][1])+F(x[i+1],y[i+1][n+1])+2*sum)
        println(i+1,"----",g1[i+1])
        sum=0.0;
    end
    return g1
end


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

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

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

sum=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


g2=gy(n)
sum1=0.0
println("\n\n")
for i=1:n-1
    sum1=sum1+g2[i+1];
    #println(g2[i+1])
end 

ts= (hx/2) * (g1[1] + g1[n+1] + 2*sum1)
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:20 ]
g1=[0.0 for i=1:20 ]
g2=[0.0 for i=1:20 ]
hy=[0.0 for i=1:20 ]

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

sum=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


g2=gy(n)
sum1=0.0
println("\n\n")
for i=1:n-1
    sum1=sum1+g2[i+1];
    #println(g2[i+1])
end 

ts= (hx/2) * (g1[1] + g1[n+1] + 2*sum1)
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.3260482565047257
3----1.7539310924648253



梯形積分計算雙重積分結果 T2=0.601507

實際值=0.500000

誤差值=0.101507

梯形積分計算雙重積分 n=4 
1----0.0
2----0.07102946671483652
3----0.3247828699826771
4----0.8402029213086307
5----1.7272219045575166



梯形積分計算雙重積分結果 T4=0.524907

實際值=0.500000

誤差值=0.024907