Thayer Watkins
Silicon Valley
& Tornado Alley

The Nature of the Instability Problem of
Centered Difference Approximations
of Derivatives in the Numerical
Solution of Differential Equations

The centered difference approximation to the first derivative has the advantage of greater accuracy over the forward and the backward difference approximations but it necessarily introduces an additional solution, called the computational mode, which does not correspond to a solution of the differential equation. If this computational mode grows in magnitude as the numerical solution is computed it may make the numerical solution worthless as an approximation to the solution of the differential equation. The analysis below points out that the growth of the computational mode requires an initial deviation from the true solution and this makes the computation at the intial steps of the solution crucial. Since the centered difference introduces a higher order for the difference equation than that of the differential equation its solution requires more initial conditions and those additional initial conditions have to be computed using some approximation other than the centered difference. The errors introduced by those initial condition approximation may grow and propagate to the detriment of the accuracy of the later computations.

Consider the simple differential equation

dy/dt = cy.

This is a linear first order differential equation and has the solution

y(t) = y(0)ect

The centered difference approximation is

(yi+1 - yi-1)/2h = cyi
which is equivalent to
yi+1 = yi-1 + (2hcy)yi

where h is the time step Δt and yi=y(ih). This is a second order difference equation. Since it is linear and homogenous it has two solutions of the form yi = λi. The parameter λ must satisfy the condition

λ2 -2hcλ -1 = 0

The two solutions to this quadratic equation are

λ = hc + ((hc)2 + 1)1/2
λ = hc - ((hc)2 + 1)1/2

When h is sufficiently small such that (1+(hc)2)1/2=1+(hc)2/2 these solutions are closely approximated by

λ = hc + 1 + (hc)2/2 = 1 + hc(1+hc/2)
λ = hc - (1 + (hc)2/2) = -1 + hc(1-hc/2)


λi = ekih
k = ln(λ)/h

When z is sufficiently small ln(1+z) is closely approximated by z. Thus one solution corresponds to

k = ln(1 + hc(1+hc/2))/h = c(1+hc/2)

In the limit as h goes to zero this goes to c and corresponds to the true solution to the differential equation. The other solution λ = -1 + hc(1-hc/2) is a cyclic solution of period 2. When c is positive the magnitude of this solution is less than unity and the difference scheme would be stable. If c is negative the magnitude is greater than unity and the solution would blow up; i.e., the deviation from the true solution would growl exponentially.

Now let us consider the process for getting the solution to the difference scheme. Although the numerical values are computed iteratively the process is equivalent to fitting the general solution of the difference equation to the initial conditions. The general solution is

yi = q1λ1i + q2λ2i

where q1 and q2 are constants to be determined from the initial conditions; i.e.,

q1λ10 + q2λ20 = y0
q1λ11 + q2λ21 = y1

Let λ1 denote the physical mode and λ2 the computational mode. The computational mode is no problem if q2=0.

While y0 is known and equal to the value of the true solution at time equal to zero, y(0), the value of y1 is not known because the true solution at t=h, y(h), is generally not known. Thus y1 usually estimated by some method outside of the iteration scheme.

The values of q1 and q2 are found as the solution to the equations:

q1 + q2 = y0
q1λ1 + q2λ2 = y1

The solutions are

q1 = (λ2y0 - y1)/(λ2 - λ1)
q2 = (y1 - λ1y0)/(λ2 - λ1)

From this we see that q2 will be zero if

y1 = λ1y0

Thus the ideal value for y1 is not the true solution value at t=h, y0ehc, but instead y0λ1. The deviation from the physical mode is directly proportional to the deviation of the value used for y1 from y0λ1. To the extent that λ1 is a good approximation of ehc the true value of the solution at t=h is a good approximation of the ideal value of y0λ1.

The points above can be nicely illustrated. First consider the case of c=1. This is the stable case. The numeric solutions based upon centered difference but using two slightly different estimates of y1 for h = 0.1 and for 10 time steps are so close it is difficult to distinguish one from the other and both lie on the line of the true solution. This would also hold true if the number of time steps were increased to 100.

For the unstable case of c=-1 the match up of the centered difference solution for the two estimates of y1 (one computed using a forward difference and the other computed using a backward difference) is quite close to each other and the true solution for the first 10 time steps, as shown below.

But when the number of time steps is increased to 50 the picture is quite different. The numberic solutions begin to oscillate above and below the true solution to a drastic extent.

But the numeric solution based upon y1 being estimated by by y0λ1 is right on the true solution throughout the 50 time steps.

When the number of time steps is increased to 100 the numeric solutions based upon centered differences and y1 from forward and backward differences blow up and run off the scale of the graph. But the numeric solution based upon y1=y0λ1 is right on target through the 100 steps, as shown below.

A note of caution: If the value of y1 is not exactly equal to the exact value of λ1 times the exact value of y0 the solution will eventually blow up no matter how impressively close it is during the initial time steps.

HOME PAGE OF applet-magic
HOME PAGE OF Thayer Watkins