Week 5 - Literature review: ODE Stability

When does a ODE become unstable?

An ordinary differential equation (ODE) is a differential equation with an independent variable and the derivatives of that variable. In computational mechanics, ODE is solved using a numerical integration method and the solution is said to converge when the absolute error between the iterations is within the tolerance defined. However, it is note that the error between the numerical solution and the analytical solution is never zero. It is this error which determines the stability of the solution.

For a stable solution, the error term is bound to decrease with every iteration. However, in certain cases, it is observed that solution diverges and this can be noted at spikes or divergence of the values with iterations. It is to be noted that the factors influencing the stability of the solution in CFD are mainly grid size and time step dt.

To illustrate the stability analysis, we use Python programming to compare the results of explicit, implicit and analytical solutions.

The explicit method uses forward difference scheme whereas the implicit scheme is based on the back difference scheme. The analytical solution is determined using computer algorithms available at Wolfram Alpha.

The python program is as follows:

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

def explicit(total_t,dt):

	t=np.arange(0,total_t,dt)
	y=np.zeros(len(t))

	for i in range(0,len(t)-1):

		y[i+1]=y[i]+ dt*(-1000*y[i]-math.exp(-t[i]))

	return t,y


def implicit(total_t,dt):

	t=np.arange(0,total_t,dt)
	y=np.zeros(len(t))

	for i in range(1,len(t)):

		y[i]=(y[i-1] - dt*math.exp(-t[i-1]))/(1+1000*dt)

	return t,y


def analytical(total_t,dt):

	t=np.arange(0,total_t,dt)
	y=np.zeros(len(t))

	for i in range(1,len(t)):

		y[i]=(math.exp(-1000*t[i])*(1-math.exp(999*t[i])))/999

	return t,y	


total_t=0.1
dt=1e-3
t,y_explicit = explicit(total_t,dt)		  	
t,y_implicit = implicit(total_t,dt)
t,y_analytical = analytical(total_t,dt)
plt.plot(t,y_explicit,color='blue')
plt.plot(t,y_implicit,color='red')
plt.plot(t,y_analytical,color='black')
plt.title('Solution for dt = %f'%dt)
plt.legend(['explicit','implicit','analytical'])
plt.show()

For a time step dt=1e-3, we find the below plot

Note that the results are very close for the numerical solutions obtained using explicit and implicit schemes when compared to that of analytical solution. However, it can be noted that implicit solution is more accurate as compared to explicit solution.

When we increase the time step size dt to 1.1e-3, we find that explicit solution starts to develop instability.

This instability is more significant for higher values of dt.

It is noted that implicit solution gives stable solution even at higher values of dt. However, a closer look would suggest an increase in the error value when compared to that of the analytical solution.

The time step is critical in determining the stability of solution when using an explicit scheme. Thus, we define constraints using CFL (Courant–Friedrichs–Lewy) number which is dependent on the grid size and the time step. For example, for a heat convection equation, solution is found to be consistent when CFL number is less than 0.5 and the solution diverges when CFL number is more than 1.

How does stability affect engineering simulations?

For engineering simulations, there needs to be a balance between the processing expenditure and the accuracy of the solution. For explicit schemes, it is observed that less processing power is consumed at the expense of accuracy of solution as there is instability in solutions defined using explicit scheme (as constraining the solution with appropriate boundary conditions is complex). On the other hand, implicit schemes tend to give stable solution but consume tremendous processing power. Thus, the effect of stabilty on engineering simulations plays a critical role in determing the appropriate solution algorithm.

Can stability condition be derived for all types of problems?

The stability condition of a particular problem is derived using the CFL number. When the condition number is less than maximum allowable value, the solution tends to be stable. The maximum value, in turn, depends on the method used. For implicit method, larger values of CFL number can be tolerated. However, determining the exact value of the allowable value is very difficult considering the complex set of equations solved for a CFD case. Thus, it is suggested to give preference to implicit method of solution to improve stabilty of solutions.