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

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

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

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

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

function D(x::Float64) #//  d(x)= 1.0
    return  (1.0)
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("梯形積分計算雙重積分")
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=10
a=1.0
b=2.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",(7/6))
println(s)
tn=abs( 7/6 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)


輸出畫面

梯形積分計算雙重積分
1----0.5000000000000001
2----0.6050000000000002
3----0.7200000000000001
4----0.8450000000000002
5----0.98
6----1.125
7----1.2800000000000005
8----1.4450000000000005
9----1.6200000000000003
10----1.8050000000000002
11----2.0000000000000004



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

實際值=1.166667

誤差值=0.000833