Week 5.2 - Literature review: ODE Stability

AIM: Literature review on stability of ordinary differential equation

OBJECTIVE:

• When does a ODE become unstable?
• How does stability affect engineering simulations?
• Can stability condition be derived for all types of problems?

import numpy as np
import matplotlib.pyplot as plt

total_t=0.001
dt=9e-4

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.3
#dt=1.2e-3
t,y_explicit=explicit_euler(total_t,dt)
t,y_implicit=implicit_euler(total_t,dt)
y_analytic=analytical(total_t,dt)
plt.plot(t,y_explicit)
plt.plot(t,y_implicit)
plt.plot(y_analytic)
plt.legend(['explicit','implicit','analytic'])
plt.show()

• In this code ,we have one array of input provided for varying time step.
• For loop is executed to store the 'y' values.
• Then Euler explicit,Euler implcit and analytic values are calculated by providing input as time step and total time.

stability of an ordinary differential equation:

A solution of ODE `y'=f(t,y)`is stable if for every `epsilon`>0,there is`delta`>0 subject to , if `y'(t)` satisfies the ODE and `abs(y'(t_0)-y(t_0))<=delta`, then `abs(y'(t)-y(t))<=epsilon` for all `t>=t_0`.

`let consider `y'=lambday`,`y(t)=y_0*e^(lambda*t)`

if `lambda>0`:exponentia divergence ,every solution is unstable.

if `lambda<0`:every solution is asymptotically stable.

• Eulers method:

`y_(k+1)=y_k+hf(t_k,y_k)`,`dy/dt=lambda*y` with `y(0)=y_0`

we get `y_(k+1)=y_k+hlambday_k`

`y_(k+1)=(1+h*lambda)y_k`

so, Eulers method is stable if `abs(1+hlambda)<=1`

• Forward Euler: `abs(1+hlambda)>1` it is unstable

`abs(1+hlambda)<=1` it is stable.

• Backward Euler: `abs(1/(1-hlambda)<=1` it is stable.

Stability affects engineering simulations

• In many engineering problems we encounter diffusivity in fluid flow problems.So in order to solve these types of problems we have differential equations to define the problems.

• while solving,we approximate and round off the errors to get the solution.But, these errors get acumulated and shows instability to solve in computer.
• One such aspect is numerical diffusion.We need to spread out these round off errors so that the solution will not get blow off.
• In von-Neumann stability condition, is used for stability analysis for finite difference schemes and also for linear partial differential equation.

stability condition for all types of problems can be derived or can be analysed by relative error concept.