Week 5.2 - Literature review: ODE Stability

ODE STABILITY - 

Stability for the ODEs can be defined as the error in the computation of ODEs getting minimised as more and more iterations and timestep are calculated. On the other hand for an unstable calcuation this error grows woth successive calculations.

METHODS FOR SOLVING ODEs -

There are generally two methods to solve the ordinary differential functions, Explicit method and Implicit method.

Explicit method :- Explicit methods are where the calculation of the future timestep is done using the values of the present timestep. it takes less computation power and calculates fast. but here stability criteria are present. examples - jacobi and gauss-seidel methods

Implicit method :- Implicit method, the calculation for the future timestep values are done using the values at the current and future time step values so this is a stable solution. and here no stability criteria is present . example - Crank nicolson method.

When does an ODE Become unstable ?

For an explicit method, there is a criteria called a the courant number with is given by

`sf"Courant number (C)" prop (Delta t)/(Delta x)`

`sf"courant number (C)" = c . (Delta t)/(Delta x)` 

where c is the proportionality constant

for a given condition if the courant number exceed a specific limit then the soltuion becomes unstable so it is necessary in case of the explicit methods that the courant number is below the limit . So here if the value of `Delta t` is very big or the if `Delta x` is too small then the solution is going to be unstable so both this values of time-step and element size must be in proportion.

For the implicit method the solution of the ODE is always stable so there is no such criteria as courant number for this methods. But as this methods uses a lot of computional it is costly so it cannot be used everywhere.

Let us now consitder an example of the a simple ODE and try to see the difference between the implicit and explicit solution.

 The ODE used for this example is shown below first then we will be deriving the solution to the ODE and compare both the methods.

`sf"ODE =>" dy/dt = -1000.y - e^(-t)`

Explicit Eulers method:

Explicit method the value at the future timstep will be calculated using the value at the present time-step

let 'n' denote present time step and 'n+1' denote the future time step 

using explicit approximation

`(y^(n+1) - y^n)/(Delta t) = - 1000 - e^(-t)`

`y^(n+1) = y^n + Delta t(-1000 y^n -e^(-t))`

so here the value of the future can be easilty calculted as we will be knowing the values of y at each node at the present time step.

Implicit  Eulers method:

Implicit method , we use the backward differencing method where the value of past time step is know at the value of the present timestep is to be calculated using both the past and present time steps.

we will be using 'n-1' to denote the values at the past time step and 'n' the denote the value of the present time steps.

`(y^(n) - y^(n-1))/(Delta t) = -1000 y^n -e^(-t)`

`y^n = y^(n-1) + Deltat (-1000y^n -e^(-t))`

`y^n + 1000.y^n.Delta t = y^(n-1) -Delta t.e^(-t)`

`y^n = ( y^(n-1) - Delta t. e^(-t))/(1+1000 Delta t)`

in this way the value at the current time step will be calculted.

Analytical method:

For the analytical method we will be integrating the ODE for y(0) = 0  and get the result as 

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

using this as the exact analytical solution we will compare both the implicit and explcit method results.

Writing the Python Code for the numerical method and comparing the results:

the python code for the comparison for both the implicit and explicit methods is given below

"""
Comparison for explicit and implicit methods
ODE => dy/dy = -1000y -e^(-t)
"""
import numpy as np 
import matplotlib.pyplot as plt

#function to solve the ODE using explicit method 
def explicit_fun(t_total, dt):
	t = np.arange(0,t_total,dt) # creating time array
	y = np.zeros(len(t)) # creating space array
	# for loop to calculate the value of y at ecah time step
	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] # return the value of t and y as output

# function to solve the ODE using implicit method
def implicit_fun(t_total, dt):
	t = np.arange(0,t_total,dt)# creating time array
	y = np.zeros(len(t)) # creatign space array
	#for loo to calculate the vaule of y for each time step
	for i in range(1,len(t)):
		#y[i] = y[i-1] - 1000*y[n-1]*dt - np.exp(-t[i-1])*dt
		y[i] = (y[i-1] - np.exp(-t[i-1])*dt)/(1+1000*dt) 
	return [t,y] # return the value of t and y as output

# function to calculate the ODE analytically
def analytical_fun(t_total,dt):
	t = np.arange(0,t_total,dt) # creatign time array
	y = np.zeros(len(t)) # creating space array

	# for loop to calculate the value of y at each node
	for i in range(1,len(t)):
		y[i] = np.exp(-1000*t[i])*(1-np.exp(999*t[i]))/999
	return y # returning the value of Y as output 




t_total = 0.1 # total time of calculation	
dt = 1e-4 # time step size

# solving the ODE using  explicit, implicit and analytical function 
t,y_explicit = explicit_fun(t_total,dt)
t,y_implicit = implicit_fun(t_total,dt)
y_analytical = analytical_fun(t_total,dt)

# plotting
plt.plot(t,y_explicit, color='red')
plt.plot(t,y_implicit,color = 'blue')
plt.plot(t,y_analytical,color = 'green')
plt.legend(["explicit", " implicit","analytical"])
plt.show()

 now here in the above code we can see that the time step size is 1e-4 and usign this timestep we get the following plot as the output

here we can see that all the calculation of implicit , explicit and analytical are almost identical . 

Now we will change the timpstep value and see the effect.

For dt = 1.5e-4

python code 

"""
Comparison for explicit and implicit methods
ODE => dy/dy = -1000y -e^(-t)
"""
import numpy as np 
import matplotlib.pyplot as plt

#function to solve the ODE using explicit method 
def explicit_fun(t_total, dt):
	t = np.arange(0,t_total,dt) # creating time array
	y = np.zeros(len(t)) # creating space array
	# for loop to calculate the value of y at ecah time step
	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] # return the value of t and y as output

# function to solve the ODE using implicit method
def implicit_fun(t_total, dt):
	t = np.arange(0,t_total,dt)# creating time array
	y = np.zeros(len(t)) # creatign space array
	#for loo to calculate the vaule of y for each time step
	for i in range(1,len(t)):
		#y[i] = y[i-1] - 1000*y[n-1]*dt - np.exp(-t[i-1])*dt
		y[i] = (y[i-1] - np.exp(-t[i-1])*dt)/(1+1000*dt) 
	return [t,y] # return the value of t and y as output

# function to calculate the ODE analytically
def analytical_fun(t_total,dt):
	t = np.arange(0,t_total,dt) # creatign time array
	y = np.zeros(len(t)) # creating space array

	# for loop to calculate the value of y at each node
	for i in range(1,len(t)):
		y[i] = np.exp(-1000*t[i])*(1-np.exp(999*t[i]))/999
	return y # returning the value of Y as output 




t_total = 0.1 # total time of calculation	
dt = 1.5e-4 # time step size

# solving the ODE using  explicit, implicit and analytical function 
t,y_explicit = explicit_fun(t_total,dt)
t,y_implicit = implicit_fun(t_total,dt)
y_analytical = analytical_fun(t_total,dt)

# plotting
plt.plot(t,y_explicit, color='red')
plt.plot(t,y_implicit,color = 'blue')
plt.plot(t,y_analytical,color = 'green')
plt.legend(["explicit", " implicit","analytical"])
plt.show()

the output of this code with timestep value of 1..5e-4 is 

here we can see that there is slight change in the plot for the explicit method near the curve .

now we will increase the timestep further to see more effect

For dt = 5e-4

python code

"""
Comparison for explicit and implicit methods
ODE => dy/dy = -1000y -e^(-t)
"""
import numpy as np 
import matplotlib.pyplot as plt

#function to solve the ODE using explicit method 
def explicit_fun(t_total, dt):
	t = np.arange(0,t_total,dt) # creating time array
	y = np.zeros(len(t)) # creating space array
	# for loop to calculate the value of y at ecah time step
	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] # return the value of t and y as output

# function to solve the ODE using implicit method
def implicit_fun(t_total, dt):
	t = np.arange(0,t_total,dt)# creating time array
	y = np.zeros(len(t)) # creatign space array
	#for loo to calculate the vaule of y for each time step
	for i in range(1,len(t)):
		#y[i] = y[i-1] - 1000*y[n-1]*dt - np.exp(-t[i-1])*dt
		y[i] = (y[i-1] - np.exp(-t[i-1])*dt)/(1+1000*dt) 
	return [t,y] # return the value of t and y as output

# function to calculate the ODE analytically
def analytical_fun(t_total,dt):
	t = np.arange(0,t_total,dt) # creatign time array
	y = np.zeros(len(t)) # creating space array

	# for loop to calculate the value of y at each node
	for i in range(1,len(t)):
		y[i] = np.exp(-1000*t[i])*(1-np.exp(999*t[i]))/999
	return y # returning the value of Y as output 




t_total = 0.1 # total time of calculation	
dt = 5e-4 # time step size

# solving the ODE using  explicit, implicit and analytical function 
t,y_explicit = explicit_fun(t_total,dt)
t,y_implicit = implicit_fun(t_total,dt)
y_analytical = analytical_fun(t_total,dt)

# plotting
plt.plot(t,y_explicit, color='red')
plt.plot(t,y_implicit,color = 'blue')
plt.plot(t,y_analytical,color = 'green')
plt.legend(["explicit", " implicit","analytical"])
plt.show()

the output for this code with timestep of 5e-4 is shown below

we can see now the explicit method is deflecting from the original results 

lets see the change if we increase the tims step size to 1e-3

For dt = 1e-3

python code

"""
Comparison for explicit and implicit methods
ODE => dy/dy = -1000y -e^(-t)
"""
import numpy as np 
import matplotlib.pyplot as plt

#function to solve the ODE using explicit method 
def explicit_fun(t_total, dt):
	t = np.arange(0,t_total,dt) # creating time array
	y = np.zeros(len(t)) # creating space array
	# for loop to calculate the value of y at ecah time step
	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] # return the value of t and y as output

# function to solve the ODE using implicit method
def implicit_fun(t_total, dt):
	t = np.arange(0,t_total,dt)# creating time array
	y = np.zeros(len(t)) # creatign space array
	#for loo to calculate the vaule of y for each time step
	for i in range(1,len(t)):
		#y[i] = y[i-1] - 1000*y[n-1]*dt - np.exp(-t[i-1])*dt
		y[i] = (y[i-1] - np.exp(-t[i-1])*dt)/(1+1000*dt) 
	return [t,y] # return the value of t and y as output

# function to calculate the ODE analytically
def analytical_fun(t_total,dt):
	t = np.arange(0,t_total,dt) # creatign time array
	y = np.zeros(len(t)) # creating space array

	# for loop to calculate the value of y at each node
	for i in range(1,len(t)):
		y[i] = np.exp(-1000*t[i])*(1-np.exp(999*t[i]))/999
	return y # returning the value of Y as output 




t_total = 0.1 # total time of calculation	
dt = 1e-3 # time step size

# solving the ODE using  explicit, implicit and analytical function 
t,y_explicit = explicit_fun(t_total,dt)
t,y_implicit = implicit_fun(t_total,dt)
y_analytical = analytical_fun(t_total,dt)

# plotting
plt.plot(t,y_explicit, color='red')
plt.plot(t,y_implicit,color = 'blue')
plt.plot(t,y_analytical,color = 'green')
plt.legend(["explicit", " implicit","analytical"])
plt.show()

the resulting lot for this value is shown below

the error has increased a lot and the curve for the explicit solution has sharp edges 

For dt = 1.2e-3

For dt = 1.5 e-3

it is clear from the above 2 images the the solution for the explicit method becomes unstabe as the timstpe value is increased beyond a certain value but for the implicit method the solution is stable for all these values of timestep.

For dt = 2e-3

for this value the solution for the explicit method is totally blown up but the implicit method is still stable .The implicit method is therefore called as unconditionally stable because what ever the value of the timestep we get a stable solution.

How Does stability affect engineering solutons?

For the engineering solution stability is very important. Here the small error can lead to untability that can grow exonentially as we can see in the example for dt = 1.5e-3 the stability was in a very small location but as we increase the timestep to dt = 2e-4 the unstability can be see throught the solution . it has incresed expoentially and if the engineering solution has such instablities it can result in huge loss in the form of expense and also my even have situation where it can dangerous to life.

So it is very important to have stability in the engineering solutions. we can see that the instability in the above equation starts at dt = 1.1e-4 . 

Can stability conditions be derived for all types of problems ?

Stability can be derived for most of all types of problems. the explicit method will require more criteria for the approximation solutions and for the implicit methods as it is more stable it can be derived easily but will require more computation power and will be costlier.