Example 5.1.1

Here are some points that we could consider to be observations of an unknown function on $[-1,1]$.

The polynomial interpolant, as computed using polyfit, looks very sensible. It's the kind of function you'd take home to meet your parents. The interpolant plot is created using @glsbegin@fplot@glsend@.

But now consider a different set of points generated in almost exactly the same way.

The points themselves are unremarkable. But take a look at what happens to the polynomial interpolant.

Surely there must be functions that are more intuitively representative of those points!

Example 5.1.3

The scipy.interpolate module makes it easy to construct some standard interpolating functions. Here is one that is linear between each pair of consecutive nodes.

We may instead request a smoother interpolant that is piecewise cubic.

Example 5.1.4

The piecewise cubic cardinal function is nowhere greater than one in absolute value. This happens to be true for all the cardinal functions, ensuring a good condition number for the interpolation. But the story for global polynomials is very different.

From the figure we can see that the condition number for polynomial interpolation on these nodes is at least 500.

Example 5.2.1

Let's define a set of 6 nodes (i.e., $n=5$ in our formulas).

We plot the hat functions $H_0,\ldots,H_5$.

Example 5.2.2

We generate a piecewise linear interpolant of $f(x)=e^{\sin 7x}$.

First we sample the function to create the data.

Now we create a callable function that will evaluate the piecewise linear interpolant at any $x$.

We measure the convergence rate for piecewise linear interpolation of $e^{\sin 7x}$.

Example 5.2.3

Since we expect convergence that is $O(h^2)=O(n^{-2})$, we use a log-log graph of error and expect a straight line of slope $-2$.

Example 5.3.1

For illustration, here is a spline interpolant using just a few nodes.

Now we look at the convergence rate as the number of nodes increases.

Since we expect convergence that is $O(h^4)=O(n^{-4})$, we use a log-log graph of error and expect a straight line of slope $-4$.

Example 5.4.2

We try to estimate the derivative of $\cos(x^2)$ at $x=0.5$ using five nodes.

We have to shift the nodes so that the point of estimation for the derivative is at $x=0$.

We can reproduce the weights in the finite difference tables by using equally spaced nodes with $h=1$. For example, here is a one-sided formula.

Example 5.5.2

Let's observe the convergence of the forward-difference formula applied to the function $\sin(e^{x+1})$ at $x=0$.

It's not clear that the sequence is converging. As predicted, however, the errors are cut approximately by a factor of 2 when $h$ is divided by 2.

Asymptotically as $h\to0$, the error is proportional to $h$.

Example 5.5.4

We'll run both formulas in parallel for a sequence of $h$ values.

In each case $h$ is decreased by a factor of 4, so that the error is reduced by a factor of 4 in the first-order method and 16 in the second-order method.

A graphical comparison can be clearer. On a log-log scale, the error should (roughly) be a straight line whose slope is the order of accuracy. However, it's conventional in convergence plots, to show $h$ decreasing from left to right, which negates the slopes.

Example 5.5.5

Let $f(x)=e^{-1.3x}$. We apply finite difference formulas of first, second, and fourth order to estimate $f'(0)=-1.3$.

Again the graph is made so that $h$ decreases from left to right. The errors are dominated at first by truncation error, which decreases most rapidly for the 4th order formula. However, increasing roundoff error eventually equals and then dominates the truncation error as $h$ continues to decrease. As the order of accuracy increases, the crossover point moves to the left (greater efficiency) and down (greater accuracy).

Example 5.6.1

The antiderivative of $e^x$ is, of course, itself. That makes evaluation of $\int_0^1 e^x\,dx$ by the Fundamental Theorem trivial.

The module scipy.integrate has multiple functions that estimate the value of an integral numerically without finding the antiderivative first. As you can see here, it's often just as accurate.

The numerical approach is far more robust. For example, $e^{\sin x}$ has no useful antiderivative. But numerically it's no more difficult.

When you look at the graphs of these functions, what's remarkable is that one of these areas is the most basic calculus while the other is almost impenetrable analytically. From a numerical standpoint, they are practically the same problem.

Example 5.6.2

We approximate the integral of the function $f(x)=e^{\sin 7x}$ over the interval $[0,2]$.

In lieu of the exact value, we will use the quadrature function to find an accurate result.

Here is the error at $n=40$.

In order to check the order of accuracy, we double $n$ a few times and observe how the error decreases.

Each doubling of $n$ cuts the error by a factor of about 4, which is consistent with second-order convergence. Another check: the slope on a log-log graph should be $-2$.

Example 5.6.3

We estimate $\displaystyle\int_0^2 x^2 e^{-2x}\, dx$ using extrapolation.

We start with the trapezoid formula on $n=N$ nodes.

We can now apply weights to get the estimate $T_f(N)$.

Now we double to $n=2N$, but we only need to evaluate $f$ at every other interior node.

As expected for a second-order estimate, the error went down by a factor of about 4. We can repeat the same code to double $n$ again.

Let us now do the first level of extrapolation to get results from Simpson's formula. We combine the elements T[i] and T[i+1] the same way for $i=1$ and $i=2$.

With the two Simpson values $S_f(N)$ and $S_f(2N)$ in hand, we can do one more level of extrapolation to get a 6th-order accurate result.

If we consider the computational time to be dominated by evaluations of $f$, then we have obtained a result with twice as many accurate digits as the best trapezoid result, at virtually no extra cost.

Example 5.7.1

This function gets increasingly oscillatory near the right endpoint.

Accordingly, the trapezoid rule is more accurate on the left half of the interval than on the right half.

Both the picture and the numbers suggest that more nodes should be used on the right half of the interval than on the left half.

Example 5.7.2

We perform the integration and show the nodes selected underneath the curve.

The error turns out to be a bit more than we requested. It's only an estimate, not a guarantee.

Let's see how the number of integrand evaluations and the error vary with the requested tolerance.

As you can see, even though the errors are not less than the estimates, the two columns decrease in tandem. If we consider now the convergence not in $h$ (which is poorly defined) but in the number of nodes actually chosen, we come close to the fourth order accuracy of the underlying Simpson scheme.