# 例題 4-7 利用 梯形 (trapezoidal Rule ) 計算 雙重積分 f(x,y)= x*exp(y) 在[0 , x] dy 與 [0,1] dx 的定積分

# 例題 4-7
# 利用 梯形 (trapezoidal Rule ) 計算 雙重積分
# f(x,y)= x*exp(y)  在[0 , x] dy  與 [0,1] dx 的定積分

n=2     T{2}={0.601507}
n=29   T{29}={0.500471}
真實值 S=0.5

#=========================================
# Based on Trapezoidal Rule to  compute the double integral.
import math
def F(x,y):
    return (x*math.exp(y) + 0.0 )
def C(x):
    return (0.0)
def D(x):
    return (x + 0.0)

def gy(n):
    sum=0.0;
    for i in range (0,n+1):
        for j in range (1,n):
            sum=sum+F(x[i],y[i][j])

        gg[i]=(hy[i]/2.0)*(F(x[i],y[i][0])+F(x[i],y[i][n])+2*sum)
        sum=0.0
     
#======================================
x=[0.0 for i in range(30)]
y=[ [0.0 for i in range(30)] for j in range (30)]
#print(y)
gg= [0.0 for i in range(30) ]
hy= [0.0 for i in range(30) ]

sum1=0.0
n=2
a=0.0
b=1.0
hx=(b-a)/n

for i in range (0 ,n+1) :
    x[i]=a+i*hx
    #print('x[%2d]={%.6f}' %(i,x[i]))
    hy[i]=(D(x[i])-C(x[i]))/n
    #print('hy[%2d]={%.6f}' %(i,hy[i]))
    for j in range (0,n+1):
        y[i][j]=C(x[i]) + j*hy[i] + 0.0
        #print('y[%2d][%2d]={%.6f}' %(i,j,y[i][j]),end='' )
    #print()
gy(n)

print('i={%2d} gg[%2d]={%.6f}' %(0,0,gg[0]))
for i in range (1,n):
    sum1=sum1+gg[i]
    print('i={%2d} gg[%2d]={%.6f}' %(i,i,gg[i]))

print('i={%2d} gg[%2d]={%.6f}' %(n,n,gg[n]))
ts=(hx/2.0)*(gg[0]+gg[n]+2*sum1)
print("T{%d}={%.6f}\n" %(n,ts))


#=============================
# n= 10 --> n=20
#=============================
sum1=0.0
n=29
a=0.0
b=1.0
hx=(b-a)/n

for i in range (0 ,n+1) :
    x[i]=a+i*hx
    #print('x[%2d]={%.6f}' %(i,x[i]))
    hy[i]=(D(x[i])-C(x[i]))/n
    #print('hy[%2d]={%.6f}' %(i,hy[i]))
    for j in range (0,n+1):
        y[i][j]=C(x[i]) + j*hy[i] + 0.0
        #print('y[%2d][%2d]={%.6f}' %(i,j,y[i][j]),end='' )
    #print()
gy(n)

print('i={%2d} gg[%2d]={%.6f}' %(0,0,gg[0]))
for i in range (1,n):
    sum1=sum1+gg[i]
    print('i={%2d} gg[%2d]={%.6f}' %(i,i,gg[i]))

print('i={%2d} gg[%2d]={%.6f}' %(n,n,gg[n]))
ts=(hx/2.0)*(gg[0]+gg[n]+2*sum1)
print("T{%d}={%.6f}\n" %(n,ts))

========= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX4-7.py ===========
i={ 0} gg[ 0]={0.000000}
i={ 1} gg[ 1]={0.326048}
i={ 2} gg[ 2]={1.753931}
T{2}={0.601507}

i={ 0} gg[ 0]={0.000000}
i={ 1} gg[ 1]={0.001210}
i={ 2} gg[ 2]={0.004924}
i={ 3} gg[ 3]={0.011275}
i={ 4} gg[ 4]={0.020400}
i={ 5} gg[ 5]={0.032443}
i={ 6} gg[ 6]={0.047556}
i={ 7} gg[ 7]={0.065898}
i={ 8} gg[ 8]={0.087633}
i={ 9} gg[ 9]={0.112934}
i={10} gg[10]={0.141983}
i={11} gg[11]={0.174969}
i={12} gg[12]={0.212091}
i={13} gg[13]={0.253555}
i={14} gg[14]={0.299577}
i={15} gg[15]={0.350385}
i={16} gg[16]={0.406216}
i={17} gg[17]={0.467316}
i={18} gg[18]={0.533944}
i={19} gg[19]={0.606371}
i={20} gg[20]={0.684879}
i={21} gg[21]={0.769764}
i={22} gg[22]={0.861333}
i={23} gg[23]={0.959909}
i={24} gg[24]={1.065829}
i={25} gg[25]={1.179445}
i={26} gg[26]={1.301122}
i={27} gg[27]={1.431246}
i={28} gg[28]={1.570217}
i={29} gg[29]={1.718452}
T{29}={0.500471}

>>>