Numerical Integration and the Fundamental Theorem of Calculus

Sean

I occasionally get to think about fun math problems while I’m doing work in applied science. The most recent example of this came up in a conference call. We are concerned about how a satellite measurement will respond to surface pressure changes. Theoretically, the measurement is an integral of the amount of a trace gas in the column, multiplied by an instrument specific kernel:

\int_0^{p^*} q(p)\Delta\sigma(p) dp

Note that the Fundamental Theorem of Calculus tells us that the sensitivity of this quantity with respect to surface pressure p^* should exactly be the integrand evaluated at p^*. One of my collaborators asked me to verify that this was also true for our model approximation, which is a simple numerical calculation of this integral. I thought it would be a piece of cake. In our model coordinate system, the vertical coordinate is specified as a function of the surface pressure so that the model levels follow the terrain. This means that the knots on which we compute the numerical approximation also depend on the surface pressure. More specifically, we have vectors \vec{A} and \vec{B} , so that the vertical coordinate is \vec{p} = \vec{A}+\vec{B}p^*. For simplicity, let f(p) = q(p) \Delta \sigma(p) A simple right hand approximation (since the surface pressure occurs at the right end of the integration domain) would be

\int_0^{p^*} f(p) dp \approx \sum_{i=1}^{N} f(p_i) (p_i-p_{i-1})
= \sum_{i=1}^{N} f(a_i + b_i p^*) \left((a_i-a_{i-1}) -(b_i-b_{i-1})p^*\right)

Now, the right hand side is a function of p^* only, and is easily differentiable using the product rule. After taking a derivative, one would expect that some nice math tricks would lead us to canceling of terms, so that we end up with only q(p^*) \Delta\sigma(p^*), but try as I might, I couldn’t get it to come out. This could easily be due to my own lack of cleverness, but a bit of time searching via Google turned up no results for the problem either. I’m hopeful that a clever mathematician out there can figure out some conditions on \vec{A} and \vec{B} that would force this to be true, but I’ve had to put it aside for now, after trying multiple approximations (left hand, trapezoid, Gaussian quadrature).

Numerical evidence suggests that in fact the statement is not true. I ran the (Enthought Distribution) Python script:

x = linspace(0,1,11)
y = x**2
for i in arange(10):
xn = linspace(0,1.+10.**(-i),11)
yn = xn**2
dquad = (trapz(xn,yn)-trapz(x,y))/(10.0**(-i))
print dquad

In this case, the surface pressure would be 1. An important facet of the problem is keeping the number of knots (11) fixed as we take the limit (as i\rightarrow\infty) . Running the script, you’ll note that the limit appears to be in the neighborhood of 2 (before roundoff kills the computation), even though it ought to be 1. Further experimentation with different degrees of polynomials points to a what may be an interesting math result, which I’ll let you discover for yourself.

I am posting this in the hope that someone will think the problem is interesting enough to have a whack at it. It seems unsatisfying that our tried and true numerical methods would fail such a simple test.

Unknown's avatar

Published by Sean

Leave a comment