College Math Teaching

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$

May 18, 2010

A Couple of Posts From Cosmic Variance

Sometimes you’ll be happily calculating along, and wind up with an equation you wouldn’t want to meet at night in a dark alley. But, a quick flip through Abramowitz & Stegun frequently turns up your nemesis, along with handy tricks for disarming it. Additional satisfaction comes when you write the paper, and get to throw off lines like “The solutions to equation 4 are confluent hypergeometric functions (of the first kind)”.

However, it will be hard top my amusement when I once discovered that the solutions to my problem were closely related to Anger functions.

Those who have taught differential equations will note the similarity to the Bessel Function of the first kind. In fact, this is a particular solution to: $d^2y/dz^2 + (1/z)dy/dz + (1 - v^2/z^2)y = ((z-v)/({\pi}z^2))sin({\pi}v)$

Note: there is a whole library of interesting functions here.

An interesting issue in elementary probability

An eccentric benefactor holds two envelopes, and explains to you that they each contain money; one has two times as much cash as the other one. You are encouraged to open one, and you find $4,000 inside. Now your benefactor — who is a bit eccentric, remember — offers you a deal: you can either keep the$4,000, or you can trade for the other envelope. Which do you choose?

If you do the usual “expected value” calculation, then your strategy should not matter no matter how much money you see in the envelope and how big of the factor there is between the two envelopes (e. g., one envelope might contain, say, 1,000 times more money than the other).

But now suppose the factor is that one envelope has, say, 1000 times more money than the other and you open your envelope and see 2 dollars. Of course, you’d switch. On the other hand, if your envelope had 10,000 dollars in it, then of course you wouldn’t. The real life complication is that there is a realistic upper bound in how much could be in the envelope and if you discover that your envelope contains something near that upper bound, then you won’t switch.

In short, the “real life” problem is really a conditional probability problem: “what is the probability that my envelope contains the larger amount given that the upper bound is roughly X dollars”.

March 9, 2010

The Importance of Integrals and Standards

One of the challenges of teaching lots of “service” courses is that one sometimes comes under heat from client departments if one flunks too many of their prospective students (especially in the engineering/math/science calculus sequence)

Sometimes, we are told that we are too hard on them or teach students what they don’t need to know.
So, it was “art to my eyes” to read the following post at Cosmic Variance:

Having recently slogged through grading an enormous pile of graduate-level problem sets, I am compelled to share one of the most useful tricks I learned in graduate school.

Make your integrals dimensionless.

This probably seems silly to the theoretical physicists in the audience, who have a habit of changing variables and units to the point where everything is dimensionless and equals one. However, in astrophysics, you frequently are integrating over real physical quantities (numbers of photons, masses of stars, luminosities of galaxies, etc) that still have units attached. While students typically do an admirable job of setting up the necessary integrals, they frequently go off the rails when actually evaluating the integrals, as they valiantly try to propagate all those extra factors.

Here’s an example of what I mean. Suppose you want to calculate some sort of rate constant for photoionization, that when multiplied by the density of atoms, will give you the rate of photo-ionizations per volume. These sorts of rates are always density times velocity times cross section: […]

the integral reduces to something that you can start to wrap your brain around: Basically, they were talking about a change of variables. Of course, the integral is NOT elementary and one would have to use some sort of technique (residue?) to evaluate it.

But the point is that the people in physics EXPECT their students to be able to handle the mathematics.

But about the heat we catch for our flunk-out rates:

to be honest, not everyone is down on us for that:

We make rules we think will help our students–you fail if you don’t do the reading, you fail if your paper isn’t turned in on time, you can rewrite anything you fail, ad infinitum–thinking it will help. Then I come to RYS and see the bodies dropping all over the damned place.

That’s why, this semester, I started to bend instead of break. Kid wants to turn it in late? Okay. Kid can’t be in class. Who cares? I go in every day, try to start a discussion, give an impromptu lecture on days they won’t bite, let them out early once I’ve told them what I guess they probably have to know. If I can make out what the paper is about, I give it at least a B-minus. I mark the hell out of them–I write comments in the margin till there ain’t no margin left, and no ink to write in it with. But the grade is always a B-minus or higher, because if it isn’t, they’ll come to my office requesting a checklist of things they can do to write more effective essays, by which they mean essays that will get better grades. […]

So I inflate grades. So should you, unless you’re teaching your students math or anything related to keeping buildings or airplanes or economic systems from falling apart.

Emphasis mine. What can I say? 🙂