It is a fact of life that numerical approximations to differential equations may exhibit unstable behavior. That is, computational results may include exponentially growing and sometimes oscillating features that bear no relation to the solution of the original differential equation. This type of behavior is referred to as a computational instability.

미분 방정식의 수치 근사에서 불안정한 거동을 보이는 것은 드문 일이 아니다.  즉, 계산 결과에서 기하 급수적으로 증가하고 때로는 부호가 자주 반전하는 등 원래의 미분 방정식의 솔루션과 아무런 관련이없는 특성을 볼 수 있습니다.  이러한 동작을 “계산 불안정성”이라고합니다.

It is important to distinguish computational instability from physical instabilities, which may occur in some physical problems. In practice, the distinction between these types of unstable behavior is not always obvious. Certainly when a computed solution has features that change sign with each step of the finite-difference increments (in time or space) the instability is probably associated with the numerical approximation. On the one hand, changes in sign during time advancement can often be cured by decreasing the size of the time increment. On the other hand, changes in sign with space increments might indicate a lack of computational resolution, which a reduction in the size of the space increments might eliminate.

어떤 물리적인 문제가 발생할 수 있는 물리적 불안정성으로 부터 게산 불안정성을 구별하는 것이 중요합니다.  그러나 실제로는 이러한 불안정한 거동의 차이는 항상 명확하지 않습니다.  계산된 솔루션에서 유한 차등 (시간 또는 공간) 증분 단계마다 부호가 반전하는 특성이 나타나는 경우 이 불안정은 틀림없이 수치 근사에 의한 것으로 생각됩니다.  한편으로는, 시간 진행에 따른 부호 반전은 종종 시간 증가를 작게하여 해결할 수 있습니다.  반면 공간 증가에 따른 부호 반전은 계산 해상도의 부족을 나타낼 수 있습니다.  이것도 공간 증가를 작게하여 문제를 제거 할 수 있는 경우가 있습니다.

In general, it would be desirable to have a means of predicting when computational as opposed to physical instabilities can occur. Unfortunately, no general method has yet been devised. There are, however, some analysis techniques that often provide enough guidance to avoid the majority of frequently encountered computational instabilities. Some of these techniques will be described in this introductory article.

일반적으로 물리적이 아니라, 계산에 의한 불안정성이 생기는 것을 예측 할 수 있는 방법이 바람직 할 것입니다.  불행하게도 그러한 일반적인 방법은 아직 고안되어 있지 않습니다.  그러나 자주 발생하는 계산 불안정성을 많이하지 않도록 하는데 충분한 지침이되는 분석 방법은 여러 가지가 있습니다.  기본 개념을 나타내는 문서에서는 이러한 방법 중 일부에 대해 설명합니다.

A Simple Example of Computational Instability

The following simple example provides a useful introduction to the subject of computational stability. Consider the simple rate equation,

여기에서는 계산 안정성의 기본 개념을 설명하기에 적합한 간단한 예를 보여줍니다.  다음과 같은 간단한 반응 속도 식을 생각합니다.

(1)    

where Q=Q₀ at t=0. The exact solution of this equation is an exponential decay Q=Q0exp(-At).

For a numerical approximation, suppose time is discretized into equal steps of size δt and that t=n δt, where n is an integer. Numerically computed values of Q will the denoted by a lower case q to distinguish them from the exact solution and an exponent n will be used to denote the time step, so that qn is the approximation to Q(n δt). A finite-difference approximation to Eq.1 can be written as,

여기서 Q = Q 0, t = 0으로하면이 식의 정확한 솔루션은 지수 적 감쇠 Q = Q 0 exp (-At)입니다.

수치 근사치를 얻기 위해 일정한 단계 크기 δt 시간을 이산화 t = n δt합니다 (n은 정수).  수치적으로 계산 된 Q 값은 정확한 솔루션과 구별하기 위해 소문자 q로 나타내며 지수 n을 사용하여 시간 단계를 나타냅니다.  그러면 Q (n δt)의 근사치는 q n입니다.  식 1의 유한 차분 근사는 다음과 같이 작성할 수 있습니다.

(2)    

where q₀=Q_0. Because of its simplicity there is an exact solution to Eq.2,

여기서, q 0 = Q 0하면 매우 간결하게 식 2의 정확한 솔루션을 얻을 수 있습니다.

(3)    

Now the interesting thing is to observe how the values of qn evolve with n as the value of Aδt is changed. The values depend on the magnitude of the bracket (1-Aδt), which is raised to the power of n, Fig. 1.

여기서 주목해야 할 것은 Aδt 값이 변화함에 따라 값 q n이 n에 따라 어떻게 변화 해 나갈 것인가 하는 점입니다.  값은 n 승되는 괄호 (1-Aδt) 값의 크기에 따라 달라집니다 (그림 1).

When Aδt is greater than 2.0, the bracket is negative and has a magnitude that is greater than 1.0, so qn changes sign alternately with increasing n and its magnitude is rapidly growing. This is a clear case of computational instability in which the numerical results provide no useful information about the evolution of q. In contrast, if the value of Aδt is less than 1.0, the bracket is positive and has a magnitude less than 1.0. In this case the evolution of qn is a decent approximation of the exponential decay in the exact solution, in which the values lie a little below the exact values. Intermediate values of Aδt (between 1.0 and 2.0) result in qn values that decay, but do so with an oscillation in sign each time step. Based on these observations we can say that the finite-difference approximation Eq. 2 is a stable and reasonable approximation to the exact equation provide the time step δt is limited by the “stability condition” Aδt<1.0.

Aδt가 2.0보다 큰 경우 괄호 안에 음수가되고, 그 절대 값은 1.0보다 커집니다.  따라서 n이 증가 할 때마다 q n의 부호는 반전 값은 크게 증가하고 있습니다.  이것은 계산 결과를 봐도 q가 어떻게 변화 해 나갈 것인가를 파악할 수없는 계산 불안정성이 보이는 전형적인 케이스입니다.  이에 대해 Aδt 값이 1.0보다 작을 때, 괄호 안은 정되고, 그 절대 값은 1.0보다 작아집니다.  이 경우 q n의 변화는 근사치가 정확한 값을 약간 밑도는 엄격 해의 지수 적 감쇠의 우수한 근사치입니다.  Aδt에 중간 값 (1.0과 2.0 사이)을 대입하면 q n 값은 감소하지만, 시간 단계마다 부호가 반전됩니다.  이러한 점에서 유한 차분 근사 식 2는 시간 단계 δt가 “안정 조건”Aδt <1.0에 의해 엄격하게 제한 될 때 식에 안정적이고 합리적인 근사치를 나타내는 것을 알 수 있습니다.

Linear Stability Analysis of von Neumann

One of the 20th century’s most famous mathematicians, John von Neumann played an important role in the development of the modern computer. So it may be no surprise that he also pioneered analytical techniques for studying the properties of finite-difference equations. His approach to evaluating the computational stability of a difference equation employs a Fourier series method and is best described in References 1 and 2.

20 세기의 가장 유명한 수학자의 한사람 인 John von Neumann은 현대 컴퓨터의 기초를 세운 업적들로 알려져 있습니다 만, 유한 차분 방정식의 특성을 조사하기위한 분석 기법도 개발하고 있습니다 .  그 방법은 차분 방정식의 계산 안정성을 평가하기 위해 푸리에 급수를 이용한 접근 방식을 취하며, 참고 문헌 1 및 2 에 대한 세부 사항을 설명하고 있습니다.

The approach von Neumann used is based on two assumptions. First, that the difference equation can be linearized with respect to a small perturbation in the solution. Second, that there exists a set of overlapping regions in each of which the coefficients of terms in the difference equation for the perturbation can be considered locally constant. When these conditions are satisfied a Fourier series can be used in each region to compute the local behavior of the perturbation.

von Neumann의 접근 방식은 두 가지 가정에 근거하고 있습니다.  차분 방정식은 솔루션에 포함 된 작은 섭동에 대해 선형화 할 수있는, 그리고 섭동에 대한 차분 방정식의 항의 계수를 국소적으로 상수로 간주 될 수 겹치는 영역이 존재하는 것입니다.  이러한 조건이 충족 될 때, 각 영역에서 푸리에 급수를 이용하여 섭동의 국소 적 거동을 계산할 수 있습니다.

The solution for the perturbation Pjn, at discrete location j and time step n, is sought as a sum of Fourier terms having the form

시간 단계 n과 공간의 이산 위치 j의 섭동 P j n은 다음의 형식을 취해 푸리에 항의 합으로 구해집니다.

(4)    

where xj is the location of the jth finite-difference point, k is the wave number of the mode and r is a time amplitude. Because of the linear nature of the difference equation, it is only necessary to consider the behavior of a typical mode. Substitution of the typical mode into the linearized difference equation leads to an algebraic equation for the time amplitude r that must be satisfied if there is to be solution of the form postulated. A full solution for the perturbation could then be constructed from a sum of Fourier modes over a range of wave numbers.

여기서 xj는 j 번째의 유한 차분 점의 위치, k는 모드의 주파수, r은 시간 진폭입니다.  차분 방정식의 선형 특성에 따라 전형적인 형태의 거동을 고찰하기 만하면 됩니다.  전형적인 모드를 선형화 된 차분 방정식에 대입하여 추정 된 형식의 솔루션을 얻기 위해 충족되어야 시간 진폭 r의 대수 방정식을 얻습니다.  그러면 범위의 주파수에 걸쳐 푸리에 모드의 합에서 섭동의 전체 솔루션을 얻을 수 있습니다.

Numerical stability implies that as time increases (i.e., with increasing n) the magnitude of each mode must not grow unboundedly, and this means that the magnitude of r must be less than or equal to unity.

수치적으로 안정되어 있다는 것은 시간이 진행하는 (즉 n이 증가)에 따라 각 모드의 크기가 제한없이 증가하고 안되어, 그것은 r 값의 크기가 1 이하 이어야한다는 것을 의미합니다.

An Example of Linear Stability Analysis

An example will make von Neumann’s technique clear. For illustration purposes an advection-diffusion equation is a good test case,

von Neumann의 방법을 명확하게 하기 위해, 예를 이용하여 설명합니다.  여기에서는 테스트 케이스에 적합 이류 확산 방정식을 예로 들어 있습니다.

(5)    

which describes the convection and diffusion of a function u(x,t). The convection velocity c and diffusion coefficient ν are constant. Solutions of this equation are known to be well-behaved (i.e., physically stable).

이것은 함수 u (x, t)의 대류와 확산을 나타내는 식입니다.  대류 속도 c와 확산 계수 ν는 정수입니다.  이 방정식의 해는 양호한 거동을 나타내는 것을 알 수 있습니다 (즉, 물리적으로 안정되어있다).

A simple finite-difference approximation using constant time steps δt and spatial steps δx for ujn= u(jδx,nδt) is

상수 인 타임 단계 δt 및 공간 단계 δx를 사용하여 u j n = u (jδx, nδt) 한 다음 간단한 유한 차분 근사 식을 얻습니다.

(6)    

Inserting von Neumann’s mode Eq. 4 to this equation, with xj=jδx, and solving for r gives

이 식은 von Neumann 모드 식 4를 대입하고 xj = jδx로 r에 대해 해결합니다.

(7)    

Because this involves real and imaginary terms the square of the magnitude of r is given by the sum of the squares of the real and imaginary parts,

이 식에는 실수 항과 허수 부분이 포함되어 있으며, r의 절대 값의 제곱은 실수 부와 허수 부를 각각 제곱 한 합계에서 구할 수 있습니다.

(8)    

From this result, limits on ν, c, and δt can be determined that ensure that r remains less than or equal to 1.0 the requirement for not having an unbounded solution of the form Eq. 4.

In general, the extreme values for r occur when the sine and cosine factors are at their extreme values of 0.0, 1.0 or -1.0. For instance, if kδx=0 then r=1 and this mode will be neutrally stable, neither increasing nor decreasing with increasing values of n (i.e., rn=1). When kδx=π the value of r2 is (1-4νδt/δx2)2 and this implies the stability condition

따라서 식 4 형식으로 답을 비 묶여 있어서는 안된다는 요구 사항, 즉 r이 1.0 이하로하는 조건을 만족하는 ν, c δt에 대한 제한이 결정됩니다.

기본적으로, r이 극값이되는 것은, sin과 cos의 인수가 각각의 극치 인 0.0,1.0 또는 -1.0 때입니다.  예를 들어, kδx = 0 일 때 r = 1입니다.  이 모드는 중립 안정이며, n의 증가에 따른 값의 증가 / 감소되지 않습니다 (즉, r n = 1).  kδx = π 일 때 r 2 값은 (1-4νδt / δx 2) 2이며 다음의 안정 조건이 도출됩니다.

(9)   

And finally, when kδx=π/2 the stability condition r2

또한 kδx = π / 2의 경우 안정적인 조건 r 2 <1은 다음과 같이됩니다.

(10)    

Since the left side of the inequality is a sum of squares, neither term can exceed 1 without the left side exceeding 1.0, thus we have two conditions for stability solutions

부등식의 좌변은 제곱의 합이기 때문에 어느 부분이 1을 초과하면 왼쪽이 1.0을 초과하게되어, 안정적인 솔루션을 얻기위한 두 가지 조건이 인도됩니다.

(11)        and  

The first condition is a new stability condition, but the second condition has already been made more restrictive in Eq. 9. In addition to these conditions, Eq. 10 admits of one more condition that can be evaluated by expanding the first term and rearranging the inequality,

첫 번째 조건은 새로운 안정 조건이지만, 두 번째 조건은 더 까다로운 것이 이미 식 9로 표시되어 있습니다.  이러한 조건 이외에 식 (10)의 제 1 항을 전개하고 불평등을 정렬하여 구해지는 또 하나의 조건이 부여됩니다.

(12)    

The minimum value on the right side, according to Eq. 9, is 1.0, which finally leads to a condition for the minimum size of diffusion ν. In total, there are three stability conditions to be satisfied,

식 9에서 우변의 최소값은 1.0입니다.  따라서 확산 계수 ν의 최소 크기에 대한 조건이 부여됩니다.  마지막으로, 충족되어야 할 3 가지 안정 조건이 인도됩니다.

(13)              

When these three conditions are satisfied no Fourier modes will grow as n increases and so no unbounded instability develops. The first condition, in particular, is interesting because it says there will be unstable solutions for any non-zero time step unless there is a sufficient amount of diffusion. In contrast, the remaining conditions are simply restrictions on the maximum time step size.

이러한 세 가지 조건이 충족 될 때 n이 증가함에 따라 푸리에 모드가 증가하는 것은 아니라 제한없는 불안정성이 발생하지 않습니다.  특히 첫 번째 조건은 충분히 확산되지 않는 한 0이 아닌 임의 시간 단계에 해가 불안정해질 수를 나타내고 있다는 점에서 매우 흥미 롭다고 할 수 있습니다.  반면 나머지 두 조건은 최대 시간 단계에 대한 단순한 제한에 지나지 않습니다.

At this point it’s worthwhile observing that a violation of the first condition in Eq. 13 leads to a solution with purely exponential growth, while violation of either of the second two conditions exhibits a growth that not only grows exponentially, but also changes sign with each increment in time. The reason for this behavior and other useful information about finite-difference equations is explained in the second CFD-101 article on stability, Heuristic Analysis.

여기서, 식 (13)의 첫 번째 조건이 충족되지 않은 경우 해석은 기하 급수적으로 증가 할뿐입니다 만, 나머지 2 개의 조건의 어느 한쪽이라도 충족되지 않은 경우에는 지수 함수 적인 증가뿐만 아니라 시간 증분마다 부호가 반전하는 것에 주목하십시오.  이 동작을 가져다 원인 및 유한 차분 방정식에 관한 기타 유용한 정보는 안정성을 설명하는 “CFD-101″의 다음 문서 “휴리스틱 분석”을 참조하십시오.

In summary, the von Neumann type of Fourier analysis of finite-difference equations is quite useful provided the equation(s) are linear and have constant coefficients within a set of overlapping regions. For more general situations other methods must be sought to analyze computational stability.

이처럼 대상 유한 차분 방정식이 선형이며, 겹치는 영역 세트에서 상수 계수가 사용되는 경우에는 von Neumann 유형의 푸리에 분석이 매우 유용하다고 할 수 있습니다.  보다 일반적인 상황에 대해 계산 안정성을 분석하는 경우에는 다른 방법을 고려해야합니다.

End Note

Readers who may be a bit rusty with exponentials and imaginary numbers may find the following set of identities useful for performing a von Neumann type analysis:

지수와 허수에 관한 노트로 von Neumann 유형의 분석을 수행하는 데 도움이되는 항등식은 다음과 같습니다.

References

1. B.G. O’Brien, M.A. Hyman and S. Kaplan, “A Study of the Numerical Solution of Partial Differential Equations,” J. Math. Phys. 29, 223 (1951).
2. J. von Neumann, Collected Works, IV, Macmillan Co., New York, NY (1963).