Week 5 - Literature review: ODE Stability

INTRODUCTION

An ODE is also known as the ordinary differential equation is a form differential equation that contains one or more functions of one independent variable. The solution of an ODE is a function that satisfies the equation. Finding the solution by using analytical methods is not always possible in the case of differential equations, so to find the solutions numerical methods are employed. The numerical algorithms that are used for the calculation determine the accuracy of the solution.

One of the important characteristics of the solution algorithm is its stability. An algorithm or numerical scheme is said to be stable if, for a small change in the input, there must only be a small change in the output. In the numerical calculation process, small errors are introduced which can amplify as the solution marches in time and the result may blow up resulting in unphysical solutions.

An important aspect of ODEs is the condition of the expressions. If there is only a slight variation in the output for a slight variation in the input the expression is said to be well-conditioned. If there is a large variation in the output for a small variation in the input, the expression is said to be ill-conditioned.

It is also important to note that for a well-conditioned expression there is a stable way of solving and if the problem is ill-conditioned there is no reliable way to solve it.

STUDY OF A SAMPLE ODE

An ODE is selected for the study and the expression is kept as general as possible. The ODE used in the solution is gien below:

\begin{equation}
\frac{dy}{dt} + cy +e^{-t} =0
\end{equation}

It can also be written in the format :

\begin{equation}
\frac{dy}{dt} = f(t)
\end{equation}

where `f(t) = -cy-e^{-t}`

\begin{equation}
\frac{dy}{dt} = -cy-e^{-t}
\end{equation}

The numerical solution can be done in two different ways, the explicit Eulers method and the implicit Eulers method.

The explicit Euler method uses the values from the previous time step to calculate the value in the current time step.

\begin{equation}
\frac{y^{n+1}-y^{n}}{\Delta t}= -cy^{n} -e^{-t}
\end{equation}

\begin{equation}
y^{n+1}= y^{n} + \left(-cy^{n} -e^{-t}\right)\Delta t
\end{equation}

Where `y^{n+1}` is the current time step and `y^n` is the previous time step. It is interesting to note that the explicit method is not unconditionally stable. The stability of this method is highly dependent on the value of the time-step.

The implicit Euler method is shown below:

\begin{equation}
\frac{y^{n}-y^{n-1}}{\Delta t}= \left(-cy^{n} -e^{-t}\right)
\end{equation}

\begin{equation}
y^{n}= \left(-cy^{n} -e^{-t}\right)\Delta t + y^{n-1}
\end{equation}

\begin{equation}
y^{n}\left(1+ c\Delta t\right) + e^{-t}\Delta t = y^{n-1}
\end{equation}

Therefore we can say

\begin{equation}
y^{n} = \frac{y^{n-1}-e^{-t}\Delta t}{\left(1+ c\Delta t\right)}
\end{equation}

For an initial value of `y(0) = 0` we can derive the analytical solution as

`y = {e^{-ct}(1-e^{(c-1)t})}/{c-1}`

The stability of the equation can be controlled using two parameters, the constant `c` and the time step `\Delta t`.

CODE AND EXPLANATION

The code used to solve the above ODE is given below:

"""
Created on Thu Dec 17 20:01:35 2020

@author: Abhijith S Pai
"""

import matplotlib.pyplot as plt
import numpy as np


def explicit_euler(t_total,dt):
    t=np.arange(0,t_total,dt)
    y=np.zeros(len(t))
    for i in range(1,len(t)):
        y[i]=y[i-1]+dt*(-c*y[i-1]-np.exp(-t[i-1]))
    return [t,y]
    
def implicit_euler(t_total,dt):
    t=np.arange(0,t_total,dt)
    y=np.zeros(len(t))
    for i in range(1,len(t)):
        y[i]=(y[i-1]-np.exp(-t[i-1])*dt)/(1+c*dt)
    return [t,y]


        
def analytical_soln(t_total,dt):
    t=np.arange(0,t_total,dt)
    y=np.zeros(len(t))
    for i in range(1,len(t)):
        y[i]=(np.exp(-c*t[i])*(1-np.exp((c-1)*t[i])))/(c-1)
    return [t,y]

# Time known constrain
t_total=0.1
c=100
dt=1e-2

t,y_explicit=explicit_euler(t_total,dt)
t,y_implicit=implicit_euler(t_total,dt)
t,y_analytical=analytical_soln(t_total,dt)


plt.figure(figsize=(10,5))
plt.plot(t,y_explicit,color='blue')
plt.plot(t,y_implicit,color='red')
plt.plot(t,y_analytical,color='orange')
plt.xlabel('Time (s)')
plt.ylabel('y')
plt.legend(['Explicit','Implicit','Analytical'])
plt.savefig('odestability.png')

The different functions namely explicit_euler(t_total,dt), implicit_euler(t_total,dt)) and analytical_soln(t_total,dt)) is used to calculate the explicit, implicit and analytical solution of the ODE. The total time of integration is t_total and the time-step is depicted using dt. The variable c represents the constant value.

RESULTS AND DISCUSSION

From the execution of the program, it is found that as the value of c increases, the time step needs to be decreased to get a stable solution.

It is found that when `c\times dt = 0.1` there is an exact match between the explicit, implicit and analytical solution.

For c = 100 and dt = 1e-3 the solution plot is given below:

For c = 1000 and dt = 1e-4 the solution plot is given below:

For c = 10000 and dt = 1e-5 the solution plot is given below:

The instability arises when `c\times dt` is greater than 1.1. When the value further increases the vibration in the solution also increases.

For c = 100 and dt = 1.1x1e-2 the explicit solution starts to become unstable. The output is given below:

When the value of dt is increased to 1.8e-2 the explicit solutions exhibit more instability as shown below:

For c = 1000 and dt = 1.1x1e-3 there is a point where there is an instability in the explicit scheme, which is shown below:

The value of dt is increased to 1.8e-3 the instability is much more profound which is depicted in the following image:

For c = 10000 and dt = 1.1x1e-4 the instability starts and the results is shown below:

When the value of dt is increased to 1.8e-4 the instability oscillations also increases as shown below:

CONCLUSION

From the detailed analysis of the ODE system, it is found that the constant c and time-step dt. It is found that if `c\times dt = 1.1`, the solution starts to become unstable.

It is also interesting to note that when `c \times dt = 0.1` the error is very less and the implicit and explicit schemes give almost the same result.

It is also found that the implicit scheme is stable irrespective of the time step in the differential equation considered in the above study. A generalization is not possible since in some cases there is a chance that the implicit scheme becomes unstable.

For different types of ODEs, the stability criterion will be different. If the criterion is not satisfied the solution becomes unstable.

An unstable solution implies that the solution is not physical and the results out of such simulation will not be realistic which cannot be used in engineering calculations.

From the literature review, it can be found out that the stability criteria can only be developed for well-conditioned equations. If the equation is ill-conditioned there is no reliable way for a stable solution.

REFERENCE

https://www.cs.princeton.edu/courses/archive/fall11/cos323/notes/cos323_f11_lecture19_pde2.pdf

https://www.math.utk.edu/~ccollins/M577/Handouts/cond_stab.pdf