Week 5.2 - Literature review: ODE Stability
import matplotlib.pyplot as plt
import numpy as np
# For explicit method
def explicit(T,dt):
t = np.arange(0,T,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]));
return t,y
# For explicit method using backward differencing
def implicit(T,dt):
t = np.arange(0,T,dt)
y = np.zeros(len(t))
for i in range(0,len(t)):
y[i] = (y[i-1] - np.exp(-t[i-1])*dt)/(1+1000*dt);
return[t,y]
T = 0.1
dt = 0.01
t,y_exp = explicit(T,dt)
t,y_imp = implicit(T,dt)
plt.plot(t,y_exp)
plt.plot(t,y_imp)
plt.legend(('Explicit','implicit'))
plt.grid("on")
plt.show()
Aim : To review the literature about Ordinary Differential Equation and to write the Python program on Stability.
Let's start by defining a simple order differentiate equation.
`dy/dt = e^t`
if i have to solve this particular ODE analyticaly then i knew that the solution
`y = e^t + c`(constant)
this constant can be determined if it is definate interval, how would i solve it numerically so i will basically say that
`((y^{New})- (y^{old)}) /(Deltat)= e^t`
e^t is nothing but original function = f(t)
`y^{New}= f(t) *DeltaT + y^{old}`
= `y^{old} + f(t) *deltat`
this is what we going to doing
once we write the program in Python, we are going to compare the answer with Analytical expression.
Let's take a look at what a backward difference method means and how its different from explicit scheme. let's basically take the ODE
`dy/dt = -1000y- e^(-t)`
first i will going to deal with explicit method
`y^{New} - y^{old}` or we can say that ;
`(y^{n+1} - (y^{n}))/ (Deltat) = -1000*y^n - e^-1`
`y^{n+1} = y^n - (1000y^n*Delta t) - (e^-t*Deltat)` , this is what we called as an explicit function
and as we demonstattion this can give us unstable results.
For, implicit Method, we are using called BDF , which stands for Backward differenecing. Only difference here is we basically write
`y^n - y^{n-1}` here we assume that`y^{n-1}` is our old value. In our explicit case we assume that y^n is the old value.~
we can use n+1 and n and get an answer but having expression like it make lot of sense. right beacuse tell us understand for which term we are using the old value and for which term we are using the current value.
so i can say that , `(y^n - y^(n+1))/(Deltat) = -1000*y^n - e^(-t)`
this is what we called backward differencing because the value of y^n is not known so that is what make its implicit. so this small change makes this scheme stable.
so let's make so rearrangement
`y^n = (-1000y^n - e^-t ) * Deltat + y^(n-1)`
similarly, `y^n - (-1000y^n - e^-t ) *Deltat = y^(n-1)`
`y^n(1 + 1000Deltat) + e^t Deltat = y^(n-1)`
`y^n = (y^(n-1) - e^(-t)Deltat)/( 1+1000Deltat)`
Ordinary Differential Equation
The differential equation containing one or more functions of an independent variable and derivatives of those functions. But in the case ODE, the word ordinary is used for derivative of the functions for the single independent variable.
In mathematics, the term 'Ordinary Differential Equation' is also known as ODE is an equation contains only one independent variable and one or more of its derivatives with respect to the variable. In other words, the ODE is represented as the relation having one independent variable x, the real dependent variable y, with some of its derivatives.
Applications of ODE
1. Modelling the growth of diseases
2. Describes the movement of electricity
3. Describes the motion of the pendulum,waves
4. Used in Newton's second law of motion and law of cooling.
When does a ODE become unstable?
Error in the solution is proportional to 'delta t'. Larger the value of delta t, larger is the number of grid points and larger is the error. when there is more number of species and reactions we cannot calculate analytical expression easily, in that case implicit method is the only way to obtain solution. At higher value of time step the solution starts to blow up, that is becoming unstable for explicit method but stable for implicit method.
How does stability affect engineering simulations?
If the solution becomes unstable the required value will deviate instead of converging. Stability gurantees that numerical solution produces bounded results. If the solution is not bounded the results starts to diverge and we don't get the proper output.
Can stability condition be derived for all types of problems?
Stability of solution can be determined by below 2 ways
1. Eigen value analysis: Real eigne values less than or equal to zero
2. Von-Neumann Stability analysis : This can be applied when domain of the solution is periodic.
PLOTS : Both Implicit and explicit schemes are performed on a function. Here the plots show how they score against each other.
When we increase the number of grid points the simulation is going to slower because we are solving more equations when we are trying to run with large number of grid points the simulation is going to run really slow not because of increase number of grid points but also because of stability criterion and decrease in time step. From the given attached plot it can be observed that implicit method gives stable solution at higher value of 'dt'. The time step plays a crucial role in determining the stability of the solution while using explicit method.
CODING PROCEDURE -
1. Import matplotlib command and calling pylpot as plt
2. Import numpy and call it as np
3. defining a function as a explicit method in which 2 variable are calling T and dt
4. Define a variable T as a np.arange with initial value 0 and final value T with dt timestep.
5. Define a variable y as a np.zeros value with length of t
6. Apply for loop with range of 0 and length of t and defining the function y, then return the values of y and t
7. In case of Implicit same things are going to happend but the key difference is the expression of the y is different .
8. Giving the value of T = 0.1 and dt = 0.02,0.002,0.01
and after this make a plot for different time steps and defining the command of legend, grid and make plots for explicit and implicit function for different timesteps.