Week 5.2 - Literature review: ODE Stability

AIM:  To understand the stability of an ODE and its effect on engineering simulation.

Ordinary Differential Equation: It contains derivatives of one or more functions of an Independent variable. We need to determine the  function or set of functions that satisfy the equation to solve an ODE.

When does a ODE become unstable?

In engineering applications numerical integrations schemes are used for getting the approximation for the ODE's and the absolute error between iterations gives us an idea about the stability criteria of the Solution. For time integrations these numerical integrations can be classified as explicit or implicit scheme. Error between the iterations sholud decrease for stable solution. Discrete point density in the domain of iteration and time step size in the iterative method are also important for the Stability of the Solution that why time step size is less and number of discrete point should be larger. If the solution remains stable for the arbitrarily large step size then solution is called unconditionally stable. If either of the properties doesn't match then solution is called unstable.

Implicit Scheme: Euler Method or forward differencing

Explicit Scheme: Backward Euler method of backward differencing

For system stabilty, Courant number(non-dimensional number) also used to evaluate the time step requirements of a transient simulation for a given mesh size and flow velocity. 

How does stability affect engineering simulations?

In a unstable system, output values does not satisfied the system and create the difference with the actual experimental results. For validation and verification purpose an unstable ODE's give unrealistic results which are not acceptable.

Can stability condition be derived for all types of problems?

Stability of an ODE can be illustrated by calculating the analytical solutions and approximations using implicit and explicit technique. ODE used is:

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

Solving the differential equation

analytical solution:

`y(t) = frac{e^(-1000t)(1-e^(999t))}{999}`

Implicit Solution:

`frac{y^n-y^(n-1)}{dt} = -1000y^n-e^(-t)`

`y^n = frac{y^(n-1)-e^(-t)dt}{1+1000dt}`

Explicit Solution:

`frac{y^(n+1)-y^(n)}{dt} = -1000y^n-e^(-t)`

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

Python Code 

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

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

Analytical solution
y(t) =  (e^(-1000t)(1-e^(999t)))/999

Explicit Method
y^(n+1) = y^n - (1000y^n + e^(-t))dt

Implicit Method 
y^(n) = (y^(n-1) - e^(-t))dt/(1+1000dt)
'''

total_t = 0.1
dt = 1.2e-3
t = np.arange(0,total_t,dt)
y = np.zeros(len(t))

def analytical(total_t,t):
    t = np.arange(0,total_t,dt)
    y = np.zeros(len(t))
    for i in range(1,len(t)):
                 y[i] = (mt.exp(-1000*t[i])*(1-mt.exp(999*t[i])))/999
    return[t,y]

def explicit(total_t,t):
    t = np.arange(0,total_t,dt)
    y = np.zeros(len(t))
    for i in range(1,len(t)):
        y[i] = y[i-1]-(1000*y[i-1] + mt.exp(-t[i-1]))*dt
    return[t,y]

def implicit(total_t,t):
    t = np.arange(0,total_t,dt)
    y = np.zeros(len(t))
    for i in range(1,len(t)):
        y[i] = (y[i-1]- (mt.exp(-t[i-1]))*dt)/(1+1000*dt)
    return[t,y]

t,y_explicit = explicit(total_t,t)
t,y_implicit = implicit(total_t,t)
t,y_analytical = analytical(total_t,t)
plt.plot(t,y_implicit,'b', label = 'Implicit')
plt.plot(t,y_explicit,'g', label = 'Explicit')
plt.plot(t,y_analytical,'r', label = 'Analytical')
plt.legend()
plt.show()

 Results 

Step Size: 1.5e-4

Step Size: 1e-3

Step Size: 1.2e-3

Step Size: 1.9e-3

It can be seen here as the time step size increased from 1e-4 to 1.9e-3,

In explicit solution- a small error generated and starts to amplify with successive increase in the step size.

In implicit solution: Plot can be seen smooth here as compared to analytical solution. Even though some gap can be seen at certain areas of inflection which can be justified by numerical approximation. 

An expression has been used always to solve the problem that is continous also in its inputs. And slight change in input can lead slight change in output also but if the expression gives large change in output with slight change in input, that considered as ill conditioned. The exact cut off between ill-conditioned and well-conditioned expression depends upon context of problem and use of the results.