# geotecha - A software suite for geotechncial engineering
# Copyright (C) 2013 Rohan T. Walker (rtrwalker@gmail.com)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see http://www.gnu.org/licenses/gpl.html.
"""Numerical integration by quadrature"""
from __future__ import division, print_function
import matplotlib.pyplot
import numpy as np
from scipy import integrate
from scipy.special import jn_zeros
from scipy.special import jn
from matplotlib import pyplot as plt
import functools
import unittest
from numpy.testing import assert_allclose
from numpy.polynomial.polynomial import Polynomial
[docs]def gauss_kronrod_abscissae_and_weights(n):
"""Gauss-Kronrod quadrature abscissae and weights
Coarse integral = Sum(f(xi) * wi1)
Fine integral = Sum(f(xi) * wi2)
For the coarse integral the unused weights are set to zero
Parameters
----------
n : [2-20, 32, 64, 100]
number of integration points for the Gauss points. Number of Kronrod
points will automatically be 2 * n + 1.
Returns
-------
xi : 1d array
Abscissae for the quadrature points.
wi1 : 1d array
Weights for the coarse integral.
wi2 : 1d array
Weights for the fine integral
References
----------
.. [2] Holoborodko, Pavel. 2011. 'Gauss-Kronrod Quadrature Nodes and
Weights. November 7.
http://www.advanpix.com/2011/11/07/gauss-kronrod-quadrature-nodes-weights/#Tabulated_Gauss-Kronrod_weights_and_abscissae
"""
if n not in [7,10,15,20,25,30]:
raise ValueError('n must be 2-20, 32, 64, or 100')
weights = {
7: {
'g': np.array(
[[-0.9491079123427585245261897, 0.1294849661688696932706114],
[ -0.7415311855993944398638648, 0.2797053914892766679014678],
[ -0.4058451513773971669066064, 0.3818300505051189449503698],
[ 0.0000000000000000000000000, 0.4179591836734693877551020],
[ 0.4058451513773971669066064, 0.3818300505051189449503698],
[ 0.7415311855993944398638648, 0.2797053914892766679014678],
[ 0.9491079123427585245261897, 0.1294849661688696932706114]],
dtype=float),
'k': np.array(
[[-0.9914553711208126392068547, 0.0229353220105292249637320],
[ -0.9491079123427585245261897, 0.0630920926299785532907007],
[ -0.8648644233597690727897128, 0.1047900103222501838398763],
[ -0.7415311855993944398638648, 0.1406532597155259187451896],
[ -0.5860872354676911302941448, 0.1690047266392679028265834],
[ -0.4058451513773971669066064, 0.1903505780647854099132564],
[ -0.2077849550078984676006894, 0.2044329400752988924141620],
[ 0.0000000000000000000000000, 0.2094821410847278280129992],
[ 0.2077849550078984676006894, 0.2044329400752988924141620],
[ 0.4058451513773971669066064, 0.1903505780647854099132564],
[ 0.5860872354676911302941448, 0.1690047266392679028265834],
[ 0.7415311855993944398638648, 0.1406532597155259187451896],
[ 0.8648644233597690727897128, 0.1047900103222501838398763],
[ 0.9491079123427585245261897, 0.0630920926299785532907007],
[ 0.9914553711208126392068547, 0.0229353220105292249637320]],
dtype=float),
'dup': np.array(
[False, True, False, True, False, True, False, True, False,
True, False, True, False, True, False], dtype=bool)
},
10: {
'g': np.array(
[[-0.9739065285171717200779640, 0.0666713443086881375935688],
[ -0.8650633666889845107320967, 0.1494513491505805931457763],
[ -0.6794095682990244062343274, 0.2190863625159820439955349],
[ -0.4333953941292471907992659, 0.2692667193099963550912269],
[ -0.1488743389816312108848260, 0.2955242247147528701738930],
[ 0.1488743389816312108848260, 0.2955242247147528701738930],
[ 0.4333953941292471907992659, 0.2692667193099963550912269],
[ 0.6794095682990244062343274, 0.2190863625159820439955349],
[ 0.8650633666889845107320967, 0.1494513491505805931457763],
[ 0.9739065285171717200779640, 0.0666713443086881375935688]],
dtype=float),
'k': np.array(
[[-0.9956571630258080807355273, 0.0116946388673718742780644],
[ -0.9739065285171717200779640, 0.0325581623079647274788190],
[ -0.9301574913557082260012072, 0.0547558965743519960313813],
[ -0.8650633666889845107320967, 0.0750396748109199527670431],
[ -0.7808177265864168970637176, 0.0931254545836976055350655],
[ -0.6794095682990244062343274, 0.1093871588022976418992106],
[ -0.5627571346686046833390001, 0.1234919762620658510779581],
[ -0.4333953941292471907992659, 0.1347092173114733259280540],
[ -0.2943928627014601981311266, 0.1427759385770600807970943],
[ -0.1488743389816312108848260, 0.1477391049013384913748415],
[ 0.0000000000000000000000000, 0.1494455540029169056649365],
[ 0.1488743389816312108848260, 0.1477391049013384913748415],
[ 0.2943928627014601981311266, 0.1427759385770600807970943],
[ 0.4333953941292471907992659, 0.1347092173114733259280540],
[ 0.5627571346686046833390001, 0.1234919762620658510779581],
[ 0.6794095682990244062343274, 0.1093871588022976418992106],
[ 0.7808177265864168970637176, 0.0931254545836976055350655],
[ 0.8650633666889845107320967, 0.0750396748109199527670431],
[ 0.9301574913557082260012072, 0.0547558965743519960313813],
[ 0.9739065285171717200779640, 0.0325581623079647274788190],
[ 0.9956571630258080807355273, 0.0116946388673718742780644]],
dtype=float),
'dup': np.array(
[False, True, False, True, False, True, False, True, False,
True, False, True, False, True, False, True, False, True,
False, True, False], dtype=bool)
},
15: {
'g': np.array(
[[-0.9879925180204854284895657, 0.0307532419961172683546284],
[ -0.9372733924007059043077589, 0.0703660474881081247092674],
[ -0.8482065834104272162006483, 0.1071592204671719350118695],
[ -0.7244177313601700474161861, 0.1395706779261543144478048],
[ -0.5709721726085388475372267, 0.1662692058169939335532009],
[ -0.3941513470775633698972074, 0.1861610000155622110268006],
[ -0.2011940939974345223006283, 0.1984314853271115764561183],
[ 0.0000000000000000000000000, 0.2025782419255612728806202],
[ 0.2011940939974345223006283, 0.1984314853271115764561183],
[ 0.3941513470775633698972074, 0.1861610000155622110268006],
[ 0.5709721726085388475372267, 0.1662692058169939335532009],
[ 0.7244177313601700474161861, 0.1395706779261543144478048],
[ 0.8482065834104272162006483, 0.1071592204671719350118695],
[ 0.9372733924007059043077589, 0.0703660474881081247092674],
[ 0.9879925180204854284895657, 0.0307532419961172683546284]],
dtype=float),
'k': np.array(
[[-0.9980022986933970602851728, 0.0053774798729233489877921],
[ -0.9879925180204854284895657, 0.0150079473293161225383748],
[ -0.9677390756791391342573480, 0.0254608473267153201868740],
[ -0.9372733924007059043077589, 0.0353463607913758462220379],
[ -0.8972645323440819008825097, 0.0445897513247648766082273],
[ -0.8482065834104272162006483, 0.0534815246909280872653431],
[ -0.7904185014424659329676493, 0.0620095678006706402851392],
[ -0.7244177313601700474161861, 0.0698541213187282587095201],
[ -0.6509967412974169705337359, 0.0768496807577203788944328],
[ -0.5709721726085388475372267, 0.0830805028231330210382892],
[ -0.4850818636402396806936557, 0.0885644430562117706472754],
[ -0.3941513470775633698972074, 0.0931265981708253212254869],
[ -0.2991800071531688121667800, 0.0966427269836236785051799],
[ -0.2011940939974345223006283, 0.0991735987217919593323932],
[ -0.1011420669187174990270742, 0.1007698455238755950449467],
[ 0.0000000000000000000000000, 0.1013300070147915490173748],
[ 0.1011420669187174990270742, 0.1007698455238755950449467],
[ 0.2011940939974345223006283, 0.0991735987217919593323932],
[ 0.2991800071531688121667800, 0.0966427269836236785051799],
[ 0.3941513470775633698972074, 0.0931265981708253212254869],
[ 0.4850818636402396806936557, 0.0885644430562117706472754],
[ 0.5709721726085388475372267, 0.0830805028231330210382892],
[ 0.6509967412974169705337359, 0.0768496807577203788944328],
[ 0.7244177313601700474161861, 0.0698541213187282587095201],
[ 0.7904185014424659329676493, 0.0620095678006706402851392],
[ 0.8482065834104272162006483, 0.0534815246909280872653431],
[ 0.8972645323440819008825097, 0.0445897513247648766082273],
[ 0.9372733924007059043077589, 0.0353463607913758462220379],
[ 0.9677390756791391342573480, 0.0254608473267153201868740],
[ 0.9879925180204854284895657, 0.0150079473293161225383748],
[ 0.9980022986933970602851728, 0.0053774798729233489877921]],
dtype=float),
'dup': np.array(
[False, True, False, True, False, True, False, True, False,
True, False, True, False, True, False, True, False, True,
False, True, False, True, False, True, False, True, False,
True, False, True, False], dtype=bool)
},
20: {
'g': np.array(
[[-0.9931285991850949247861224, 0.0176140071391521183118620],
[ -0.9639719272779137912676661, 0.0406014298003869413310400],
[ -0.9122344282513259058677524, 0.0626720483341090635695065],
[ -0.8391169718222188233945291, 0.0832767415767047487247581],
[ -0.7463319064601507926143051, 0.1019301198172404350367501],
[ -0.6360536807265150254528367, 0.1181945319615184173123774],
[ -0.5108670019508270980043641, 0.1316886384491766268984945],
[ -0.3737060887154195606725482, 0.1420961093183820513292983],
[ -0.2277858511416450780804962, 0.1491729864726037467878287],
[ -0.0765265211334973337546404, 0.1527533871307258506980843],
[ 0.0765265211334973337546404, 0.1527533871307258506980843],
[ 0.2277858511416450780804962, 0.1491729864726037467878287],
[ 0.3737060887154195606725482, 0.1420961093183820513292983],
[ 0.5108670019508270980043641, 0.1316886384491766268984945],
[ 0.6360536807265150254528367, 0.1181945319615184173123774],
[ 0.7463319064601507926143051, 0.1019301198172404350367501],
[ 0.8391169718222188233945291, 0.0832767415767047487247581],
[ 0.9122344282513259058677524, 0.0626720483341090635695065],
[ 0.9639719272779137912676661, 0.0406014298003869413310400],
[ 0.9931285991850949247861224, 0.0176140071391521183118620]],
dtype=float),
'k': np.array(
[[-0.9988590315882776638383156, 0.0030735837185205315012183],
[ -0.9931285991850949247861224, 0.0086002698556429421986618],
[ -0.9815078774502502591933430, 0.0146261692569712529837880],
[ -0.9639719272779137912676661, 0.0203883734612665235980102],
[ -0.9408226338317547535199827, 0.0258821336049511588345051],
[ -0.9122344282513259058677524, 0.0312873067770327989585431],
[ -0.8782768112522819760774430, 0.0366001697582007980305572],
[ -0.8391169718222188233945291, 0.0416688733279736862637883],
[ -0.7950414288375511983506388, 0.0464348218674976747202319],
[ -0.7463319064601507926143051, 0.0509445739237286919327077],
[ -0.6932376563347513848054907, 0.0551951053482859947448324],
[ -0.6360536807265150254528367, 0.0591114008806395723749672],
[ -0.5751404468197103153429460, 0.0626532375547811680258701],
[ -0.5108670019508270980043641, 0.0658345971336184221115636],
[ -0.4435931752387251031999922, 0.0686486729285216193456234],
[ -0.3737060887154195606725482, 0.0710544235534440683057904],
[ -0.3016278681149130043205554, 0.0730306903327866674951894],
[ -0.2277858511416450780804962, 0.0745828754004991889865814],
[ -0.1526054652409226755052202, 0.0757044976845566746595428],
[ -0.0765265211334973337546404, 0.0763778676720807367055028],
[ 0.0000000000000000000000000, 0.0766007119179996564450499],
[ 0.0765265211334973337546404, 0.0763778676720807367055028],
[ 0.1526054652409226755052202, 0.0757044976845566746595428],
[ 0.2277858511416450780804962, 0.0745828754004991889865814],
[ 0.3016278681149130043205554, 0.0730306903327866674951894],
[ 0.3737060887154195606725482, 0.0710544235534440683057904],
[ 0.4435931752387251031999922, 0.0686486729285216193456234],
[ 0.5108670019508270980043641, 0.0658345971336184221115636],
[ 0.5751404468197103153429460, 0.0626532375547811680258701],
[ 0.6360536807265150254528367, 0.0591114008806395723749672],
[ 0.6932376563347513848054907, 0.0551951053482859947448324],
[ 0.7463319064601507926143051, 0.0509445739237286919327077],
[ 0.7950414288375511983506388, 0.0464348218674976747202319],
[ 0.8391169718222188233945291, 0.0416688733279736862637883],
[ 0.8782768112522819760774430, 0.0366001697582007980305572],
[ 0.9122344282513259058677524, 0.0312873067770327989585431],
[ 0.9408226338317547535199827, 0.0258821336049511588345051],
[ 0.9639719272779137912676661, 0.0203883734612665235980102],
[ 0.9815078774502502591933430, 0.0146261692569712529837880],
[ 0.9931285991850949247861224, 0.0086002698556429421986618],
[ 0.9988590315882776638383156, 0.0030735837185205315012183]],
dtype=float),
'dup': np.array(
[False, True, False, True, False, True, False, True, False,
True, False, True, False, True, False, True, False, True,
False, True, False, True, False, True, False, True, False,
True, False, True, False, True, False, True, False, True,
False, True, False, True, False], dtype=bool)
},
25: {
'g': np.array(
[[-0.9955569697904980979087849, 0.0113937985010262879479030],
[ -0.9766639214595175114983154, 0.0263549866150321372619018],
[ -0.9429745712289743394140112, 0.0409391567013063126556235],
[ -0.8949919978782753688510420, 0.0549046959758351919259369],
[ -0.8334426287608340014210211, 0.0680383338123569172071872],
[ -0.7592592630373576305772829, 0.0801407003350010180132350],
[ -0.6735663684734683644851206, 0.0910282619829636498114972],
[ -0.5776629302412229677236898, 0.1005359490670506442022069],
[ -0.4730027314457149605221821, 0.1085196244742636531160940],
[ -0.3611723058093878377358217, 0.1148582591457116483393255],
[ -0.2438668837209884320451904, 0.1194557635357847722281781],
[ -0.1228646926107103963873598, 0.1222424429903100416889595],
[ 0.0000000000000000000000000, 0.1231760537267154512039029],
[ 0.1228646926107103963873598, 0.1222424429903100416889595],
[ 0.2438668837209884320451904, 0.1194557635357847722281781],
[ 0.3611723058093878377358217, 0.1148582591457116483393255],
[ 0.4730027314457149605221821, 0.1085196244742636531160940],
[ 0.5776629302412229677236898, 0.1005359490670506442022069],
[ 0.6735663684734683644851206, 0.0910282619829636498114972],
[ 0.7592592630373576305772829, 0.0801407003350010180132350],
[ 0.8334426287608340014210211, 0.0680383338123569172071872],
[ 0.8949919978782753688510420, 0.0549046959758351919259369],
[ 0.9429745712289743394140112, 0.0409391567013063126556235],
[ 0.9766639214595175114983154, 0.0263549866150321372619018],
[ 0.9955569697904980979087849, 0.0113937985010262879479030]],
dtype=float),
'k': np.array(
[[-0.9992621049926098341934575, 0.0019873838923303159265079],
[ -0.9955569697904980979087849, 0.0055619321353567137580402],
[ -0.9880357945340772476373310, 0.0094739733861741516072077],
[ -0.9766639214595175114983154, 0.0132362291955716748136564],
[ -0.9616149864258425124181300, 0.0168478177091282982315167],
[ -0.9429745712289743394140112, 0.0204353711458828354565683],
[ -0.9207471152817015617463461, 0.0240099456069532162200925],
[ -0.8949919978782753688510420, 0.0274753175878517378029485],
[ -0.8658470652932755954489970, 0.0307923001673874888911090],
[ -0.8334426287608340014210211, 0.0340021302743293378367488],
[ -0.7978737979985000594104109, 0.0371162714834155435603306],
[ -0.7592592630373576305772829, 0.0400838255040323820748393],
[ -0.7177664068130843881866541, 0.0428728450201700494768958],
[ -0.6735663684734683644851206, 0.0455029130499217889098706],
[ -0.6268100990103174127881227, 0.0479825371388367139063923],
[ -0.5776629302412229677236898, 0.0502776790807156719633253],
[ -0.5263252843347191825996238, 0.0523628858064074758643667],
[ -0.4730027314457149605221821, 0.0542511298885454901445434],
[ -0.4178853821930377488518144, 0.0559508112204123173082407],
[ -0.3611723058093878377358217, 0.0574371163615678328535827],
[ -0.3030895389311078301674789, 0.0586896800223942079619742],
[ -0.2438668837209884320451904, 0.0597203403241740599790993],
[ -0.1837189394210488920159699, 0.0605394553760458629453603],
[ -0.1228646926107103963873598, 0.0611285097170530483058590],
[ -0.0615444830056850788865464, 0.0614711898714253166615441],
[ 0.0000000000000000000000000, 0.0615808180678329350787598],
[ 0.0615444830056850788865464, 0.0614711898714253166615441],
[ 0.1228646926107103963873598, 0.0611285097170530483058590],
[ 0.1837189394210488920159699, 0.0605394553760458629453603],
[ 0.2438668837209884320451904, 0.0597203403241740599790993],
[ 0.3030895389311078301674789, 0.0586896800223942079619742],
[ 0.3611723058093878377358217, 0.0574371163615678328535827],
[ 0.4178853821930377488518144, 0.0559508112204123173082407],
[ 0.4730027314457149605221821, 0.0542511298885454901445434],
[ 0.5263252843347191825996238, 0.0523628858064074758643667],
[ 0.5776629302412229677236898, 0.0502776790807156719633253],
[ 0.6268100990103174127881227, 0.0479825371388367139063923],
[ 0.6735663684734683644851206, 0.0455029130499217889098706],
[ 0.7177664068130843881866541, 0.0428728450201700494768958],
[ 0.7592592630373576305772829, 0.0400838255040323820748393],
[ 0.7978737979985000594104109, 0.0371162714834155435603306],
[ 0.8334426287608340014210211, 0.0340021302743293378367488],
[ 0.8658470652932755954489970, 0.0307923001673874888911090],
[ 0.8949919978782753688510420, 0.0274753175878517378029485],
[ 0.9207471152817015617463461, 0.0240099456069532162200925],
[ 0.9429745712289743394140112, 0.0204353711458828354565683],
[ 0.9616149864258425124181300, 0.0168478177091282982315167],
[ 0.9766639214595175114983154, 0.0132362291955716748136564],
[ 0.9880357945340772476373310, 0.0094739733861741516072077],
[ 0.9955569697904980979087849, 0.0055619321353567137580402],
[ 0.9992621049926098341934575, 0.0019873838923303159265079]],
dtype=float),
'dup': np.array(
[False, True, False, True, False, True, False, True, False,
True, False, True, False, True, False, True, False, True,
False, True, False, True, False, True, False, True, False,
True, False, True, False, True, False, True, False, True,
False, True, False, True, False, True, False, True, False,
True, False, True, False, True, False], dtype=bool)
},
30: {
'g': np.array(
[[-0.9968934840746495402716301, 0.0079681924961666056154659],
[ -0.9836681232797472099700326, 0.0184664683110909591423021],
[ -0.9600218649683075122168710, 0.0287847078833233693497192],
[ -0.9262000474292743258793243, 0.0387991925696270495968019],
[ -0.8825605357920526815431165, 0.0484026728305940529029381],
[ -0.8295657623827683974428981, 0.0574931562176190664817217],
[ -0.7677774321048261949179773, 0.0659742298821804951281285],
[ -0.6978504947933157969322924, 0.0737559747377052062682439],
[ -0.6205261829892428611404776, 0.0807558952294202153546949],
[ -0.5366241481420198992641698, 0.0868997872010829798023875],
[ -0.4470337695380891767806099, 0.0921225222377861287176327],
[ -0.3527047255308781134710372, 0.0963687371746442596394686],
[ -0.2546369261678898464398051, 0.0995934205867952670627803],
[ -0.1538699136085835469637947, 0.1017623897484055045964290],
[ -0.0514718425553176958330252, 0.1028526528935588403412856],
[ 0.0514718425553176958330252, 0.1028526528935588403412856],
[ 0.1538699136085835469637947, 0.1017623897484055045964290],
[ 0.2546369261678898464398051, 0.0995934205867952670627803],
[ 0.3527047255308781134710372, 0.0963687371746442596394686],
[ 0.4470337695380891767806099, 0.0921225222377861287176327],
[ 0.5366241481420198992641698, 0.0868997872010829798023875],
[ 0.6205261829892428611404776, 0.0807558952294202153546949],
[ 0.6978504947933157969322924, 0.0737559747377052062682439],
[ 0.7677774321048261949179773, 0.0659742298821804951281285],
[ 0.8295657623827683974428981, 0.0574931562176190664817217],
[ 0.8825605357920526815431165, 0.0484026728305940529029381],
[ 0.9262000474292743258793243, 0.0387991925696270495968019],
[ 0.9600218649683075122168710, 0.0287847078833233693497192],
[ 0.9836681232797472099700326, 0.0184664683110909591423021],
[ 0.9968934840746495402716301, 0.0079681924961666056154659]],
dtype=float),
'k': np.array(
[[-0.9994844100504906375713259, 0.0013890136986770076245516],
[ -0.9968934840746495402716301, 0.0038904611270998840512672],
[ -0.9916309968704045948586284, 0.0066307039159312921733198],
[ -0.9836681232797472099700326, 0.0092732796595177634284411],
[ -0.9731163225011262683746939, 0.0118230152534963417422329],
[ -0.9600218649683075122168710, 0.0143697295070458048124514],
[ -0.9443744447485599794158313, 0.0169208891890532726275723],
[ -0.9262000474292743258793243, 0.0194141411939423811734090],
[ -0.9055733076999077985465226, 0.0218280358216091922971675],
[ -0.8825605357920526815431165, 0.0241911620780806013656864],
[ -0.8572052335460610989586585, 0.0265099548823331016106017],
[ -0.8295657623827683974428981, 0.0287540487650412928439788],
[ -0.7997278358218390830136689, 0.0309072575623877624728843],
[ -0.7677774321048261949179773, 0.0329814470574837260318142],
[ -0.7337900624532268047261711, 0.0349793380280600241374997],
[ -0.6978504947933157969322924, 0.0368823646518212292239111],
[ -0.6600610641266269613700537, 0.0386789456247275929503487],
[ -0.6205261829892428611404776, 0.0403745389515359591119953],
[ -0.5793452358263616917560249, 0.0419698102151642461471475],
[ -0.5366241481420198992641698, 0.0434525397013560693168317],
[ -0.4924804678617785749936931, 0.0448148001331626631923556],
[ -0.4470337695380891767806099, 0.0460592382710069881162717],
[ -0.4004012548303943925354762, 0.0471855465692991539452615],
[ -0.3527047255308781134710372, 0.0481858617570871291407795],
[ -0.3040732022736250773726771, 0.0490554345550297788875282],
[ -0.2546369261678898464398051, 0.0497956834270742063578116],
[ -0.2045251166823098914389577, 0.0504059214027823468408931],
[ -0.1538699136085835469637947, 0.0508817958987496064922975],
[ -0.1028069379667370301470968, 0.0512215478492587721706563],
[ -0.0514718425553176958330252, 0.0514261285374590259338629],
[ 0.0000000000000000000000000, 0.0514947294294515675583404],
[ 0.0514718425553176958330252, 0.0514261285374590259338629],
[ 0.1028069379667370301470968, 0.0512215478492587721706563],
[ 0.1538699136085835469637947, 0.0508817958987496064922975],
[ 0.2045251166823098914389577, 0.0504059214027823468408931],
[ 0.2546369261678898464398051, 0.0497956834270742063578116],
[ 0.3040732022736250773726771, 0.0490554345550297788875282],
[ 0.3527047255308781134710372, 0.0481858617570871291407795],
[ 0.4004012548303943925354762, 0.0471855465692991539452615],
[ 0.4470337695380891767806099, 0.0460592382710069881162717],
[ 0.4924804678617785749936931, 0.0448148001331626631923556],
[ 0.5366241481420198992641698, 0.0434525397013560693168317],
[ 0.5793452358263616917560249, 0.0419698102151642461471475],
[ 0.6205261829892428611404776, 0.0403745389515359591119953],
[ 0.6600610641266269613700537, 0.0386789456247275929503487],
[ 0.6978504947933157969322924, 0.0368823646518212292239111],
[ 0.7337900624532268047261711, 0.0349793380280600241374997],
[ 0.7677774321048261949179773, 0.0329814470574837260318142],
[ 0.7997278358218390830136689, 0.0309072575623877624728843],
[ 0.8295657623827683974428981, 0.0287540487650412928439788],
[ 0.8572052335460610989586585, 0.0265099548823331016106017],
[ 0.8825605357920526815431165, 0.0241911620780806013656864],
[ 0.9055733076999077985465226, 0.0218280358216091922971675],
[ 0.9262000474292743258793243, 0.0194141411939423811734090],
[ 0.9443744447485599794158313, 0.0169208891890532726275723],
[ 0.9600218649683075122168710, 0.0143697295070458048124514],
[ 0.9731163225011262683746939, 0.0118230152534963417422329],
[ 0.9836681232797472099700326, 0.0092732796595177634284411],
[ 0.9916309968704045948586284, 0.0066307039159312921733198],
[ 0.9968934840746495402716301, 0.0038904611270998840512672],
[ 0.9994844100504906375713259, 0.0013890136986770076245516]],
dtype=float),
'dup': np.array(
[False, True, False, True, False, True, False, True, False,
True, False, True, False, True, False, True, False, True,
False, True, False, True, False, True, False, True, False,
True, False, True, False, True, False, True, False, True,
False, True, False, True, False, True, False, True, False,
True, False, True, False, True, False, True, False, True,
False, True, False, True, False, True, False], dtype=bool)
},
}
w = weights[n]
dup=w['dup']
xi = w['k'][:,0]
wi1 = np.zeros_like(xi)
wi1[dup] = w['g'][:, 1]
wi2 = w['k'][:,1]
return xi, wi1, wi2
[docs]def gauss_legendre_abscissae_and_weights(n):
"""Gauss-Legendre quadrature abscissae and weights
Integral = Sum(f(xi) * wi)
Parameters
----------
n : [2-20, 32, 64, 100]
Number of integration points
Returns
-------
xi, wi : 1d array of len(n)
Abscissae and weights for numericla integration
References
----------
.. [1] Holoborodko, Pavel. 2014. 'Numerical Integration'. Accessed
April 24.
http://www.holoborodko.com/pavel/numerical-methods/numerical-integration/.
"""
if n not in [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,
32, 64, 100]:
raise ValueError('n must be 2-20, 32, 64, or 100')
weights = {
2: np.array(
[[-0.5773502691896257645091488, 1.0000000000000000000000000],
[ 0.5773502691896257645091488, 1.0000000000000000000000000]],
dtype=float),
3: np.array(
[[-0.7745966692414833770358531, 0.5555555555555555555555556],
[ 0, 0.8888888888888888888888889],
[ 0.7745966692414833770358531, 0.5555555555555555555555556]],
dtype=float),
4: np.array(
[[-0.8611363115940525752239465, 0.3478548451374538573730639],
[ -0.3399810435848562648026658, 0.6521451548625461426269361],
[ 0.3399810435848562648026658, 0.6521451548625461426269361],
[ 0.8611363115940525752239465, 0.3478548451374538573730639]],
dtype=float),
5: np.array(
[[-0.9061798459386639927976269, 0.2369268850561890875142640],
[ -0.5384693101056830910363144, 0.4786286704993664680412915],
[ 0, 0.5688888888888888888888889],
[ 0.5384693101056830910363144, 0.4786286704993664680412915],
[ 0.9061798459386639927976269, 0.2369268850561890875142640]],
dtype=float),
6: np.array(
[[-0.9324695142031520278123016, 0.1713244923791703450402961],
[ -0.6612093864662645136613996, 0.3607615730481386075698335],
[ -0.2386191860831969086305017, 0.4679139345726910473898703],
[ 0.2386191860831969086305017, 0.4679139345726910473898703],
[ 0.6612093864662645136613996, 0.3607615730481386075698335],
[ 0.9324695142031520278123016, 0.1713244923791703450402961]],
dtype=float),
7: np.array(
[[-0.9491079123427585245261897, 0.1294849661688696932706114],
[ -0.7415311855993944398638648, 0.2797053914892766679014678],
[ -0.4058451513773971669066064, 0.3818300505051189449503698],
[ 0, 0.4179591836734693877551020],
[ 0.4058451513773971669066064, 0.3818300505051189449503698],
[ 0.7415311855993944398638648, 0.2797053914892766679014678],
[ 0.9491079123427585245261897, 0.1294849661688696932706114]],
dtype=float),
8: np.array(
[[-0.9602898564975362316835609, 0.1012285362903762591525314],
[ -0.7966664774136267395915539, 0.2223810344533744705443560],
[ -0.5255324099163289858177390, 0.3137066458778872873379622],
[ -0.1834346424956498049394761, 0.3626837833783619829651504],
[ 0.1834346424956498049394761, 0.3626837833783619829651504],
[ 0.5255324099163289858177390, 0.3137066458778872873379622],
[ 0.7966664774136267395915539, 0.2223810344533744705443560],
[ 0.9602898564975362316835609, 0.1012285362903762591525314]],
dtype=float),
9: np.array(
[[-0.9681602395076260898355762, 0.0812743883615744119718922],
[ -0.8360311073266357942994298, 0.1806481606948574040584720],
[ -0.6133714327005903973087020, 0.2606106964029354623187429],
[ -0.3242534234038089290385380, 0.3123470770400028400686304],
[ 0, 0.3302393550012597631645251],
[ 0.3242534234038089290385380, 0.3123470770400028400686304],
[ 0.6133714327005903973087020, 0.2606106964029354623187429],
[ 0.8360311073266357942994298, 0.1806481606948574040584720],
[ 0.9681602395076260898355762, 0.0812743883615744119718922]],
dtype=float),
10: np.array(
[[-0.9739065285171717200779640, 0.0666713443086881375935688],
[ -0.8650633666889845107320967, 0.1494513491505805931457763],
[ -0.6794095682990244062343274, 0.2190863625159820439955349],
[ -0.4333953941292471907992659, 0.2692667193099963550912269],
[ -0.1488743389816312108848260, 0.2955242247147528701738930],
[ 0.1488743389816312108848260, 0.2955242247147528701738930],
[ 0.4333953941292471907992659, 0.2692667193099963550912269],
[ 0.6794095682990244062343274, 0.2190863625159820439955349],
[ 0.8650633666889845107320967, 0.1494513491505805931457763],
[ 0.9739065285171717200779640, 0.0666713443086881375935688]],
dtype=float),
11: np.array(
[[-0.9782286581460569928039380, 0.0556685671161736664827537],
[ -0.8870625997680952990751578, 0.1255803694649046246346943],
[ -0.7301520055740493240934163, 0.1862902109277342514260976],
[ -0.5190961292068118159257257, 0.2331937645919904799185237],
[ -0.2695431559523449723315320, 0.2628045445102466621806889],
[ 0, 0.2729250867779006307144835],
[ 0.2695431559523449723315320, 0.2628045445102466621806889],
[ 0.5190961292068118159257257, 0.2331937645919904799185237],
[ 0.7301520055740493240934163, 0.1862902109277342514260976],
[ 0.8870625997680952990751578, 0.1255803694649046246346943],
[ 0.9782286581460569928039380, 0.0556685671161736664827537]],
dtype=float),
12: np.array(
[[-0.9815606342467192506905491, 0.0471753363865118271946160],
[ -0.9041172563704748566784659, 0.1069393259953184309602547],
[ -0.7699026741943046870368938, 0.1600783285433462263346525],
[ -0.5873179542866174472967024, 0.2031674267230659217490645],
[ -0.3678314989981801937526915, 0.2334925365383548087608499],
[ -0.1252334085114689154724414, 0.2491470458134027850005624],
[ 0.1252334085114689154724414, 0.2491470458134027850005624],
[ 0.3678314989981801937526915, 0.2334925365383548087608499],
[ 0.5873179542866174472967024, 0.2031674267230659217490645],
[ 0.7699026741943046870368938, 0.1600783285433462263346525],
[ 0.9041172563704748566784659, 0.1069393259953184309602547],
[ 0.9815606342467192506905491, 0.0471753363865118271946160]],
dtype=float),
13: np.array(
[[-0.9841830547185881494728294, 0.0404840047653158795200216],
[ -0.9175983992229779652065478, 0.0921214998377284479144218],
[ -0.8015780907333099127942065, 0.1388735102197872384636018],
[ -0.6423493394403402206439846, 0.1781459807619457382800467],
[ -0.4484927510364468528779129, 0.2078160475368885023125232],
[ -0.2304583159551347940655281, 0.2262831802628972384120902],
[ 0, 0.2325515532308739101945895],
[ 0.2304583159551347940655281, 0.2262831802628972384120902],
[ 0.4484927510364468528779129, 0.2078160475368885023125232],
[ 0.6423493394403402206439846, 0.1781459807619457382800467],
[ 0.8015780907333099127942065, 0.1388735102197872384636018],
[ 0.9175983992229779652065478, 0.0921214998377284479144218],
[ 0.9841830547185881494728294, 0.0404840047653158795200216]],
dtype=float),
14: np.array(
[[-0.9862838086968123388415973, 0.0351194603317518630318329],
[ -0.9284348836635735173363911, 0.0801580871597602098056333],
[ -0.8272013150697649931897947, 0.1215185706879031846894148],
[ -0.6872929048116854701480198, 0.1572031671581935345696019],
[ -0.5152486363581540919652907, 0.1855383974779378137417166],
[ -0.3191123689278897604356718, 0.2051984637212956039659241],
[ -0.1080549487073436620662447, 0.2152638534631577901958764],
[ 0.1080549487073436620662447, 0.2152638534631577901958764],
[ 0.3191123689278897604356718, 0.2051984637212956039659241],
[ 0.5152486363581540919652907, 0.1855383974779378137417166],
[ 0.6872929048116854701480198, 0.1572031671581935345696019],
[ 0.8272013150697649931897947, 0.1215185706879031846894148],
[ 0.9284348836635735173363911, 0.0801580871597602098056333],
[ 0.9862838086968123388415973, 0.0351194603317518630318329]],
dtype=float),
15: np.array(
[[-0.9879925180204854284895657, 0.0307532419961172683546284],
[ -0.9372733924007059043077589, 0.0703660474881081247092674],
[ -0.8482065834104272162006483, 0.1071592204671719350118695],
[ -0.7244177313601700474161861, 0.1395706779261543144478048],
[ -0.5709721726085388475372267, 0.1662692058169939335532009],
[ -0.3941513470775633698972074, 0.1861610000155622110268006],
[ -0.2011940939974345223006283, 0.1984314853271115764561183],
[ 0, 0.2025782419255612728806202],
[ 0.2011940939974345223006283, 0.1984314853271115764561183],
[ 0.3941513470775633698972074, 0.1861610000155622110268006],
[ 0.5709721726085388475372267, 0.1662692058169939335532009],
[ 0.7244177313601700474161861, 0.1395706779261543144478048],
[ 0.8482065834104272162006483, 0.1071592204671719350118695],
[ 0.9372733924007059043077589, 0.0703660474881081247092674],
[ 0.9879925180204854284895657, 0.0307532419961172683546284]],
dtype=float),
16: np.array(
[[-0.9894009349916499325961542, 0.0271524594117540948517806],
[ -0.9445750230732325760779884, 0.0622535239386478928628438],
[ -0.8656312023878317438804679, 0.0951585116824927848099251],
[ -0.7554044083550030338951012, 0.1246289712555338720524763],
[ -0.6178762444026437484466718, 0.1495959888165767320815017],
[ -0.4580167776572273863424194, 0.1691565193950025381893121],
[ -0.2816035507792589132304605, 0.1826034150449235888667637],
[ -0.0950125098376374401853193, 0.1894506104550684962853967],
[ 0.0950125098376374401853193, 0.1894506104550684962853967],
[ 0.2816035507792589132304605, 0.1826034150449235888667637],
[ 0.4580167776572273863424194, 0.1691565193950025381893121],
[ 0.6178762444026437484466718, 0.1495959888165767320815017],
[ 0.7554044083550030338951012, 0.1246289712555338720524763],
[ 0.8656312023878317438804679, 0.0951585116824927848099251],
[ 0.9445750230732325760779884, 0.0622535239386478928628438],
[ 0.9894009349916499325961542, 0.0271524594117540948517806]],
dtype=float),
17: np.array(
[[-0.9905754753144173356754340, 0.0241483028685479319601100],
[ -0.9506755217687677612227170, 0.0554595293739872011294402],
[ -0.8802391537269859021229557, 0.0850361483171791808835354],
[ -0.7815140038968014069252301, 0.1118838471934039710947884],
[ -0.6576711592166907658503022, 0.1351363684685254732863200],
[ -0.5126905370864769678862466, 0.1540457610768102880814316],
[ -0.3512317634538763152971855, 0.1680041021564500445099707],
[ -0.1784841814958478558506775, 0.1765627053669926463252710],
[ 0, 0.1794464703562065254582656],
[ 0.1784841814958478558506775, 0.1765627053669926463252710],
[ 0.3512317634538763152971855, 0.1680041021564500445099707],
[ 0.5126905370864769678862466, 0.1540457610768102880814316],
[ 0.6576711592166907658503022, 0.1351363684685254732863200],
[ 0.7815140038968014069252301, 0.1118838471934039710947884],
[ 0.8802391537269859021229557, 0.0850361483171791808835354],
[ 0.9506755217687677612227170, 0.0554595293739872011294402],
[ 0.9905754753144173356754340, 0.0241483028685479319601100]],
dtype=float),
18: np.array(
[[-0.9915651684209309467300160, 0.0216160135264833103133427],
[ -0.9558239495713977551811959, 0.0497145488949697964533349],
[ -0.8926024664975557392060606, 0.0764257302548890565291297],
[ -0.8037049589725231156824175, 0.1009420441062871655628140],
[ -0.6916870430603532078748911, 0.1225552067114784601845191],
[ -0.5597708310739475346078715, 0.1406429146706506512047313],
[ -0.4117511614628426460359318, 0.1546846751262652449254180],
[ -0.2518862256915055095889729, 0.1642764837458327229860538],
[ -0.0847750130417353012422619, 0.1691423829631435918406565],
[ 0.0847750130417353012422619, 0.1691423829631435918406565],
[ 0.2518862256915055095889729, 0.1642764837458327229860538],
[ 0.4117511614628426460359318, 0.1546846751262652449254180],
[ 0.5597708310739475346078715, 0.1406429146706506512047313],
[ 0.6916870430603532078748911, 0.1225552067114784601845191],
[ 0.8037049589725231156824175, 0.1009420441062871655628140],
[ 0.8926024664975557392060606, 0.0764257302548890565291297],
[ 0.9558239495713977551811959, 0.0497145488949697964533349],
[ 0.9915651684209309467300160, 0.0216160135264833103133427]],
dtype=float),
19: np.array(
[[-0.9924068438435844031890177, 0.0194617882297264770363120],
[ -0.9602081521348300308527788, 0.0448142267656996003328382],
[ -0.9031559036148179016426609, 0.0690445427376412265807083],
[ -0.8227146565371428249789225, 0.0914900216224499994644621],
[ -0.7209661773352293786170959, 0.1115666455473339947160239],
[ -0.6005453046616810234696382, 0.1287539625393362276755158],
[ -0.4645707413759609457172671, 0.1426067021736066117757461],
[ -0.3165640999636298319901173, 0.1527660420658596667788554],
[ -0.1603586456402253758680961, 0.1589688433939543476499564],
[ 0, 0.1610544498487836959791636],
[ 0.1603586456402253758680961, 0.1589688433939543476499564],
[ 0.3165640999636298319901173, 0.1527660420658596667788554],
[ 0.4645707413759609457172671, 0.1426067021736066117757461],
[ 0.6005453046616810234696382, 0.1287539625393362276755158],
[ 0.7209661773352293786170959, 0.1115666455473339947160239],
[ 0.8227146565371428249789225, 0.0914900216224499994644621],
[ 0.9031559036148179016426609, 0.0690445427376412265807083],
[ 0.9602081521348300308527788, 0.0448142267656996003328382],
[ 0.9924068438435844031890177, 0.0194617882297264770363120]],
dtype=float),
20: np.array(
[[-0.9931285991850949247861224, 0.0176140071391521183118620],
[ -0.9639719272779137912676661, 0.0406014298003869413310400],
[ -0.9122344282513259058677524, 0.0626720483341090635695065],
[ -0.8391169718222188233945291, 0.0832767415767047487247581],
[ -0.7463319064601507926143051, 0.1019301198172404350367501],
[ -0.6360536807265150254528367, 0.1181945319615184173123774],
[ -0.5108670019508270980043641, 0.1316886384491766268984945],
[ -0.3737060887154195606725482, 0.1420961093183820513292983],
[ -0.2277858511416450780804962, 0.1491729864726037467878287],
[ -0.0765265211334973337546404, 0.1527533871307258506980843],
[ 0.0765265211334973337546404, 0.1527533871307258506980843],
[ 0.2277858511416450780804962, 0.1491729864726037467878287],
[ 0.3737060887154195606725482, 0.1420961093183820513292983],
[ 0.5108670019508270980043641, 0.1316886384491766268984945],
[ 0.6360536807265150254528367, 0.1181945319615184173123774],
[ 0.7463319064601507926143051, 0.1019301198172404350367501],
[ 0.8391169718222188233945291, 0.0832767415767047487247581],
[ 0.9122344282513259058677524, 0.0626720483341090635695065],
[ 0.9639719272779137912676661, 0.0406014298003869413310400],
[ 0.9931285991850949247861224, 0.0176140071391521183118620]],
dtype=float),
32: np.array(
[[-0.9972638618494815635449811, 0.0070186100094700966004071],
[ -0.9856115115452683354001750, 0.0162743947309056706051706],
[ -0.9647622555875064307738119, 0.0253920653092620594557526],
[ -0.9349060759377396891709191, 0.0342738629130214331026877],
[ -0.8963211557660521239653072, 0.0428358980222266806568786],
[ -0.8493676137325699701336930, 0.0509980592623761761961632],
[ -0.7944837959679424069630973, 0.0586840934785355471452836],
[ -0.7321821187402896803874267, 0.0658222227763618468376501],
[ -0.6630442669302152009751152, 0.0723457941088485062253994],
[ -0.5877157572407623290407455, 0.0781938957870703064717409],
[ -0.5068999089322293900237475, 0.0833119242269467552221991],
[ -0.4213512761306353453641194, 0.0876520930044038111427715],
[ -0.3318686022821276497799168, 0.0911738786957638847128686],
[ -0.2392873622521370745446032, 0.0938443990808045656391802],
[ -0.1444719615827964934851864, 0.0956387200792748594190820],
[ -0.0483076656877383162348126, 0.0965400885147278005667648],
[ 0.0483076656877383162348126, 0.0965400885147278005667648],
[ 0.1444719615827964934851864, 0.0956387200792748594190820],
[ 0.2392873622521370745446032, 0.0938443990808045656391802],
[ 0.3318686022821276497799168, 0.0911738786957638847128686],
[ 0.4213512761306353453641194, 0.0876520930044038111427715],
[ 0.5068999089322293900237475, 0.0833119242269467552221991],
[ 0.5877157572407623290407455, 0.0781938957870703064717409],
[ 0.6630442669302152009751152, 0.0723457941088485062253994],
[ 0.7321821187402896803874267, 0.0658222227763618468376501],
[ 0.7944837959679424069630973, 0.0586840934785355471452836],
[ 0.8493676137325699701336930, 0.0509980592623761761961632],
[ 0.8963211557660521239653072, 0.0428358980222266806568786],
[ 0.9349060759377396891709191, 0.0342738629130214331026877],
[ 0.9647622555875064307738119, 0.0253920653092620594557526],
[ 0.9856115115452683354001750, 0.0162743947309056706051706],
[ 0.9972638618494815635449811, 0.0070186100094700966004071]],
dtype=float),
64: np.array(
[[-0.9993050417357721394569056, 0.0017832807216964329472961],
[ -0.9963401167719552793469245, 0.0041470332605624676352875],
[ -0.9910133714767443207393824, 0.0065044579689783628561174],
[ -0.9833362538846259569312993, 0.0088467598263639477230309],
[ -0.9733268277899109637418535, 0.0111681394601311288185905],
[ -0.9610087996520537189186141, 0.0134630478967186425980608],
[ -0.9464113748584028160624815, 0.0157260304760247193219660],
[ -0.9295691721319395758214902, 0.0179517157756973430850453],
[ -0.9105221370785028057563807, 0.0201348231535302093723403],
[ -0.8893154459951141058534040, 0.0222701738083832541592983],
[ -0.8659993981540928197607834, 0.0243527025687108733381776],
[ -0.8406292962525803627516915, 0.0263774697150546586716918],
[ -0.8132653151227975597419233, 0.0283396726142594832275113],
[ -0.7839723589433414076102205, 0.0302346570724024788679741],
[ -0.7528199072605318966118638, 0.0320579283548515535854675],
[ -0.7198818501716108268489402, 0.0338051618371416093915655],
[ -0.6852363130542332425635584, 0.0354722132568823838106931],
[ -0.6489654712546573398577612, 0.0370551285402400460404151],
[ -0.6111553551723932502488530, 0.0385501531786156291289625],
[ -0.5718956462026340342838781, 0.0399537411327203413866569],
[ -0.5312794640198945456580139, 0.0412625632426235286101563],
[ -0.4894031457070529574785263, 0.0424735151236535890073398],
[ -0.4463660172534640879849477, 0.0435837245293234533768279],
[ -0.4022701579639916036957668, 0.0445905581637565630601347],
[ -0.3572201583376681159504426, 0.0454916279274181444797710],
[ -0.3113228719902109561575127, 0.0462847965813144172959532],
[ -0.2646871622087674163739642, 0.0469681828162100173253263],
[ -0.2174236437400070841496487, 0.0475401657148303086622822],
[ -0.1696444204239928180373136, 0.0479993885964583077281262],
[ -0.1214628192961205544703765, 0.0483447622348029571697695],
[ -0.0729931217877990394495429, 0.0485754674415034269347991],
[ -0.0243502926634244325089558, 0.0486909570091397203833654],
[ 0.0243502926634244325089558, 0.0486909570091397203833654],
[ 0.0729931217877990394495429, 0.0485754674415034269347991],
[ 0.1214628192961205544703765, 0.0483447622348029571697695],
[ 0.1696444204239928180373136, 0.0479993885964583077281262],
[ 0.2174236437400070841496487, 0.0475401657148303086622822],
[ 0.2646871622087674163739642, 0.0469681828162100173253263],
[ 0.3113228719902109561575127, 0.0462847965813144172959532],
[ 0.3572201583376681159504426, 0.0454916279274181444797710],
[ 0.4022701579639916036957668, 0.0445905581637565630601347],
[ 0.4463660172534640879849477, 0.0435837245293234533768279],
[ 0.4894031457070529574785263, 0.0424735151236535890073398],
[ 0.5312794640198945456580139, 0.0412625632426235286101563],
[ 0.5718956462026340342838781, 0.0399537411327203413866569],
[ 0.6111553551723932502488530, 0.0385501531786156291289625],
[ 0.6489654712546573398577612, 0.0370551285402400460404151],
[ 0.6852363130542332425635584, 0.0354722132568823838106931],
[ 0.7198818501716108268489402, 0.0338051618371416093915655],
[ 0.7528199072605318966118638, 0.0320579283548515535854675],
[ 0.7839723589433414076102205, 0.0302346570724024788679741],
[ 0.8132653151227975597419233, 0.0283396726142594832275113],
[ 0.8406292962525803627516915, 0.0263774697150546586716918],
[ 0.8659993981540928197607834, 0.0243527025687108733381776],
[ 0.8893154459951141058534040, 0.0222701738083832541592983],
[ 0.9105221370785028057563807, 0.0201348231535302093723403],
[ 0.9295691721319395758214902, 0.0179517157756973430850453],
[ 0.9464113748584028160624815, 0.0157260304760247193219660],
[ 0.9610087996520537189186141, 0.0134630478967186425980608],
[ 0.9733268277899109637418535, 0.0111681394601311288185905],
[ 0.9833362538846259569312993, 0.0088467598263639477230309],
[ 0.9910133714767443207393824, 0.0065044579689783628561174],
[ 0.9963401167719552793469245, 0.0041470332605624676352875],
[ 0.9993050417357721394569056, 0.0017832807216964329472961]],
dtype=float),
100: np.array(
[[-0.9997137267734412336782285, 0.0007346344905056717304063],
[ -0.9984919506395958184001634, 0.0017093926535181052395294],
[ -0.9962951347331251491861317, 0.0026839253715534824194396],
[ -0.9931249370374434596520099, 0.0036559612013263751823425],
[ -0.9889843952429917480044187, 0.0046244500634221193510958],
[ -0.9838775407060570154961002, 0.0055884280038655151572119],
[ -0.9778093584869182885537811, 0.0065469484508453227641521],
[ -0.9707857757637063319308979, 0.0074990732554647115788287],
[ -0.9628136542558155272936593, 0.0084438714696689714026208],
[ -0.9539007829254917428493369, 0.0093804196536944579514182],
[ -0.9440558701362559779627747, 0.0103078025748689695857821],
[ -0.9332885350430795459243337, 0.0112251140231859771172216],
[ -0.9216092981453339526669513, 0.0121314576629794974077448],
[ -0.9090295709825296904671263, 0.0130259478929715422855586],
[ -0.8955616449707269866985210, 0.0139077107037187726879541],
[ -0.8812186793850184155733168, 0.0147758845274413017688800],
[ -0.8660146884971646234107400, 0.0156296210775460027239369],
[ -0.8499645278795912842933626, 0.0164680861761452126431050],
[ -0.8330838798884008235429158, 0.0172904605683235824393442],
[ -0.8153892383391762543939888, 0.0180959407221281166643908],
[ -0.7968978923903144763895729, 0.0188837396133749045529412],
[ -0.7776279096494954756275514, 0.0196530874944353058653815],
[ -0.7575981185197071760356680, 0.0204032326462094327668389],
[ -0.7368280898020207055124277, 0.0211334421125276415426723],
[ -0.7153381175730564464599671, 0.0218430024162473863139537],
[ -0.6931491993558019659486479, 0.0225312202563362727017970],
[ -0.6702830156031410158025870, 0.0231974231852541216224889],
[ -0.6467619085141292798326303, 0.0238409602659682059625604],
[ -0.6226088602037077716041908, 0.0244612027079570527199750],
[ -0.5978474702471787212648065, 0.0250575444815795897037642],
[ -0.5725019326213811913168704, 0.0256294029102081160756420],
[ -0.5465970120650941674679943, 0.0261762192395456763423087],
[ -0.5201580198817630566468157, 0.0266974591835709626603847],
[ -0.4932107892081909335693088, 0.0271926134465768801364916],
[ -0.4657816497733580422492166, 0.0276611982207923882942042],
[ -0.4378974021720315131089780, 0.0281027556591011733176483],
[ -0.4095852916783015425288684, 0.0285168543223950979909368],
[ -0.3808729816246299567633625, 0.0289030896011252031348762],
[ -0.3517885263724217209723438, 0.0292610841106382766201190],
[ -0.3223603439005291517224766, 0.0295904880599126425117545],
[ -0.2926171880384719647375559, 0.0298909795933328309168368],
[ -0.2625881203715034791689293, 0.0301622651051691449190687],
[ -0.2323024818449739696495100, 0.0304040795264548200165079],
[ -0.2017898640957359972360489, 0.0306161865839804484964594],
[ -0.1710800805386032748875324, 0.0307983790311525904277139],
[ -0.1402031372361139732075146, 0.0309504788504909882340635],
[ -0.1091892035800611150034260, 0.0310723374275665165878102],
[ -0.0780685828134366366948174, 0.0311638356962099067838183],
[ -0.0468716824215916316149239, 0.0312248842548493577323765],
[ -0.0156289844215430828722167, 0.0312554234538633569476425],
[ 0.0156289844215430828722167, 0.0312554234538633569476425],
[ 0.0468716824215916316149239, 0.0312248842548493577323765],
[ 0.0780685828134366366948174, 0.0311638356962099067838183],
[ 0.1091892035800611150034260, 0.0310723374275665165878102],
[ 0.1402031372361139732075146, 0.0309504788504909882340635],
[ 0.1710800805386032748875324, 0.0307983790311525904277139],
[ 0.2017898640957359972360489, 0.0306161865839804484964594],
[ 0.2323024818449739696495100, 0.0304040795264548200165079],
[ 0.2625881203715034791689293, 0.0301622651051691449190687],
[ 0.2926171880384719647375559, 0.0298909795933328309168368],
[ 0.3223603439005291517224766, 0.0295904880599126425117545],
[ 0.3517885263724217209723438, 0.0292610841106382766201190],
[ 0.3808729816246299567633625, 0.0289030896011252031348762],
[ 0.4095852916783015425288684, 0.0285168543223950979909368],
[ 0.4378974021720315131089780, 0.0281027556591011733176483],
[ 0.4657816497733580422492166, 0.0276611982207923882942042],
[ 0.4932107892081909335693088, 0.0271926134465768801364916],
[ 0.5201580198817630566468157, 0.0266974591835709626603847],
[ 0.5465970120650941674679943, 0.0261762192395456763423087],
[ 0.5725019326213811913168704, 0.0256294029102081160756420],
[ 0.5978474702471787212648065, 0.0250575444815795897037642],
[ 0.6226088602037077716041908, 0.0244612027079570527199750],
[ 0.6467619085141292798326303, 0.0238409602659682059625604],
[ 0.6702830156031410158025870, 0.0231974231852541216224889],
[ 0.6931491993558019659486479, 0.0225312202563362727017970],
[ 0.7153381175730564464599671, 0.0218430024162473863139537],
[ 0.7368280898020207055124277, 0.0211334421125276415426723],
[ 0.7575981185197071760356680, 0.0204032326462094327668389],
[ 0.7776279096494954756275514, 0.0196530874944353058653815],
[ 0.7968978923903144763895729, 0.0188837396133749045529412],
[ 0.8153892383391762543939888, 0.0180959407221281166643908],
[ 0.8330838798884008235429158, 0.0172904605683235824393442],
[ 0.8499645278795912842933626, 0.0164680861761452126431050],
[ 0.8660146884971646234107400, 0.0156296210775460027239369],
[ 0.8812186793850184155733168, 0.0147758845274413017688800],
[ 0.8955616449707269866985210, 0.0139077107037187726879541],
[ 0.9090295709825296904671263, 0.0130259478929715422855586],
[ 0.9216092981453339526669513, 0.0121314576629794974077448],
[ 0.9332885350430795459243337, 0.0112251140231859771172216],
[ 0.9440558701362559779627747, 0.0103078025748689695857821],
[ 0.9539007829254917428493369, 0.0093804196536944579514182],
[ 0.9628136542558155272936593, 0.0084438714696689714026208],
[ 0.9707857757637063319308979, 0.0074990732554647115788287],
[ 0.9778093584869182885537811, 0.0065469484508453227641521],
[ 0.9838775407060570154961002, 0.0055884280038655151572119],
[ 0.9889843952429917480044187, 0.0046244500634221193510958],
[ 0.9931249370374434596520099, 0.0036559612013263751823425],
[ 0.9962951347331251491861317, 0.0026839253715534824194396],
[ 0.9984919506395958184001634, 0.0017093926535181052395294],
[ 0.9997137267734412336782285, 0.0007346344905056717304063]],
dtype=float),
}
return weights[n][:,0], weights[n][:,1]
[docs]def shanks_table(seq, table=None, randomized=False):
r'''Copied from sympy.mpmath.mpmath.calculus.extrapolation.py
This shanks function is taken almost verbatim (minus an initial ctx
argument???) from sympy.mpmath.mpmath.calculus.extrapolation.py:
- http://docs.sympy.org/dev/modules/mpmath/calculus/sums_limits.html#mpmath.shanks
- https://github.com/sympy/sympy/blob/master/sympy/mpmath/calculus/extrapolation.py
mpmath is BSD license
Notes
-----
Given a list ``seq`` of the first `N` elements of a slowly
convergent infinite sequence `(A_k)`, :func:`~mpmath.shanks` computes the iterated
Shanks transformation `S(A), S(S(A)), \ldots, S^{N/2}(A)`. The Shanks
transformation often provides strong convergence acceleration,
especially if the sequence is oscillating.
The iterated Shanks transformation is computed using the Wynn
epsilon algorithm (see [1]). :func:`~mpmath.shanks` returns the full
epsilon table generated by Wynn's algorithm, which can be read
off as follows:
- The table is a list of lists forming a lower triangular matrix,
where higher row and column indices correspond to more accurate
values.
- The columns with even index hold dummy entries (required for the
computation) and the columns with odd index hold the actual
extrapolates.
- The last element in the last row is typically the most
accurate estimate of the limit.
- The difference to the third last element in the last row
provides an estimate of the approximation error.
- The magnitude of the second last element provides an estimate
of the numerical accuracy lost to cancellation.
For convenience, so the extrapolation is stopped at an odd index
so that ``shanks(seq)[-1][-1]`` always gives an estimate of the
limit.
Optionally, an existing table can be passed to :func:`~mpmath.shanks`.
This can be used to efficiently extend a previous computation after
new elements have been appended to the sequence. The table will
then be updated in-place.
The Shanks transformation:
The Shanks transformation is defined as follows (see [2]): given
the input sequence `(A_0, A_1, \ldots)`, the transformed sequence is
given by
.. math ::
S(A_k) = \frac{A_{k+1}A_{k-1}-A_k^2}{A_{k+1}+A_{k-1}-2 A_k}
The Shanks transformation gives the exact limit `A_{\infty}` in a
single step if `A_k = A + a q^k`. Note in particular that it
extrapolates the exact sum of a geometric series in a single step.
Applying the Shanks transformation once often improves convergence
substantially for an arbitrary sequence, but the optimal effect is
obtained by applying it iteratively:
`S(S(A_k)), S(S(S(A_k))), \ldots`.
Wynn's epsilon algorithm provides an efficient way to generate
the table of iterated Shanks transformations. It reduces the
computation of each element to essentially a single division, at
the cost of requiring dummy elements in the table. See [1] for
details.
Precision issues:
Due to cancellation effects, the sequence must be typically be
computed at a much higher precision than the target accuracy
of the extrapolation.
If the Shanks transformation converges to the exact limit (such
as if the sequence is a geometric series), then a division by
zero occurs. By default, :func:`~mpmath.shanks` handles this case by
terminating the iteration and returning the table it has
generated so far. With *randomized=True*, it will instead
replace the zero by a pseudorandom number close to zero.
(TODO: find a better solution to this problem.)
Examples (truncated from original)
We illustrate by applying Shanks transformation to the Leibniz
series for `\pi`:
>>> S = [4*sum((-1)**n/(2*n+1) for n in range(m))
... for m in range(1,30)]
>>>
>>> T = shanks_table(S[:7])
>>> for row in T:
... print('['+', '.join(['{:.6g}'.format(v) for v in row])+']')
...
[-0.75]
[1.25, 3.16667]
[-1.75, 3.13333, -28.75]
[2.25, 3.14524, 82.25, 3.14234]
[-2.75, 3.13968, -177.75, 3.14139, -969.938]
[3.25, 3.14271, 327.25, 3.14166, 3515.06, 3.14161]
'''
if len(seq) < 2:
raise ValueError("seq should be of minimum length 2")
if table:
START = len(table)
else:
START = 0
table = []
STOP = len(seq) - 1
if STOP & 1:
STOP -= 1
one = 1.0#ctx.one
eps = np.spacing(1)#+ctx.eps
if randomized:
from random import Random
rnd = Random()
rnd.seed(START)
for i in range(START, STOP):
row = []
for j in range(i+1):
if j == 0:
a, b = 0, seq[i+1]-seq[i]
else:
if j == 1:
a = seq[i]
else:
a = table[i-1][j-2]
b = row[j-1] - table[i-1][j-1]
if not b:
if randomized:
b = rnd.getrandbits(10)*eps
elif i & 1:
return table[:-1]
else:
return table
row.append(a + one/b)
table.append(row)
return table
[docs]def shanks(seq, ind=0):
"""Iterated Shanks transformation to accelerate series convergence
Though normally applied to a 1d array, `shanks` will actually operate on
the last dimension of seq which allows for multi-dimensional arrays. e.g.
for 2d data each row of sequence whould be a separate sequence
Parameters
----------
seq : list or array
If seq is a numpy array then it's elements will be modified in-place.
If seq is a list then seq will not be modified.
ind : int, optional
Start index for extrapolation. Can be negative, e.g. ind=-5
will extrapolate based on the last 5 elements of the `seq`.
default ind=0 i.e. use all elements.
Returns
-------
out : array with 1 dim less than `seq`, or float if seq is only 1d.
Extrapolated value. If `seq` is a numpy array then due to in-place
modification the result will also be in seq[..., -1].
See also
--------
shanks_table : Copy of sympy.mpmath.calculus.extrapolation.shanks
Provides the whole epsilon table and error estimates.
numpy.apply_along_axis : If your sequence is not in the last dimension of
an array then use np.apply_along_axis to apply it along a specific
axis.
Notes
-----
I think this will also work on multi-dimensional data. The shanks
extrapolation will be performed on the last dimension of the data.
So for 2d data each row is a separate sequence.
For sequence:
.. math A=\\sum_{m=0}^{\\infty}a_m
The partial sum is first defined as:
.. math:: A_n=\\sum_{m=0}^{n}a_m
This forms a new sequence, the convergence of which can be sped up by
repeated use of:
.. math:: S(A_n)=\\frac{A_{n+1}A_{n-1}-A_n^2}{A_{n+1}-2A_n+A_{n-1}}
"""
seq = np.atleast_1d(seq)
if ind is None:
return +seq[..., -1]
if ind < 0:
ind = seq.shape[-1] + ind
ind = max(ind, 0)
for i in range(ind, seq.shape[-1] - 2, 2):
denom = (seq[..., i + 2:] - 2 * seq[..., i + 1: -1] + seq[..., i:-2])
if np.any(denom==0):
return +seq[..., -1]
seq[..., i + 2:] = (
(seq[..., i + 2:] * seq[..., i:-2] - seq[..., i + 1:-1]**2) /
denom)
return +seq[...,-1]
[docs]def gk_quad(f, a, b, args=(), n=10):
"""Integration by Gauss-Kronrod quadrature between intervals
Parameters
----------
f : function or method
Function to integrate.
a, b : 1d array
Limits of integration. Must have len(a)==len(b).
args : tuple, optional
`args` will be passed to f using f(x, *args). Default args=().
n : [7,10,15,20,25,30], optional
Number of gauss quadrature evaluation points. Default n=10. There will
be 2*n+1 Kronrod quadrature points.
Returns
-------
igral : 1d float array of size len(a)
Integral of f between a and b. Each value in igral corresponds to
the corresponding a-b interval. If you are just using a and b
to subdivide a larger interval then simply sum igral for the entire
integral.
err_estimate : 1d float array of size len(a)
Estimate of the error in the integral. i.e. fine integral minus
coarse integral.
"""
ai = np.atleast_1d(a)
bi = np.atleast_1d(b)
xj_, wj1, wj2 = gauss_kronrod_abscissae_and_weights(n)
# dim1 = each integration limits, a and b
# dim2 = each quadrature point
ai = ai[:, np.newaxis]
bi = bi[:, np.newaxis]
xj_ = xj_[np.newaxis, :]
wj1 = wj1[np.newaxis, :]
wj2 = wj2[np.newaxis, :]
bma = (bi - ai) / 2 # b minus a
bpa = (ai + bi) /2 # b plus a
xij = bma * xj_ + bpa # xj_ are in [-1, 1] so need to transform to [a, b]
fij = f(xij, *args)
# igral1 = np.ravel(bma) * np.sum(fij * wj1, axis=1)
# igral2 = np.ravel(bma) * np.sum(fij * wj2, axis=1)
igral1 = bma[:, 0] * np.sum(fij * wj1, axis=1)
igral2 = bma[:, 0] * np.sum(fij * wj2, axis=1)
err_estimate = igral2 - igral1
return igral2, err_estimate
[docs]def gl_quad(f, a, b, args=(), n=10, shanks_ind=False):
"""Integration by Gauss-Legendre quadrature with subdivided interval
Parameters
----------
f : function or method
function to integrate.
a, b : 1d array
limits of integration
args : tuple, optional
args will be passed to f using f(x, *args). default=()
n : [2-20, 32, 64, 100], optional
number of quadrature evaluation points. default=10
Returns
-------
igral : 1d float array of size len(a)
Integral of f between a and b. Each value in igral corresponds to
the corresponding a-b interval. If you are just using a and b
to subdivide a larger interval then simply sum igral for the entire
integral.
Notes
-----
Be careful when using large values of n.There may be precision issues.
"""
ai = np.atleast_1d(a)
bi = np.atleast_1d(b)
xj_, wj = gauss_legendre_abscissae_and_weights(n)
# dim1 = each integration limits, a and b
# dim2 = each quadrature point
ai = ai[:, np.newaxis]
bi = bi[:, np.newaxis]
xj_ = xj_[np.newaxis, :]
wj = wj[np.newaxis, :]
bma = (b - a) / 2 # b minus a
bpa = (a + b) /2 # b plus a
xij = bma * xj_ + bpa # xj_ are in [-1, 1] so need to transform to [a, b]
fij = f(xij, *args)
igral = np.sum(bma * fij * wj, axis=1)
return igral
if __name__ == '__main__':
import nose
nose.runmodule(argv=['nose', '--verbosity=3', '--with-doctest'])
# nose.runmodule(argv=['nose', '--verbosity=3'])
