Importance of stability in solution for solving ODE's

Objective:- To understand importance of stability in solution for solving ODE's.

Theory:-

Differential equation:

An ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions. And a differential equation is an equation that contains a function with one or more derivatives. In the case ODE, the word ordinary is used for derivative of the functions for the single independent variable. And in case of partial differential equation (PDE) it consists more than one independent variables.

Example of ODE:

`(d^2y)/(dx^2)+(dy)/(dx)+y=c;" where y is function of x"`

Example of PDE:

`(du)/(dx)+(dv)/(dy)+(dw)/(dz)=c`

Euler method:

In mathematics and computational science, the Euler method (also called forward Euler method or an explicit method) is a first-order numerical procedure for solving ordinary differential equations (ODEs) with a given initial value. The Euler method is a first-order method, which means that the local error (error per step) is proportional to the square of the step size, and the global error (error at a given time) is proportional to the step size. The Euler method often serves as the basis to construct more complex methods like predictor–corrector method. The Euler method can also be numerically unstable, especially for stiff equations, meaning that the numerical solution grows very large for equations where the exact solution does not. A simple modification of the Euler method which eliminates the stability problems. This differs from the (standard, or forward) Euler method in that any function is evaluated at the end point of the step, instead of the starting point. The backward Euler method or simply called as an implicit method.

The stability criteria for Euler method can be understand by following simple equation-

`y_("new")=Exxy_("old")`

For absolute stability, `|E|<=1`, failing this condition results in error and unstable solution.

ODE function we are solving here is as follows-

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

By using forward scheme of FDM we get,

`(dy)/dx=(y_("new")-y_("old"))/(Deltax)`

As we know FDM is can be conducted by two type of approaches namely Explicit and implicit.

An explicit solution results from a method that is independent of other values (for the same level), a single equation is used to evaluate new variables for a single time step.

An implicit solution contains information obtained from solving simultaneous equations for previous and current time step. This is computationally more demanding but allows for larger time steps and better stability.

Explicit solution:

`y_(n+1)=y_(n)+(Deltax)xx(-1000y_(n)-e^(-x))`

Implicit solution:

`y_n=(y_(n-1)-(e^(-x)xxDeltax))/(1+1000Deltax)`

As seen in above equation denominator term is very much bigger than numerator term which results in almost stable solution. Hence implicit methods are more complicated to use but usually much more stable than explicit methods.

Analytical solution:

Please click here to get analytical derivative of above function

By using above link calculator we get analytical derivative as-

`y=(e^(-1000 x) (1 - e^(999 x)))/999`

Python program:-

'''
Program to compare results obtained by numerical methods for solving ODE's

ODE is as follows-
dy/dx = -1000*y - exp(-x)

'''

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

# Analytical solution
def analytical_f(x_last,dx):

	x=np.arange(0,x_last,dx)
	y=np.zeros(len(x))

	for i in range(1,len(x)):
		y[i]=(math.exp(-c*x[i])*(1-math.exp((c-1)*x[i])))/(c-1)
	return [x,y]

# Numerical solution using explicit method
def explicit_f(x_last,dx):
	
	x=np.arange(0,x_last,dx)
	y=np.zeros(len(x))
	
	for i in range(1,len(x)):
		y[i]=y[i-1]+dx*(-c*y[i-1]-math.exp(-x[i-1]))
	return [x,y]

# Numerical solution using implicit method
def implicit_f(x_last,dx):

	x=np.arange(0,x_last,dx)
	y=np.zeros(len(x))

	for i in range(1,len(x)):
		y[i]=(y[i-1]-math.exp(-x[i-1])*dx)/(1+c*dx)
	return [x,y]

# Boundary conditions 
c=1000

x_0 = 0 #First node
x_last = 0.1 #End node

n = 1000 #Number of nodes

dx = (x_last - x_0)/(n-1)

# Solving by using above functions
x,y_analytical=analytical_f(x_last,dx)
x,y_explicit=explicit_f(x_last,dx)
x,y_implicit=implicit_f(x_last,dx)

# Plotting
plt.grid('on')
plt.plot(y_analytical,color='red')
plt.plot(y_explicit,'-.',color='blue')
plt.plot(y_implicit,'--',color='green')
plt.xlabel('Number of nodes')
plt.ylabel('y magnitude')
plt.legend(['Analytical solution','Explicit solution','Implicit solution'])
plt.title(' Plot for values of dx= '+str(dx)+' & n = '+str(n))
plt.show()

Output plots:-

As observed from all of the above figures its clearly visible as number of nodes changes the value for dx is also changing which affects stability of simulation. The increase in dx value increases un-stability of solution. The solution obtained by implicit schemes are more stable and accurate and precise than solution obtained by using explicit scheme. Though the implicit solution and analytical solution looks same line but they do have minor difference in solution value.

In the numerical simulations (computational analysis), numerical stability is a generally desirable property. The precise definition of stability depends on the context. One is numerical linear algebra and the other is algorithms for solving ordinary and partial differential equations by discrete approximation. In numerical linear algebra the principal concern is instabilities caused by proximity to singularities of various kinds, such as very small or nearly colliding eigenvalues. On the other hand, in numerical algorithms for differential equations the concern is the growth of round-off errors and/or small fluctuations in initial data which might cause a large deviation of final answer from the exact solution.

Some numerical algorithms may damp out the small fluctuations (errors) in the input data; others might magnify such errors. Calculations that can be proven not to magnify approximation errors are called numerically stable. One of the common tasks of numerical analysis is to try to select algorithms which are robust – that is to say, do not produce a wildly different result for very small change in the input data.

Engineering problems are real life problems and stability plays important role in solution, so if solution for such problem is unstable it will lead in failure of product. So it is undesirable and hazardous to have unstable solution in engineering problems.

Stability in engineering simulations depends on some factors like upwind schemes, solving methods (implicit or explicit), time step size, mesh grid size, relaxation factor, etc.

For example, while using coarse mesh with fist order upwind scheme, we get converged solution fast (with few number of iterations) but it would not be accurate and stable. So for coarse mesh choosing second order scheme is a good choice. Second order scheme use a greater number of neighboring cells to computes flow variables at a point, therefore it gives accurate results. Using under relaxation factor also helps in giving solution greater stability, though it takes more time in convergence as compared to relaxed or over relaxed system. Under relaxation factor prevents oscillation in our results. Another condition is CFL condition or ‘Courant number’, it also affects simulation results. This term gives unstable results if size of our time step is larger than the distance between two neighboring cells in mesh grid. So using correct time step size is very important.

It is possible to derive stability criteria using Von Neumann stability and/or linear stability analysis. The stability conditions (criteria) depend upon ODE, boundary conditions and solver type.