11  Differential equations

Exercise 11.1 - Euler’s methods

The simplest idea to approximate an ordinary differential equation is to discretize time with a step \(h\) and approximate the time derivative on each interval of length \(h\). There are two simple ways to approximate the time derivative. The first is the forward finite difference approximation:

\[ \dot{\boldsymbol{x}}(t) \approx \frac{\boldsymbol{x}(t+h)-\boldsymbol{x}(t)}{h} \,, \]

the second, the backward finite difference:

\[ \dot{\boldsymbol{x}}(t) \approx \frac{\boldsymbol{x}(t)-\boldsymbol{x}(t-h)}{h} \,. \]

The unknowns being the evaluations of the solution \(\boldsymbol{x}\) at times \(t_i = ih\) for \(i\geq0\), i.e., \(\boldsymbol{x}_i = \boldsymbol{x}(t_i)\). The differential equation can thus be approximated using forward finite differences by:

\[ \frac{\boldsymbol{x}_{i+1}-\boldsymbol{x}_i}{t_{i+1}-t_i} = f(t_i, \boldsymbol{x}_i) \,, \]

which gives the explicit Euler formula:

\[ \boldsymbol{x}_{i+1} = \boldsymbol{x}_i + (t_{i+1}-t_i) f(t_i, \boldsymbol{x}_i) \,. \]

With the backward finite difference approximation, we obtain the implicit Euler method (also called backward Euler method):

\[ \boldsymbol{x}_{i} = \boldsymbol{x}_{i-1} + (t_i - t_{i-1}) f(t_i, \boldsymbol{x}_i) \,. \]

On the one hand, the explicit Euler formula allows to compute directly all \(\boldsymbol{x}_i\) by recurrence knowing \(\boldsymbol{x}_0\). On the other hand, the implicit Euler formula requires at each time step the resolution of a nonlinear equation for \(\boldsymbol{x}_i\), for example, with the Newton’s method.

a. Write a function euler_explicit(f,x0,t) that given an initial data x0 returns the values \(\boldsymbol{x}_0, \boldsymbol{x}_1, \dots, \boldsymbol{x}_m\) computed with the explicit Euler method at times \((t_i)_ {i=0}^m\) represented by the vector t.

b. Use the explicit Euler method to solve the differential equation:

\[ \dot{x}(t) + x(t) = \sin(t) \,, \quad x(0)=1 \,, \]

for \(t\in[0,10]\). Compare the results with the exact solution:

\[ x(t) = \frac{1}{2}\big( \sin(t) - \cos(t) + 3e^{-t} \big) \,, \]

for different time discretizations.

c. Solve the previous problem with the implicit Euler method.

d. !! Define a function euler_implicit(f, Dxf, x0, t) implementing the implicit Euler method for nonlinear equations. Note that to solve the nonlinear problem with Newton’s method, the derivative of \(f\) according to \(\boldsymbol{x}\) is required.

e. ! Use the previous methods to find an approximate solution of the system:

\[ \begin{align*} \dot{x}(t)+\cos(y(t)) & =\sin(t)\,, & x(0) & =1\,,\\ \dot{y}(t)+\cos(x(t)) & =0\,, & y(0) & =0\,. \end{align*} \]

Exercise 11.2 - Runge-Kutta methods

The purpose of this exercise is to introduce a class of methods more accurate than Euler’s methods for solving ordinary differential equations. Instead of doing a first-order approximation in \(h\) the idea is to do a higher-order approximation of the derivative.

The basic idea is to construct a sequence \(\boldsymbol{x}_i\) giving an approximation of the solution of \(\dot{\boldsymbol{x}}(t)=f(t,\boldsymbol{x})\) at time \(t_i\) for \(i\in\mathbb{N}\). This sequence is defined by:

\[ \boldsymbol{x}_{i+1} = \boldsymbol{x}_i + M(t_i, \boldsymbol{x}_i, t_{i+1}-t_i) \,, \]

for a certain function \(M\) called method. For example, for the explicit Euler method, the function \(M\) is given by:

\[ M(t, \boldsymbol{x}, h) = h f(t,\boldsymbol{x}) \,. \]

A Runge-Kutta method of order two is given by:

\[ M(t,\boldsymbol{x}, h) = h f\bigg(t+\frac{h}{2}, \boldsymbol{x} + \frac{h}{2} f(t,\boldsymbol{x}) \bigg) \,. \]

A Runge-Kutta method of order four is given by:

\[ M(t,\boldsymbol{x}, h) = \frac{h}{6} \big(\boldsymbol{k}_1 + 2\boldsymbol{k}_2 + 2\boldsymbol{k}_3 + \boldsymbol{k}_4\big) \,, \]

where

\[ \begin{align*} \boldsymbol{k}_{1} & =f(t,\boldsymbol{x})\,,\\ \boldsymbol{k}_{2} & =f\bigg(t+\frac{h}{2},\boldsymbol{x}+\frac{h}{2}\boldsymbol{k}_{1}\bigg)\,,\\ \boldsymbol{k}_{3} & =f\bigg(t+\frac{h}{2},\boldsymbol{x}+\frac{h}{2}\boldsymbol{k}_{2}\bigg)\,,\\ \boldsymbol{k}_{4} & =f(t+h,\boldsymbol{x}+h\boldsymbol{k}_{3})\,. \end{align*} \]

Note that more generally, a Runge-Kutta method of order \(s\) is given by:

\[ M(t,\boldsymbol{x},h) = h \sum_{i=1}^s b_i \boldsymbol{k}_i \,, \]

where

\[ \begin{align*} \boldsymbol{k}_{1} & =f(t,\boldsymbol{x})\,,\\ \boldsymbol{k}_{2} & =f(t+c_{2}h,\boldsymbol{x}+ha_{21}\boldsymbol{k}_{1})\,,\\ \boldsymbol{k}_{3} & =f(t+c_{3}h,\boldsymbol{x}+h(a_{31}\boldsymbol{k}_{1}+a_{32}\boldsymbol{k}_{2}))\,,\\ & \vdots\\ \boldsymbol{k}_{s} & =f(t+c_{s}h,\boldsymbol{x}+h(a_{s1}\boldsymbol{k}_{1}+a_{s2}\boldsymbol{k}_{2}+\cdots+a_{s,s-1}\boldsymbol{k}_{s-1}))\,. \end{align*} \]

The coefficients \(a_{ij}\) (for \(1 \leq j < i \leq s\)), \(c_i\) (for \(2 \leq i \leq s\)), and \(b_i\) (for \(1 \leq i \leq s\)), are often represented in a Butcher table:

\[\def\arraystretch{1.4} \begin{array}{c|ccccc} 0\\ c_{2} & a_{21}\\ c_{3} & a_{31} & a_{32}\\ \vdots & \vdots & & \ddots\\ c_{s} & a_{s1} & a_{s2} & \cdots & a_{s,s-1}\\ \hline & b_{1} & b_{2} & \cdots & b_{s-1} & b_{s} \end{array}\]

For example, the Butcher array from the previous method of order two is:

\[\def\arraystretch{1.4} \begin{array}{c|cc} 0\\ \frac{1}{2} & \frac{1}{2}\\[2pt] \hline & 0 & 1 \end{array}\]

and that of the fourth-order method:

\[\def\arraystretch{1.4} \begin{array}{c|cccc} 0\\ \frac{1}{2} & \frac{1}{2}\\ \frac{1}{2} & 0 & \frac{1}{2}\\ 1 & 0 & 0 & 1\\ \hline & \frac{1}{6} & \frac{1}{3} & \frac{1}{3} & \frac{1}{6} \end{array}\]

a. Define a function integrate(f, x0, t, M) which for a given list of times \((t_i)_ {i=0}^N\) returns the corresponding values \(\boldsymbol{x}_0,\boldsymbol{x}_1,\dots,\boldsymbol{x}_N\) with method \(M\).

b. Implement the functions M(f,t,x,h) for the explicit Euler method and the Runge-Kutta method of order two. Compare the two methods.

c. Implement the function M(f,t,x,h) for the Runge-Kutta method of order four. Compare with the second-order method.

Exercise 11.3 - Movement of a planet

The goal is to simulate the two-dimensional motion of a planet orbiting around a fixed star. The star is supposed to be fixed at the origin and the position of the planet in the plane is described by the vector \(\boldsymbol{x}\in\mathbb{R}^2\). The star is supposed to interact with the planet with the potential:

\[ V(\boldsymbol{x}) = \frac{1}{\alpha} |\boldsymbol{x}|^\alpha \,, \]

for a certain \(\alpha\in\mathbb{R}\), where \(|\boldsymbol{x}|\) denotes the euclidean norm of the vector \(\boldsymbol{x}\). Note that the gravitational potential corresponds to \(\alpha=-1\). The equation of the planet in this force field is given by:

\[ \ddot{\boldsymbol{x}} = -\boldsymbol{\nabla}V(\boldsymbol{x}) = -\boldsymbol{x} |\boldsymbol{x}|^{\alpha-2} \,. \]

a. Rewrite the second-order differential equation as a first-order differential equation for \(\boldsymbol{x}\) and \(\boldsymbol{p}=\dot{\boldsymbol{x}}\).

b. Implement the function f(t,xp) corresponding to the equation found in the previous point.

c. Using the Runge-Kutta method of fourth-order, solve the differential equation for different initial data and different values of \(\alpha\) and plot the \(\boldsymbol{x}(t)\) trajectories in the plane. Interpret the results and explain in particular why the cases \(\alpha=-1\) and \(\alpha=2\) are different from the others.

Exercise 11.4 - Lorenz attractor

The Lorenz model is a system of three coupled differential equations of the form

\[ \begin{align*} \dot{x} & =\sigma(y-x) \,,\\ \dot{y} & =x(\rho-z)-y \,,\\ \dot{z} & =xy-\beta z \,, \end{align*} \]

where \(rho,\sigma,\beta\) are three real parameters. This is a very simplified model of coupling between the atmosphere and the ocean proposed in 1963 by Edward Lorenz.

a. Write mathematically the expression of the function \(f:\mathbb{R}\times\mathbb{R}^3\) allowing to put the Lorenz system in the form

\[ \dot{\boldsymbol{x}} = f(t,\boldsymbol{x}) \,, \]

where \(\boldsymbol{x}\) is the vector \((x,y,z)\). Implement a function f(t,x,rho,sigma,beta) corresponding to the function \(f\).

b. Write a function plot_lorenz(rho,sigma,beta) which for given parameters \(\rho,\sigma,\beta\), plots the trajectory \((x(t),z(t))\) for \(t\in[0,20]\) from the initial data \(\boldsymbol{x}_0=(x_0,y_0,z_0)=(1,1,1)\). Use, for example, the Runge-Kutta method of order four with a time step \(\Delta t = 0.001\). Test with \(\sigma=10\), \(\beta=8/3\), and the values \(\rho=10,15,20,25\), and describe what is observed.

c. Using SymPy, determine the stationary solutions according to the parameters \(\rho,\sigma,\beta\), i.e., the solutions of \(f(t,\boldsymbol{x})=\boldsymbol{0}\) for all \(t>0\). Interpret the previous graphs in the light of this.

Exercise 11.5 - Cubic wave equation !!

The goal is to solve numerically the nonlinear wave equation on \(\mathbb{R}\):

\[ \begin{align*}-\frac{\partial^{2}u}{\partial t^{2}}+\frac{\partial^{2}u}{\partial x^{2}} & =u^{3}\,, & u(0,\cdot) & =u_{0}\,, & \frac{\partial u}{\partial t}(0,\cdot) & =v_{0}\,,\end{align*} \]

for \(u:\mathbb{R}^+ \times \mathbb{R}\to\mathbb{R}\) with \(u_0, v_0:\mathbb{R}\to\mathbb{R}\) two given functions.

a. Rewrite the previous equation as two first-order equations in time for \(u\) and \(v=\frac{\partial u}{\partial t}\).

b. By approximating the second derivative in space by finite differences as in Exercise 9.4 show that the equation can be approximated as follows:

\[ \begin{align*} \frac{\partial u_{n}}{\partial t} & =v_{n} \,, & u_{n}(0) & =u_{0}(x_{n})\,,\\ \frac{\partial v_{n}}{\partial t} & =\frac{u_{n-1}-2u_{n}+u_{n+1}}{h^{2}}-u_{n}^{3} \,, & v_{n}(0) & =v_{0}(x_{n})\,, \end{align*} \]

where \((x_n)_ {n=0}^{N}\) denotes \(N+1\) evenly spaced points from \(h\) in the interval \([-L,L]\) and \(u_n(t) = u(t,x_n)\) and \(v_n(t) = v(t,x_n)\). For the conditions at the boundary of the domain, i.e., when \(n=0\) or \(n=N\), we take:

\[ \frac{\partial v_{0}}{\partial t} = 0 \,, \qquad \frac{\partial v_{N}}{\partial t} = 0 \,. \]

c. Determine the function \(f:\mathbb{R}^{2N+2} \to \mathbb{R}^{2N+2}\) allowing to put the previous approximation in the form \(\dot{\boldsymbol{u}} = f(t,\boldsymbol{u})\) for \(\boldsymbol{u}=(u,v)\) and implement this function.

d. Solve the differential equation given by \(\dot{\boldsymbol{u}} = f(t,\boldsymbol{u})\), for example, with the fourth-order Runge-Kutta method. A good choice of parameters is \(L=100\), \(N=1 000\) and for the initial data \(u_0(x)=e^{-x^2}\) and \(v_0(x)=0\). The speed of propagation of the wave is one and, therefore, after a time greater than \(L\), the wave leaves the box \([-L,L]\) and no longer corresponds to a good approximation of the initial equation.

e. Using the animation module of Matplotlib, make a video showing the evolution of the wave as a function of time.

Exercise 11.6 - Bogacki-Shampine methods !!!

By combining two Runge-Kutta methods of different orders (for example, (2,3) or (4,5)), one will obtain an empirical estimate of the error over a time step. Using this error estimate, it is possible to adapt the time step, either by increasing or decreasing it, and thus adapt to the equation.

For a Runge-Kutta method of order \(s\), an internal method of lower order (usually \(s-1\)) is given by:

\[ M^*(t,\boldsymbol{x},h) = h \sum_{i=1}^s b_i^* \boldsymbol{k}_i \,, \]

where the \(\boldsymbol{k}_i\) are identical to those of the \(s\) order method. An estimate of the error is then given by:

\[ E(t,\boldsymbol{x},h) = M(t,\boldsymbol{x},h) - M^*(t,\boldsymbol{x},h) =h \sum_{i=1}^s (b_i - b_i^*) \boldsymbol{k}_i \,. \]

Such a method is given by an extended Butcher table:

\[\def\arraystretch{1.4} \begin{array}{c|ccccc} 0\\ c_{2} & a_{21}\\ c_{3} & a_{31} & a_{32}\\ \vdots & \vdots & & \ddots\\ c_{s} & a_{s1} & a_{s2} & \cdots & a_{s,s-1}\\ \hline & b_{1} & b_{2} & \cdots & b_{s-1} & b_{s}\\ & b_{1}^* & b_{2}^* & \cdots & b_{s-1}^* & b_{s}^* \end{array}\]

a. Implement the Bogacki-Shampine method of order (4,5). The original article is available at doi:10.1016/0898-1221(96)00141-1.