Julia語言 例題4-6 雙重定積分 a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=2 求積分值

Julia語言 例題4-6 雙重定積分

           b  d

        a   c
a=0.0  , b=1.0  , c=1.0 ,d=2.0  , f(x,y)= x^2 y 取n=2 求積分值

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

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


x=[0.0 for i=1:20]
g1=[0.0 for i=1:20 ]
y=[0.0 for i=1:20 ]

sum=0.0
n=2
a=1.0
b=2.0
hx=(b-a)/n
ax=a
for i=1:n+1
    x[i]=ax
    ax=ax+hx
end

c=0.0
d=1.0
hy=(d-c)/n
ay=c
for i=1:n+1
    y[i]=ay
    ay+=hy
end

for i=1:n+1
    s=@sprintf("i=%2d , x[%1d]=%0.3f , y[%1d]=%0.3f ",i,i,x[i],i,y[i])
    println(s)
end 

println("\n梯形理則 (trapezoidal Rule) ")
Gx0=0.5*hy*(F(x[1],y[1]) + 2*F(x[1],y[2]) +F(x[1],y[3]))
s=@sprintf("Gx0=%0.5f",Gx0)
println(s)

Gx1=0.5*hy*(F(x[2],y[1])+2* F(x[2],y[2])+ F(x[2],y[3]) )
s=@sprintf("Gx1=%0.5f",Gx1)
println(s)

Gx2=0.5*hy*(F(x[3],y[1])+2* F(x[3],y[2])+ F(x[3],y[3]) )
s=@sprintf("Gx2=%0.5f",Gx2)
println(s)

Tn=0.5*hx*(Gx0+2*Gx1+Gx2)
s=@sprintf("0.5*hy*(Gx0+2*Gx1+Gx2)= %0.5f",Tn)
println(s)

#===============================================#
println("\n辛普森理則（Simpson's Rule)")
Gx0=(1/3)*hy*( F(x[1],y[1])+4* F(x[1],y[2])+ F(x[1],y[3]) )
s=@sprintf("Gx0=%0.5f",Gx0)
println(s)

Gx1=(1/3)*hy*( F(x[2],y[1])+4* F(x[2],y[2])+ F(x[2],y[3]) )
s=@sprintf("Gx1=%0.5f",Gx1)
println(s)

Gx2=(1/3)*hy*( F(x[3],y[1])+4* F(x[3],y[2])+ F(x[3],y[3]) )
s=@sprintf("Gx2=%0.5f",Gx2)
println(s)
Tn=(1/3)*hy*(Gx0+4*Gx1+Gx2)
s=@sprintf("(1/3)*hy*(Gx0+4*Gx1+Gx2)= %0.5f", Tn)
println(s)
#===============================================#

輸出畫面

i= 1 , x[1]=1.000 , y[1]=0.000 
i= 2 , x[2]=1.500 , y[2]=0.500 
i= 3 , x[3]=2.000 , y[3]=1.000 

梯形理則 (trapezoidal Rule) 
Gx0=0.50000
Gx1=1.12500
Gx2=2.00000
0.5*hy*(Gx0+2*Gx1+Gx2)= 1.18750

辛普森理則（Simpson's Rule)
Gx0=0.50000
Gx1=1.12500
Gx2=2.00000
(1/3)*hy*(Gx0+4*Gx1+Gx2)= 1.16667