Literature review: ODE Stability

INTRODUCTION:

Ordinary Differential Equations (ODE) can be solved through various numerical methods, including numerical discretization and numerical integration. Some Odes can be solved directly through mathematical integration. However, most practical phenomena contain a system of ODEs and PDEs which may be coupled or uncoupled and cannot be solved directly. For this purpose, numerical schemes are implemented to find the solution of the system of equations. different numerical schemes have different levels of accuracy and can be stable or unstable. Numerical schemes can mainly be classified into explicit and implicit schemes. For time-explicit schemes, the value of the variable is calculated using the known values from the previous time steps. For time-implicit schemes, the value of a variable is defined using a coupled set of equations including the value at that time-step and is solved using iterative techniques.

Instability of ODEs

Implicit schemes are inherently stable as they use iterative techniques for a given time-step. The same cannot be said about explicit schemes. Explicit scheme stability is highly dependant on the values of the time-step chosen. The ideal value of the time step depends on the given ODE. The amount of instability depends on the time-step. If the time-step value is higher than this ideal value, the solution becomes unstable by producing oscillating results. The higher the value of the time-step, the more oscillating the solution. Foy very high time-step values, the solution completely diverges from the actual solution. A lower time-step ensures stability but is computationally more expensive. To demonstrate the effect of the time-step on the stability of a solution, consider the following ODE:

`(dy)/(dt)=-1000y-e^(-t)`

The differential `(dy)/(dt)`can be approximated through discretization:

 `(dy)/(dt)=(y^n-y^(n-1))/(Deltat)`

where the superscript n denotes the time-step and `Deltat` is the value of the time-step chosen. The ODE can be written in explicit and implicit schemes as:

Explicit :

`(y^n-y^(n-1))/(Deltat)=-100y^(n-1)-e^(-t^n)rArry^n=(1-1000Deltat)y^(n-1)-e^(-t^n)Deltat`

Implicit:

`(y^n-y^(n-1))/(Deltat)=-1000y^n-e^(-t^n)rArry^n=(y^(n-1)-e^(-t^n)Deltat)/(1+1000Deltat)`

The analytical solution for this ODE for an initial condition of `y(t=0)=0`is given as:

`y=(e^(-1000t))/999-(e^(-t))/999`

The explicit and implicit solution is compared to the analytical solution using various time steps:

import numpy as np
import matplotlib.pyplot as plt
import math

def yprime(y,t):
	return -1000*y - math.exp(-t)

t = np.arange(0, 0.1, step=0.0001)
dt = t[1] - t[0]
y_expl = [0]*len(t)
y_impl = [0]*len(t)
y_analytical = [(math.exp(-1000*tp)/999)-(math.exp(-tp)/999) for tp in t]

for i in range(1,len(t)):
	y_expl[i] = y_expl[i-1] + yprime(y_expl[i-1],t[i])*dt
	y_impl[i] = (y_impl[i-1] - math.exp(-t[i])*dt)/(1+1000*dt)

fig, ax = plt.subplots()
plt.plot(t,y_expl,color='red',linewidth=2,label='Explicit Solution')
plt.plot(t,y_impl,color='blue',linewidth=2,label='Implicit Solution')	 
plt.plot(t,y_analytical,'--',color='black',linewidth=2,label='Analytical Solution')	
plt.grid('on')

plt.xlabel('Time (sec)',fontsize=14)
plt.ylabel('Y',fontsize=16)
plt.legend(fontsize=14)
plt.title('Time step = ' + str(dt) , fontsize=16, fontweight='bold')
plt.show()

we can see that `Deltat=0.0001` is the minimum value required for the stability of the explicit solution, whereas the implicit solution is inherently stable for any value of `Deltat`.

How does stability affect engineering simulations?

Stability is a very important aspect when performing engineering simulations. To obtain a solution, the numerical schemes chosen must be stable. Sometimes, even unstable schemes provide a solution, but that solution would not be accurate and might even be physically impossible. Hence, a proper study of the stability of the scheme is essential before the solution can be taken as nearly accurate.

One such measure of stability is to study mesh convergence. A Mesh convergence study involves performing the simulations for different mesh sizes and comparing the solution. The solution should return consistent values with a lesser and lesser difference with respect to the mesh sizes. If a mesh convergence study produces inconsistent results, then the numerical scheme might be unstable even though we are obtaining a solution.

Can stability conditions be derived for all types of problems?

Yes, stability conditions can be derived for all types of problems. The stability condition generally depends on the nature of ODE given. For a linear equation, as shown above, the stability conditions can be met by changing the value of the time-step. However, the stability of some equations might be uncontrolled such as second-order ODEs whose stability depends on eigenvalues of a particular matrix. This instability is reflected in the physical system as well and can be controlled by adding some control measures in the system itself. The two main types of stability analysis that are applied to numerical schemes are Von Neumann and Eigenvalue stability conditions.