Exercise 1



In [1]:

    
import numpy
from matplotlib import pyplot
% matplotlib inline



In [2]:

    
import os, sys
sys.path.append(os.path.split(os.path.split(os.getcwd())[0])[0])



In [3]:

    
import utils.quadrature as quad

(a)

Let $f = x^6$. And $\int_{-1}^{1} f dx = \frac{2}{7}$.

Use Gauss-Lobatto-Legendre quadrature to numerically evaluate $\int_{-1}^{1} f dx$ with number of quadrature points $Q=4,\ 5,\ 6$.



In [4]:

    
def f(x):
    """the integrand: x**6"""
    
    return x**6



In [5]:

    
print("The exact solution is: ", 2./7.)

for Qi in [4, 5, 6]:
    qd = quad.GaussLobattoJacobi(Qi)
    ans = qd(f, -1, 1)
    print("Q = {0}, ans = {1}, absolute error is {2}".format(Qi, ans, abs(ans-2./7.)))









    



The exact solution is:  0.2857142857142857
Q = 4, ans = 0.3466666666666666, absolute error is 0.060952380952380925
Q = 5, ans = 0.2857142857142857, absolute error is 0.0
Q = 6, ans = 0.2857142857142856, absolute error is 1.1102230246251565e-16

(b)

Let $f = x^6$. And $\int_{-1}^{2} f dx = \frac{129}{7}$.

Use Gauss-Lobatto-Legendre quadrature to numerically evaluate $\int_{-1}^{2} f dx$ with number of quadrature points $Q=4,\ 5,\ 6$.



In [6]:

    
print("The exact solution is: ", 129./7.)

for Qi in [4, 5, 6]:
    qd = quad.GaussLobattoJacobi(Qi)
    ans = qd(f, -1, 2)
    print("Q = {0}, ans = {1}, absolute error is {2}".format(Qi, ans, abs(ans-129./7.)))









    



The exact solution is:  18.428571428571427
Q = 4, ans = 19.469999999999995, absolute error is 1.0414285714285683
Q = 5, ans = 18.428571428571438, absolute error is 1.0658141036401503e-14
Q = 6, ans = 18.428571428571416, absolute error is 1.0658141036401503e-14

(c)

Let $f = \sin{x}$. And $\int_{0}^{\pi/2} f dx = 1$.

Use Gauss-Lobatto-Legendre quadrature to numerically evaluate $\int_{0}^{\pi/2} f dx$ with number of quadrature points $2 \le Q \le 8$ and plot the error to check the convergence order.



In [7]:

    
def f(x):
    """the integrand: x**6"""
    
    return numpy.sin(x)



In [8]:

    
print("The exact solution is: ", 1.)

Q = range(2, 9)
err = numpy.zeros_like(Q, dtype=numpy.float64)
for i, Qi in enumerate(Q):
    qd = quad.GaussLobattoJacobi(Qi)
    ans = qd(f, 0., numpy.pi/2.)
    err[i] = numpy.abs(ans - 1.)
    print("Q = {0}, ans = {1}, absolute error is {2}".format(Qi, ans, err[i]))









    



The exact solution is:  1.0
Q = 2, ans = 0.7853981633974483, absolute error is 0.21460183660255172
Q = 3, ans = 1.0022798774922104, absolute error is 0.0022798774922103693
Q = 4, ans = 0.9999891898309784, absolute error is 1.0810169021602256e-05
Q = 5, ans = 1.0000000284807518, absolute error is 2.8480751756987388e-08
Q = 6, ans = 0.9999999999525426, absolute error is 4.7457371366022016e-11
Q = 7, ans = 1.0000000000000555, absolute error is 5.551115123125783e-14
Q = 8, ans = 0.9999999999999978, absolute error is 2.220446049250313e-15



In [9]:

    
pyplot.semilogy(Q, err, 'k.-', lw=2, markersize=15)
pyplot.title(r"Error of numerical integration of $\int_{0}^{\pi/2}\sin{(x)} dx$" + 
             "\n with Gauss-Lobatto-Legendre quadrature", y=1.08)
pyplot.xlabel(r"$Q$" + ", number of quadratiure points")
pyplot.ylabel("absolute error")
pyplot.grid();