Week 5.2 - Literature review: ODE Stability

AIM: To carry out a literature review on ordinary differential equations stability.

Introduction

A stable problem is the one for which small changes in the initial conditions elicit only small changes in the solution. In other words, stability in mathematics, is the state in which a little change in a system does not have a significantly disruptive impact on that system. A function f(x) is referred to as stable in terms of the solution of a differential equation if any alternative solution of the equation that begins sufficiently near to it at x = 0 continues to be close to it for subsequent values of x. The solution is said to be asymptotically stable if the difference between the solutions decreases as x rises. A solution is referred to as unstable if it lacks one of these characteristics. Physical problems must have stable solutions because otherwise, the mathematical equations describing the problem will not be able to correctly predict the future. Slight deviations from the mathematical model brought on by unavoidable measurement errors must have an equivalently slight effect on the solution. Physical problems must have stable solutions because otherwise, the mathematical equations describing the problem will not be able to correctly predict the future. Slight deviations from the mathematical model brought on by unavoidable measurement errors must have an equivalently slight effect on the solution. In this work, a literature review on ODE stability is carried out and with the same occasion three main questions are answered as follows: when does an ODE become unstable? How does stability affect engineering simulations? Can stability condition be derived for all types of problems?

Theory

An equation that contains some ordinary derivatives of a function (as opposed to partial derivatives) is known as an ordinary differential equation (ODE). Solving an ODE, or figuring out which function or functions fulfil the equation.

How do you locate the function itself if you know what a function's derivative is? You must integrate in order to find the antiderivative. For instance, if you receive

                                                                    `dx/dt (t)=cos⁡(t) (1)`

Then what is function x(t)? Cos(t) antiderivative is sin(t), thus x(t) must also be sin(t). But there's one thing we overlooked: if we only knew the derivative, there was always some arbitrary constant we couldn't figure out. As a result, the above equation only allows us to determine that.

x(t)=sin(t) +C

for some arbitrary constant C. You can verify that indeed x(t) satisfies the equation (1).

Applications of Ordinary Differential Equations

Ordinary differential equations are used to compute the flow or movement of electricity, the motion of an object like a pendulum, and to clarify thermodynamics ideas. Additionally, they are employed in medicine to graphically monitor the progression of diseases.

Stability analysis Methods for ODEs

Generally, different methods are used to carry out the stability analysis of ordinary differential equations. These methods include Euler’s method, Trapezoid method, Explicit and Implicit methods. In order to get a clear overview about the ODE’s stability, Euler’s method will be exploited through examples.

According to Euler's approach, tangents to the curve are drawn along the length of each subinterval when a curve is plotted with the independent variable at the abscissa and the dependent variable at the ordinate. The curve is therefore regarded as linear along the length of these subintervals. The expression that results from the aforesaid procedure is,  `y_(i+1)= y_i+h.f(x_i,y_i)`

   where  `f(x_i,y_i)`  is any function and h is the length of the subinterval.

Considering the procedure described above, this phrase is pretty simple. This is how Euler solves problems. Considering the procedure described above, this phrase is pretty simple. Since Yi+1 only appears on the left side, this is Euler's explicit approach for solving ODEs.

The aforementioned formula, however, is only stable for low values of h. Since we assume that the profile is linear along the length of this subinterval, as the length of the subinterval increases, so does the error, causing the system to become unstable. To address this issue, Euler's approach is slightly adjusted to make the system inherently stable for any value of the length of the subinterval.

The updated Euler's method's expression is,

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

Since y i+1 appears in both the LHS and the RHS, the modified Euler's technique is an implicit scheme. Using a general single step method equation, the stability conditions may be understood.

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

If |E| <1 for absolute stability, the error shoots up to a very high value and produces illogical results. The following example will help us understand the concept of stability.

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

1. Using Euler’s Method (Explicit scheme)

In order to illustrate the stability concept, a Python code is written to solve the above equation.

"""

Explicit Euler demo

"""
import math
import matplotlib.pyplot as plt



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


t_start = 0
t_end = 0.1
n = 500

dt = (t_end - t_start)/(n-1)

t = 0;
y_old = 0;
y_anal = []
y_new = []

for i in range(0, n):
	
	y_new.append(y_old + dt*f(t,y_old))
	y_anal.append(f(t,y_old))
	y_old = y_new[i]
	t = t + dt

	


plt.plot(y_new, color = 'green')
#plt.plot(y_anal, color = 'red')
plt.show()

Code Explanation

In the above code, the “math” module as well as “matplotlib.pyplot” module is imported. The function     is defined and the time step is taken as 0.1. The umber of grids is set to 500. Different variables are defined with initial values which further affect the solution. A “for loop” is then defined to compute the values of “y_new” for “n” number of grid points using the equation in the example. The initial values of “y_old” and “t” are set to zero. Then these values are plotted using “matplotlib. pyplot” library as shown below.

Unstable Condition

In this case, the time step is kept unchanged as previously but the number of grid points is changed to 50. Other parts of the program remaining the same.

"""

Explicit Euler demo

"""
import math
import matplotlib.pyplot as plt



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


t_start = 0
t_end = 0.1
n = 50

dt = (t_end - t_start)/(n-1)

t = 0;
y_old = 0;
y_anal = []
y_new = []

for i in range(0, n):
	
	y_new.append(y_old + dt*f(t,y_old))
	y_anal.append(f(t,y_old))
	y_old = y_new[i]
	t = t + dt

	


plt.plot(y_new, color = 'green')
#plt.plot(y_anal, color = 'red')
plt.show()

The result obtained is as follows:

We can see through the graph that the value of “y_new” starts fluctuating and the solution keeps diverging from the actual solution with increasing frequency instead of converging. Clearly we can say that a small change in input parameters leads to large change in output values which is undesirable for the convergence of the solution referred to as unstable solution.

2. Using Modified Euler’s Method (Implicit Scheme)

The above equation is now solved implicitly using modified Euler’s method. After simplification, the expression is as follows:

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

The code used is displayed below.

import numpy as np
import matplotlib.pyplot as plt

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] + dt*(-1000*y[i-1]-np.exp(-t[i-1]));

	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] - np.exp(-t[i-1])*dt)/(1 + 1000*dt);


	return [t, y]

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] = np.exp(-1000*t[i])*(1 - np.exp(999*t[i]))/999

	return y


total_t = 0.1
dt = 1.0e-4
t, y_explicit = explicit_euler(total_t, dt)
t, y_implicit = implicit_euler(total_t, dt)
y_analytical = analytical(total_t, dt)
plt.plot(t, y_explicit, color = 'blue')
plt.plot(t, y_implicit, color = 'red')
plt.plot(t, y_analytical, color = 'orange')
plt.legend(['Explicit', 'Implicit', 'Analytical'])
plt.show()

Explanation of the program

In this program, three functions are defined: the first function (explicit_euler) is employed to solve the given equation using Euler’s explicit method; the second function (implicit_euler) is employed to solve the same equation using Euler’s implicit method and the last function (analytical) is employed to compute the analytical solution of the given equation. The time span is taken as total_t equal tom 0.1s. The time step dt is taken as 1e-04. The function values are computed using “for loops” under the three functions or methods. After that, these values are plotted in the same graph to give a comparative view of the explicit and implicit methods.

 Output

We can see that, for the specific input parameters values, the explicit solution, implicit solution and the analytical solution are converging.

When does an ODE become unstable?

Information received by solving simultaneous equations for the entire grid for each time step is contained in an implicit solution. Although this requires more computer power, it offers better stability and larger time steps.

The behavior of the solution as the time-step dt is raised is what is meant by numerical stability. The method is referred to as unconditionally stable if the solution continues to behave well for arbitrary big values of the time step. Explicit approaches, which are always conditionally stable and become unstable as dt (step) value increases, never lead to this circumstance. Therefore, when dt grows, ODE becomes unstable for this solution.

How does stability affect engineering simulations?

When the exact differential equation's solution is bounded, stability ensures that the numerical approach yields a bounded solution. A solution cannot be found for a specific simulation if the solution is not bounded because it diverges.

Can Stability condition be derived for all types of problems?

An implicit method is one in which the differential equation is discretized in such a way that there are multiple unknowns at n+1 time level on the LHS of the equation and the terms on the RHS are known ones at n time level. A neutrally stable scheme is one in which errors remain constant as the computations are carried forward. The advantage of implicit schemes is that no stability requirements need to be met. They are permanently steady. However, because you require more storage, they are computationally expensive. So long as the time step size is chosen properly, an implicit scheme can be used to solve any kind of problem. In essence, two stability requirements must be met:

When dealing with linear partial differential equations, von Neumann stability analysis, also referred to as Fourier stability analysis, is a process used to examine the stability of finite difference schemes. The Fourier decomposition of numerical inaccuracy or error serves as the foundation for the study.

If the eigenvalue of the recursive system represented by the matrix is greater than 1, the associated eigenvector will be amplified repeatedly without any limitations, making the system (matrix) unstable. Eigenvalues less or equal to 1 are also stable.