Week 5 - Literature review: ODE Stability

Aim: Analysis of the stability of ODE solvers

Objective:

• Understanding how and when an ODE becomes unstable
• Understanding the effect of ODE stability on engineering simulations
• Checking stability conditions for a variety of problems

Theory:

In mathematics, Ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions. The term 'ordinary' is used in contrast with the term partial differential equation, which may have more than one independant variable.

ODEs can be solved using various numerical methods usch as numerical integration and numerical discretization. Some ODEs can be solved through mathematical procedures such as (separation methods, Jacobian method, ..). However, in solving engineering problems, one may arrive at a system of ODEs that may be coupled or uncoupled and cannot be solved directly, hence, we need to implement numerical methods. Different numerical methods have different levels of accuracy depending on the way it converges to the solution, and hence can have their own stability or unstability factors. These numerical schemes can be broadly classified into explicit and implicit schemes.

For time-explicit schemes, the value of a 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.

When does ODE become unstable?

The explicit schemes faces the problem of unstability because of its dependance on the value or size of time step chosen. i.e. the degree of unstability depends on the time step being chosen. The ideal value of time-step depends on the given ODE. When time step value is higher than the ideal value, oscillations in solution are observed. For extremely high time step, the solution gets completely diverged from the required solution. Implicit schemes are inherently stable as they use iterative techniques for a given time-step.

To observe this effect, we take the following example:

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

The DE can be approximated through discretization as:

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

where, `n` denotes the time steps and `deltat` is the size (value) of time step. Hence, we have

For Explicit:

`(y^n - y^(n-1))/(Deltat) = - 1000y^(n-1) - e^{-t^n}`

`=> y^n = (1 - 1000 Delta t) y^(n-1) - e^{-t^n} Delta t`

For Implicit:

`(y^n - y^(n-1))/(Deltat) = - 1000y^(n-1) - e^{-t^n}`

`=> y^n = (y^(n-1) - e^{-t^n})/(1 + 1000 Delta t)`

The analytical solution for the ODE for initial condition `y(t=0)` is

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

No, we compare explicit and implicit solution at different time steps using various time steps:

"""
Program to compare explicit and implicit solution at various time steps
"""

from numpy import *
from matplotlib.pyplot import *
"""
Program to compare explicit and implicit solution at various time steps
"""

from numpy import *
from matplotlib.pyplot import *
from math import *

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

for i in range(1,20):
    t = arange(0, 0.1, step=(i*1e-4))
    dt = t[1] - t[0]
    y_expl = [0]*len(t)
    y_impl = [0]*len(t)
    y_anal = [(exp(-1000*tp)/999) - (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] - exp(-t[i])*dt)/(1+1000*dt)

    fig, ax = subplots()
    plot(t,y_expl,color='red',linewidth=2,label='Explicit Solution')
    plot(t,y_impl,color='blue',linewidth=2,label='Implicit Solution')
    plot(t,y_anal,'--',color='black',linewidth=2,label='Analytical Solution')
    grid('on')
    ticklabel_format(axis='y',style='sci',scilimits=(0,0))
    xlabel('Time (sec)',fontsize=14)
    ylabel('Y',fontsize=16)
    legend(fontsize=14)
    title('Time step = ' + str(dt),fontsize=16, fontweight='bold')
    savefig('Solutions_for_' + str(dt) + ".png")
    show()

Results:

beginning of oscillations

oscillations increase with increase in deviation from ideal time step

Hence, we observe that a time step in range of 1e-4 `(1 * 10^-4)`. And a deviation from ideal time step leads to oscillations in solution in explicit method.

How does stability affect engineering simulations?

Engineering simulations needs stability in function and inputs. Hence, to obtain a solution, a stable metod must be chosen. Or else, the soluion would be correct mathematically, but may not feasible or can be impractical (for example a function returning a negative mass or a function calculates speed greater than speed of light). Hence, a proper scheme needs to be adopted for obtaining accurate results.

Can stability conditions ber derived for all types of problems?

Stability condition can be derived for all problems and it depends on the given ODE. For a linear equation, (our example) stability can be obtained by changing time step. However, for second order ODEs, stability is harder to control because it is dependent on eigen values. Hence, stability can be induced by adding a suitable contol method.