Section 5.2 showed how to obtain exact values for definite integrals of some simple functions (low-degree polynomials) by using areas of rectangles. For functions with no closed-form antiderivative, the rectangle method typically produces an approximate value of the definite integral—the more rectangles, the better the approximation.
For example, suppose you wanted to evaluate
\[\int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx \nonumber \]
with the rectangle method. This means finding the area of the shaded region in the figure on the right. Computers have eliminated the need to do these sorts of calculations by hand. Though the rectangle method is simple to implement in a traditional programming language (e.g. via a looping construct), there are easier ways in a domain-specific language (DSL) geared toward scientific computing. One such DSL is MATLAB, or its free open-source clone Octave.10
Implementing the rectangle method from scratch in Octave is a one-liner. For example, suppose you divide the interval \(\ival{0}{\sqrt{\pi}}\) into 100,000 (\(10^5\)) subintervals of equal length, producing 100,001 (1 + 1e5 in scientific notation) equally spaced points in \(\ival{0}{\sqrt{\pi}}\) (including \(0\) and \(\sqrt{\pi}\,\)). First use the left endpoints of the \(10^5\) subintervals:
octave> sum(sin(linspace(0,sqrt(pi),1+1e5)(1:end-1).^2)*sqrt(pi)/1e5)
ans = 0.8948314693913354
Now use the right endpoints:
octave> sum(sin(linspace(0,sqrt(pi),1+1e5)(2:end).^2)*sqrt(pi)/1e5)
ans = 0.8948314693913354
Finally, use the midpoints of the subintervals:
octave> sum(sin((linspace(0,sqrt(pi),1+1e5)(1:end-1)+sqrt(pi)/2e5).^2)*
sqrt(pi)/1e5)
ans = 0.8948314695305527
The true value of the integral up to 15 decimal places is 0.894831469484145, so all three approximations are accurate to 9 decimal places. The syntax in the above commands can be explained with some examples. The following command creates 4 equally spaced points in the interval \(\ival{1}{7}\) (including \(x=1\) and \(x=7\)), thus dividing \(\ival{1}{7}\) into 3 subintervals each of length \((7-1)/3 = 2\):
octave> linspace(1,7,4)
ans =
1 3 5 7
Get all but the last number in the above list:11
octave> linspace(1,7,4)(1:end-1)
ans =
1 3 5
Now square each number in that list of numbers (the dot before the exponentiation operator applies the squaring operation 2 element-wise in the list):
octave> linspace(1,7,4)(1:end-1).^2
ans =
1 9 25
Now take the sine of each of those squared numbers (measured in radians):
octave> sin(linspace(1,7,4)(1:end-1).^2)
ans =
0.8414709848078965 0.4121184852417566 -0.132351750097773
Now multiply each of those numbers (the heights of the rectangles) by the width (2) of each rectangle, then add up those areas:
octave> sum(sin(linspace(1,7,4)(1:end-1).^2)*2)
ans = 2.24247543990376
Before the advent of modern computing, the rectangle method was considered inefficient, and so alternative methods were created. Two such methods are the trapezoid rule and Simpson’s rule. The idea behind both methods is to take advantage of a nonlinear function’s changing slope by using nonrectangular regions. For the trapezoid rule those regions are trapezoids, while Simpson’s rule uses quasi-rectangular regions whose top edges are parabolas, as shown in Figure [fig:nummethods]:
For a partition \(P= \lbrace a=x_0 < x_1 < \cdots < x_{n-1} < x_n = b \rbrace\) of an interval \(\ival{a}{b}\) into \(n \ge 1\) subintervals of equal width \(h = (b-a)/n\), let \(y_i = f(x_i)\) for \(i = 0\), \(1\), \(\ldots\), \(n\). The trapezoid rule adds up the areas of trapezoids on each subinterval \(\ival{x_i}{x_{i+1}}\), with the top edge being the line segment joining the points \((x_i,y_i)\) and \((x_{i+1},y_{i+1})\). The approximation formula is straightforward to derive, based on areas of trapezoids:
Simpson’s rule depends on pairs of neighboring subintervals: \(\ival{x_0}{x_1}\) and \(\ival{x_1}{x_2}\), \(\ival{x_2}{x_3}\) and \(\ival{x_3}{x_4}\), \(\ldots\), \(\ival{x_{n-2}}{x_{n-1}}\) and \(\ival{x_{n-1}}{x_n}\). Thus, \(n \ge 2\) must be even. The top edge of the region over each pair \(\ival{x_{i}}{x_{i+1}}\) and \(\ival{x_{i+1}}{x_{i+2}}\) is the unique parabola joining the 3 points \((x_i,y_i)\), \((x_{i+1},y_{i+1})\), and \((x_{i+2},y_{i+2})\). The approximation formula is then:12
Example \(\PageIndex{1}\): trapsimp
Add text here.
Solution
Approximate the value of \(~\displaystyle\int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx~\) by using the trapezoid rule and Simpson’s rule with \(n=10^5\) subintervals.
Solution: Since \(x_0 = 0\) and \(x_n = \sqrt{\pi}\), then \(y_0 = \sin\,(x_0^2) = \sin\,0 = 0\) and \(y_n = \sin\,(x_n^2) = \sin\,\pi = 0\). Thus, \(y_0\) and \(y_n\) contribute nothing to the summation formulas for both rules. In particular the trapezoid rule approximation becomes
\[\int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx ~\approx~ \frac{h}{2}\,(0 ~+~ 2y_1 ~+~ 2y_2 ~+~ \cdots ~+~ 2y_{n-1} ~+~ 0) ~=~ h\,\cdot\,\left(\sum_{k=1}^{n-1} y_k\right) \nonumber \]
which is simple to implement in Octave:
octave> x = linspace(0,sqrt(pi),1+1e5);
octave> h = sqrt(pi)/1e5;
octave> h*sum(sin(x(2:end-1).^2))
ans = 0.8948314693913405
Likewise, the Simpson’s rule approximation becomes
\[\int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx ~\approx~ \frac{h}{3}\,(4y_1 ~+~ 2y_2 ~+~ 4y_3 ~+~ 2y_4 ~+~ \cdots ~+~ ~+~ 2y_{n-2} ~+~ 4y_{n-1}) \nonumber \]
which can be implemented easily by using Octave’s powerful indexing features:
octave> (h/3)*(4*sum(sin(x(2:2:end-1).^2)) + 2*sum(sin(x(3:2:end-1).^2)))
ans = 0.8948314694841457
In the above command, the statement x(3:2:end-1) allows you to skip every other element in the list x after position 3, by moving up the list in increments of 2 positions all the way to the next-to-last position in the list (end-1). Similarly for x(2:2:end-1), which starts at position 2 and then moves up in increments of 2.
Note that Simpson’s rule gives essentially the true value in this case, and the value from the trapezoid rule is virtually the same as the value produced by the built-in trapz function in Octave/MATLAB:
octave> trapz(x,sin(x.^2))
ans = 0.8948314693913402
In general you are better off using these sorts of built-in functions instead of implementing your own.
Typically Simpson’s rule is slightly more efficient than the trapezoid rule, which is slightly more efficient than the rectangle method. However, in the above examples all the approximations were accurate to at least 9 decimal places (equivalent to getting the distance between Detroit and Chicago correct within the thickness of a toothpick). The running time of each calculation was only a few thousandths of a second. Modern computing has generally made the efficiency differences negligible. Notice that the approximations in the rectangle method, the trapezoid rule and Simpson’s rule can all be written as linear combinations of function values \(f(a_i)\) multiplied by “weights” \(w_i\):
\[\int_a^b f(x)~\dx ~\approx~ \sum_{i=0}^{n} w_i\,f(a_i) \nonumber \]
For example, the weights in Simpson’s rule are \(w_i = \frac{h}{3}\), \(\frac{2h}{3}\), or \(\frac{4h}{3}\), depending on the points \(a_i\) in the interval \(\ival{a}{b}\). The method of Gaussian quadrature transforms an integral over any interval \(\ival{a}{b}\) into an integral over the specific interval \(\ival{-1}{1}\) and then uses a standard set of points in \(\ival{-1}{1}\) and known weights for those points:13
[tbl:gaussquad]
Example \(\PageIndex{1}\): gaussquad
Add text here.
Solution
Approximate the value of \(~\displaystyle\int_0^2 \dfrac{\dx}{1 + x^3}~\) by using Gaussian quadrature with \(n=4\) points.
Solution: For \(a=0\) and \(b=2\), use the substitution \(u = \frac{1}{b-a}(2x - a - b) = x-1\), so that \(x=u+1\) and \(\dx = \du\). Thus, \(g(u) = f(u+1) = \frac{1}{1 + (u+1)^3}\). Using \(n=4\) in Table [tbl:gaussquad], the points \(a_i\) and weights \(w_i\) are
\[\begin{aligned} {3} a_1 ~&=~ -0.339981 \quad\quad & w_1 ~&=~ 0.652145\\ a_2 ~&=~ 0.339981 \quad\quad & w_2 ~&=~ 0.652145\\ a_3 ~&=~ -0.861136 \quad\quad & w_3 ~&=~ 0.347855\\ a_4 ~&=~ 0.861136 \quad\quad & w_4 ~&=~ 0.347855\end{aligned} \nonumber \]
and so
\[\begin{aligned} \int_0^2 \frac{\dx}{1 + x^3} ~&=~ \frac{2-0}{2}\,\int_{-1}^1 g(u)~\du ~=~ \int_{-1}^1 \frac{\du}{1 + (u+1)^3} ~\approx~ \sum_{i=1}^4 w_i\,g(a_i) ~=~ \sum_{i=1}^4 w_i\,\cdot\,\frac{1}{1 + (a_i+1)^3}\
\[6pt] &\approx~ \frac{0.652145}{1 + (-0.339981+1)^3} ~+~ \frac{0.652145}{1 + (0.339981+1)^3} ~+~ \frac{0.347855}{1 + (-0.861136+1)^3} ~+~ \frac{0.347855}{1 + (0.861136+1)^3}\
\[4pt] &\approx~ 1.091621\end{aligned} \nonumber \]
The true value of the integral to six decimal places is 1.090002.
Using more points (e.g. \(n=7\)) is easy to implement in Octave, using element-wise operations on arrays:
octave> a = [ 0 -0.405845 0.405845 -0.741531 0.741531 -0.949108 0.949108 ];
octave> w = [ 0.417959 0.381830 0.381830 0.279705 0.279705 0.129485 0.129485 ];
octave> sum(w./(1 + (a+1).^3))
ans = 1.090016688064804
Gaussian quadrature can be applied to improper integrals. For example,
\[\int_0^{\infty} f(x)\,e^{-x}\,\dx ~\approx~ \sum_{i=1}^n w_i\,f(a_i) \nonumber \]
using the points \(a_i\) and weights \(w_i\) in Table [tbl:gaussquadimp]14 for \(n=3, 4\), or 5 points in \(\lival{0}{\infty}\):
[tbl:gaussquadimp]
Example \(\PageIndex{1}\): gqinf
Add text here.
Solution
Approximate the value of \(~\displaystyle\int_0^{\infty} x^5\,e^{-x}\,\dx~\) by using Gaussian quadrature with \(n=3\) points in Table [tbl:gaussquadimp].
Solution: For \(n=3\), Table [tbl:gaussquadimp] gives \(a_1=0.415775\), \(a_2=2.294280\), \(a_3=6.289945\), and \(w_1=0.711093\), \(w_2=0.278518\), \(w_3=0.010389\). Then for \(f(x)=x^5\),
\[\begin{aligned} \int_0^{\infty} x^5\,e^{-x}\,\dx ~&\approx~ \sum_{i=1}^n w_i\,f(a_i) ~=~ \sum_{i=1}^n w_i\,a_i^5\
\[4pt] &\approx~ 0.711093\,(0.415775)^5 ~+~ 0.278518\,(2.294280)^5 ~+~ 0.010389\,(6.289945)^5\\ &\approx~ 119.9974709727211\end{aligned} \nonumber \]
The true value is \(\Gamma\,(6) = 5! = 120\).
Note: The points \(a_i\) in Table [tbl:gaussquadimp] are the roots of the Laguerre polynomials of degree \(n\).
[sec6dot6]
[exer:pendper] A simple pendulum of length \(l\) swings through an angle of \(90\Degrees\) on either side of the vertical with period \(P\), given by
\[P ~=~ 4\,\sqrt{\dfrac{l}{g}}\,\int_0^{\pi/2} \frac{\dtheta}{\sqrt{1 \;-\; 0.5\,\sin^2 \theta}} \nonumber \]
where \(g = 9.8\) m/s2 is the acceleration due to gravity. Use the rectangle method (with left endpoints), the trapezoid rule, and Simpson’s rule to write \(P\) as a constant multiple of \(\sqrt{l/g}\). Preferably, use a computer and \(n= 10^5\) subintervals of equal width (or \(n=10\) subintervals if calculating by hand).
Repeat Exercise [exer:pendper] using Gaussian quadrature with \(n=5\) points.
[exer:gqsinx2] Approximate the value of \(~\displaystyle\int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx~\) by using Gaussian quadrature with \(n=7\) points.
2
Repeat Exercise [exer:gqsinx2] with \(n=9\) points.
Repeat Exercise [exer:gqsinx2] with \(n=10\) points.
The points \(a_i\) in Table [tbl:gaussquad] for Gaussian quadrature are the roots of the Legendre polynomials \(P_n(x)\), with \(P_0(x) = 1\), \(P_1(x) = x\), and \(P_n(x)\) defined for integers \(n \ge 2\) by the recursion formula
\[n\,P_{n}(x) ~=~ (2n-1)\,x\,P_{n-1}(x) ~-~ (n-1)\,P_{n-2}\,(x) \nonumber \]
- Write out \(P_n(x)\) explicitly in standard polynomial form for \(n = 2, 3, 4, 5\).
- Verify that the roots of \(P_n(x)\) match the \(n\) points \(a_1\), \(\ldots\), \(a_n\) in Table [tbl:gaussquad] for \(n = 2, 3, 4, 5\).
- With no calculations, explain why \(~\int_{-1}^1 P_2(x)\,P_3(x)\,\dx = \int_{-1}^1 P_3(x)\,P_4(x)\,\dx = \int_{-1}^1 P_4(x)\,P_5(x)\,\dx = 0\).
- For \(n=0, 1, 2, 3\), verify that
\[\int_{-1}^1 P_n^2(x)~\dx ~=~ \frac{2}{2n+1} ~. \nonumber \]
Repeat Example
Example \(\PageIndex{1}\): gqinf
Add text here.
Solution
with \(n=4\) points.
Repeat Example
Example \(\PageIndex{1}\): gqinf
Add text here.
Solution
with \(n=5\) points.
Use Table [tbl:gaussquadimp] to approximate the value of \(~\displaystyle\int_0^{\infty} \ln\,(1+x)\,e^{-x}\,\dx~\) with \(n=5\) points.
