San José State University
Thayer Watkins
Silicon Valley
& TornadoAlley

Nonlinear Computational Instability

In 1959 Norman A. Phillips introduced the computational meteorology world to the concept of nonlinear computational instability. Phillips had created an experimental general circulation model, one of the first, and it appeared to work, but when the computation was extended beyond a certain time the computed values simply exploded to astronomical values. Computational instability was not unususal in meteorological models but what Phillips found was different from previously observed instabilities. For one thing, the instability was unaffected by a reduction in the time step. For another, the magnitude of the solution values did not simply grow exponentially as in a computational mode solution but instead remained relatively small and nearly constant before exploding into astronomical values. The solution values even decreased a bit before the computational explosion, as shown below in the graph of Phillips' results.

Phillips explained nonlinear computational instability in terms of aliasing. He felt that the effects for wavelengths shorter than twice the gridlength were accumulating in in the two-gridlength cycles and creating enormous values. The meteorological profession seems to have generally accepted Phillips' explanation. Since the nonlinear computational instability could not be eliminated by shortening the time step Phillips looked to averaging or filtering to moderate the instability. These remedies were not really the solution to the problem. The solution to the problem of nonlinear computational instability came from Akio Arakawa of UCLA. Arakawa's solution was staggered grids. The explanation of how one of Arakawa's staggered grids could eliminate the problem of nonlinear computational instability that is presented in the literature is in terms of aliasing and the prevention of Moire-type phenomena.

Although Phillips' explanation of the source of nonlinear computational instability may have some validity as an explanation of the phenomenon there is a more mundane explanation. The nonlinearity of the difference equations creates super-exponential growth of the form

xt = Knt

where K is a constant and n is the degree of the nonlinearity. The explosion in values arises because there is an exponential raised to an exponential power.

The Single Variable Case

Consider the differential equation dy/dt = cy2. The solution to this equation is

1/y(t) - 1/y(0) = ct
which upon rearrangement becomes
y(t) = y(0)/(1+cy(0)t).

This is such a tame solution compared to what comes out of the centered difference approximation. That approximation is

yi+1 = yi-1 + (2Δtc)yi2

where yi represents y(iΔt). The solution depends upon the first two points y0 and y1. It needs to be noted here that forward and backward difference approximations of the derivative do not have the extreme problems of instability that prevail for the centered differences.

For some initial values the solution is bounded but for other initial values the computed results very quickly reach magnitudes of 10100 and beyond.

When the computed values reach such levels the value of the yi-1 term in the difference equation is insignificant compared to the term involving yi2. In this case the solution to the above difference equation is approximated by the solution to the difference equation

zi+1 = (2Δtc)zi2

The solution to this difference equation is readily apparent from the first few values:

z1 = (2Δtc)z02
z2 = (2Δtc)3z04
z3 = (2Δtc)7z08
z4 = (2Δtc)15z016

That is to say, the solution is

zi = (2Δtc)2i-1 z02i
or, more conveniently
zi = (2Δtcz0)2i/(2Δtc)

Thus the solution depends upon the magnitude of (2Δtcz0). Once the difference equation creates a value Y such that |2ΔtcY| > 1 then the explosive sequence of values follows.

The difference equation

yi+1 = yi-1 + (2Δtc)yi2

may have a range of i such that the values of y are nearly constant, say Y. In this case the above difference equation's solution would be approximated by the solution to

zi+1 = zi-1 + (2ΔtcY)yi

which has a solution of the form

zi = d1λ1i + d2λ2i

The values of λ1 and λ2 are the roots of the quadratic equation

λ2 - (2ΔtcY)λ - 1 = 0

For small values of (ΔtcY) these roots are approximately

1 + (ΔtcY) and -1+(ΔtcY)

Thus the solution for zi would be, for positive c,

zi = d1(1 + ΔtcY)i + d2(-1)i(1 - ΔtcY)i

The negative root above explains the source of the 2Δt cycle.

Although for small (ΔtcY) the growth rate of the solution may be slow it will eventually reach a point where |2Δtczi|>1 and explosive growth will follow.

If the nonlinear differential equation involves a polynomial of degree n then the base of the power function will be n rather than 2.

The Multiple Variable Case

The two-variable case illustrates what is involved for the m-variable case. The general two variable case is:

dy1/dt = c12,0y12 + c11,1y1y2 + c10,2y22
dy2/dt = c22,0y12 + c21,1y1y2 + c20,2y22

where a special notation has been introduced. For the coefficient chj,k the superscript h refers to the variable number, yh, to which it applies and the subscripts (j,k) refer to the powers of y1 and y2, respectively, to which it applies. Since j+k=2, one of the subscripts could be dispensed with here but later it is useful to have both so that the sum can supply some information.

The centered finite difference approximation to the above differential equations is, in terms of yh(iΔt) = yh,i

y1,i+1 = y1,i-1 + (2Δt)[c12,0y1,i2 + c11,1y1,iy2,i + c10,2y2,i2]
y2,i+1 = y2,i-1 + (2Δt)[c22,0y1,i2 + c21,1y1,iy2,i + c20,2y2,i2]

When the computation gets into the explosive phase the (i-1) terms can be ignored and the difference equation system becomes

z1,i+1 = (2Δt)[c12,0z1,i2 + c11,1z1,iz2,i + c10,2z2,i2]
z2,i+1 = (2Δt)[c22,0z1,i2 + c21,1z1,iz2,i + c20,2z2,i2]


z1,1 = (2Δt)[c12,0z1,02 + c11,1z1,0z2,0 + c10,2z2,02]
z2,1 = (2Δt)[c22,0z1,02 + c21,1z1,0z2,0 + c20,2z2,02]

When these expressions are substituted into the expressions for z1,2 and z2,2 the result is of the form

z1,2 = (2Δt)3[d14,0z1,04 + d13,1z1,03z2,0 + d12,2z1,02z2,02 + d11,3z1,0z2,03 + d10,4z2,04]
z2,2 = (2Δt)3[d24,0z1,04 + d23,1z1,03z2,0 + d22,2z1,02z2,02 + d21,3z1,0z2,03 + d20,4z2,04]

where the d-coefficients are determined from the c-coefficients. Generally the coefficient dhj,k is the sum of all the products of c-coefficients that correspond to terms of the form z1,0jz2,0k. In particular, the dh4,0 coefficient can only arise from products of the terms involving z1,02 and is thus c12,03. This carries forward so that

d12i,0 = c12,02i-1

This means that z1,i will include a term of the form


The other d-coefficient terms for the i-th values of z1,i and z2,i will likewise be polynomials of degree (2i-1) in the c-coefficents. Therefore zh,i is the sum of terms of the form


As in the single variable case, the variables for the i-th step are polynomials of degree 2i in the initial values and degree (2i-1) in the c-coefficients and the time step term (2Δt). If some term such as (2Δtch2,0zh,0) achieves a magnitude greater than unity then the values increase explosively. It may take a period of time for such a term to arise and during that time the computations may appear to be stable. However once such terms occur the instability is immediately apparent. Thus the two-variable case involves the possibility that the variables will grow on the order of K2i, a super-exponential type of growth that might more aptly be called exponential-exponential and which can justifiably be characterized as explosive. Clearly the same conclusion applies for the m-variable case.


Phillips, Norman A. (1959). "An Example of Non-Linear Computational Instability." in The Atmosphere and the Sea in Motion, edited by Bert Bolin, pp. 501-04. New York: Rockefeller Institute Press.

Phillips, Norman A. (1956). "The General Circulation of the Atmosphere: A Numerical Experiment." Quarterly Journal of the Royal Meteorological Society Vol. 82: pp. 123-64.

HOME PAGE OF applet-magic
HOME PAGE OF Thayer Watkins