Julia語言 範例EX4-8 梯行積分法 雙重積分計算

Julia語言 範例EX4-8 梯行積分法 雙重積分計算
n=10

           b d(x)

         a   c(x)

a=0  ,  b=1 , c(x)=0 , d(x)=1-x^2   f(x,y)=4xy


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

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

function C(x::Float64) #// c(x)=0
    return  (0.0)
end
function D(x::Float64) #// d(x)= 1-x^2
    return  (1-x*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]=(hy[i+1]/2.0)*(F(x[i+1],y[i+1][1])+F(x[i+1],y[i+1][n+1])+2*sum)
        sum=0.0;
    end
    return g1
end

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=10
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)
println("\n\n")
for i=1:n-1
    sum=sum+g1[i+1];
end 

ts=(hx/2.0)*(g1[1]+g1[n+1]+2*sum)
s=@sprintf("梯形積分計算雙重積分結果 T%d=%.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%.6lf\n",(1/3))
println(s)

s=@sprintf("誤差值=%.6lf\n",(abs(1/3)- ts))
println(s)


輸出畫面

梯形積分計算雙重積分結果 T10=0.331650

實際值=0.333333

誤差值=0.001683