#!/usr/bin/env python
# A. Pletzer Tue Mar 20 11:42:05 EST 2001
# G. Genellina 2009-09-10: Minor syntax changes, compatibility with Python 2.4 and above
"""
Gauss Integration
"""
from itertools import izip
_nodes =(
(0.,),
(-0.5773502691896257,
0.5773502691896257,),
(-0.7745966692414834,
0.,
0.7745966692414834,),
(-0.861136311594053,
-0.3399810435848562,
0.3399810435848562,
0.861136311594053,),
(-0.906179845938664,
-0.5384693101056829,
0.,
0.5384693101056829,
0.906179845938664,),
(-0.932469514203152,
-0.6612093864662646,
-0.2386191860831968,
0.2386191860831968,
0.6612093864662646,
0.932469514203152,),
(-0.949107912342759,
-0.7415311855993937,
-0.4058451513773972,
0.,
0.4058451513773971,
0.7415311855993945,
0.949107912342759,),
(-0.960289856497537,
-0.7966664774136262,
-0.5255324099163289,
-0.1834346424956498,
0.1834346424956498,
0.5255324099163289,
0.7966664774136262,
0.960289856497537,),
(-0.968160239507626,
-0.836031107326637,
-0.6133714327005903,
-0.3242534234038088,
0.,
0.3242534234038088,
0.6133714327005908,
0.836031107326635,
0.968160239507627,),
(-0.973906528517172,
-0.865063366688984,
-0.6794095682990246,
-0.433395394129247,
-0.1488743389816312,
0.1488743389816312,
0.433395394129247,
0.6794095682990246,
0.865063366688984,
0.973906528517172,),
(-0.97822865814604,
-0.88706259976812,
-0.7301520055740422,
-0.5190961292068116,
-0.2695431559523449,
0.,
0.2695431559523449,
0.5190961292068117,
0.73015200557405,
0.887062599768093,
0.978228658146058,),
(-0.981560634246732,
-0.904117256370452,
-0.7699026741943177,
-0.5873179542866143,
-0.3678314989981804,
-0.1252334085114688,
0.1252334085114688,
0.3678314989981804,
0.5873179542866143,
0.7699026741943177,
0.904117256370452,
0.981560634246732,),
)
_weights=(
(2.,),
(1.,
1.,),
(0.5555555555555553,
0.888888888888889,
0.5555555555555553,),
(0.3478548451374539,
0.6521451548625462,
0.6521451548625462,
0.3478548451374539,),
(0.2369268850561887,
0.4786286704993665,
0.5688888888888889,
0.4786286704993665,
0.2369268850561887,),
(0.1713244923791709,
0.3607615730481379,
0.4679139345726913,
0.4679139345726913,
0.3607615730481379,
0.1713244923791709,),
(0.129484966168868,
0.2797053914892783,
0.3818300505051186,
0.4179591836734694,
0.3818300505051188,
0.279705391489276,
0.1294849661688697,),
(0.1012285362903738,
0.2223810344533786,
0.3137066458778874,
0.3626837833783619,
0.3626837833783619,
0.3137066458778874,
0.2223810344533786,
0.1012285362903738,),
(0.0812743883615759,
0.1806481606948543,
0.2606106964029356,
0.3123470770400029,
0.3302393550012597,
0.3123470770400025,
0.2606106964029353,
0.1806481606948577,
0.0812743883615721,),
(0.06667134430868681,
0.149451349150573,
0.2190863625159832,
0.2692667193099968,
0.2955242247147529,
0.2955242247147529,
0.2692667193099968,
0.2190863625159832,
0.149451349150573,
0.06667134430868681,),
(0.05566856711621584,
0.1255803694648743,
0.1862902109277404,
0.2331937645919927,
0.2628045445102466,
0.2729250867779006,
0.2628045445102466,
0.2331937645919933,
0.1862902109277339,
0.1255803694649132,
0.05566856711616958,),
(0.04717533638647547,
0.1069393259953637,
0.1600783285433586,
0.2031674267230672,
0.2334925365383534,
0.2491470458134027,
0.2491470458134027,
0.2334925365383534,
0.2031674267230672,
0.1600783285433586,
0.1069393259953637,
0.04717533638647547,),
)
_nodesLog =(
(0.3333333333333333,),
(0.1120088061669761,
0.6022769081187381,),
(0.06389079308732544,
0.3689970637156184,
0.766880303938942,),
(0.04144848019938324,
0.2452749143206022,
0.5561654535602751,
0.848982394532986,),
(0.02913447215197205,
0.1739772133208974,
0.4117025202849029,
0.6773141745828183,
0.89477136103101,),
(0.02163400584411693,
0.1295833911549506,
0.3140204499147661,
0.5386572173517997,
0.7569153373774084,
0.922668851372116,),
(0.0167193554082585,
0.100185677915675,
0.2462942462079286,
0.4334634932570557,
0.6323509880476823,
0.81111862674023,
0.940848166743287,),
(0.01332024416089244,
0.07975042901389491,
0.1978710293261864,
0.354153994351925,
0.5294585752348643,
0.7018145299391673,
0.849379320441094,
0.953326450056343,),
(0.01086933608417545,
0.06498366633800794,
0.1622293980238825,
0.2937499039716641,
0.4466318819056009,
0.6054816627755208,
0.7541101371585467,
0.877265828834263,
0.96225055941096,),
(0.00904263096219963,
0.05397126622250072,
0.1353118246392511,
0.2470524162871565,
0.3802125396092744,
0.5237923179723384,
0.6657752055148032,
0.7941904160147613,
0.898161091216429,
0.9688479887196,),
(0.007643941174637681,
0.04554182825657903,
0.1145222974551244,
0.2103785812270227,
0.3266955532217897,
0.4554532469286375,
0.5876483563573721,
0.7139638500230458,
0.825453217777127,
0.914193921640008,
0.973860256264123,),
(0.006548722279080035,
0.03894680956045022,
0.0981502631060046,
0.1811385815906331,
0.2832200676673157,
0.398434435164983,
0.5199526267791299,
0.6405109167754819,
0.7528650118926111,
0.850240024421055,
0.926749682988251,
0.977756129778486,),
)
_weightsLog=(
(-1.,),
(-0.7185393190303845,
-0.2814606809696154,),
(-0.5134045522323634,
-0.3919800412014877,
-0.0946154065661483,),
(-0.3834640681451353,
-0.3868753177747627,
-0.1904351269501432,
-0.03922548712995894,),
(-0.2978934717828955,
-0.3497762265132236,
-0.234488290044052,
-0.0989304595166356,
-0.01891155214319462,),
(-0.2387636625785478,
-0.3082865732739458,
-0.2453174265632108,
-0.1420087565664786,
-0.05545462232488041,
-0.01016895869293513,),
(-0.1961693894252476,
-0.2703026442472726,
-0.239681873007687,
-0.1657757748104267,
-0.0889432271377365,
-0.03319430435645653,
-0.005932787015162054,),
(-0.164416604728002,
-0.2375256100233057,
-0.2268419844319134,
-0.1757540790060772,
-0.1129240302467932,
-0.05787221071771947,
-0.02097907374214317,
-0.003686407104036044,),
(-0.1400684387481339,
-0.2097722052010308,
-0.211427149896601,
-0.1771562339380667,
-0.1277992280331758,
-0.07847890261203835,
-0.0390225049841783,
-0.01386729555074604,
-0.002408041036090773,),
(-0.12095513195457,
-0.1863635425640733,
-0.1956608732777627,
-0.1735771421828997,
-0.135695672995467,
-0.0936467585378491,
-0.05578772735275126,
-0.02715981089692378,
-0.00951518260454442,
-0.001638157633217673,),
(-0.1056522560990997,
-0.1665716806006314,
-0.1805632182877528,
-0.1672787367737502,
-0.1386970574017174,
-0.1038334333650771,
-0.06953669788988512,
-0.04054160079499477,
-0.01943540249522013,
-0.006737429326043388,
-0.001152486965101561,),
(-0.09319269144393,
-0.1497518275763289,
-0.166557454364573,
-0.1596335594369941,
-0.1384248318647479,
-0.1100165706360573,
-0.07996182177673273,
-0.0524069547809709,
-0.03007108900074863,
-0.01424924540252916,
-0.004899924710875609,
-0.000834029009809656,),
)
_NGMAX = len(_nodes)
_NGMIN = 1
def gauss(xmin, xmax, funct, ng=10):
"""
Gauss quadature (weight function = 1.0):
xmin, xmax: boundaries of integration domain
funct: integrand function
ng: Gauss integration order
"""
ng = max(min(ng, _NGMAX), _NGMIN)
ns = _nodes[ng-1]
ws = _weights[ng-1]
dx = xmax - xmin
xs = [(dx*y + xmin + xmax)/2. for y in ns]
return 0.5*dx*sum(funct(x)*w for x,w in izip(xs,ws))
def gaussLog(xmin, xmax, funct, ng=10):
"""
Gauss quadature with Log singularity at x=xmin:
xmin, xmax: boundaries of integration domain
funct: integrand function
ng: Gauss integration order
"""
ng = max(min(ng, _NGMAX), _NGMIN)
ns = _nodesLog[ng-1]
ws = _weightsLog[ng-1];
dx = xmax - xmin
xs = [(dx*y + xmin) for y in ns]
return dx*sum(funct(x)*w for x,w in izip(xs,ws))
####
if __name__ == '__main__':
from math import *
def f2(x):
return x**2
def f3(x):
return x**4
def f4(x):
return cos(2.*pi*(x-0.128726465))
def f5(x):
return 2.*cos(2.*pi*(x-0.128726465))**2
print '-'*80
print 'Gauss (weight function = 1)'
print '-'*80
# simple tests
print 'gauss(0., 1., f3, 1)=', gauss(0., 1., f3, 1)
print 'gauss(0., 1., f3, 2)=', gauss(0., 1., f3, 2)
print 'gauss(0., 1., f4, 3)=', gauss(0., 1., f3, 3)
print 'gauss(0., 1., f3 )=', gauss(0., 1., f3 )
# convergence test
ng = range(_NGMIN, _NGMAX+1)
print """\n
Integrate[Cos[2.*Pi*(x-0.128726465)], {x, 0, 1}]
\n"""
error = [gauss(0., 10.0, f4, n) for n in ng]
print ' n = ', '%8d'*len(ng) % tuple(ng)
print 'error = ', '%8.0e'*len(error) % tuple(error)
print """\n
Integrate[2.*Cos[2.*Pi*(x-0.128726465)]^2, {x, 0, 1}]
\n"""
error = [gauss(0., 1.0, f5, n)-1.0 for n in ng]
print ' n = ', '%8d'*len(ng) % tuple(ng)
print 'error = ', '%8.0e'*len(error) % tuple(error)
print '-'*80
print 'Gauss with Log singularity at left boundary'
print '-'*80
a, b = 0., 1.
print """\n
Integrate[Log[x]*x^2, {x, 0, 1}]
\n"""
exact = -1./9.
error = [gaussLog(a, b, f2, n) - exact for n in ng]
print ' n = ', '%8d'*len(ng) % tuple(ng)
print 'error = ', '%8.0e'*len(error) % tuple(error)
print """\n
Integrate[Log[x]*2.*Cos[2.*Pi*(x-0.128726465)]^2, {x, 0, 1}]
\n"""
exact = -1.242002481967963
error = [gaussLog(a, b, f5, n) - exact for n in ng]
print ' n = ', '%8d'*len(ng) % tuple(ng)
print 'error = ', '%8.0e'*len(error) % tuple(error)
