San José State University 

appletmagic.com Thayer Watkins Silicon Valley & TornadoAlley USA 


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 twogridlength 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 Moiretype 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 superexponential growth of the form
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.
Consider the differential equation dy/dt = cy^{2}. The solution to this equation is
This is such a tame solution compared to what comes out of the centered difference approximation. That approximation is
where y_{i} represents y(iΔt). The solution depends upon the first two points y_{0} and y_{1}. 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 10^{100} and beyond.
When the computed values reach such levels the value of the y_{i1} term in the difference equation is insignificant compared to the term involving y_{i}^{2}. In this case the solution to the above difference equation is approximated by the solution to the difference equation
The solution to this difference equation is readily apparent from the first few values:
That is to say, the solution is
Thus the solution depends upon the magnitude of (2Δtcz_{0}). Once the difference equation creates a value Y such that 2ΔtcY > 1 then the explosive sequence of values follows.
The difference equation
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
which has a solution of the form
The values of λ_{1} and λ_{2} are the roots of the quadratic equation
For small values of (ΔtcY) these roots are approximately
Thus the solution for z_{i} would be, for positive c,
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Δtcz_{i}>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 twovariable case illustrates what is involved for the mvariable case. The general two variable case is:
where a special notation has been introduced. For the coefficient c^{h}_{j,k} the superscript h refers to the variable number, y_{h}, to which it applies and the subscripts (j,k) refer to the powers of y_{1} and y_{2}, 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 y_{h}(iΔt) = y_{h,i}
When the computation gets into the explosive phase the (i1) terms can be ignored and the difference equation system becomes
Thus
When these expressions are substituted into the expressions for z_{1,2} and z_{2,2} the result is of the form
where the dcoefficients are determined from the ccoefficients. Generally the coefficient d^{h}_{j,k} is the sum of all the products of ccoefficients that correspond to terms of the form z_{1,0}^{j}z_{2,0}^{k}. In particular, the d^{h}_{4,0} coefficient can only arise from products of the terms involving z_{1,0}^{2} and is thus c^{1}_{2,0}^{3}. This carries forward so that
This means that z_{1,i} will include a term of the form
The other dcoefficient terms for the ith values of z_{1,i} and z_{2,i} will likewise be polynomials of degree (2^{i}1) in the ccoefficents. Therefore z_{h,i} is the sum of terms of the form
As in the single variable case, the variables for the ith step are polynomials of degree 2^{i} in the initial values and degree (2^{i}1) in the ccoefficients and the time step term (2Δt). If some term such as (2Δtc^{h}_{2,0}z_{h,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 twovariable case involves the possibility that the variables will grow on the order of K^{2i}, a superexponential type of growth that might more aptly be called exponentialexponential and which can justifiably be characterized as explosive. Clearly the same conclusion applies for the mvariable case.
References:
Phillips, Norman A. (1959). "An Example of NonLinear Computational Instability." in The Atmosphere and the Sea in Motion, edited by Bert Bolin, pp. 50104. 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. 12364.
HOME PAGE OF Thayer Watkins 