Week 5 - Literature review: ODE Stability

A numerical method is said to be unstable if small perturbations cause the solution to oscillate and diverge.

Consider the ODE:

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

The analytical solution used for comparing the solutions obtained with the help of wolfram alpha is given below:

y = exp( − 1000 × t) × (1 − exp(999 × t)) / 999

An ODE or PDE system can be solved numerically. Numerical solution can be by two method i.e explicit and implicit schemes. In explicit the values of the previous time step are used for computation and in implicit the values from the current time step are used for computation.

Euler's forward Difference (Explicit):

`y_(n + 1) = y_n + h f (y_n, t_n )`

`y_(n + 1) = y_n + Δt ( − 1000y_n − e^(-t_n) )`

y_{n + 1} = (1- 1000y_n Δt ) y_n − Δt e^-tn

In above equation, (1 − 1000Δt) is a growth factor.

If |1 − 1000Δt| ≤ 1 , then the solution is stable.

|1 − 1000Δt| > 1 , then the solution is unstable.

Euler's Backward Difference(Implicit):

y_n = y_n-1 + h f (y_n, t_n )

y_n = y_n-1 + Δt ( − 1000y_n − e^-tn )

From above (1/ (1000 Δt +1 ) ) the value of the time step can be determined.

From this relation the value of the time step can be determined.

"""
A numerical method is said to be unstable if small perturbations cause the solution to oscillate and diverge.
"""



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

def analytical(t_total,dt):
	t=np.arange(0,t_total,dt) 
	y=np.zeros(len(t)) 

	for i in range(0,len(t)):
		y[i] = math.exp(-1000*t[i])*(1-math.exp(999*t[i]))/999;   # analytical solutions obtained with the help of  wolfram alpha 
	return[t,y]	


def explicit(t_total,dt):
	t=np.arange(0,t_total,dt) 
	y=np.zeros(len(t)) 

	for i in range(0,len(t)):
		y[i]=y[i-1]+dt*(-1000*y[i-1] -np.exp(-t[i-1]));    # yn + 1 = yn + h f (yn, tn )
	return[t,y]	

def implicit(t_total,dt):
	t=np.arange(0,t_total,dt)
	y=np.zeros(len(t))

	for i in range(0,len(t)):
		y[i]=(y[i-1]-np.exp(-t[i])*dt)/(1+1000*dt);      #  yn = yn-1 + h f (yn, tn )
	return[t,y]	

t_total=0.5
dt=0.003
	
t,y_ana=analytical(t_total,dt)
t,y_exp=explicit(t_total,dt)
t,y_imp=implicit(t_total,dt)
plt.plot(t,y_ana)
plt.plot(t,y_exp)
plt.plot(t,y_imp)
plt.legend(('Analytical','Explicit','Implicit'))
plt.grid(True)
title_1=('Timestep:',dt)
plt.title(title_1)
plt.show()

From the above two graphs its is clear that implicit schemes are stable and explicit method becomes unstable at the larger time steps and spatial steps.whereas i.e the explicit schemes are conditionally stable.

From the above images we can note that implicit is stable and nearly accurate to the analytical solution but the explicit method solution depends on the time step dt used, giving unstable results when values of dt is higher. To avoid this, the time step should be reduced.

When does a ODE become unstable?

Stability, in mathematics, condition in which a slight disturbance in a system does not produce too disrupting an effect on that system. In terms of the solution of a differential equation, a function f(x) is said to be stable if any other solution of the equation that starts out sufficiently close to it when x = 0 remains close to it for succeeding values of x.If the difference between the solutions approaches zero as x increases, the solution is called asymptotically stable. If a solution does not have either of these properties, it is called unstable.(* https://www.britannica.com/science/stability-solution-of-equations) .

Stability has to do with the behavior of the solution with an increse in the time-step dt. In explicit methods, which is conditionally stable and is unstable as dt (step) value increases. Hence ODE becomes unstable for this solution as dt increases.

How does stability affect engineering simulations?

When ODE solution becomes unstable results in unphysical solutions as a result error increses in the system, these unphysical solutions will represrnts unaccurate results when doing an engineering simulation.

Can stability condition be derived for all types of problems?

In numerical analysis, von Neumann stability analysis I.e Fourier stability analysis is a procedure used to check the stability of finite difference schemes as applied to linear partial differential equations. The analysis is based on the Fourier decomposition of numerical error. But detection of stability of models with several variables is complex as compared to one-variable models. To check the stability of model , Eigenvalues are required to calculate. Eigenvalues are generally complex numbers. If real parts of all eigenvalues are negative, then the solution is stable. If at least one eigenvalue has a positive real part, then the solution is unstable.