A study on Numerical Stability Analysis for ODE systems

Objective:

To study and understand the concept of stability of ODE systems and answer the following questions with the backing of suitable computer programs.

1. When does an ODE become unstable?

2. How does stability affect engineering simulations?

3. Can stability condition be derived for all types of problems?

Solution:

We know that in CFD, we numerically simulate a physical system by solving its governing equations. These governing equations are either partial differential equations or ordinary differential equations. To solve these differential equations, we use various kinds of discretization schemes. However, not all discretization schemes may result in the correct or desired solution. Therefore, it is essential to understand why and when a discretization scheme may give an undesired and sometimes a totally unphysical solution. The two factors that determine this are stability and consistency. Let us define each factor.

Stability:

A scheme is said to be stable if the errors do not grow as the solution advances. The errors directly correspond to the residuals in the system. If the residuals do not grow, the solution will eventually converge a steady state solution.

Consistency:

A scheme is said to be consistent if the solution approaches the actual analytical solution as the grid size approaches to zero. The error in a solution is directly proportional to the grid size. When the grid size tends to 0, the errors also should tend to zero and the solution also should tend to the analytical solution.

If a scheme satisfies both stability and consistency, then the solution is said to be convergent i.e., it will reach a steady state, desired physical result. Though at most times, a stable solution is also a consistent one, that is not always true.

So, when does an ODE become unstable?

Let us understand the stability criteria for ODE using this general form of ODE

`dy/dt = -ay`

The analytical solution for this ODE is be given by

`y = y_0 e^(-at)`

The numerical solution can be computed using the explicit Euler as

`(y^(n+1) - y^n)/h = -ay^n`

`y^(n+1) = y^n - ay^nh`

`y^(n+1) = y^n(1 - ah)`

We know that the instability occurs only when the value of the function at (n+1)th iteration, y^(n+1) is greater in magnitude than the value of the function at nth iteration, y^n as the solution advances. This means that the value of |1-ah| < 1 i.e., h < 2/a for the solution to be stable.

This implies that whenever h > 2/a, the solution is unstable.

Let us demonstrate this using a simple ODE

`dy/dt = -1000y + 3000 - 2000e^-t`

The analytical solution for the initial condition y(0) = 0, is

`y = 3 - 0.998 e^(-1000 t) - 2.002 e^-t`

The numerical solution using explicit Euler or forward difference can be given by

`(y^(n+1) - y^n)/h = -1000 y^n + 3000 - 2000 e^-t`

`y^(n+1) = y^n + h(-1000 y^n + 3000 - 2000 e^-t)`

We know that for stability, h < 2/a, where a is the constant that multiplies the function on the right hand side.

Therefore, for the ODE to be stable the time step size, h < 2/1000, or h < 0.002. Let us corroborate the stability criteria graphically using a program.

(i) h < 0.002

For h < 0.002, the solution was stable as shown below;

(ii) h = 0.002

When h = 0.002, the solution does not have exponentially growing errors i.e., the errors remain constant but it does not converge as well.

(iii) h > 0.002

For h > 0.002, the errors start growing and the solution starts to diverge rather than to converge.

Effect of Stability on Engineering Simulations

We know that an unstable solution is characterized by exponentially growing errors. These errors or residuals are a direct function of the derivatives of the quantities being solved in the solution. These residuals directly correspond to the spurious oscillations of the solution and would result in a totally unphysical solution. Such a solution is doomed to fail and cannot be used in real-time.

On the contrary, a stable solution is characterized by decreasing residuals. This means that the oscillations in the solution decrease as the solution advances and over time settle to give a steady state value. Such a solution is what is desired while performing numerical simulations.

It is important to note that a barely stable solution is not enough for engineering simulation which can be achieved by satisfying the stability condition. The solution also needs to be accurate which can be achieved by further minimizing the time step size.

Can the stability condition be derived for all types of problems?

The stability condition can only be derived for linear ODE systems. However, for non-linear ODEs, we can predict the stability. This is done by creating a Jacobian matrix(i.e., a matrix containing the first-order derivatives of all function w.r.t. all independent variables) and computing all its eigenvalues. If all the eigenvalues are negative, it means that the solution is stable and even if one of the eigenvalues is positive, it means that the solution will be unstable. This method of predicting stability using eigenvalues of the Jacobian matrix can also be used for linear systems.

Python Script

\"\"\"
Stability analysis for the ODE
dy/dt = -1000*y + 3000 - 2000*e^-t

The analytical solution is
y = 3 - 0.998*e^(-1000*t) - 2.002*e^-t

The numerical explicit solution is
y_(n+1) = y_n + h*(-1000*y_n + 3000 - 2000*e^-t)
\"\"\"

import numpy as np
import matplotlib.pyplot as plt

t_end = 0.5	# simulation time
h = 0.002	# step size
t_array = np.arange(h, t_end, h)
time = np.arange(0, t_end, h)
y0 = 0	# initial value
y_numerical = [y0]
y_analytical = [y0]
y_old = y0

# loop for solution advance
for t in t_array:
	
	# numerical solution using explicit Euler or forward difference
	y_new = y_old + h*(-1000*y_old + 3000 - 2000*np.exp(-t))
	y_numerical.append(y_new)
	y_old = y_new

	# analytical solution
	y_analytical.append(3 - 0.998*np.exp(-1000*t) - 2.002*np.exp(-t))

# plot
fig1 = plt.figure(1)
plt.plot(time, y_numerical, alpha = 0.7, label = \'Numerical Sol\')
plt.plot(time, y_analytical, alpha = 0.7, label = \'Analytical Sol\')
plt.legend(loc = \'best\')
plt.xlabel(\'time\')
plt.ylabel(\'y = f(t)\')
fig1.suptitle(\'Stability analysis for a simple ODE for time step h = {}\'.format(h))
plt.show()