例題4-6 # 利用 梯形與 辛普森 理則計算 雙重積分 f(x,y)= x^2 y 在[1,2] dx 與 [0,1] dy 的定積分

例題4-6
# 利用 梯形 (trapezoidal Rule )與  辛普森 理則 （Simpson's Rule)
#計算 雙重積分   f(x,y)= x^2 y  在[1,2] dx  與 [0,1] dy 的定積分

from math import *

def F(x,y):
    temp= x*x*y+0.0
    return temp
#===============================================
#在[1,2] dx  與 [0,1] dy 的定積分
#===============================================
a1=1.0
b1=2.0
n=2
hx=(b1-a1)/n
x=[0.0 for i in range(10)]
x[0]=a1
for i in range (1,n+1):
    x[i]=(x[i-1]+hx)
#print(x)
#===============================================
#在[1,2] dx  與 [0,1] dy 的定積分
#=============================================== 
a2=0.0
b2=1.0
n=2
hy=(b2-a2)/n
y=[0.0 for i in range(10)]
y[0]=a2
for i in range (1,n+1):
    y[i]=(y[i-1]+hy)
#print(y)

print("梯形理則 (trapezoidal Rule )")
Gx0=0.5*hy*(   F(x[0],y[0])+2* F(x[0],y[1])+ F(x[0],y[2]) )
print ('Gx0=',Gx0,'\t\t',end='')
Gx1=0.5*hy*(   F(x[1],y[0])+2* F(x[1],y[1])+ F(x[1],y[2]) )
print ('Gx1=',Gx1,'\t',end='')
Gx2=0.5*hy*(   F(x[2],y[0])+2* F(x[2],y[1])+ F(x[2],y[2]) )
print ('Gx2=',Gx2)
Tn=0.5*hy*(Gx0+2*Gx1+Gx2)
print('0.5*hy*(Gx0+2*Gx1+Gx2)=',end='')
print ('T[%2d] = {%8.6f}' %(n,Tn))
print('')
#===============================================
print("辛普森理則（Simpson's Rule)")
Gx0=(1/3)*hy*( F(x[0],y[0])+4* F(x[0],y[1])+ F(x[0],y[2]) )
print ('Gx0=',Gx0,'\t\t',end='')
Gx1=(1/3)*hy*( F(x[1],y[0])+4* F(x[1],y[1])+ F(x[1],y[2]) )
print ('Gx1=',Gx1,'\t',end='')
Gx2=(1/3)*hy*( F(x[2],y[0])+4* F(x[2],y[1])+ F(x[2],y[2]) )
print ('Gx2=',Gx2)
Tn=(1/3)*hy*(Gx0+4*Gx1+Gx2)
print('(1/3)*hy*(Gx0+4*Gx1+Gx2)=',end='')
print ('T[%2d] = {%8.6f}' %(n,Tn))
print('')
#===============================================
a1=1.0
b1=2.0
n=10
hx=(b1-a1)/n
x=[0.0 for i in range(15)]
x[0]=a1
for i in range (1,n+1):
    x[i]=(x[i-1]+hx)
#print(hx)

a2=0.0
b2=1.0
n=10 
hy=(b2-a2)/n
y=[0.0 for i in range(15)]
y[0]=a2
for i in range (1,n+1):
    y[i]=(y[i-1]+hy)
#print(hy,x,y)

print("辛普森理則（Simpson's Rule) ,  n=10")
Gx0=(hx/3)*( F(x[0],y[0])+4* F(x[0],y[1])+ 2*F(x[0],y[2]) +4* F(x[0],y[3])+ 2*F(x[0],y[4]) + 4* F(x[0],y[5])+ 2*F(x[0],y[6]) + 4* F(x[0],y[7])+ 2*F(x[0],y[8])  +  4*F(x[0],y[9])  + F(x[0],y[10] ))
print ('Gx0=',Gx0)
Gx1=(hx/3)*( F(x[1],y[0])+4* F(x[1],y[1])+ 2*F(x[1],y[2]) +4* F(x[1],y[3])+ 2*F(x[1],y[4])  + 4* F(x[1],y[5])+ 2*F(x[1],y[6]) + 4* F(x[1],y[7])+ 2*F(x[1],y[8])  +   4*F(x[1],y[9])  + F(x[2],y[10] ))
print ('Gx1=',Gx1)
Gx2=(hx/3)*( F(x[2],y[0])+4* F(x[2],y[1])+ 2*F(x[2],y[2]) +4* F(x[2],y[3])+ 2*F(x[2],y[4]) + 4* F(x[2],y[5])+ 2*F(x[2],y[6]) + 4* F(x[2],y[7])+ 2*F(x[2],y[8])  +   4*F(x[2],y[9] ) + F(x[2],y[10] ))
print ('Gx2=',Gx2)
Gx3=(hx/3)*( F(x[3],y[0])+4* F(x[3],y[1])+ 2*F(x[3],y[2]) +4* F(x[3],y[3])+ 2*F(x[3],y[4])  + 4* F(x[3],y[5])+ 2*F(x[3],y[6]) + 4* F(x[3],y[7])+ 2*F(x[3],y[8])  +   4*F(x[3],y[9] ) + F(x[3],y[10]))
print ('Gx3=',Gx3)
Gx4=(hx/3)*( F(x[4],y[0])+4* F(x[4],y[1])+ 2*F(x[4],y[2]) +4* F(x[4],y[3])+ 2*F(x[4],y[4])  + 4* F(x[4],y[5])+ 2*F(x[4],y[6]) + 4* F(x[4],y[7])+ 2*F(x[4],y[8])  +   4*F(x[4],y[9]  )+ F(x[4],y[10] ))
print ('Gx4=',Gx4)
Gx5=(hx/3)*( F(x[5],y[0])+4* F(x[5],y[1])+ 2*F(x[5],y[2]) +4* F(x[5],y[3])+ 2*F(x[5],y[4])  + 4* F(x[5],y[5])+ 2*F(x[5],y[6]) + 4* F(x[5],y[7])+ 2*F(x[5],y[8])  +   4*F(x[4],y[9] ) + F(x[5],y[10] ))
print ('Gx5=',Gx5)
Gx6=(hx/3)*( F(x[6],y[0])+4* F(x[6],y[1])+ 2*F(x[6],y[2]) +4* F(x[6],y[3])+ 2*F(x[6],y[4])  + 4* F(x[6],y[5])+ 2*F(x[6],y[6]) + 4* F(x[6],y[7])+ 2*F(x[6],y[8])  +   4*F(x[6],y[9] ) + F(x[6],y[10] ))
print ('Gx6=',Gx6)
Gx7=(hx/3)*( F(x[7],y[0])+4* F(x[7],y[1])+ 2*F(x[7],y[2]) +4* F(x[7],y[3])+ 2*F(x[7],y[4])  + 4* F(x[7],y[5])+ 2*F(x[7],y[6]) + 4* F(x[7],y[7])+ 2*F(x[7],y[8])  +   4*F(x[7],y[9] ) + F(x[7],y[10] ))
print ('Gx7=',Gx7)
Gx8=(hx/3)*( F(x[8],y[0])+4* F(x[8],y[1])+ 2*F(x[8],y[2]) +4* F(x[8],y[3])+ 2*F(x[8],y[4])  + 4* F(x[8],y[5])+ 2*F(x[8],y[6]) + 4* F(x[8],y[7])+ 2*F(x[8],y[8])  +   4*F(x[8],y[9] ) +F(x[8],y[10] ))
print ('Gx8=',Gx8)
Gx9=(hx/3)*( F(x[9],y[0])+4* F(x[9],y[1])+ 2*F(x[9],y[2]) +4* F(x[9],y[3])+ 2*F(x[9],y[4])  + 4* F(x[9],y[5])+ 2*F(x[9],y[6]) + 4* F(x[9],y[7])+ 2*F(x[9],y[8])  +   4*F(x[9],y[9] ) + F(x[9],y[10] ))
print ('Gx9=',Gx9)
Gx10=(hx/3)*( F(x[10],y[0])+4* F(x[10],y[1])+ 2*F(x[10],y[2]) +4* F(x[10],y[3])+ 2*F(x[10],y[4])  + 4* F(x[10],y[5])+ 2*F(x[10],y[6]) + 4* F(x[10],y[7])+ 2*F(x[10],y[8])  +   4*F(x[10],y[9] ) + F(x[10],y[10] ))
print ('Gx10=',Gx10)

Tn=(hy/3)*(Gx0+4*Gx1+2*Gx2 +4*Gx3+2*Gx4 + 4*Gx5+2*Gx6+ 4*Gx7+2*Gx8 + 4*Gx9 + Gx10   )
print('(hy/3)*(Gx0+4*Gx1+2*Gx2 +4*Gx3+2Gx4 + 4*Gx5+2*Gx6+ 4*Gx7+2*Gx8 + 4*Gx9  +Gx10)=',end='')
print ('T[%2d] = {%8.6f}' %(n,Tn))


======= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX4-6-1.py ===========
梯形理則 (trapezoidal Rule )
Gx0= 0.5 Gx1= 1.125 Gx2= 2.0
0.5*hy*(Gx0+2*Gx1+Gx2)=T[ 2] = {1.187500}

辛普森理則（Simpson's Rule)
Gx0= 0.5 Gx1= 1.125 Gx2= 2.0
(1/3)*hy*(Gx0+4*Gx1+Gx2)=T[ 2] = {1.166667}

辛普森理則（Simpson's Rule) ,  n=10
Gx0= 0.5
Gx1= 0.6126666666666667
Gx2= 0.7200000000000002
Gx3= 0.8450000000000004
Gx4= 0.9800000000000005
Gx5= 1.0902000000000005
Gx6= 1.2800000000000011
Gx7= 1.4450000000000007
Gx8= 1.6200000000000012
Gx9= 1.8050000000000013
Gx10= 2.0000000000000018
(hy/3)*(Gx0+4*Gx1+2*Gx2 +4*Gx3+2Gx4 + 4*Gx5+2*Gx6+ 4*Gx7+2*Gx8 + 4*Gx9  +Gx10)=T[10] = {1.163049}
>>>