Week 5 - Literature review: ODE Stability

STABILITY :
A solution is going to be a little different when solved using numerical approach than that using analytical approach. The difference between the two can be called as error .A solution is called stable if the error in it decreases with time,i.e. it should diminish and the answer should reach more close to the actual answer. If that does not happen , the solution becomes unstable.Explicit way of solving the ordinary differential equation has a limit of giving stable solution whereas the implicit takes the stablity of a solution as an in buit part of it, i.e an implicit solution is always going to be stable although it can take some complicated number of steps, but it will be stable whereas explicit solution needs to have some constraints.

In explicit solution, the value of dt i.e. increment in time step can only be increased upto a certain value.

Implicit solution will give stable answer everytime but.( is not so that implicit solution will give accurate value for any value of dt, the value of dt can be big so that solution is achieved in less steps but it is also less accurate when dt increases because truncation error increases.)

If the value of dt is increased very much i.e lesser no. of grid points, the answer will become inaccuarte, it may be still stable but in accurate irrespective of the type of solution used.

Example- dy/dt= -200y-4e^(-t)

Explicit solution:

Backward differencing, (y(n)-y(n-1))/dt=-200y(n-1)-4e^t

y(n)=y(n-1)-(200y(n-1)+4e^t)dt

Implicit Solution:

Backward differencing, y(n)-y(n-1)/dt=-200y(n)-4e^t

y(n)(1+200dt)=y(n-1)-(4e^t)*dt

Analytical solution: (-4/201)(e^(-200*t)(e^(201*t)-1)

Keeping the time=0.1 sec and the dt=1e-4

Performing the code:

import matplotlib.pyplot as plt 
import numpy as np
A=200
def explicit_euler(t_total,dt):
	t=np.arange(0,t_total,dt)
	y=np.zeros(len(t))
	for i in range(1,len(t)):
		y[i]=y[i-1]-dt*(A*y[i-1]+4*np.exp(t[i-1]))
	return[t,y]
def implicit_euler(t_total,dt):
	t=np.arange(0,t_total,dt)
	y=np.zeros(len(t))

	for i in range(1,len(t)):
		y[i]=(y[i-1]-(dt*4*np.exp(t[i-1])))/(1+A*dt)
	return[t,y]

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

	for i in range(1,len(t)):
		y[i]=(-4/201)*(np.exp(-200*t[i])*(np.exp(201*t[i])-1))
		
	return[t,y]


dt=5e-4
t_total=0.1

t,y_explicit=explicit_euler(t_total,dt)
t,y_implicit=implicit_euler(t_total,dt)
t,y_analytical=analytical_euler(t_total,dt)
plt.plot(t,y_analytical,color='blue')
plt.plot(t,y_implicit,color='black')
plt.plot(t,y_explicit,color='red')

plt.show()

Observation:

when dt=5e-4, i.e. very less, the explicit, implicit and the analytical solution are same and coincide.

red-explicit

blue-analytical

black-implicit

If i increase the value of time increment, keeping the equation same, then we can see that implicit solution remains almost close to the analytical answer, but the explicit solution will become unstable after some time.

The code:

import matplotlib.pyplot as plt 
import numpy as np

A=200
for dt in np.arange(5e-4,60e-4,0.001):
	def explicit_euler(t_total,dt):
		t=np.arange(0,t_total,dt)
		y=np.zeros(len(t))
		for i in range(1,len(t)):
			y[i]=y[i-1]-dt*(A*y[i-1]+4*np.exp(t[i-1]))
		return[t,y]
		
	def implicit_euler(t_total,dt):
		t=np.arange(0,t_total,dt)
		y=np.zeros(len(t))

		for i in range(1,len(t)):
			y[i]=(y[i-1]-(dt*4*np.exp(t[i-1])))/(1+A*dt)
		return[t,y]
	
	def analytical_euler(t_total,dt):
		t=np.arange(0,t_total,dt)
		y=np.zeros(len(t))

		for i in range(1,len(t)):
			y[i]=(-4/201)*(np.exp(-200*t[i])*(np.exp(201*t[i])-1))
			
		return[t,y]

	
	
	
	t_total=0.1
	t,y_explicit=explicit_euler(t_total,dt)
	t,y_implicit=implicit_euler(t_total,dt)
	t,y_analytical=analytical_euler(t_total,dt)
	plt.plot(t,y_analytical,color='blue')
	plt.plot(t,y_implicit,color='black')
	plt.plot(t,y_explicit,color='red')
	plt.xlabel('time')
	plt.ylabel('value of y')
plt.show()

dt is varied from 5e-4 to 60e-4, and the curve looks as:

Blue line denotes the analytical answer, the red line denotes the explicit solution which goes on becoming unstable after some time steps and the black line denotes the implicit solution, which remains close to the analytical answer.

Observation: We can see that with increase in dt, the explicit solution becomes unstable and after a certain amount of time , it will become totally unstable.Thus, the ODE stability depends on the time increment and the no.of node points , when it is performed explicitly. Whereas here , implicit solution wil always give stable answer.

Keeping time step same but changing the value of coffiecient of y in equation

Also keeping the time step same and changing the value of cofficient of y( given as A in the code) in the equation will change the no. of time steps it take to reach the solution. As the value of A increases, the value of dt should be decreased in order to keep the solution stable, mostly such as value of dt x value of constant (given as A in the code)= same always for all variations.

The code:

import matplotlib.pyplot as plt 
import numpy as np


for A in range(200,800,100):
	def explicit_euler(t_total,dt):
		t=np.arange(0,t_total,dt)
		y=np.zeros(len(t))
		for i in range(1,len(t)):
			y[i]=y[i-1]-dt*(A*y[i-1]+4*np.exp(t[i-1]))
		return[t,y]
		
	def implicit_euler(t_total,dt):
		t=np.arange(0,t_total,dt)
		y=np.zeros(len(t))

		for i in range(1,len(t)):
			y[i]=(y[i-1]-(dt*4*np.exp(t[i-1])))/(1+A*dt)
		return[t,y]
	'''
	def analytical_euler(t_total,dt):
		t=np.arange(0,t_total,dt)
		y=np.zeros(len(t))

		for i in range(1,len(t)):
			y[i]=(-4/201)*(np.exp(-200*t[i])*(np.exp(201*t[i])-1))
			
		return[t,y]

	'''
	
	t_total=0.1
	dt=20e-4
	
	t,y_explicit=explicit_euler(t_total,dt)
	t,y_implicit=implicit_euler(t_total,dt)
	#t,y_analytical=analytical_euler(t_total,dt)
	#plt.plot(t,y_analytical,color='blue')
	plt.plot(t,y_implicit,color='black')
	plt.plot(t,y_explicit,color='red')
	plt.xlabel('time')
	plt.ylabel('value of y')
plt.show()

The graph below shows that by changing the value of A, and keeping the time same as any value , eg here 20e-4 , we can see that the same time step value can be regarded as a unfit value when A value gets increased, with increase in A , the implicit curve still is stable but the explicit cutve becomes unstable after some value of A for the same time step.

So,anODE becomes unstable depending on its time step increment value according to that equation. We need to check this stability only for explicit solutions.

The stability criteria can not always be easily checked and hence we tend to perform implicit solution more when we are not sure that explicit will give the stable answer.

Stability affect engineering problems in crucial ways, it can totally change the answer( as shown by the curve above that how red line of explicit is giving totally destorted results after a certain time) and destroy the results of the solution giving totally wrong answers which are of no use, that is why it is more safe to use implicit solutions although it takes a little more time to give output.