Solving Ordinary Non Linear Differential Equations

Introduction

Let $[a,b]$ be in $\mathbb{R}$ and $y_0$ be a real value. Then we define function $f$: $$\begin{array}{c}f:[a,b]\times \mathbb{R}\to \mathbb{R}\\ (t,y)\mapsto f(t,y)\end{array}$$ We wish to solve the following Cauchy problem, for all: $$\begin{array}{rl}& \{\begin{array}{c}{y}^{\prime }\left(t\right)=f(t,y)t\in [a,b]\\ y\left(a\right)=y_0\end{array}\end{array}$$ We have to find a continuous function $y$, derivable on $[a,b]$ such that it yields: $$forallt\in [a,b],y^{\prime }\left(t\right)=f(t,y\left(t\right))$$ with $y\left(a\right)=y_0.$

Function $y\left(t\right)={\displaystyle \frac{t^2}{4}}$ is a solution of the following problem: $$\{\begin{array}{l}{y}^{\prime }\left(t\right)=\sqrt{y(t)}\text{ pour }t\in [0,4]\\ y\left(0\right)=0\end{array}$$ Nevertheless, we can show that for any $\alpha >0$ functions: $$\{\begin{array}{l}y\left(t\right)=0\text{ for }t\in [0,\alpha ]\\ y\left(t\right)={\displaystyle \frac{{(t-\alpha )}^2}{4}}\text{ si }t>\alpha \end{array}$$ are also solutions to this Cauchy problem... The solution is not unique.

Euler Algorithm (1768)

In general, we can't find the exact solution $y$ of the Cauchy problem $\mathcal{C}$. It's why we try to find an approximation $y_1$ of $y(t_0+h)=y\left(t_1\right)$ where $h>0$ is the step of the numerical method. We set $t_i=t_{i-1}+h$ for $i=1,\dots ,n$ with $t_0=a$ and $t_n=b.$ The Euler algorithm try to approximate the exact solution of $y(t_0+h)$ at $t_0+h$ by using the tangent value of $y\left(t\right)$ at $(t_0,y_0).$ The tangent equation can be written by: $$y-y_0=y^{\prime }\left(t_0\right)(x-t_0)$$ Hence, for $x=t_0+h$ we have: $$y_1=y_0+f(t_0,y_0).h$$ We then solve this new Cauchy problem: $$\left({\mathcal{C}}_1\right)\{\begin{array}{l}{y}^{\prime }\left(t\right)=f(t,y)\\ y\left(t_1\right)=y_1\end{array}$$ which allows us to approximate $y(t_2)$ by $y_2$ at $y(t_1+h):$ $$y_2=y_1+h.f(t_1,y_1).$$ Iteratively, we obtain the $n$ values $y_i$ which approximate function $y$ at each $t_i$.

Exercise [Non Linear Equation of order $1$]
Let us given the following Cauchy problem: $$\{\begin{array}{lll}{y}^{\prime }(t)& =& 1+(y(t){)}^2,t\in [0,1],\\ y(0)& =& 0.\end{array}$$

• Solve this Cauchy problem with the scalar Euler algorithm. Draw the solution curve with Python. For this you have to define the following function:

import numpy as np
import matplotlib.pyplot as plt

""" ODE function """
def f(t,y):
    return 1+y*y

""" Euler Algorithm"""
def Euler(y0,a,b,nb):
    h=(b-a)/nb
    t=a
    X=[a]
    Y=[y0]
    y1=y0
    for i in range(1,nb):
        #TO DO
        Y.append(y2)
        X.append(t)
        y1=y2
    return X,Y

y0= 0.0
a = 0.0
b = 1.0
niter = 150
x,y = Euler(y0,a,b,niter)
plt.plot(x,y,'r', label=r"Euler method", linestyle='--')
plt.legend(loc='upper left', bbox_to_anchor=(1.1, 0.95),fancybox=True, shadow=True)
plt.show()

• Solve this Cauchy problem with the scalar Runge-Kutta of order $2$ algorithm (see Runge Kutta Algorithms ).

• Solve this Cauchy problem with the scalar Runge-Kutta of order $4$ algorithm.

• Solve this Cauchy problem with the second order Adams-Moulton algorithm.

Exercise [Ricatti (1712) Non Linear Equation of order $1$]
Let us given the following Cauchy problem: $$\{\begin{array}{lll}{y}^{\prime }(t)& =& t^2+(y(t){)}^2,t\in [0,{\displaystyle \frac{1}{2}}],\\ y(0)& =& 0.\end{array}$$

• Solve this Cauchy problem with the scalar Euler and Runge-Kutta algorithm.

• Draw the solution curves with Python.

• Compare the approximted values of $y_{nb}$ with the reference solution
  $y(\frac{1}{2})=$ 0.041 791 146 154 681 863 220 76

Exercise [Numerical instability]
Let us given the following Cauchy problem: $$\{\begin{array}{lll}{y}^{\prime }(t)& =& 3y(t)-4e^{-t},t\in [0,10],\\ y(0)& =& 1.\end{array}$$

1. Use the fourth-order Runge-Kutta method with step $h=0.1$ to solve this problem.
2. Compare the result with the analytical solution $y(t)=e^{-t}$.

Let us give $$(D)\{\begin{array}{l}{y}^{\prime \prime }=f(t,y,y^{\prime })\\ y\left(t_0\right)=y_0\\ y^{\prime }\left(t_0\right)=y_0^{\prime }\end{array}$$ where function $f$ is of class $C^1$. By setting $$V\left(t\right)=\left(\begin{array}{c}y\left(t\right)\\ y^{\prime }\left(t\right)\end{array}\right)$$ we obtain the following first order Cauchy problem: $$\{\begin{array}{l}{\displaystyle \frac{dV\left(t\right)}{dt}}=\left(\begin{array}{c}{y}^{\prime }\\ y^{\prime \prime }\end{array}\right)=\left(\begin{array}{c}{y}^{\prime }\\ f(t,y,y^{\prime })\end{array}\right)=F(t,V)\\ V\left(t_0\right)=\left(\begin{array}{c}{y}_0\\ y_0^{\prime }\end{array}\right)\end{array}$$ Thus, we can solve this problem by using the Euler algorithm where the first iteration is given by: $$V_1=V_0+h.F(t_0,V_0)$$

Exercise [Non Linear Equation of order $2$]
We want to solve this second order non linear differential equation: $$\{\begin{array}{l}m{l}^2{\theta }^{\prime \prime }\left(t\right)=-mgl\sin \left(\theta \left(t\right)\right)+c\left(t\right)\text{ pour }t\in [0,T],\\ \theta \left(0\right)=0\\ {\theta }^{\prime }\left(0\right)=1.\end{array}$$ with torque $c\left(t\right)=A\sin \left(t\right)$, $m=0.1$kg, $l=30$cm, $T=5$s and $A=0.1$N.m.

1. Solve this Cauchy problem with the Euler vector algorithm.

import numpy as np

""" Euler Vector Algorithm"""

def F(t,V):
    L=[V[1], # TO DO]
    return np.array(L)

def NEuler(y0,a,b,nb):
    h=(b-a)/nb
    t=a
    X=[a]
    Y=[V0[0]]
    V1=V0
    for i in range(1,nb):
        #TO DO
        t=t+h
        Y.append(V2[0])
        X.append(t)
        V1=V2
    return X,Y

2. Solve this Cauchy problem with the Runge-Kutta of order $2$ algorithm.

3. Solve this Cauchy problem with the Runge-Kutta of order $4$ algorithm.

4. Solve this Cauchy problem with the second order Adams-Moulton algorithm.

Exercise [Mass Spring]

Consider the mass–spring system where dry friction is present between the block and the horizontal surface. The frictional force has a constant magnitude $\mu .m.g$ ($\mu$ is the coefficient of friction) and always opposes the motion. We denote by $y$ the measure from the position where the spring is unstretched.

1. Write the differential equation of the motion of the block
2. If the block is released from rest at $y=y_0$, verify by numerical integration that the next positive peak value of $y$ is $y_0-4\mu .m.\frac{g}{k}$. This relationship can be derived analytically.

Exercise [Iron block]

The magnetized iron block of mass $m$ is attached to a spring of stiffness $k$ and free length $L$. The block is at rest at $x=L$. When the electromagnet is turned on, it exerts a repulsive force $F=c/x^2$ on the block.

1. Write the differential equation of the motion of the mass
2. Determine the amplitude and the period of the motion by numerical integration

Exercise [Magnus Effect applied to a ball]
On June 3, 1997, Roberto Carlos scored against France football team a free kick goal with an improbable effect. We want to simulate in this exercice what happened that day. For this, we set:

• $\rho =1.2Kg.m^{-3}$ : the air density
• $R=0.11m$ the radius of the ball
• $m=0.430Kg$ : the mass of the ball
• $M(t)=(x(t),y(t),z(t))$: the location of the ball at $t\ge 0$
• $v_0=36m.s^{-1}$ : the initial speed of the ball
• $v(t)$ : the speed of the ball at $t\ge 0$.
• $\overrightarrow{v}(t)$ : the speed vector of the ball at $t\ge 0$.
• $\overrightarrow{\omega }=(0,0,{\omega }_0)$ : the angular velocity of the ball which is assumed to be constant during the trajectory
• $\alpha$ : the angle between $\overrightarrow{\omega }$ and $\overrightarrow{v}$
• $\beta$, $\gamma$ : are the initial angles of the force applied to the ball

• its weight
• a friction force: $\overrightarrow{f}$, $\Vert \overrightarrow{f}\Vert =\frac{1}{2}{C}_f\pi R^2\rho v^2,$ where $C_f\approx 0.45$
• a lift force $\overrightarrow{F}=\frac{1}{2}\pi \rho R^3\sin (\alpha )\overrightarrow{\omega }\wedge \overrightarrow{v}$.

1. Define the differential equation apply to the ball.

2. Define the initial conditions.

3. Solve this Cauchy problem with the Runge-Kutta of order $4$ algorithm.

4. Solve this Cauchy problem with the second order Adams-Moulton algorithm.

Runge-Kutta Algorithms

We want to solve the following Cauchy problem: $$\{\begin{array}{l}{y}^{\prime }\left(t\right)=f(t,y)\\ y\left(t_0\right)=y_0\end{array}\text{ for }t\in [t_0,t_n]$$ where $f$ is of class $C^1$. For this, we set $$\begin{array}{rl}{k}_1& =h.f(t_0,y_0)\\ k_2& =h.f(t_0+ah,y_0+\alpha k_1)\end{array}$$ and $$y_1=y_0+R_1{k}_1+R_2{k}_2$$ Consequently, we have to calculate the unknows $a,\alpha ,R_1,R_2$ such that we minimize the following error: $$e_1=y\left(t_1\right)-y_1$$ To do so, we identify the Taylor expansion of $y\left(t_1\right)$ and $y_1.$ We assume that $f$ is of class $C^2$. As $$\begin{array}{rl}\frac{dy}{dt}& =f(t,y)\\ \frac{d^{2}y}{d{t}^2}& =\frac{\partial f(t,y)}{\partial t}+\frac{\partial f(t,y)}{\partial y}f(t,y)\end{array}$$ it yields, $$\begin{array}{rl}y\left(t_1\right)& =y(t_0+h)\\ & \\ y\left(t_1\right)& =y_0+h.\frac{dy\left(t_0\right)}{dt}+\frac{h^2}{2}\frac{d^{2}y\left(t_0\right)}{d{t}^2}+O\left(h^3\right)\\ & \\ y\left(t_1\right)& =y_0+h.f(t_0,y_0)+\frac{h^2}{2}[\frac{\partial f(t_0,y_0)}{\partial t}+\frac{\partial f(t_0,y_0)}{\partial y}f(t_0,y_0)]+O\left(h^3\right)\end{array}$$ As $\phi (t,y,h)=R_1{k}_1+R_2{k}_2$ is of class $C^2$, we can apply the Taylor expansion at $h=0$. Then we have, $$\begin{array}{rl}\phi (t,y,h)& =\phi (t,y,0)+h\frac{\partial \phi (t,y,0)}{\partial h}+\frac{h^2}{2}\frac{{\partial }^2\phi (t,y,0)}{\partial h^2}+O\left(h^3\right)\\ \phi (t,y,h)& =0+h[R_1.f(t_0,y_0)]+h.[R_2.f(t_0,y_0)+R_{2}ha\frac{\partial f(t_0,y_0)}{\partial t}]\\ & +h\left[R_{2}h\alpha \frac{\partial f(t_0,y_0)}{\partial y}f(t_0,y_0)\right]+O\left(h^3\right)\end{array}$$ It comes $$\begin{gathered} y_1=y_0+(R_1+R_2)hf(t_0,y_0) \\ +h^2[a{R}_2\frac{\partial f(t_0,y_0)}{\partial t}+\alpha R_2\frac{\partial f(t_0,y_0)}{\partial y}f(t_0,y_0)]+O\left(h^3\right). \end{gathered}$$ Then, the difference $$\begin{array}{rl}y\left(t_1\right)-y_1& =(1-R_1-R_2)hf(t_0,y_0)+h^2\left[(\frac{1}{2}-a{R}_2)\frac{\partial f(t_0,y_0)}{\partial t}\right]\\ & +h^2\left[(\frac{1}{2}-\alpha R_2)\frac{\partial f(t_0,y_0)}{\partial y}f(t_0,y_0)\right]+O\left(h^3\right)\end{array}$$ is of order $O\left(h^3\right)$ if $R_1+R_2=1$ and $a{R}_2=\alpha R_2={\displaystyle \frac{1}{2}}.$ It follows that:

• if $R_2=0$ then we obtain the previous Euler method
• if $R_1=0,R_2=1$ and $a=\alpha ={\displaystyle \frac{1}{2}}$ we obtain the modified Euler method where $y_1=y_0+h.f(t_0+\frac{h}{2},y_0+\frac{h}{2}f(t_0,y_0))$
• if $R_1=R_2={\displaystyle \frac{1}{2}}$ and $a=\alpha =1$ then we obtain the Euler-Cauchy method where $y_1=y_0+{\displaystyle \frac{1}{2}}hf(t_0,y_0)+{\displaystyle \frac{1}{2}}h.f(t_0+h,y_0+hf(t_0,y_0))$.

The fourth-order Runge–Kutta method is obtained from the Taylor series along the samelines as the second-ordermethod. Since the derivation is rather long and not very instructive, we skip it. The final form of the integration formula again depends on the choice of the parameters; that is, there is no unique Runge–Kutta fourth-order formula. The most popular version, which is known simply as the Runge–Kutta method, entails the following sequence of operations: $$\begin{array}{rl}{k}_1& =h.f(t_0,y_0)\\ k_2& =h.f(t_0+\frac{h}{2},y_0+\frac{k_1}{2})\\ k_3& =h.f(t_0+\frac{h}{2},y_0+\frac{k_2}{2})\\ k_4& =h.f(t_0+h,y_0+k_3)\\ y_1& =y_0+\frac{1}{6}(k_1+2k_2+2k_3+k_4)\end{array}$$ The main drawback of this method is that it does not lend itself to an estimate of the truncation error. Therefore, we must guess the integration step size $h$, or determine it by trial and error. In contrast, the so-called adaptive methods can evaluate the truncation error in each integration step and adjust the value of $h$ accordingly (but at a higher cost of computation). The adaptive method is introduced in the next Section.

Multi-step Algorithms

For all $n$, we set $f(t_n,y_n)=f_n.$ We call $k$-step method an algorithm of the kind: $$\{\begin{array}{c}{y}_n=\sum _{i=1}^k{\alpha }_i{y}_{n-i}+h\sum _{i=0}^k{\beta }_i{f}_{n-i}\text{ for }n\ge k\\ y_{0,}\dots y_{k-1}\text{ are given}\end{array}$$ where the ${\alpha }_i$ and ${\beta }_i$ are real values which are not dependant of $h$.

• If ${\beta }_0=0$ then we have an explicit method.
• If ${\beta }_0\ne 0$ then we have an implicit method.

If we integrate the differential equation we have: $$y\left(t_n\right)=y\left(t_{n-1}\right)+{\int }_{t_{n-1}}^{t_n}f(t,y\left(t\right))dt$$ In order to calculte this integral, we substitute $f(t,y\left(t\right))$ by an interpolation polynomial $p\left(t\right)$ of degree $k$ which satisfies $p\left(t_j\right)=f(t_j,y_j)=f_j$ for $j=n-k,\dots ,n$. For this, we need to take into account the $k$ values $y_{n-k},\dots ,y_{n-1}$ that were calculated previously and $y_n$ which is known. If we write this polynomial in the Newton'polynomial basis we obtain: $$y_n=y_{n-1}+\sum _{j=0}^k\left(\prod _{i=0}^{j-1}(t-t_{n-i})f[t_n,\dots ,t_{n-j}]\right).$$ If the abscissa $t_j$ are uniform then it yields: $$y_n=y_{n-1}+h\underset{i=0}{\sum ^k}{\beta }_i{f}_{n-i}\text{.}$$ Coefficients ${\beta }_i$ have to be defined such that the formula are exacts for maximum degree polynomials. If $k=2$, $$\begin{array}{rl}\text{for }y(t)& =1\text{ then }{y}^{\prime }(t)=0\\ \text{and }y\left(t_n\right)& =1=y\left(t_{n-1}\right)+h.(f(t_n,y_n).{\beta }_0+f(t_{n-1},y_{n-1}).{\beta }_1+f(t_{n-2},y_{n-2}).{\beta }_2)\\ \text{It yields }1& =1+h.[0.{\beta }_0+0.{\beta }_1+0.{\beta }_2]\end{array}$$ $$\begin{array}{rl}\text{for }y(t)& =t-t_n\text{ then }{y}^{\prime }(t)=1\\ \text{and }y\left(t_n\right)& =0=-h+h.({\beta }_0+{\beta }_1+{\beta }_2)\end{array}$$ $$\begin{array}{rl}\text{for }y(t)& ={(t-t_n)}^2\text{ then }{y}^{\prime }(t)=2(t-t_n)\\ \text{and }y\left(t_n\right)& =0=h^2+h.(0.{\beta }_0-2h.{\beta }_1-4h.{\beta }_2)\end{array}$$ From $$\begin{array}{rl}{\beta }_0+{\beta }_1+{\beta }_2& =1\\ 2.{\beta }_1-4.{\beta }_2& =1\end{array}$$ we obtain that $$\begin{array}{rl}{\beta }_0& ={\beta }_1=\frac{1}{2}\\ {\beta }_2& =0\end{array}$$ Consequently, the Adams-Moulton of order $2$ can be written: $$y_n=y_{n-1}+\frac{h}{2}(f_n+f_{n-1})$$ Order 3: $$y_n=y_{n-1}+\frac{h}{12}(5f_n+8f_{n-1}-f_{n-2})$$

• To calculate $f_n$ we need to know the value $y_n$ then we have implicit methods. To solve then, we can approximate $y_n$ by $y_n^{\ast }$ by using an explicit method and insert this value in the implicit schema (correction step).
• Usually the implicit methods are more robust than the explicit ones.

Adams-Bashforth Algorithm
Order 2 : $$y_n=y_{n-1}+\frac{h}{2}(3f_{n-1}-f_{n-2})$$ Order 3 : $$y_n=y_{n-1}+\frac{h}{12}(23f_{n-1}-16f_{n-2}+5f_{n-3})$$