# 例題 4-7 用辛普森理則 計算 雙重積分 # f(x,y)= x*exp(y) 在[0 , x] dy 與 [0,1] dx 的定積分

# 例題 4-7
# 利用 辛普森理則 (Simpson's Rule ) 計算 雙重積分
# f(x,y)= x*exp(y)  在[0 , x] dy  與 [0,1] dx 的定積分
# Based on Simpson's  Rule to  compute the double integral.
#======================================
T{2}={0.574057}
T{4}={0.854994}
T{20}={3.479674}
真實值 =0.5
?????

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):
    sum1=0.0
    sum2=0.0
    for i in range (0,n+1):
        for j in range (1,n):
            if(j%2==0):
                sum2=sum2+F(x[i],y[i][j])
            else:
                sum1=sum1+F(x[i],y[i][j])
        gg[i]=(hy[i]/3.0)*(F(x[i],y[i][0])+F(x[i],y[i][n])+2.0*sum2+4.0*sum1)
        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
    hy[i]=(D(x[i])-C(x[i]))/n
    for j in range (0,n+1):
        y[i][j]=C(x[i]) + j*hy[i] + 0.0

gy(n)
sum1=0.0
sum2=0.0
print('i={%2d} gg[%2d]={%.6f}' %(0,0,gg[0]))
for i in range (1,n):
    if(i%2==0):
        sum2=sum2+gg[i]
    else:
        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/3.0)*(gg[0]+gg[n]+2.0*sum2+4.0*sum1)
print("T{%d}={%.6f}\n" %(n,ts))

#=============================
# n= 2 ---> 10 --> n=20
#=============================
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=4
a=0.0
b=1.0
hx=(b-a)/n

for i in range (0 ,n+1) :
    x[i]=a+i*hx
    hy[i]=(D(x[i])-C(x[i]))/n
    for j in range (0,n+1):
        y[i][j]=C(x[i]) + j*hy[i] + 0.0

gy(n)
sum1=0.0
sum2=0.0
print('i={%2d} gg[%2d]={%.6f}' %(0,0,gg[0]))
for i in range (1,n):
    if(i%2==0):
        sum2=sum2+gg[i]
    else:
        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/3.0)*(gg[0]+gg[n]+2.0*sum2+4.0*sum1)
print("T{%d}={%.6f}\n" %(n,ts))

#=============================
# n= 10 --> n=20
#=============================
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=20
a=0.0
b=1.0
hx=(b-a)/n

for i in range (0 ,n+1) :
    x[i]=a+i*hx
    hy[i]=(D(x[i])-C(x[i]))/n
    for j in range (0,n+1):
        y[i][j]=C(x[i]) + j*hy[i] + 0.0

gy(n)
sum1=0.0
sum2=0.0
print('i={%2d} gg[%2d]={%.6f}' %(0,0,gg[0]))
for i in range (1,n):
    if(i%2==0):
        sum2=sum2+gg[i]
    else:
        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/3.0)*(gg[0]+gg[n]+2.0*sum2+4.0*sum1)
print("T{%d}={%.6f}\n" %(n,ts))


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

i={ 0} gg[ 0]={0.000000}
i={ 1} gg[ 1]={0.071006}
i={ 2} gg[ 2]={0.442582}
i={ 3} gg[ 3]={1.418856}
i={ 4} gg[ 4]={3.415314}
T{4}={0.854994}

i={ 0} gg[ 0]={0.000000}
i={ 1} gg[ 1]={0.002564}
i={ 2} gg[ 2]={0.015473}
i={ 3} gg[ 3]={0.046959}
i={ 4} gg[ 4]={0.105811}
i={ 5} gg[ 5]={0.201419}
i={ 6} gg[ 6]={0.343806}
i={ 7} gg[ 7]={0.543668}
i={ 8} gg[ 8]={0.812420}
i={ 9} gg[ 9]={1.162238}
i={10} gg[10]={1.606105}
i={11} gg[11]={2.157865}
i={12} gg[12]={2.832273}
i={13} gg[13]={3.645053}
i={14} gg[14]={4.612959}
i={15} gg[15]={5.753836}
i={16} gg[16]={7.086688}
i={17} gg[17]={8.631749}
i={18} gg[18]={10.410562}
i={19} gg[19]={12.446051}
i={20} gg[20]={14.762615}
T{20}={3.479674}

>>>