This site is being phased out.

Stability Analysis

From Mathematics Is A Science
Jump to navigationJump to search

Perhaps the most important thing that we need to know about these various numerical schemes is WHEN THEY WILL CONVERGE TO THE CORRECT SOLUTION OF THE GIVEN CONTINUOUS Initial Value Problem from which they were obtained, when the values of both the time and space steps are allowed to all tend to zero. Due to known results from Numerical Analysis which will be explicitly cited later, we know that we get CONVERGENCE as long as we have used the correct numerical scheme to discretize the equation, and we also have STABILITY. Of course, stability will turn out to depend on the parameter values delta-t and delta-x RELATIVE TO ONE ANOTHER. We prove that, in-general, for the scheme that we have given for the WAVE-EQUATION ON A SQUARE GRID, the scheme will be STABLE if

$$2 - 4*\frac{c^2*(\Delta t)^2}{(\Delta x)^2} > 0$$

which is equivalent to saying that

$$(\frac{c\Delta t}{\Delta x})^2 \leq \frac{1}{2}$$


We obtain this condition using the methods of VON NEUMANN STABILITY ANALYSIS, whereby we separate variables in the discrete recurrence-relation for the wave-height (note that this is simply a discrete time-dependent differential two-form on our cell-complex of squares). We argue via the discretized case of Sturm-Liouville Theory, which will be discussed in more detail later, that since our recurrence relation is LINEAR in both space variables and in time, all solutions of the recurrence can be written as linear combinations of solutions obtained via Separation of Variables. Thus, in order to obtain a stability condition for our algorithm, we merely need to obtain a stability-condition for the SEPARATED solutions. We THEN assume that our time and space dependent components of the separated solution are each PERIODIC, with a period LARGER THAN THE OVERALL GRID-SIZE. . .this is always valid, clearly, since we have chosen our period large-enough to accommodate any possible set of values for our discrete two-form (wave height) within the grid. By doing so, we use discrete Fourier Analysis, TOGETHER WITH OUR BOUNDARY CONDITIONS, to find expressions for the time and space components of our separated solution; using this, and a clever trick from algebra, we can finally prove our stability condition given above!

After this, we analyze two cases of problems on a NON-SQUARE grid. The first case is of the discretized wave-equation on a grid of EQUILATERAL TRIANGLES, arranged into the overall shape of a Rhombus. Just like in the case of the square grid, we set the value of our differential two-form (wave height) equal to zero at all boundary-triangles. Although, of course, the actual numerical algorithm could be carried out with other boundary-conditions, as well. In this case, I obtain the correct form of the discrete recurrence-relation for the wave-height, and I HYPOTHESIZE a condition on the parameters which will guarantee stability of the scheme, which is obtained intuitively using a simple analogy with the earlier case of the square-grid.

The final analysis is carried-out for the Heat Equation on a REGULAR-HEXAGONAL grid.


Stability Analysis for the Wave-Equation on a Square Grid With Boundary Held at Zero

If we take our recurrence relation for the Wave-Equation on a (N+1)x(N+1) Square Grid with arbitrary time and space steps, but rewrite F as u and our variables within the function as subscripts, with i denoting x, j denoting y, and n denoting our time-step t, we obtain the recurrence-relation

$$u_{i,j}^{(n+1)} = (2 - 4(\frac{c\Delta t}{\Delta x})^2)*u_{i,j}^{(n)} - u_{i,j}^{(n-1)} + (\frac{c\Delta t}{\Delta x})^2*(u_{i-1,j}^{(n)} + u_{i+1,j}^{(n)} + u_{i,j-1}^{(n)} + u_{i,j+1}^{(n)})$$

By writing $$u_{i,j}^{(n)} = a_i b_j c_n$$, which is simply Separation of Variables, and performing some simple algebra, we obtain,

$$\frac{c_{n} + c_{n-2}}{c_{n-1}} = 2 - (\frac{c\Delta t}{\Delta x})^2 * (4 - \frac{a_{i-1} + a_{i+1}}{a_{i}} - \frac{b_{j-1} + b_{j+1}}{b_{j}}) $$

Now, in the last expression above, THE LEFT HAND SIDE ONLY DEPENDS ON n, and THE RIGHT HAND SIDE ONLY DEPENDS ON i AND j. So, BOTH SIDES MUST BE EQUAL TO SOME CONSTANT; let's denote this constant by the Greek letter ZETA from here on. So,

$$c_n = \zeta*c_{n-1} - c_{n-2}$$ and $$\zeta = 2 - (\frac{c\Delta t}{\Delta x})^2 * (4 - \frac{a_{i-1} + a_{i+1}}{a_{i}} - \frac{b_{j-1} + b_{j+1}}{b_{j}})$$

Now, by the theory of LINEAR HOMOGENOUS RECURRENCE RELATIONS WITH CONSTANT COEFFICIENTS ([1]), we can obtain the value of c_n at time-step n as

$$c_n = A*(\frac{\zeta - \sqrt{\zeta^2-4}}{2})^n + B*(\frac{\zeta + \sqrt{\zeta^2-4}}{2})^n$$

where A and B are some constants. Now, note that IF ZETA IS BETWEEN -2 AND 2, then the values being taken to the n-th power in the above expression both have absolute-value of UNITY, and therefore, the expression STAYS BOUNDED. Moreover, IF ZETA IS NOT BETWEEN -2 AND 2, the expression becomes unbounded! Since this is the only time-dependent part of the expression for the separated solution, this means that AS TIME GOES TO INFINITY, THE SOLUTION STAYS BOUNDED AS LONG AS ZETA IS BETWEEN -2 AND 2. Now, since our formulas are all LINEAR, the error terms (difference between solutions with two different sets of initial conditions) grow via the same recurrence as the original terms, and therefore, if the original terms stay bounded, so do the error terms! Now, the definition of stability is that the error terms stay bounded over time; therefore, our separated solution is stable if zeta is between -2 and 2! We use this fact, along with the boundary-conditions and certain assumptions about a_i and b_j, to derive our desired stability condition, as follows.

We can assume that a_i and b_j are each 2N-PERIODIC in i and j, respectively! The fundamental realization here is that it is ALWAYS VALID TO MAKE THIS ASSUMPTION, as we are only working on a (N+1)x(N+1) grid, so no matter what the values are within the grid, they can always be COVERED by a 2N-periodic solution! However, we already know FROM FOURIER ANALYSIS that all 2N-periodic functions can be written as linear combinations of sines and cosines!

BY LINEARITY, it is acceptable without any loss of generality to only look at terms of the form

$$a_i = C*\cos{(\frac{\pi*i*p}{N})} + D*\sin{(\frac{\pi*i*q}{N})}$$ and $$b_j = E*\cos{(\frac{\pi*j*r}{N})} + F*\sin{(\frac{\pi*j*s}{N})}$$ where p, q, r, and s are some integers modulo N and C, D, E, and F are some constants.

Now, our BOUNDARY CONDITIONS say that

$$a_0 = a_N = b_0 = b_N = 0$$

so, without loss of generality, we may say that

$$C=E=0 \ and \ D=F=1$$

(Note: In many applications of this method, the numbers (q,s) would be called the WAVE-NUMBER for the particular solution being analyzed).

Now, WE ALREADY HAVE A FORMULA FOR THE VALUE OF ZETA IN TERMS OF a_i, b_j, etc.; so, for terms of the type we are now analyzing, zeta will be dependent on both the integer q and the integer s, and for stability, we require that ALL VALUES OF ZETA(q,s) BE BETWEEN -2 AND 2!

Our formula for the value of Zeta in terms of a_i and b_j given above may now be rewritten, using trigonometric identities, as

$$\zeta_{q,s} = 2 - 2 (\frac{c\Delta t}{\Delta x})^2(2 - \cos{(\frac{iq\pi}{N})} - \cos{(\frac{js\pi}{N})})$$

Of course, Zeta clearly also depends upon the values of i and j, despite being written only with subscripts q and s. Now we require, for stability, that Zeta be between -2 and 2 for all values of (i,j,q,s).

Now, we already know that, FOR REAL VALUES OF x, 1-cos(x) is between 0 and 2, and this inequality is precise; therefore, our above expression becomes

$$2 - 8(\frac{c\Delta t}{\Delta x})^2 \leq \zeta_{q,s} \leq 2$$

For stability, we merely require this inequality to not violate our requirement that Zeta always be BETWEEN -2 AND 2! Therefore, we have that:

$$The \ scheme \ is \ stable \ if \ \ (\frac{c\Delta t}{\Delta x})^2 \leq \frac{1}{2} $$

Proposed Stability Analysis for the Wave-Equation on a Rhomboidal Grid of EQUILATERAL TRIANGLES

Say we have a 2NxN Rhomboidal grid of equilateral triangles, with a wave-height u defined for the interior of each cell (thus, u is a discrete 2-form on our cell complex of triangles), illustrated here for the case N=3 to show what we mean:

Triangular Grid.jpg


Meanwhile, the LINEARLY DAMPENED version of the equation with damping-constant k will take the form:

$$u_{i,j}^{(n+1)} = (2 -3 \frac{c^2\Delta t ^2}{\Delta x ^2} - k\Delta t)*u_{i,j}^{(n)} - (1 - k\Delta t)u_{i,j}^{(n-1)} + \frac{c^2\Delta t ^2}{\Delta x ^2}*(u_{i,j+1}^{(n)} + u_{i,j-1}^{(n)} + u_{i+(-1)^{(j-1)},j}^{(n)})$$

The term

$$(u_{i,j+1}^{(n)} + u_{i,j-1}^{(n)} + u_{i+(-1)^{(j-1)},j}^{(n)})$$

is due to taking the spatial Laplacian at (i,j), which involves the average of the values of u at the three cells adjacent to (i,j). Due to the alternating orientations of the cells, we end up with our subscript of the form (-1)^(j-1) in this formula!


I HYPOTHESIZE THAT THIS SCHEME IS STABLE WHEN

$$(2 -3 \frac{c^2\Delta t ^2}{\Delta x ^2} - k\Delta t) \geq 0$$

Proposed Stability Analysis for the Heat-Equation on an approximately-Rhomboidal grid of REGULAR HEXAGONS

Say we have an NxN grid of regular hexagons arranged into the shape of a rhombus (illustrated below for the case N=3 to show what we mean here) where in each 2-cell, we have a temperature T which is thus a discrete 2-form.

Hexagonal Grid Annotated.jpg


In this case, if we write k for the constant of heat-conductivity through the walls, then by Newton's Law of Cooling, we obtain for the value T in any INTERIOR cell:

$$T_{i,j}^{(n+1)} = (1 - \frac{k\Delta t}{(\Delta x)^2})*T_{i,j}^{(n)} + (\frac{k\Delta t}{6(\Delta x)^2})*(T_{i-1,j}^{(n)} + T_{i-1,j+1}^{(n)} + T_{i,j-1}^{(n)} + T_{i,j+1}^{(n)} + T_{i+1,j}^{(n)} + T_{i+1,j-1}^{(n)})$$

I HYPOTHESIZE THAT THIS SCHEME IS STABLE WHEN

$$ 0 \leq \frac{k\Delta t}{(\Delta x)^2} \leq \frac{6}{5}$$

I believe this stability condition is sufficient, but not necessary, as I suspect that the given constant could range all the way up to 2 and still give stability; however, 6/5 is the best constant that I was able to obtain for a proof so far.

Sketch Proof of this Theorem:

On Certain Matters of Rigor: The Value of the Space-Steps, and the Relationship Between Stability and Convergence, as well as some background on Sturm-Liouville Theory

Clearly, the value of the space-steps "Delta x" must be defined in each case which has been considered; for the square-grid, the space steps are clearly equal in value to the length of the side of any given grid-square. For the Hexagonal-Grid, meanwhile, the space-steps should be the same size AS THE DISTANCE BETWEEN THE CENTERS OF TWO ADJACENT HEXAGONS! And, for the triangular grid, it would seem that the space-steps should be equal to the side-length of any given grid-triangle.

It must also be made clear why, in the first place, we wish to consider Stability of our numerical schemes rather than the question of convergence to the continuous solutions! The reason for this is the Lax Equivalence Theorem Lax Equivalence Theorem which states that such a Numerical Scheme converges to the continuous IVP from which it was obtained if and only if it is stable!

The fundamental concepts from Sturm-Liouville theory which must be utilized in the above analysis are that, for any equation defined in terms of a linear-operator applied to a function being equated to a zero on the other side, all possible solutions of the equation can be written as linear combinations of solutions with separated-variables. It is immediately apparent where, in the above analyses, this concept has been applied.