College Math Teaching

April 17, 2014

Rant: misconceptions we give to our students

Filed under: editorial, pedagogy — Tags: , , , — collegemathteaching @ 1:53 pm

I gave a take home test to my numerical analysis students. This was one of the problems (a “warm up”):


I wanted them to choose a technique (say, a Lagrange polynomial or a cubic spline) to approximate the function at those values.

Most did fine. However one student’s answer bothered me: the student used a proper method but then rounded…to an integer each time.

Yes, in this special case, by pure chance..this turned out to be correct as this data comes from the gamma function.

But in general, that is terrible intuition. There is no reason a function should take integers to integers, even if that always happens with polynomials with integer coefficients.

Unfortunately, too many elementary calculus textbooks and too many math educators reinforce such bad intuition.

April 15, 2014

Excel’s weird way of dealing with powers

Filed under: pedagogy — Tags: — collegemathteaching @ 6:47 pm

I noticed something a bit weird about how Excel handles powers. For example, the expression e^{\frac{-(x+1)^2}{8}} isn’t ambiguous at all. But in Excel, you get something every different if you enter EXP(-(a2+1)^2/8) or if you enter EXP(-((a2+1)^2/8). Evidently, EXP(-(a2+1)^2/8) is read as EXP((-(a2+1))^2/8) which strikes me as strange.


April 7, 2014

Numerical integration: why the brain is still required…

Filed under: class room experiment, integrals, numerical methods, pedagogy — Tags: — collegemathteaching @ 4:59 pm

I gave the following demonstration in class today: \int^1_0 sin^2(512 \pi x) dx =

Now, of course, even a C student in calculus II would be able to solve this exactly using sin^2(u) = \frac{1}{2} - \frac{1}{2}cos(2u) to obtain: \int^1_0 sin^2(512 \pi x) dx=\frac{1}{2}

But what about the “just bully” numerical methods we’ve learned?

Romberg integration fails miserably, at least at first:


(for those who don’t know about Romberg integration: the first column gives trapezoid rule approximations, the second gives Simpson’s rule approximations and the third gives Boole’s rule; the value of \Delta x gets cut in half as the rows go down).

I said “at first” as if one goes to, say, 20 rows, one can start to get near the correct answer.

Adaptive quadrature: is even a bigger fail:


The problem here is that this routine quits when the refined Simpson’s rule approximation agrees with the less refined approximation (to within a certain tolerance), and here, the approximations are both zero, hence there is perfect agreement, very early in the process.

So, what to do?

One should note, of course, that the integrand is positive except for a finite number of points where it is zero. Hence one knows right away that the results are bogus.

One quick way to get closer: just tweak the limits of integration by a tiny amount and calculate, say, \int^{.999}_{.001} sin(512*\pi *x) dx and do some mathematics!


The point: the integration routines cannot replace thinking.

April 4, 2014

Gaussian Quadrature and Legendre polynomials

In a previous post, I discussed the Legendre polynomials. We now know that if P_n is the n’th Legendre polynomial then:

1. If k < n, \int^{1}_{-1}P_k(x)P_n(x)dx = 0 .
2. If Q(x) has degree k < n then \int^{1}_{-1}Q(x)P_n(x)dx = 0 .
3. P_n(x) has n simple roots in (-1,1) which are symmetric about 0.
4. P_n(x) is an odd function if n is odd and an even function if n is even.
5. P_0, P_1, P_2,....P_n form an orthogonal basis for the vector pace of polynomials of degree n or less.

To make these notes complete, I’d like to review the Lagrange Interpolating polynomial:

Given a set of points (x_0, y_0), (x_1, y_1), .....(x_n, y_n) where the x_0, x_1, ...x_n are all distinct, the Lagrange polynomial through these points (called “nodes”) is defined to be:
L_n(x) = y_0 \frac{(x-x_1)(x-x_2)....(x-x_n)}{(x_0 -x_1)(x_0-x_2)...(x_0-x_n)} +y_1 \frac{(x-x_0)(x-x_2)...(x-x_n)}{(x_1-x_0)(x_1-x_2)....(x_1-x_n)} + ....y_n\frac{(x-x_0)(x-x_2)...(x-x_{n-1})}{(x_n-x_0)(x_n-x_2)....(x_n-x_{n-1})}

It is easy to see that L_n(x_i) = y_i because the associated coefficient for the y_i term is 1 and the other coefficients are zero. We call the coefficient polynomials L_{n,i} = \Pi^n_{j=0, j \ne i}\frac{(x-x_j)}{x_j-x_i}

It can be shown by differentiation that if f is a function that has n+1 continuous derivatives on some open interval containing (a,b) and if L_n(x) is the Lagrange polynomial running through the nodes (x_i, f(x_i)) where all of the x_i \in [a,b] then the maximum of |f(x) - L_n(x)| is bounded by |K f^{n+1}(\omega) | where \omega \in [a,b] . So if the n+1’st derivative of f is zero, then the interpolation is exact.

So, we introduce some notation: let P_m(x) be the m'th Legendre polynomial, let r_1, r_2,...r_m be the roots, and c_i = L_(m,i) = \Pi^n_{j=0, j \ne i}\frac{(x-x_j)}{x_j-x_i} , the j’th coefficient polynomial for the Lagrange polynomial through the nodes (r_i, 1) .

Now let P(x) be any polynomial of degree less than m. First: the Lagrange polynomial L_m through the nodes (r_1, P(r_1)), (r_2, P(r_2)), ...(r_m, P(r_m)) represents P(x) exactly as the error term is a multiple of the m+1'st derivative of P(x) which is zero.

Therefore \int^{1}_{-1} P(x) dx = \int^{1}_{-1} L_m (x)dx = \int^{1}_{-1} P(r_1 )L_{m,1} +P(r_2 )L_{m,2} + ....P(r_m )L_{m,m} dx = \sum^{m}_{i=1} P(r_i )\int^1_{-1} L_{m,i}(x) dx

Now here is the key: \int^1_{-1} L_{m,i}(x) dx only depends on the Legendre polynomial being used and NOT on P(x) .

Hence we have a quadrature formula that is exact for all polynomials of degree less than m .

But, thanks to the division algorithm, we can do even better. We show that this quadrature scheme is exact for all polynomials P whose degree are less than 2m . This is where the orthogonality of the Legendre polynomials will be used.

Let P be a polynomial of degree less than 2m . Then by the division algorithm, P(x) = Q(x)P_m(x) + R(x) where BOTH Q(x) and R(x) have degree less than m.

Now note the following: P(r_i) = Q(r_i)P_m(r_i) + R(r_i) = R(r_i) ; that is where the fact that the r_i are the roots of the Legendre polynomials matters.
Then if we integrate: \int^1_{-1} P(x) dx = \int^{1}_{-1} Q(x)P_m(x) dx + \int^{1}_{-1} R(x)dx . Now because of property 2 above, \int^{1}_{-1} Q(x)P_m(x)dx =0 .

Now because the degree of R is less than m , \int^1_{-1}R(x) dx = \sum^{m}_{i=1} R(r_i )\int^1_{-1} L_{m,i}(x) dx =\sum^{m}_{i=1} P(r_i )\int^1_{-1} L_{m,i}(x) dx

This means: this quadrature scheme using the Legendre polynomials of degree m is EXACT for all polynomials of degree less than 2m.

Of course, one must know the roots r_i and the corresponding \int^1_{-1} L_{m,i}(x) dx but these have been tabulated.

Why this is exciting:

Let’s try this technique on: \int^1_{-1} e^x dx with, say, n = 5. We know that the exact answer is e-e^{-1} . But using the quadrature method with a spreadsheet:

The roots r_i are in column A, rows 5-9. The weights (coefficients; the \int^1_{-1} L_{m,i}(x) dx ) are in column B, rows 5-9. In column F we apply the function e^x to the roots in column A; the products of the weights with e^{r_i} are in column G, and then we do a simple sum to get a very accurate approximation to the integral.

The accuracy: this is due to the quadrature method being EXACT on the first 9 terms of the series for e^{x} = \sum^{\infty}_{j=0} \frac{x^j}{j!} . The inexactness comes from the higher order Taylor terms.

Note: to convert \int^{b}_{a} f(x) dx , one uses the change of variable: u = \frac{2}{b-a} x - \frac{a+b}{b-a} to convert this to \frac{(b-a)}{2}\int^{1}_{-1} f(u \frac{(b-a)}{2}+\frac{a+b}{2})dx

April 1, 2014

Legendre Polynomials: elementary linear algebra proof of orthogonality

In our numerical analysis class, we are coming up on Gaussian Quadrature (a way of finding a numerical estimate for integrals). Here is the idea: given an interval [a,b] and a positive integer n we’d like to select numbers x_i \in [a,b], i \in \{1,2,3,...n\} and weights c_i so that \int^b_a f(x) dx is estimated by \sum^n_{i=1} c_i f(x_i) and that this estimate is exact for polynomials of degree n or less.

You’ve seen this in calculus classes: for example, Simpson’s rule uses x_1 =a, x_2 = \frac{a+b}{2}, x_3 = b and uses c_1 = \frac{b-a}{6}, c_2 =\frac{2(b-a)}{3}, c_3 =\frac{b-a}{6} and is exact for polynomials of degree 3 or less.

So, Gaussian quadrature is a way of finding such a formula that is exact for polynomials of degree less than or equal to a given fixed degree.

I might discuss this process in detail in a later post, but the purpose of this post is to discuss a tool used in developing Gaussian quadrature formulas: the Legendre polynomials.

First of all: what are these things? You can find a couple of good references here and here; note that one can often “normalize” these polynomials by multiplying by various constants.

One way these come up: they are polynomial solutions to the following differential equation: \frac{d}{dx}((1-x^2)\frac{d}{dx} P_n(x)) + n(n+1)P_n(x) = 0 . To see that these solutions are indeed polynomials (for integer values of n ). To see this: try the power series method expanded about x = 0 ; the singular points (regular singular points) occur at x = \pm 1 .

Though the Legendre differential equation is very interesting, it isn’t the reason we are interested in these polynomials. What interests us is that these polynomials have the following properties:

1. If one uses the inner product f \cdot g = \int^1_{-1} f(x) g(x) dx for the vector space of all polynomials (real coefficients) of finite degree, these polynomials are mutually orthogonal; that is, if n \ne m, P_m(x) \cdot P_n (x) = \int^1_{-1} P_n(x)P_m(x) dx = 0 .

2. deg(P_n(x)) = n .

Properties 1 and 2 imply that for all integers n , \{P_0(x), P_1(x), P_2(x), ....P_n(x) \} form an orthogonal basis for the vector subspace of all polynomials of degree n or less. If follows immediately that if Q(x) is any polynomial of degree k < m , then Q(x) \cdot P_m(x) = 0 (Q(x) is a linear combination of P_j(x) where each j < m )

Now these properties can be proved from the very definitions of the Legendre polynomials (see the two references; for example one can note that P_n is an eigenfunction for the Hermitian operator \frac{d}{dx}((1-x^2)\frac{d}{dx} P_n(x)) with associated eigenvalue n(n+1) and such eigenfunctions are orthogonal.

This little result is fairly easy to see: call the Hermitian operator A and let m \ne n, A(P_m) =\lambda_m P_m, A(P_n) =\lambda_n = A(P_n) and \lambda_n \ne \lambda_m .

Then consider: (A(P_m) \cdot P_n) = (\lambda_m P_m \cdot P_n) = \lambda_m (P_m \cdot P_n ) . But because A is Hermitian, (A(P_m) \cdot P_n) = (P_m \cdot A(P_n)) = (P_m \cdot \lambda_n P_n) = \lambda_n (P_m \cdot P_n) . Therefore, \lambda_m (P_m \cdot P_n ) = \lambda_n(P_m \cdot P_n) which means that P_m \cdot P_n = 0 .

Of course, one still has to show that this operator is Hermitian and this is what the second reference does (in effect).

The proof that the operator is Hermitian isn’t hard: assume that f, g both meet an appropriate condition (say, twice differentiable on some interval containing [-1,1] ).
Then use integration by parts with dv =\frac{d}{dx} ((1-x^2) \frac{d}{dx}f(x)), u =g(x) : \int^1_{-1} \frac{d}{dx} ((1-x^2) \frac{d}{dx}f(x))g(x) = ((1-x^2) \frac{d}{dx}f(x))g(x)|^1_{-1}-\int^1_{-1}(1-x^2)\frac{d}{dx} f(x) \frac{d}{dx}g(x) dx . But ((1-x^2) \frac{d}{dx}f(x))g(x)|^1_{-1} =0 and the result follows by symmetry.

But not every student in my class has had the appropriate applied mathematics background (say, a course in partial differential equations).

So, we will take a more basic, elementary linear algebra approach to these. For our purposed, we’d like to normalize these polynomials to be monic (have leading coefficient 1).

Our approach

Use the Gram–Schmidt process from linear algebra on the basis: 1, x, x^2, x^3, x^4.....

Start with P_0 = 1 and let U_0 = \frac{1}{\sqrt{2}} ; here the U_i are the polynomials normalized to unit length (that is, \int^{1}_{-1} (U_k(x))^2 dx = 1 . That is, U_i(x) = \sqrt{\frac{1}{\int^1_{-1}(P_i(x))^2 dx}} P_i(x)

Next let P_1(x) =x, U_1(x) = \sqrt{\frac{2}{3}} x

Let P_2(x) = x^2 - \sqrt{\frac{2}{3}} x \int^{1}_{-1} (\sqrt{\frac{2}{3}} x)x^2 -\frac{1}{\sqrt{2}}\int^{1}_{-1} \frac{1}{\sqrt{2}}x^2 = x^2 -\frac{1}{3} Note that this is not too bad since many of the integrals are just integrals of an odd function over [-1,1] which become zero.

So the general definition:

P_{n+1}(x) = x^{n+1} - U_n \int^1_{-1}x^{n+1} U_n(x) dx - U_{n-1}\int^1_{-1} U_{n-1} x^{n+1}dx .... - \frac{1}{\sqrt{2}}\int^1_{-1} \frac{1}{\sqrt{2}}x^{n+1} dx

What about the roots?
Here we can establish that each P_m(x) has m distinct, real roots in (-1,1) . Suppose P_m(x) has only k < m distinct roots of odd multiplicity in (-1,1) , say x_1, x_2, ...x_k . Let W(x) = (x-x_1)(x-x_2)...(x-x_k) ; note that W has degree k < m . Note that P_m(x)W(x) now has all roots of even multiplicity; hence the polynomial P_m(x)W(x) cannot change sign on [-1,1] as all roots have even multiplicity. But \int^{1}_{-1} P_m(x)W(x) dx = 0 because W has degree strictly less than m . That is impossible. So P_m(x) has at least m distinct roots of odd multiplicity, but since P_m(x) has degree m, they are all simple roots.

Create a free website or blog at