Num Methods: Numerical Quadrature

30 days ago by solin

Numerical Quadrature

This worksheet illustrates basic principles of numerical quadrature.  You can do a variety of things here:

  1. In an interval (a, b), calculate quadrature weights for a given list of quadrature points.
  2. Integrate a given function f(x) in interval (a, b) given quadrature points and weights.
  3. Integrate a given function f(x) in interval (a, b) using a Gauss rule with n points.
  4. Integrate a given function f(x) in interval (a, b) using an adaptive Gauss rule.
  5. Write Gauss quadrature rule with an arbitrary number of points.
  6. Play with a symbolic integration package (SymPy). See that not every function f(x) can be integrated on paper.

The following active input text box contains the main code. You need to click on it to see an "Evaluate" link appear at its bottom. Click on this link to have the code sent to our computer and interpreted. Then you can evaluate/modify/evaluate the examples in the other text boxes in the same way. Have fun and let us know at femhub@googlegroups.com with any problems or questions.

The Code:
# Import symbolics, linear algebra, Legendre # polynomials, and numerical quadrature functions. from numpy import arange, zeros from sympy import (pi, sin, cos, log, exp, tan, atan, asin, acos, sinh, cosh, tanh, asinh, acosh, sqrt) from sympy import integrate, var, lambdify from scipy.linalg import solve from scipy.special.orthogonal import legendre from scipy.special.orthogonal import p_roots from scipy.integrate import quadrature # we will be using the symbol 'x' for independent variable var("x") # Monomial basis. def mono(i, x): return float(x**i) # Define Legendre basis. def P(i, x): return legendre(i)(x) # Calculate weights in an interval (a, b) # for a given list of quadrature points def calculate_weights(a, b, quad_x): n = len(quad_x) rhs = zeros(n) # zero vector m = zeros((n, n)) # zero matrix for i in range(0, n): # loop over monomials x**i for j in range(0, n): # loop over points m[i][j] = float(quad_x[j])**i rhs[i] = float(integrate(x**i, (x, a, b))) print "Matrix: ", m print "RHS: ", rhs # solve the matrix problem weights = solve(m, rhs) print "Weights: ", weights return weights # reference mapping from (-1, 1) onto (a, b) def refmap(a, b, x): return (a+b)/2. + x*(b-a)/2. # Integrate a function f(x) in interval (a, b) using a list of quadrature # points quad_x and list of corresponding weights quad_w def num_quad_general(f, x, a, b, quad_x, quad_w): f_num = lambdify(x, f, modules=["math"]) val = 0 for i in range(len(quad_x)): val += f_num(quad_x[i]) * quad_w[i] print "Approx integral =", val val_exact = float(integrate(f, (x, a, b))) print "Exact integral =", val_exact return val # Integrate a function f(x) in interval (a, b) using fixed-order # Gauss quadrature with n points def num_quad_gauss(f, x, a, b, n): print "Preparing Gauss quadrature in interval (%g, %g)..." %(float(a), float(b)) points, weights = p_roots(n) quad_x = [float(refmap(a, b, xx)) for xx in points] quad_w = [float(w*(b - a)/2.) for w in weights] print "Points:", quad_x print "Weights:", quad_w val = num_quad_general(f, x, a, b, quad_x, quad_w) return val # Integrate a function f(x) in interval (a, b) using adaptive # Gauss quadrature with n points def num_quad_gauss_adapt(f, x, a, b, n, eps): val = num_quad_gauss(f, x, a, b, n) c = (a+b)/2. val_left = num_quad_gauss(f, x, a, c, n) val_right = num_quad_gauss(f, x, c, b, n) if abs(val - (val_left + val_right)) < eps: return float(val_left + val_right) val_left = num_quad_gauss_adapt(f, x, a, c, n, eps/2.) val_right = num_quad_gauss_adapt(f, x, c, b, n, eps/2.) val_approx = val_left + val_right val_exact = float(integrate(f, (x, a, b))) print "-----------------------------------------------------------------------" print "Interval (%g, %g): approximate %.15f, exact %.15f" %(a, b, val_approx, val_exact) print "-----------------------------------------------------------------------" return val_approx 
       

Task 1

In an interval (a, b), calculate quadrature weights for a given list of quadrature points.
# arbitrary points quad_w = calculate_weights(0, 1, [0, 0.5, 1]) 
        
Matrix:  [[ 1.    1.    1.  ]
 [ 0.    0.5   1.  ]
 [ 0.    0.25  1.  ]]
RHS:  [ 1.          0.5         0.33333333]
Weights:  [ 0.16666667  0.66666667  0.16666667]
# equidistant points # notice negative weights for larger n a = 0 b = 1. n = 11 h = (b-a)/(n-1) pts_list = [a + h*i for i in range(0, n)] print "Point list:", pts_list quad_w = calculate_weights(a, b, pts_list) 
       
# Chebyshev points # notice much better weights than with equidistant points a = 0. b = 1. n = 11 pts_list = [(a+b)/2 + float(cos(i*pi/(n-1)))*(b-a)/2 for i in range(0, n)] print "Point list:", pts_list quad_w = calculate_weights(a, b, pts_list) 
       

Task 2

Integrate a given function f(x) in interval (a, b) given quadrature points and weights.
a = 0 b = pi quad_x = ([0, pi/2, 1]) quad_w = calculate_weights(a, b, quad_x) val = num_quad_general(sin(x), x, a, b, quad_x, quad_w) 
       

Task 3

Integrate a given function f(x) in interval (a, b) using a Gauss rule.
val = num_quad_gauss(1/(1+x*x), x, 0, 5, 4) 
       

Task 4

Integrate a given function f(x) in interval (a, b) using an adaptive Gauss rule with n points.
val = num_quad_gauss_adapt(1/(1+x*x), x, 0, 5, 1, 1e-4) 
       
val = num_quad_gauss_adapt(sqrt(x), x, 0, 1, 3, 1e-4) 
        
WARNING: Output truncated!  
full_output.txt


Preparing Gauss quadrature in interval (0, 1)...
Points: [0.11270166537925841, 0.5, 0.8872983346207417]
Weights: [0.27777777777777779, 0.44444444444444442,
0.2777777777777779]
Approx integral = 0.669179633899
Exact integral = 0.666666666667
Preparing Gauss quadrature in interval (0, 0.5)...
Points: [0.056350832689629204, 0.25, 0.44364916731037085]
Weights: [0.1388888888888889, 0.22222222222222221,
0.13888888888888895]
Approx integral = 0.236590728481
Exact integral = 0.235702260396
Preparing Gauss quadrature in interval (0.5, 1)...
Points: [0.55635083268962915, 0.75, 0.94364916731037085]
Weights: [0.1388888888888889, 0.22222222222222221,
0.13888888888888895]
Approx integral = 0.430964722058
Exact integral = 0.430964406271
Preparing Gauss quadrature in interval (0, 0.5)...
Points: [0.056350832689629204, 0.25, 0.44364916731037085]
Weights: [0.1388888888888889, 0.22222222222222221,
0.13888888888888895]
Approx integral = 0.236590728481
Exact integral = 0.235702260396
Preparing Gauss quadrature in interval (0, 0.25)...
Points: [0.028175416344814602, 0.125, 0.22182458365518543]
Weights: [0.069444444444444448, 0.1111111111111111,
0.069444444444444475]
Approx integral = 0.0836474542374
Exact integral = 0.0833333333333
Preparing Gauss quadrature in interval (0.25, 0.5)...
Points: [0.27817541634481457, 0.375, 0.47182458365518543]
Weights: [0.069444444444444448, 0.1111111111111111,
0.069444444444444475]
Approx integral = 0.15236903871
Exact integral = 0.152368927062
Preparing Gauss quadrature in interval (0, 0.25)...
Points: [0.028175416344814602, 0.125, 0.22182458365518543]
Weights: [0.069444444444444448, 0.1111111111111111,
0.069444444444444475]
Approx integral = 0.0836474542374
Exact integral = 0.0833333333333
Preparing Gauss quadrature in interval (0, 0.125)...
Points: [0.014087708172407301, 0.0625, 0.11091229182759271]
Weights: [0.034722222222222224, 0.055555555555555552,
0.034722222222222238]
Approx integral = 0.0295738410601
Exact integral = 0.0294627825494
Preparing Gauss quadrature in interval (0.125, 0.25)...
Points: [0.13908770817240729, 0.1875, 0.23591229182759271]
Weights: [0.034722222222222224, 0.055555555555555552,
0.034722222222222238]
Approx integral = 0.0538705902572
Exact integral = 0.0538705507839
Preparing Gauss quadrature in interval (0, 0.125)...
Points: [0.014087708172407301, 0.0625, 0.11091229182759271]
Weights: [0.034722222222222224, 0.055555555555555552,
0.034722222222222238]
Approx integral = 0.0295738410601
Exact integral = 0.0294627825494
Preparing Gauss quadrature in interval (0, 0.0625)...
Points: [0.0070438540862036506, 0.03125, 0.055456145913796356]
Weights: [0.017361111111111112, 0.027777777777777776,
0.017361111111111119]
Approx integral = 0.0104559317797
Exact integral = 0.0104166666667
Preparing Gauss quadrature in interval (0.0625, 0.125)...
Points: [0.069543854086203644, 0.09375, 0.11795614591379636]
Weights: [0.017361111111111112, 0.027777777777777776,
0.017361111111111119]
Approx integral = 0.0190461298387
Exact integral = 0.0190461158828

...

Weights: [0.0086805555555555559, 0.013888888888888888,
0.0086805555555555594]
Approx integral = 0.0103261439761
Exact integral = 0.0103261439339
--------------------------------------------------------------------\
---
Interval (0, 0.125): approximate 0.029462859754258, exact
0.029462782549439
--------------------------------------------------------------------\
---
Preparing Gauss quadrature in interval (0.125, 0.25)...
Points: [0.13908770817240729, 0.1875, 0.23591229182759271]
Weights: [0.034722222222222224, 0.055555555555555552,
0.034722222222222238]
Approx integral = 0.0538705902572
Exact integral = 0.0538705507839
Preparing Gauss quadrature in interval (0.125, 0.1875)...
Points: [0.13204385408620364, 0.15625, 0.18045614591379636]
Weights: [0.017361111111111112, 0.027777777777777776,
0.017361111111111119]
Approx integral = 0.0246638059618
Exact integral = 0.0246638051871
Preparing Gauss quadrature in interval (0.1875, 0.25)...
Points: [0.19454385408620364, 0.21875, 0.24295614591379636]
Weights: [0.017361111111111112, 0.027777777777777776,
0.017361111111111119]
Approx integral = 0.0292067457159
Exact integral = 0.0292067455968
--------------------------------------------------------------------\
---
Interval (0, 0.25): approximate 0.083333411431963, exact
0.083333333333333
--------------------------------------------------------------------\
---
Preparing Gauss quadrature in interval (0.25, 0.5)...
Points: [0.27817541634481457, 0.375, 0.47182458365518543]
Weights: [0.069444444444444448, 0.1111111111111111,
0.069444444444444475]
Approx integral = 0.15236903871
Exact integral = 0.152368927062
Preparing Gauss quadrature in interval (0.25, 0.375)...
Points: [0.26408770817240729, 0.3125, 0.36091229182759271]
Weights: [0.034722222222222224, 0.055555555555555552,
0.034722222222222238]
Approx integral = 0.0697597777818
Exact integral = 0.0697597755906
Preparing Gauss quadrature in interval (0.375, 0.5)...
Points: [0.38908770817240729, 0.4375, 0.48591229182759271]
Weights: [0.034722222222222224, 0.055555555555555552,
0.034722222222222238]
Approx integral = 0.0826091518085
Exact integral = 0.0826091514716
--------------------------------------------------------------------\
---
Interval (0, 0.5): approximate 0.235702341022225, exact
0.235702260395516
--------------------------------------------------------------------\
---
Preparing Gauss quadrature in interval (0.5, 1)...
Points: [0.55635083268962915, 0.75, 0.94364916731037085]
Weights: [0.1388888888888889, 0.22222222222222221,
0.13888888888888895]
Approx integral = 0.430964722058
Exact integral = 0.430964406271
Preparing Gauss quadrature in interval (0.5, 0.75)...
Points: [0.52817541634481457, 0.625, 0.72182458365518543]
Weights: [0.069444444444444448, 0.1111111111111111,
0.069444444444444475]
Approx integral = 0.197310447694
Exact integral = 0.197310441497
Preparing Gauss quadrature in interval (0.75, 1)...
Points: [0.77817541634481457, 0.875, 0.97182458365518543]
Weights: [0.069444444444444448, 0.1111111111111111,
0.069444444444444475]
Approx integral = 0.233653965727
Exact integral = 0.233653964774
--------------------------------------------------------------------\
---
Interval (0, 1): approximate 0.666666754443864, exact
0.666666666666667
--------------------------------------------------------------------\
---

Task 5

Write Gauss quadrature rule with an arbitrary number of points.
# This script prints Gauss quadrature points and # weights for a given number n of points. from scipy.special.orthogonal import p_roots n = 5 print "n =", n print "orders: %d, %d" % (2*n - 2, 2*n - 1) roots, weights = p_roots(n) print "{" for root, weight in zip(roots, weights): print " { %.15f, %.15f }," % (root, weight) print "};" 
        
n = 5
orders: 8, 9
{
  { -0.906179845938664, 0.236926885056189 },
  { -0.538469310105683, 0.478628670499367 },
  { -0.000000000000000, 0.568888888888889 },
  { 0.538469310105683, 0.478628670499367 },
  { 0.906179845938664, 0.236926885056189 },
};

Task 6

Use SymPy to integrate some functions "on paper". See where the limits of symbolic integration are.
integrate(x**2, x) 
       
var("n") integrate(x**n, x) 
       
from sympy import log integrate(cos(log(x))/x) 
       
from sympy import log integrate(cos(log(x**2))/x) 
       
from sympy import exp integrate(exp(x**2))