Academic Articles & Supplements

Mathematica in Action: Check Your Answers... But How?

Rob Knapp, Stan Wagon

Author

Rob Knapp, Stan Wagon

Title

Mathematica in Action: Check Your Answers... But How?

Description

Numerical instabilities can take a variety of forms. How can one recognize cases where such instabilities cause the numerical solution to a differential equation to be invalid?

Check Your Answers . . .
But How?

Numerical instabilities can take a variety of forms. How can one recognize cases where such instabilities cause the numerical solution to a differential equation to be invalid?

by Rob Knapp and Stan Wagon

One of the joys of using Mathematica is the ability it gives us to check the results of complicated computations in several different ways. For symbolic work such as derivatives, antiderivatives, or infinite sums, this is pretty easy: just make some numerical approximation to the problem and see if the results agree. But the situation regarding work that is numerical in the first place is very different.

NDSolve, the built-in command for finding approximate solutions to differential equations, has AccuracyGoal and PrecisionGoal options that allow you to control the error. The difficulty arises because the all of the methods NDSolve is based on are iterative, taking many steps to arive at a solution. The error control goals only apply to what is referred to as the local truncation error at each individual step. For most differential equations, this is not a problem since the overall, or global, error accumulated over the steps is typically proportional to the number of steps and the local truncation error. However, there are differential equations with sensitivities that cause the global error to be much greater than this.

When S.W. taught differential equations a couple of years ago, he tried hard to discover appropriate techniques to tell whether the solution to a complicated system, such as that governing the motion of a double pendulum or a forced pendulum, was reliable. Techniques for verifying a numerical solution are, as far as we can tell, not mentioned in textbooks. Of course, the techniques presented here are not foolproof. But they are very simple to use and we should surely be teaching them to our students. If readers know of examples that shed new light on these sort of difficulties, we would be pleased to learn about them.

Method
1
: Don’t Forget Symbolics

There are many reasons why numerical algorithms are better than symbolic ones. But don’t write off symbolic methods entirely, since they can provide an alternative approach. We begin with an example due to Courtney Coleman (Harvey Mudd College).

soln=x[t]/.First[NDSolve[ {

′

[t]==2x[t]+Cos[t],x[0]==-2/5},x[t],{t,0,7.5}]]

InterpolatingFunction[{{0.,7.5}},<>][t]

Plot[soln, {t,0,7.5}];

Figure

: An approximate solution to

′

(t)=2x(t)+cost

The plotted solution looks plausible, but in fact it is highly inaccurate. The sections that follow will present techniques for discovering that such a solution is wrong. However, it is worth noting that this example is easily handled by DSolve. That gives us the symbolic solution

(sint-2cost)/5

, and that can be checked by differentiation, so we know it is correct.

trueSoln=x[t]/.First[DSolve[ {

′

[t]==2x[t]+Cos[t],x[0]==-2/5},x[t],t]]

(-2Cos[t]+Sin[t])

Plot[trueSoln,{t,0,7.5}];

Figure

: The symbolic solution to the equation

′

(t)=2x(t)+cost

is very different from the machine-precision numerical solution in Figure

We can easily plot a family of solutions corresponding to initial conditions that differ from

by multiples of

-5

. As can be seen in Figure

, the difficulty with this equation comes from the tremendous divergence that occurs for solutions with initial values near

soln=x[t]/.FirstDSolve

′

[t]==2x[t]+Cos[t],x[0]==-25+i

-5

,x[t],t;Plot[Evaluate[Table[soln,{i,-10,10}]],{t,0,15},PlotRange{-2,10},FrameTrue,AxesFalse];

Figure

: The family of solutions corresponding to initial conditions near

; they are all unbounded except for the one corresponding to the exact starting value

Method
2
: Look Critically at the Results

Here is an innocent-looking but ultimately nasty equation found by John Polking (Rice University):

″

(t)=

′

(t)

-x(t)

. It will be more convenient to look at this as a system:

′

-x

. Small initial values cause no problem (try them), but when

x(0)

gets to about 9 something interesting happens. Here we solve the equation and plot the pair

(x,

′

)

in the phase plane.

=9;tmax=13;Polking={

′

[t]==y[t],

′

[t]==

y[t]

-x[t]};soln={x[t],y[t]}/.First[NDSolve[{Polking,x[0]==

,y[0]==0},{x[t],y[t]},{t,0,tmax}]];

trajectory=ParametricPlot[Evaluate[soln],{t,0,tmax},FrameTrue,AxesOrigin{0,0},AxesTrue,Epilog{AbsolutePointSize[4],Point[{

,0}]},FrameTicks{Range[0,9,3],Range[-2,2],None,None}];

Figure

: The solution to Polking’s equation that starts at

(9,0)

looks plausible, but is in fact very wrong.

We can use this example to illustrate two ways of checking a numerical result. First, if we just look at more of the trajectory we will discover something interesting. The code that follows colors the previously computed trajectory gray.

soln1={x[t],y[t]}/.First[NDSolve[{Polking,x[0]==

,y[0]==0},{x[t],y[t]},{t,-tmax,tmax}]];ShowBlock{$DisplayFunction=Identity}, {ParametricPlot[Evaluate[soln1],{t,0,tmax},FrameTrue,AxesOrigin{0,0},AxesTrue,PlotStyle{AbsoluteThickness[4],GrayLevel[0.6]},Epilog{AbsolutePointSize[4],Point[{

,0}]},FrameTicks{Range[0,9,3],Range[-2,2],None,None}],ParametricPlot[Evaluate[soln1],{t,-tmax,0}]};

Figure

: By coloring part of the orbit gray we see that this orbit crosses itself, which cannot happen for an autonomous system. Therefore there must be an error in the solution.

It is evident from Figure

that the trajectory crosses itself, an impossibility for an autonomous system. We could examine more orbits to get a sense of what is going on (they circulate around the equilibrium point at

(0,0)

It turns out that some algebraic cleverness due to Polking allows us to see the true trajectories. We omit the details (see [Schwalbe and Wagon 1997, §11.2]), but it is not hard to show that the solution with initial value

(

,0)

must lie on the curve defined by

+1)

2(x-

)

-2x-1+2

. Thus we can get an image of the orbits by combining several contour plots as follows. There is one parabolic solution (

) that serves as a separatrix, as it separates closed orbits from unbounded ones.

PolkingOrbit[x0_]:=Block{$DisplayFunction=Identity}, ContourPlot(2x0+1)

2(x-x0)

-2x-1+2

,{x,-3,20},{y,-5,5},ContourStyleAbsoluteThickness[0.6],Contours{0},ContourShadingFalse,PlotPoints120Show[PolkingOrbit/@-3,-2,-1,-0.7,-0.6,-0.53,-0.52,-0.51,-0.5-

-3

,-0.5-

-5

,-0.5-

-7

,-0.5-

-9

,-0.5-

-13

,-0.5,0.3,0.5,1,2,3,4,6,8,10,12,14,16,18,20,Epilog->{AbsolutePointSize[5],Point[{14,0}]}];

Figure

: The true solutions to Polking’s equation can be seen by combining 28 contour plots.

Method
3
: Work Backwards

In some cases, a reasonable approach is to time-reverse the equation and use the final value of a numerical solution as the initial value of the auxiliary system. By "time-reverse" we mean let

u(t)=x(-t)

(similarly for more variables) and solve as t goes from

. However, for such a technique to be useful it is important that the equation be invariant under time-reversal. Consider, for example,

′

(t)=-x(t)

. Nearby solutions converge exponentially towards each other (and towards zero), so nearly all numerical methods do quite well on this. But time-reversal yields

′

(t)=u(t)

, for which nearby solutions diverge from each other exponentially. So, if you started with a final value that had only a small, insignificant roundoff error, working backwards would magnify that error to the visual scale; thus the technique would have indicated a problem where none existed.

But many equations do have symmetries that can be exploited by reversing time, and the Polking equation is one of them. Introduce u(t), v(t) to start at the place where the first solution ends and be the time-reversed versions of x(t) and y(t). So

u(t)=x(-t)

and

v(t)=y(-t)

. Then

′

(t)=-

′

(-t)=-y(-t)=-v(t)

and

′

(t)=-

′

(-t)=-

y(-t)

-x(-t)=u(t)-

v(t)

. Even though this system is not identical to the Polking equation, a simple substitution gets us one that is: let

w(t)

-v(t)

. Then the system becomes

′

(t)=w(t)

and

′

(t)=

w(t)

-u(t)

, which is identical to the original x-y system! We now solve this system starting from the last point of soln. Before plotting we flip the second coordinate so that the plot looks like a time-reversed trajectory.

backSol={u[t],w[t]}/.FirstNDSolve 

′

[t]==w[t],

′

[t]==

w[t]

-u[t], u[-tmax]==soln1/.ttmax, w[-tmax]==-soln2/.ttmax,{u[t],w[t]},{t,-tmax,0};

Show[Block[{$DisplayFunction=Identity},ParametricPlot[Evaluate[backSol{1,-1}],{t,-tmax,0}, PlotStyle{GrayLevel[0.7],AbsoluteThickness[4]}],trajectory ],FrameTrue,PlotRange{{-1,11},{-4,4}}];

Figure

: When we work backward from the last point of the solution of Figure

we see that we do not get back to the initial condition (9, 0).

The reversed solution agrees with the original solution for a while, but then falls away, and does not come anywhere near

(9,0)

. This is not absolute proof that the original solution is bad, but it should raise a red flag.

Method
4
: Perturb the Initial Conditions

Although this is the last method we present, it is perhaps the easiest to use and most general. Numerical solvers proceed step by step, but in an adaptive way, attempting to make each step satisfy a certain error tolerance. Mathematica’s NDSolve tries for 6-digit precision at each step. This aspect of the error is called local truncation error (as opposed to roundoff error or propagation error). Now, since the result of the first step is accurate to only 6 digits, this inaccuracy might propagate unacceptably through the perhaps thousands of subsequent steps. An easy way to check for this is to perturb the initial conditions by 10–6.

We illustrate with the Duffing equation. The equation, which models a certain type of forced, damped pendulum, is given below. Note how perfect rationals are used, because we will want to use high-precision later and so must avoid machine-precision reals.

Duffing=

″

[t]+15/100

′

[t]-x[t]+

x[t]

==3/10Cos[t];lowPrecSoln=NDSolve[{Duffing,x[0]==-1,

′

[0]==1},x,{t,0,100},MaxSteps2000]〚1〛;perturbed=NDSolve[{Duffing,x[0]==-1+

-6

′

[0]==1},x,{t,0,100},MaxSteps2000]〚1〛;x[100]/.{lowPrecSoln,perturbed}

{-0.431506,-1.48062}

The serious disagreement indicates that a 10–6 difference can lead to a radically different solution. Of course — and this is easy to forget — we are not comparing the perturbed solution to a true solution. We are just observing that small changes of the sort that will occur in the numerical algorithm cause major changes in the computed result. Therefore we can place no confidence in either solution. A view of the graphs (Figure

) provides some evidence that the solution is all right up to t = 70.

Plot[Evaluate[x[t]/.{lowPrecSoln,perturbed}],{t,0,100},PlotStyle{{AbsoluteThickness[4],GrayLevel[0.6]},{AbsoluteThickness[1.5]}},FrameTrue];

Figure

: Two Duffing orbits starting from nearby initial conditions. They agree only up to about t = 70. The unperturturbed solution is shown as a thick gray curve and the perturbed case is in black.

Just like the other verification techniques we present, the perturbation method is by no means a certain way of catching problems. For example, this method fails to indicate any problem for the Polking equation example shown in sections

and

because the rapid convergence of nearby initial conditions to the left of the equilibrium point in Figure

eradicates the perturbation.

A perturbation of 10–10 leads to

x(100)≈-0.42

, and so has only a small effect on the computed result. This is because the perturbation is negligible compared to the local truncation error, so it is effectively lost. For these tests to be meaningful, the perturbation must be of the same order of magnitude or larger than the local truncation error.To see more of the correct solutions we can bring Mathematica’s variable precision capabilities to bear. We can change the local error of the numerical algorithm by increasing the accuracy goal to, say, 10. (In Mathematica, “precision” refers to relative error and “accuracy” refers to absolute error; we will use the latter here since these Duffing solutions have modest scales.) When we increase the accuracy goal to 10, then a 10-10 perturbation leads to the results 0.23 and

–0.52

t=100

. Thus, as before, we are still getting inaccurate solutions. But the increased accuracy goal (which causes the numerical algorithm to work harder) has improved things in that the perturbed and unperturbed solutions now agree out to about

t=80

(Figure

lowPrecSoln=NDSolve[{Duffing,x[0]==-1,

′

[0]==1},x,{t,0,110},MaxSteps2000,AccuracyGoal10]〚1〛;perturbed10=NDSolve[{Duffing,x[0]==-1+

-10

′

[0]==1},x,{t,0,110},MaxSteps2000,AccuracyGoal10]〚1〛;x[100]/.{lowPrecSoln,perturbed10};

Plot[Evaluate[x[t]/.{lowPrecSoln,perturbed10}],{t,0,100},PlotStyle{AbsoluteThickness[4],GrayLevel[0.8]}, {AbsoluteThickness[1]},FrameTrue,FrameTicks{Automatic,Automatic,None,None}];

Figure

: A higher accuracy goal causes the perturbed solution to agree with the unperturbed case out to

t=80

. The unperturturbed solution is shown as a thick, light gray curve, while the solution corresponding to a perturbation of 10–10 in a black curve.

With the more accurate solution, it is clear that the divergence of solutions with initial conditions of order 10–6 is not simply an artifact of the numerical method or the choice of size of the local truncation error. On the contrary, the divergence is a property of the differential equation. However, it is true that the smaller the perturbation (or truncation error for that matter), the longer the time interval for which solutions match before visibly diverging.

One must not yield to the temptation to simply increase the accuracy goal to, say, 20, in the hope that this would yield good results past

t=100

. That would be very misleading because the inherent imprecision of hardware arithmetic (typically 16 digits) makes such an accuracy goal unattainable. Increasing the WorkingPrecision option is required in that case.

To see that the perturbation technique does indeed discriminate, consider a different initial condition: (0.6, 1.3). Then, going back to the default accuracy goal, a 10–6 perturbation leads to results that are very close: –0.0147109 and –0.0147116. And a plot of the solutions will show that they are essentially identical and asymptotically periodic.

x[100]/.First[NDSolve[ {Duffing,x[0]==#〚1〛,

′

[0]==#〚2〛}, x,{t,0,100},MaxSteps∞]]&/@{{0.6,1.3},{0.6+

-6

,1.3}}

{-0.0147111,-0.0147126}

Of course, there are other responses to an equation as sensitive as the Duffing equation. It might be that the initial condition is only an approximation, so that it is meaningless to talk about starting at exactly

(-1,1)

. But the idea of starting exactly there makes sense, and Mathematica can help us see what the actual solution looks like. Several differential equations textbooks (see [Knapp and Wagon 1996]) compute trajectories of the Duffing equation for precise initial values with no indication of the numerical instability that exists.

For many sensitive equations there is a theory called shadowing that states that while a trajectory may not correspond to the initial conditions provided, it is the exact solution for some nearby initial condition. However, such theorems are proved only in some cases, and it is still not proved that shadowing applies to the Duffing equation. It can be important to distinguish between an inaccurate approximate solution that may be used to draw some information about the orbits and an accurate solution that can be used in a numerical proof of properties of the system. In fact, due to this very distinction Guckenheimer and Holmes say [Guckenheimer and Holmes 1993, p. 91], “There is a substantial theoretical question as to whether this ‘strange attractor’ of the Duffing equation is an artifact of the noise and is absent in the ideal deterministic system.”

Method
5
: A Conservative Approach

Sometimes you may be able to get additional information about the system without actually having a solution. Many systems have quantities that mathematically must stay constant over the entire course of the solutions. A whole class of systems that have at least one conserved quantity are Hamiltonian systems. Conserved quantities can be very useful in checking the quality of a solution since one can simply check to see how close such a quantity stays to its original value.

For example, this is effectively what we did with the Polking system at the end of section

. As it was pointed out, the quantity

c[x0_][x_,y_]:=(2x0+1)

2(x-x0)

-2x-1+2

should always be zero for a solution starting at the intial condition

x(0)=

y(0)=0

. Let’s look at how c evolves over the solution we computed.

Plot[c[

]@@soln,{t,0,tmax},FrameTrue,AxesNone];

Figure

: The conserved quantity for Polking’s equations should stay close to 0, but near

t=9

serious deterioration sets in.

It looks like this starts to deviate at about t = 9. It is easier to see when we use logarithms.

Plot[Log[10,Abs[c[

]@@soln]],{t,0,tmax}];

Figure

: A logarithmic view of the status of the conserved quantity for Polking’s equation.

Now you can see that even though this grows for a while at the beginning of the solution, it does not grow steadily until after about

t=9

, which, not surprisingly, is near the point at which the curve goes wrong on the phase plane plot.Now you can see that even though this grows for a while at the beginning of the solution, it does not grow steadily until after about

t=9

, indicating the real problem starts there. The error is not clearly visible in the phase plane plot of Figure

until

t=12

, when the error in this conserved quantity is order one and the growth has had a very large effect.For an example of a Hamiltonian system we use the simple pendulum equation

′

=v,

′

=-gsinθ

, for which the conserved quantity is the energy,

(

/2)-gcosθ

. One way of checking the energy is to look at it after the solution is computed (as in the preceding example) or you can watch the energy as the solution progresses. The latter involves some tricks, so we illustrate that. The function h just returns its second argument, which will be v, but h is designed to act on numbers only, thus avoiding some preprocessing that NDSolve would do if h were symbolic. Then it is a simple matter to add a side condition to h so that it prints out the energy. We use a counter so that the printing occurs only at every 1500th step.

g=9.8;energy=

/2-9.8Cos[θ];i=0;h[a_?NumericQ,b_,c_]:=i++; If[Mod[i,1500]==0,Print[StringForm[ "`` steps; t = ``, energy = ``", i,NumberForm[a,3],energy/.{vb,θc}]]]; bNDSolve[{

′

[t]==h[t,v[t],θ[t]],

′

[t]==-gSin[θ[t]],v[0]==0,θ[0]==0.123},{θ[t],v[t]},{t,0,250},MaxSteps∞];

1500 steps; t = 43.7, energy = -9.72599

3000 steps; t = 89.2, energy = -9.72603

4500 steps; t = 135., energy = -9.72606

6000 steps; t = 181., energy = -9.72609

7500 steps; t = 225., energy = -9.72613

This shows that the energy is pretty constant, though because the time interval is so large some deterioration, but only a little, is evident.

An interesting example to show where the failure to properly conserve a quantity indicates some real problems with the solution is a problem in partial differential equations: the nonlinear Schrödinger equation, which is uses to model a variety of things from fiber optics to water waves.

NLSE[u_[x_,t_]]:=ID[u[x,t],t]==D[u[x,t],x,x]+2

Abs[u[x,t]]

u[x,t]

This equation is completely integrable and has an infinite number of conserved quantities (including a Hamiltonian). For this example, we will track the conservation of the analogue to the mass of a wave packet,

∞

∫

-∞

u[x,t]

. We set parameters as follows.

[x_]:=a(1+εCos[

x]);ε=1/10;L=4

π;a=0.35;

:=2πn/L;T=100;

We first approximate the solution with NDSolve.

u=ϕ/.First[NDSolve[{NLSE[ϕ[x,t]],ϕ[x,0]==

[x],ϕ[0,t]==ϕ[L,t]},ϕ,{x,0,L},{t,0,T},MaxSteps20000]]

InterpolatingFunction[{{0,17.7715},{0.,100.}},<>]

Plot3D[Abs[u[x,t]],{x,0,L},{t,0,T},PlotPoints->70];

Figure

: The solution to a nonlinear Schrödinger equation using NDSolve’s defaults.

Since the conserved quantity is itself an integral, we need to approximate the integration to get an approximate value to the quantity. One way is to use NIntegrate.

M[u_][t_]:=NIntegrate[

Abs[u[x,t]]

,{x,0,L}]

but this has the drawback of being quite slow

Timing[M[u][0]]

{0.366667Second,2.18793}

An alternative method is to use a cruder approximation (essentially the Riemann sum computed based on the spatial grid points), which is still good enough to track the value for verification purposes. (The third part of an InterpolatingFunction object contains the points it is defined at, in each dimension, so the first of the third part gives the points in the first (spatial in this case) dimension.)

Clear[M];M[u_][t_]:=Module[{grid=Drop[u〚3,1〛,-1],h}, h=grid2-grid1; Plus@@

Abs[u[grid,t]]

Not surprisingly, this is much faster, and allows us to generate a plot or a table without taking hours.

Table[M[u][t],{t,0,T,T/20}]

{2.1879,2.18784,2.18668,2.18914,2.18903,2.18893,2.18624,2.18927,2.18965,2.19238,2.19233,2.19119,2.19708,2.19688,2.20043,2.2009,2.20712,2.20672,2.21374,2.21361,2.21637}

These values are changing too rapidly, and so indicate that the solution is most likely not correct.

It is worth noting that this computation is done very near an instability in the nonlinear Schrödinger equation, so it is very difficult, if not impossible, to get the solution correctly over a long time interval, even with a specialized method. Still, the result of NDSolve can be improved by specifying a small enough spatial grid step size, and properly balancing the error criteria between spatial and temporal error.

=First[ϕ/.NDSolve[ {NLSE[ϕ[x,t]],ϕ[x,0]==

[x],ϕ[0,t]==ϕ[L,t]},ϕ,{x,0,L},{t,0,T}, StartingStepSize{L/32,Automatic},MaxSteps20000,AccuracyGoal{4,6},PrecisionGoal{4,6}]];

Now the plot of the solution is much more regular.

Plot3D[Abs[

[x,t]],{x,0,L},{t,0,T},PlotPoints70];

Figure

: Controlling the starting step size yields a much better behaved solution to the nonlinear Schrödinger equation.

And the mass analogue is much better conserved, though not perfectly.

Table[M[

][t],{t,0,T,T/20}]

{2.1879,2.18788,2.18787,2.18788,2.18793,2.18794,2.18794,2.18794,2.18795,2.188,2.18803,2.18803,2.18802,2.188,2.18803,2.1881,2.1881,2.18809,2.18808,2.18811,2.18817}

Be cautioned that while a poorly conserved quantity is an indicator of a poor solution, a well conserved quantity does not necessarily mean that you have a good solution. In fact, there is a whole class of methods, known as symplectic methods, that are designed to conserve Hamiltonians to within roundoff error, but their solutions are by no means that accurate.

Method
6
: The Residual

Mathematica’s framework for solving differential equations numerically gives us the results as interpolating functions. This is an incredibly useful feature because such functions can be differentiated, integrated, plotted, and so on. So it is an easy matter to substitute the solution back into the equation to see how well it solves the equation. The residual refers to the result that we expect to be zero, that is, the difference between the left and right sides of the differential equation when a suspected solution is substituted. Here is a simple example.

soln=NDSolve[{

′

[t]==x[t],x[0]==1},x,{t,0,1}]

{{xInterpolatingFunction[{{0.,1.}},<>]}}

residual[t_]:=Evaluate[

′

[t]-x[t]/.First[soln]]

rplot=Plot[residual[t],{t,0,1}];

Figure

: A typical plot of the residual of an approximate solution to a differential equation.

This example is quite typical in that the residual appears to be larger at the beginning (with more frequent oscillations). The reason for this is that NDSolve begins by using small steps with a lower-order method, and tries to work up to larger steps with higher-order methods. Because the order of the method is lower, the interpolation error between the actual points where the derivative was evaluated tends to be greater. Since the interpolation error between the points actually hit by the method is in some sense spurious, you may want to check only the points hit by the method. You can check what these points are by looking at the first part of the third part of the returned InterpolatingFunction object.

grid=(x/.First[soln])〚3,1〛

{0.,0.000894427,0.00178885,0.0223764,0.042964,0.0635516,0.122802,0.182053,0.241303,0.300553,0.399063,0.497572,0.596081,0.69459,0.793099,0.945494,1.}

Here is a plot of the residual at only those points.

rplot1=ListPlot[Transpose[{grid,residual/@grid}],PlotStyleAbsolutePointSize[4]];

Figure

: Examining the residual only at the points NDSolve generates is a way of eliminating interpolation error and focusing on the actual error in the solution algorithm.

And when you show this with the other plot, you see that most of the error visible in the residual is interpolation error.

Show[rplot,rplot1,Ticks{{.5,1},{-1,1}0.00005}];

Figure

: Comparing the algorithm error (dots) and the overall error (curve) shows that most of the error one sees in a residual is interpoaltion error.

It might seem comforting that the residual is small but, because NDSolve tries to keep the residual down locally, one expects a small residual with just about any equation. The reader should compute the residuals for all the examples in this paper; there are differences, but basically all the residuals are small. However, residuals can be useful in a variety of ways. First, they serve as a check that a move to high precision is really causing the algorithm to do the right thing. The residual should decrease in proportion to the amount of precision used. Second, it is worth noting that a small residual does mean that the computed solution is an exact solution to a nearby equation. More precisely, suppose r(t) is the residual for an approximate solution X(t) to the equation

′

(t)=h[x(t),t]

. Then it is trivial to verify that X is an exact solution to the modified equation

′

(t)=h[x(t),t]+r(t)

. If

r(t)

-10

(<10)

, say, and we are in a situation where some of the terms in the equation may not be known to 10 significant digits, then X might be a perfectly valid solution from a physical standpoint. And third, the residual can detect singularities that the numerical method might have missed. Here’s an example where NDSolve misses a variation in the function because it occurs between time steps.

First, define an approximation to a Dirac delta function

ApproximateDiracDelta[t_,Δ_]:=Δ

-π

(tΔ)

Then try to solve the harmonic oscillator equation with an approximate impulse at t = 0.

sol=NDSolve[{

′

[t]+x[t]==ApproximateDiracDelta[t,10],x[-π]==0,

′

[-π]==0},x,{t,-π,π}]

{{xInterpolatingFunction[{{-3.14159,3.14159}},<>]}}

Plot[Evaluate[x[t]/.sol],{t,-π,π},FrameTrue,AxesNone];

Figure

: The solution to a differential equation with a large impulse in a very small neighborhood of

t=0

. The numerical algorithm misses the impulse and returns the constant function 0.

The plot, constantly equal to 0, indicates that the approximate impulse had no effect. Is this correct, or is it a numerical artifact? The residual plot in Figure

shows that we are likely to have a problem here!

Plot[Evaluate[

′

[t]-x[t]-ApproximateDiracDelta[t,10]/.sol],{t,-π,π},PlotRangeAll];

Figure

: The residual plot for the Dirac delta example yields a residual that is much too big, indicating that something is not right.

When you know there is an approximate singularity or a discontinuity in your problem that NDSolve could miss, you can tell NDSolve to be sure to evaluate the function there by including that with the range of values for the independent variable.

sol1=NDSolve[{x''[t]+x[t]==ApproximateDiracDelta[t,10],x[-π]==0,x'[-π]==0},x,{t,-π,π,0}];

Plot[Evaluate[x[t]/.sol1],{t,-π,π},FrameTrue,AxesNone];

Figure

: Forcing NDSolve to look at the equation at t = 0 brings the Dirac delta approximation to life, and the impulse at that t-value is detected.

This technique yields an approximate solution that is identical to the true solution (which you can find symbolically using DSolve).

Method
7
: Old Faithful: Increase the Precision

A traditional approach to numerical computations of any sort is to move from single precision to double precision and compare the results. If the high-precision results differ from the low-precision results, the chances are excellent that the latter is wrong. But Mathematica allows us to increase the precision in steps hopefully seeing convergence to the true solution. Of course, one cannot just change the PrecisionGoal (or AccuracyGoal) without paying attention to the working precision. For NDSolve, a setting of WorkingPrecision to wp, greater than machine precision, causes the precision goal to be set to

wp-10

(unless it is specifcally set otherwise). One should always allow the working precision several more digits than the precision goal to account for rounding errors. We illustrate by revisiting Coleman’s example from section

tmax=10;solnCC[pg_]:=x[t]/.First[NDSolve[{

′

[t]==2x[t]+Cos[t],x[0]==-2/5},x[t],{t,0,3π;tmax},PrecisionGoalpg,WorkingPrecisionpg+10,MaxSteps5000]];

Plot[Evaluate[Table[solnCC[pg],{pg,6,16,10}]],{t,0,tmax},PlotRange{-1,10},FrameTrue,PlotStyle->{{},AbsoluteThickness[2],AbsoluteThickness[3]}];

Figure

: When increasing the precision leads to different results, it is fairly certain that the latter is incorrect. The thicker curves correspond to higher precision in the numerical algorithm.

The fact that the machine-precision solution diverges from the others at

t=5

indicates that something is surely going wrong there. In fact, this example illustrates the common phenomenon of divergence of solutions, which occurs in the Polking equation as well.

One drawback of resorting to higher precision is that computations of solutions to high precision can be quite time-consuming. Not only is the arithmetic more expensive, but the tightened local truncation error criterion results in more steps.

As a further warning, we can’t resist presenting the following example, due to S. M. Rump [Rump 1988] and presented by WRI’s Mark Sofroniou in a recent lecture. Let

f(x,y)

be the innocuous function that follows.

f[x_,y_]:=

1335

(11

-121

-2)+

;

The next computation gives us several results, all of which are wrong. Note that the high-precision results agree with each other; in many situations that is taken as evidence of accuracy, but for this example, the agreement means little.

Table[f[N[77617,wp],N[33096,wp]],{wp,12,34,4}]

{-1.18059×

,-1.18059×

,0.×

,-0.82739606011`-0,-0.82739606011`-0,-0.82739606011`-0}

Getting the correct answer requires nothing more than machine precision, provided it is used at just the right time: the actual f-value is simply a modest rational number.

{ans=f[77617,33096],N[ans]}

-

54767

66192

,-0.827396

Conclusion

Numerically sensitive problems will always pose difficulties, despite increased speed and power of hardware and software. Serious applications demand that any approximate solution should be considered guilty until proven innocent! But Mathematica gives us diverse tools for proving beyond a reasonable doubt that a numerical solution is valid.

REFERENCES

GUCKENHEIMER, J., and P. HOLMES. 1993. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, New York.

KNAPP, R., and S. WAGON, 1996. Orbits worth betting on!, C•ODE•E Newsletter (Consortium for Ordinary Differential Equations Experiments), Winter, 8–13.

RUMP, S. M. 1988. Algorithm for verified inclusions — theory and practice, in Reliability in Computing, R. E. Moore ed., Academic Press, San Diego, pp. 109–126.

SCHWALBE, D., and S. WAGON, 1997. VisualDSolve: Visualizing Differential Equations with Mathematica, Springer/TELOS, New York.

ABOUT THE AUTHORS

Rob Knapp works on the development of numerical technology in Mathematica.

Rob Knapp
Wolfram Research, Inc.
100 Trade Center Drive
Champaign, IL 61820
rknapp@wolfram.com

Stan Wagon uses Mathematica extensively in his research, teaching, and exposition, and is especially appreciative of the new ways Mathematica allows us to look at mathematical objects. He is author or coauthor of several books, including Mathematica in Action, Animating Calculus, VisualDSolve, and Which Way Did the Bicycle Go?. He also teaches a high-altitude course, Rocky Mountain Mathematica, every summer in the mountains of Colorado.

Stan Wagon
Department of Mathematics and Computer Science
Macalester College
St. Paul, MN 55105
wagon@macalester.edu

Cite this as: Rob Knapp, Stan Wagon, "Mathematica in Action: Check Your Answers... But How?" from the Notebook Archive (2002), https://notebookarchive.org/2018-10-10pqlua

Check Your Answers . . . But How?

Method 1: Don’t Forget Symbolics

Method 2: Look Critically at the Results

Method 3: Work Backwards

Method 4: Perturb the Initial Conditions

Method 5: A Conservative Approach

Method 6: The Residual

Method 7: Old Faithful: Increase the Precision

Conclusion

REFERENCES

ABOUT THE AUTHORS

Check Your Answers . . .
But How?

Method
1
: Don’t Forget Symbolics

Method
2
: Look Critically at the Results

Method
3
: Work Backwards

Method
4
: Perturb the Initial Conditions

Method
5
: A Conservative Approach

Method
6
: The Residual

Method
7
: Old Faithful: Increase the Precision