GaussQuadrature.mw

Gaussian Quadrature Rules 

2-Point Rule 

We develop Gaussian quadrature rules for the integral int(f(x), x = -1 .. 1) on the standard interval [-1, 1]. For a general interval of integration [a, b] one can perform a simple transformation of variables. 

 

For the 2-point rule, we must choose evaluation points x1, x2 and corresponding weights w1, w1 such that the quadrature rule 

 

> `:=`(Q2, `+`(`*`(w1, `*`(f(x1))), `*`(w2, `*`(f(x2)))))
 

`+`(`*`(w1, `*`(f(x1))), `*`(w2, `*`(f(x2)))) (1)
 

will be exact for integrating polynomials of as high a degree as possible. 

 

Noting that {1, x, `*`(`^`(x, 2)), `*`(`^`(x, 3))} forms a basis for all polynomials of degree 3, we set up four (nonlinear) equations in the four unknowns x1, x2, w1, w2 as follows. 

 

> `:=`(f, proc (x) options operator, arrow; 1 end proc)
 

proc (x) options operator, arrow; 1 end proc (2)
 

> `:=`(eqn1, Q2 = int(f(x), x = -1 .. 1))
 

`+`(w1, w2) = 2 (3)
 

> `:=`(f, proc (x) options operator, arrow; x end proc)
 

proc (x) options operator, arrow; x end proc (4)
 

> `:=`(eqn2, Q2 = int(f(x), x = -1 .. 1))
 

`+`(`*`(w1, `*`(x1)), `*`(w2, `*`(x2))) = 0 (5)
 

> `:=`(f, proc (x) options operator, arrow; `*`(`^`(x, 2)) end proc)
 

proc (x) options operator, arrow; `*`(`^`(x, 2)) end proc (6)
 

> `:=`(eqn3, Q2 = int(f(x), x = -1 .. 1))
 

`+`(`*`(w1, `*`(`^`(x1, 2))), `*`(w2, `*`(`^`(x2, 2)))) = `/`(2, 3) (7)
 

> `:=`(f, proc (x) options operator, arrow; `*`(`^`(x, 3)) end proc)
 

proc (x) options operator, arrow; `*`(`^`(x, 3)) end proc (8)
 

> `:=`(eqn4, Q2 = int(f(x), x = -1 .. 1))
 

`+`(`*`(w1, `*`(`^`(x1, 3))), `*`(w2, `*`(`^`(x2, 3)))) = 0 (9)
 

> `:=`(_EnvExplicit, true)
 

true (10)
 

> solve({eqn1, eqn2, eqn3, eqn4}, {w1, w2, x1, x2})
 

{w1 = 1, w2 = 1, x1 = `+`(`*`(`/`(1, 3), `*`(`^`(3, `/`(1, 2))))), x2 = `+`(`-`(`*`(`/`(1, 3), `*`(`^`(3, `/`(1, 2))))))}, {w1 = 1, w2 = 1, x1 = `+`(`-`(`*`(`/`(1, 3), `*`(`^`(3, `/`(1, 2)))))), x2 = ... (11)
 

>
 

 

The solution we are interested in is the one with x1, x2 in order from left to right. 

 

Therefore the 2-point Gaussian quadrature rule for a general function f(x) is as follows. 

 

> `:=`(f, 'f')
 

f (12)
 

> `:=`(x1, x2, `+`(`-`(`*`(`/`(1, 3), `*`(sqrt(3))))), `+`(`*`(`/`(1, 3), `*`(sqrt(3)))))
 

`+`(`-`(`*`(`/`(1, 3), `*`(`^`(3, `/`(1, 2)))))), `+`(`*`(`/`(1, 3), `*`(`^`(3, `/`(1, 2))))) (13)
 

> `:=`(w1, w2, 1, 1)
 

1, 1 (14)
 

> Q2
 

`+`(f(`+`(`-`(`*`(`/`(1, 3), `*`(`^`(3, `/`(1, 2))))))), f(`+`(`*`(`/`(1, 3), `*`(`^`(3, `/`(1, 2))))))) (15)
 

>
 

3-Point Rule 

 

> restart
 

> `:=`(Q3, `+`(`*`(w1, `*`(f(x1))), `*`(w2, `*`(f(x2))), `*`(w3, `*`(f(x3)))))
 

`+`(`*`(w1, `*`(f(x1))), `*`(w2, `*`(f(x2))), `*`(w3, `*`(f(x3)))) (16)
 

> `:=`(f, proc (x) options operator, arrow; 1 end proc)
 

proc (x) options operator, arrow; 1 end proc (17)
 

> `:=`(eqn1, Q3 = int(f(x), x = -1 .. 1))
 

`+`(w1, w2, w3) = 2 (18)
 

> `:=`(f, proc (x) options operator, arrow; x end proc)
 

proc (x) options operator, arrow; x end proc (19)
 

> `:=`(eqn2, Q3 = int(f(x), x = -1 .. 1))
 

`+`(`*`(w1, `*`(x1)), `*`(w2, `*`(x2)), `*`(w3, `*`(x3))) = 0 (20)
 

> `:=`(f, proc (x) options operator, arrow; `*`(`^`(x, 2)) end proc)
 

proc (x) options operator, arrow; `*`(`^`(x, 2)) end proc (21)
 

> `:=`(eqn3, Q3 = int(f(x), x = -1 .. 1))
 

`+`(`*`(w1, `*`(`^`(x1, 2))), `*`(w2, `*`(`^`(x2, 2))), `*`(w3, `*`(`^`(x3, 2)))) = `/`(2, 3) (22)
 

> `:=`(f, proc (x) options operator, arrow; `*`(`^`(x, 3)) end proc)
 

proc (x) options operator, arrow; `*`(`^`(x, 3)) end proc (23)
 

> `:=`(eqn4, Q3 = int(f(x), x = -1 .. 1))
 

`+`(`*`(w1, `*`(`^`(x1, 3))), `*`(w2, `*`(`^`(x2, 3))), `*`(w3, `*`(`^`(x3, 3)))) = 0 (24)
 

> `:=`(f, proc (x) options operator, arrow; `*`(`^`(x, 4)) end proc)
 

proc (x) options operator, arrow; `*`(`^`(x, 4)) end proc (25)
 

> `:=`(eqn5, Q3 = int(f(x), x = -1 .. 1))
 

`+`(`*`(w1, `*`(`^`(x1, 4))), `*`(w2, `*`(`^`(x2, 4))), `*`(w3, `*`(`^`(x3, 4)))) = `/`(2, 5) (26)
 

> `:=`(f, proc (x) options operator, arrow; `*`(`^`(x, 5)) end proc)
 

proc (x) options operator, arrow; `*`(`^`(x, 5)) end proc (27)
 

> `:=`(eqn6, Q3 = int(f(x), x = -1 .. 1))
 

`+`(`*`(w1, `*`(`^`(x1, 5))), `*`(w2, `*`(`^`(x2, 5))), `*`(w3, `*`(`^`(x3, 5)))) = 0 (28)
 

> `:=`(_EnvExplicit, true)
 

true (29)
 

> `:=`(soln, solve({eqn1, eqn2, eqn3, eqn4, eqn5, eqn6}, {w1, w2, w3, x1, x2, x3}))
 

{w1 = `/`(8, 9), w2 = `/`(5, 9), w3 = `/`(5, 9), x1 = 0, x2 = `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))), x3 = `+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))))}, {w1 = `/`(8, 9), w2 = `/`(5, 9), w3...
{w1 = `/`(8, 9), w2 = `/`(5, 9), w3 = `/`(5, 9), x1 = 0, x2 = `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))), x3 = `+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))))}, {w1 = `/`(8, 9), w2 = `/`(5, 9), w3...
{w1 = `/`(8, 9), w2 = `/`(5, 9), w3 = `/`(5, 9), x1 = 0, x2 = `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))), x3 = `+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))))}, {w1 = `/`(8, 9), w2 = `/`(5, 9), w3...
{w1 = `/`(8, 9), w2 = `/`(5, 9), w3 = `/`(5, 9), x1 = 0, x2 = `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))), x3 = `+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))))}, {w1 = `/`(8, 9), w2 = `/`(5, 9), w3...
{w1 = `/`(8, 9), w2 = `/`(5, 9), w3 = `/`(5, 9), x1 = 0, x2 = `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))), x3 = `+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))))}, {w1 = `/`(8, 9), w2 = `/`(5, 9), w3...
(30)
 

> nops([soln])
 

6 (31)
 

> soln[1]
 

{w1 = `/`(8, 9), w2 = `/`(5, 9), w3 = `/`(5, 9), x1 = 0, x2 = `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))), x3 = `+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))))} (32)
 

> soln[2]
 

{w1 = `/`(8, 9), w2 = `/`(5, 9), w3 = `/`(5, 9), x1 = 0, x2 = `+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2)))))), x3 = `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2)))))} (33)
 

> soln[3]
 

{w1 = `/`(5, 9), w2 = `/`(8, 9), w3 = `/`(5, 9), x1 = `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))), x2 = 0, x3 = `+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))))} (34)
 

> soln[4]
 

{w1 = `/`(5, 9), w2 = `/`(8, 9), w3 = `/`(5, 9), x1 = `+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2)))))), x2 = 0, x3 = `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2)))))} (35)
 

> soln[5]
 

{w1 = `/`(5, 9), w2 = `/`(5, 9), w3 = `/`(8, 9), x1 = `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))), x2 = `+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2)))))), x3 = 0} (36)
 

> soln[6]
 

{w1 = `/`(5, 9), w2 = `/`(5, 9), w3 = `/`(8, 9), x1 = `+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2)))))), x2 = `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))), x3 = 0} (37)
 

>
 

 

The solution we are interested in is the one with x1, x2, x3 in order from left to right, namely 

 

> soln[4]
 

{w1 = `/`(5, 9), w2 = `/`(8, 9), w3 = `/`(5, 9), x1 = `+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2)))))), x2 = 0, x3 = `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2)))))} (38)
 

 

Therefore the 3-point Gaussian quadrature rule for a general function f(x) is as follows. 

 

> `:=`(f, 'f')
 

f (39)
 

> `:=`(x1, x2, x3, `+`(`-`(`*`(`/`(1, 5), `*`(sqrt(15))))), 0, `+`(`*`(`/`(1, 5), `*`(sqrt(15)))))
 

`+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2)))))), 0, `+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))) (40)
 

> `:=`(w1, w2, w3, `/`(5, 9), `/`(8, 9), `/`(5, 9))
 

`/`(5, 9), `/`(8, 9), `/`(5, 9) (41)
 

> Q3
 

`+`(`*`(`/`(5, 9), `*`(f(`+`(`-`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))))))), `*`(`/`(8, 9), `*`(f(0))), `*`(`/`(5, 9), `*`(f(`+`(`*`(`/`(1, 5), `*`(`^`(15, `/`(1, 2))))))))) (42)
 

>