In [1]:
from __future__ import print_function, division
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
from sympy.core import S, Dummy, pi
from sympy.functions.combinatorial.factorials import factorial
from sympy.functions.elementary.trigonometric import sin, cos
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.functions.special.gamma_functions import gamma
from sympy.polys.orthopolys import legendre_poly, laguerre_poly, hermite_poly, jacobi_poly
from sympy.polys.rootoftools import RootOf
In [2]:
def gauss_lobatto(n, n_digits):
r"""
Computes the Gauss-Legendre quadrature [1]_ points and weights.
The Gauss-Lobatto quadrature approximates the integral:
.. math::
\int_{-1}^1 f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)
The nodes `x_i` of an order `n` quadrature rule are the roots of `P'_(n-1)`
and the weights `w_i` are given by:
.. math::
w_i = \frac{2}{n(n-1) \left[P_{n-1}(x_i)\right]^2},\quad x\neq\pm 1
w_i = \frac{2}{n(n-1)},\quad x=\pm 1
Parameters
==========
n : the order of quadrature
n_digits : number of significant digits of the points and weights to return
Returns
=======
(x, w) : the ``x`` and ``w`` are lists of points and weights as Floats.
The points `x_i` and weights `w_i` are returned as ``(x, w)``
tuple of lists.
Examples
========
>>> #from sympy.integrals.quadrature import gauss_lobatto
>>> x, w = gauss_lobatto(3, 5)
>>> x
[-1, 0, 1]
>>> w
[0.33333, 1.3333, 0.33333]
>>> x, w = gauss_lobatto(4, 5)
>>> x
[-1, -0.44721, 0.44721, 1]
>>> w
[0.16667, 0.83333, 0.83333, 0.16667]
See Also
========
gauss_legendre,gauss_laguerre, gauss_gen_laguerre, gauss_hermite, gauss_chebyshev_t,\
gauss_chebyshev_u, gauss_jacobi
References
==========
.. [1] http://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
.. [2] http://people.math.sfu.ca/~cbm/aands/page_888.htm
"""
x = Dummy("x")
p = legendre_poly(n-1, x, polys=True)
pd = p.diff(x)
xi = []
w = []
for r in pd.real_roots():
if isinstance(r, RootOf):
r = r.eval_rational(S(1)/10**(n_digits+2))
xi.append(r.n(n_digits))
w.append((2/(n*(n-1) * p.subs(x, r)**2)).n(n_digits))
xi.insert(0, -1)
xi.append(1)
w.insert(0,(S(2)/(n*(n-1))).n(n_digits))
w.append((S(2)/(n*(n-1))).n(n_digits))
return xi, w
In [3]:
%matplotlib notebook
plt.style.use("ggplot")
In [4]:
y = sym.symbols('y')
Gauss-Lobatto quadrature can integrate exactly polynomials up to $2n-3$, where $n$ is the number of points used
In [5]:
poly = 4*y**4 + 3*y**3 + y - 1
sym.plot(poly, (y, -1, 1));
In [6]:
print(sym.integrate(poly,(y,-1,1)))
print(sym.integrate(poly,(y,-1,1)).n(8))
In [7]:
for k in range(2,10):
xi, w = gauss_lobatto(k, 8)
poly_int = sum(w[j]*poly.subs(y, xi[j]) for j in range(k))
print(k, poly_int)
In [8]:
poly = 2*y**5 - 2*y**4 + 3*y**3 - 4*y**2 + 5*y - 1
sym.plot(poly, (y, -1, 1));
In [9]:
print(sym.integrate(poly,(y,-1,1)))
print(sym.integrate(poly,(y,-1,1)).n(8))
In [10]:
for k in range(2,10):
xi, w = gauss_lobatto(k, 8)
poly_int = sum(w[j]*poly.subs(y, xi[j]) for j in range(k))
print(k, poly_int)
In [11]:
fun = y*sym.sin(5*y) + y**2
sym.plot(fun, (y, -1, 1));
In [12]:
fun_int = sym.integrate(fun,(y,-1,1))
print(fun_int)
print(fun_int.n(8))
In [13]:
print("{}\t{}\t{}".format('n', 'Integral', 'Relative error'))
print("===========================================")
for k in range(2,20):
xi, w = gauss_lobatto(k, 8)
poly_int = sum(w[j]*fun.subs(y, xi[j]) for j in range(k))
print("{}\t{}\t{}".format(k, float(poly_int), float((fun_int - poly_int)/poly_int)))
In [14]:
norder = 10
min_cheb = []
max_cheb = []
min_lob = []
max_lob = []
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
x_cheb = []
x_lob = []
y = []
for n in range(1, norder + 1):
x_ch = [np.cos((2*i - 1)/(2*n)*np.pi) for i in range(n, 0, -1)]
x_ch = [-1] + x_ch + [1]
x_cheb += x_ch
x_lo, _ = gauss_lobatto(n + 2, 16)
x_lob += x_lo
dx_cheb = [x_ch[k + 1] - x_ch[k] for k in range(0, n + 1)]
dx_lob = [x_lo[k + 1] - x_lo[k] for k in range(0, n + 1)]
min_cheb.append(min(dx_cheb))
max_cheb.append(max(dx_cheb))
min_lob.append(min(dx_lob))
max_lob.append(max(dx_lob))
y += [n for i in range(n + 2)]
plt.plot(x_cheb, y, marker='s', lw=0)
plt.plot(x_lob, y, marker='v', lw=0)
#%%
plt.xlim(-1.2, 1.2)
plt.ylim(0, norder)
plt.subplot(1, 2, 2)
plt.plot(range(2, norder + 2), min_cheb, label="Min distance Chebyshev",
marker="s", lw=0)
plt.plot(range(2, norder + 2), max_cheb, label="Max distance Chebyshev",
marker="o", lw=0)
plt.plot(range(2, norder + 2), min_lob, label="Min distance Lobatto",
marker="v", lw=0)
plt.plot(range(2, norder + 2), max_lob, label="Max distance Lobatto",
marker="^", lw=0)
plt.legend(loc="best")
plt.show()
In [ ]:
In [15]:
from IPython.core.display import HTML
def css_styling():
styles = open('./styles/custom_barba.css', 'r').read()
return HTML(styles)
css_styling()
Out[15]:
In [ ]: