Regular Perturbation

Motion in a Resistive Medium

A body of mass $m$ and initial velocity $V_0$ moves in a straight line in a resistive medium that offers a resistive force of magnitude $av-bv^2$, where $v=v(\tau)$ is the velocity of the object at time $\tau$. Newton's second law will lead us to the following initial value problem:

$$ m\frac{dv}{d \tau}=-av+bv^2, ~~~~~~~~~~~~~~~~~~ v(0)=V_0 \tag{1}$$

$a$ and $b$ are positive constants and we assume that $b<

$u=\frac{v}{V_0}$ and $t=\frac{\tau}{m/a}$.

This rescaling will give us the following initial value problem:

$$u'=-u+\epsilon y^2, ~~~~~~~~~~~~~~~~~~ u(0)=1. \tag{2}$$

Where $\epsilon\equiv \frac{bV_0}{a} <<1$.

The exact solution for eq. (2) is

$$ u_{ex}=\frac{e^{-t}}{1+\epsilon (e^{-t}-1)}. $$

In regular perturbation we attempt to obtain solution of the type:

$$ u=u_0+\epsilon u_1+\epsilon^2 u_2+ \dots, \tag{3}$$

We substitute the perturbation series eq. (3) into the intial value problem eq. (2) and compare terms containig the same power of \epsilon. This gives us the following set of equations:

$\epsilon^0$: $~~~~~~~$ $u_0'=-u_0, ~~~~~~~~~~~~~~~~~~~~ u_0(0)=1$

$\epsilon~~$: $~~~~~~~$ $u_1'=-u_1+{u_0}^2, ~~~~~~~~~~~ u_1(0)=0$

$\epsilon^2$: $~~~~~~~$ $u_2'=-u_2+2{u_0}u_1,\dots ~~ u_2(0)=0$

Solving these equations give us $u_0=e^{-t}$, $u_1=e^{-t}-e^{-2t}$ and $u_2=e^{-t}-2e^{-2t} +e^{-3t}$

Therefore an approximate solution for eq. (2) involving the first three terms from perturbation series will be

$u_{approx}=e^{-t}+\epsilon(e^{-t}-e^{-2t})+\epsilon^2(e^{-t}-2e^{-2t} +e^{-3t})$

In [1]:
# Regular perturbation 
# import all the packages required 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
%matplotlib inline

# Declare a global value for epsilon 
ep=0.075


#Write differential equation in form of a function 
def f1(u,t):
    du=-u+ep*u**2
    return du


plt.figure(figsize=(10,6))

# Times at which the solution is to be computed.
t = np.linspace(0, 1, 5000)

# Initial values 
uinit=np.array([1.0])


#Numerical solution of initial value problem 
usolu = odeint(f1,uinit, t)

plt.plot(t, usolu[:,0],color='grey',ls='-',lw=6,label='Numerical Solution')





#Exact solution 
t = np.linspace(0, 1, 70)
uexac=np.exp(-t)/(1+ep*(np.exp(-t)-1))
plt.plot(t, uexac,'r.',ms=10,label='Exact Solution $u_{ex}$')



# Perturabtion solution 
t = np.linspace(0, 1, 50)
# Perturbation terms 
u0=np.exp(-t)
u1=ep*(np.exp(-t)-np.exp(-2*t))
u2=ep*(np.exp(-t)-np.exp(-2*t))

plt.plot(t, u0,'k--',lw=4,label=' $u_{0}$')
plt.plot(t, u1+u0,'g.',ms=15,label='Approx. Solution $u_{approx}$')



plt.xlim(0.4,1); plt.ylim(0.38,0.7)   
plt.xlabel('time (t)',fontsize=20); 

plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.,fontsize=20)
plt.show;