Literature Review ODE Stability

The aim of this challenge is to discuss the stability of an ODE/PDE system.

Introduction:

Before we begin, it is essential to understand the different ways in which an ODE/PDE system can be solved numerically. Numerical solution schemes can be classified into two main types viz., implicit and explicit schemes. In an explicit scheme, a direct computation of the dependent variables is made in terms of known quantities whereas, when the dependent variable is defined by coupled sets of equations and either a matrix or iterative technique is needed to obtain the solution, the numerical method is said to be implicit.

Let us discuss one example.

The selected ODE is as follows:

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

The analytical solution for the above ODE is as follows:

`(e^(-1000 t) (1 - e^(999 t)))/999`

Python Code:

\'t\' varies from 0 to 0.1.

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

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 y

total_t=0.1
dt=1e-4
t=np.arange(0,total_t,dt)
y_analytical = analytical(total_t,dt)
plt.plot(t,y_analytical,color =\'orange\')
plt.xlabel(\'t\')
plt.ylabel(\'y\')
plt.legend([\'analytical\'])
plt.show()

Analytical Solution:

Now we will try to implement the implicit scheme and compare our results with the analytical solution.

Implicit Scheme:

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

Upon discretizing the above equation using a backward step, the equation is as follows:

`(y^n-y^(n-1))/(∆t)=-1000y^n -e^(-t)`

upon rearranging the terms we get

`y^n(1+1000∆t)+e^(-t)∆t=y^(n-1)`

Here the values of the previous step are dependent on the next step. Hence, it is implicit. Let us rearrange it further to make the values of y to be solvable for the present step.

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

import matplotlib.pyplot as plt
import numpy as np
import math 
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 y

def explicit_euler(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]*(1-1000*dt)-math.exp(-t[i-1])*dt
	return [t,y]

def implicit_euler(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]-math.exp(-t[i-1])*dt)/(1+1000*dt)
	return [t,y]

total_t=0.7
dt=1e-4

y_analytical = analytical(total_t,dt)
t, y_explicit = explicit_euler(total_t,dt)
t, y_implicit = implicit_euler(total_t,dt)
error = abs(y_analytical-y_explicit)

At dt=1e-4 Implicit Scheme

At dt=1e-4. the results obtained by the implicit scheme are the same as those obtained in an analytical way.

Error Graph for Implicit at dt =1e-4

Since we have used the differencing scheme of the first order, the error is of order ∆t (i.e., the error is less than 1e-4) as our dt is 1e-4

At dt=1e-1

Upon increasing the time step with the implicit scheme the error is higher but still the pattern of the solution obtained by the implicit scheme comply with the analytical solution.

Explicit Scheme:

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

Upon discretizing the above equation using a forward step, the equation is as follows:

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

upon rearranging the terms we get

`(y^(n+1) =(1-1000*∆t)y^n -e^(-t)(∆t)`

Here the values of the present step are dependent on the previous step. Hence, it is explicit.

At dt=1e-3 Explicit Scheme

At dt=1e-3. the results obtained by the explicit scheme are the same as those obtained in an analytical way. The order of error is ∆t.

Now let\'s increase the step size.

At dt = 1.5e-3

At dt=1e-2

With the increase in the step size, we find that the results overshoot and errors are of the order 1e5. This is due to the failure of the stability criterion with the explicit schemes.

As it is evident from the above case study, implicit schemes are unconditionally stable whereas, the explicit schemes are conditionally stable. ODE\'s solved by explicit method becomes unstable owing to the larger time steps and spatial steps.

The stability of ODE is important in Engineering problems because it is used to obtain the converged and desired solution. If the set of equations formed from ODE loses its stability then the results are diverged and are pretty garbage.

Stability:

Finite Difference Equations of ODE is stable if it produces a bounded solution for a stable ODE. Several methods have been devised to analyze the stability of ODE\'s. Few of them are as follows:

1. The discrete perturbation method.
2. The von Neumann method.
3. The matrix method.

Let us derive the stability criterion for the above problem by using Von Neumann Method.

The given equation is:

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

Upon discretizing the above equation using a forward step, the explicit equation is as follows:

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

upon rearranging the terms we get

`(y^(n+1) =(1-1000*∆t)y^n -e^(-t)(∆t)`

From Von-Neumann Stability Analysis:

`Gy^n=y^(n+1)` & `y^n ≈e^(-t)(∆t)`

Upon putting it in the above equation

we get,

`Gy^n=-1000*∆t*y^(n)`

The criterion for the VN-Analysis is the amplification factor G has to be:

`|G|<=1` implies

`-1<=-1000*∆t<=1`

We can divide the above equations into two equations,

`-1000*∆t<=1` which is always true since `∆t>=0` it is redundant

let us try the other equation

`-1<=-1000*∆t`implies

`1>=1000*∆t` implies

`∆t<= 1/1000 =1e-3`

Hence in the above problem, for time step greater than 1e-3 like we have seen in the above results. Oscillations, deviation and overshooting appear. If that constant was 1e4 then the criterion for the time step would be 1e-4.

Yes, strictly speaking, the criterion for all the stable ODE problems could be defined but for the complex set of equations involving unstable ODE determining the criterion would involve complex calculations.

Disadvantages of Implicit Scheme:

It is computationally expensive and requires more time to solve the set of FDE.

Disadvantages of Explicit Scheme:

It may be fast and inexpensive computationally but stimulating the components involving shocks and high-frequency change would result in diverged solutions.