Literature review: ODE Stability

Introduction:

In mathematics, when a certain small amount of disturbance is inserted into the system, if the system tend to minimize it, it is called as a stable system. If the disturbance amplifies as the system progresses in time and does not decrease in value, it is called as an unstable system.

In Numerical discretization, we insert a certain small error, mathematically, during approximation. This error needs to be checked every time its value increases, for the system to be stable.

A differential equation can be solved (numerically) by several numerical methods like Euler's method and Range-Kutta method. In this case, Euler's methods are of intrest.

Euler's method:

In Euler's method, if a curve is plotted with independent variable being at abscissa and dependent variable being at ordiante, the x-axis is divided into several sub intervals and tangents are drawn to the curve along the length of each sub interval. So, along the length of these sub intervals, the curve is considered linear.

The expression that can be derived from the process above is,

`y_(i+1) = y_i + h*f(x_i,y_i)`, where h is length of the sub interval and f(x,y) is any function.

This expression is very straightforward, considering the process mentioned above. This is Euler's method of solving ODEs explicitly (since `y_(i+1)` appears only on the left side).

However, the above expression is stable only for small values of h. When the length of sub interval increases, the error increases (as we are considering the profile to be linear along the length of this subinterval) and so the system becomes unstable.

To overcome this drawback, Euler's method is modified a little so that the system would be inherently stable, for any value of length of sub interval.

The expression for the modified Euler's method is,

`y_(i+1) = y_i + h*f(x_(i+1),y_(i+1))`

Modified Euler's method is an implicit scheme as `y_(i+1)` is appearing in both LHS and RHS.

The stability criteria can be understoor by using a general single step method equation,

`y_(i+1) = E*y_i`

For absolute stability, `|E| <= 1`, failing which the error shoots up to a very high value and gives un physical results.

Example:

Consider a differential equation,

`(dy)/(dx) = -1000*y - e^(-x)`

A. Euler's method (explicit)

Python code to solve the above equation by Euler's method (explicit scheme) is,

import matplotlib.pyplot as plt
import math

def f(x,y):
	return -1000*y-math.exp(-x)

x_0 = 0
x_end = 0.1
n = 500
h = (x_end - x_0)/(n-1)

y_initial = 1
x_initial = 0
y_final = []
for i in range(0,n):
	y_final.append(y_initial + h*f(x_initial,y_initial))
	y_initial = y_final[i]
	x_initial = x_initial + h

plt.plot(y_final)
plt.show()

When the above code is executed, the resultant plot is,

When the number of steps is changed to 50, the resultant plot is,

B. Modified Euler's method (implicit scheme):

The same above equation is solved implicitly using modified Euler's method. The expression after simplification becomes,

`y_(i+1) = (y_i - he^(-x_(i+1)))/(1+1000*h)`

Python code to solve the ODE implicitly is,

import matplotlib.pyplot as plt
import math


x_0 = 0
x_end = 0.1
n = 50
h = (x_end - x_0)/(n-1)

y_initial = 1
x_initial = 0
y_final = []
for i in range(0,n):
	x_initial = x_initial + h
	ans = (y_initial - h*math.exp(-x_initial))/(1+1000*h)
	y_initial = ans
	y_final.append(ans)
	

plt.plot(y_final)
plt.show()

Notice that the number of steps is at same value where the explicit scheme has failed i.e 50. When the above code is executed, the resultant plot is,

Based on the results above,

• The ODE becomes unstable in case of explicit scheme when the number of sub intervals are less or the value of length of sub interval is more. Implicit schemes are stable at all values of number of sub intervals
• In engineering simulations, there will always be a certain amount of error as we discretize and approximate, this error needs to be bounded troughout the simulation. So, stability criterion need to be met before simulating any engineering problem
• Stability conditions can be derived for all types of problems in which there is an approximation