{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# TD1 : Temporal Discretization\n",
"\n",
"The goal of temporal discretization methods is to numerically solve Cauchy problems such as :\n",
"$$\n",
" \\begin{cases}\n",
" u^\\prime = f(u)\\\\\n",
" u(0) = u_0\n",
" \\end{cases}\n",
"$$\n",
"where $u$ is typically a function $u : \\mathbb{R} \\rightarrow \\mathbb{R}^d$ with $d$ the dimension of the problem (for $d = 2$ we may also find $u : \\mathbb{R} \\rightarrow \\mathbb{C}$).\n",
"In this first TD we will review some of the basic methods to solve Ordinary Differential Equations (ODEs) before studying one of the most famous chaotic dynamical system : the Lorenz system."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Part 1 : The Exponential system\n",
"\n",
"The purpose of this part is to study several linear one-step methods to solve the following system on $[0,T]$ :\n",
"$$\\begin{cases}\n",
" u^\\prime=\\alpha u \\\\\n",
" u (0) = u_0\n",
"\\end{cases}\n",
"$$\n",
"where $u : [0,T] \\rightarrow \\mathbb{C}$, $\\alpha \\in \\mathbb{C}$ and $u_0 \\in \\mathbb{C}$.\n",
"We sample the interval $[0, T]$ with $N + 1$ equally spaced points $\\{0, \\Delta t, ..., N \\Delta t\\}$ defining the time step $\\Delta t$ as $\\Delta t = \\frac{T}{N}$. The solution of this Cauchy problem is analytical and we will note $u(t) = u_0 e^{\\alpha t}$ its unique solution. Our numerical attempts to solve this problem will be noted $u_{\\text{num},*}$ where the * character refers to the specified numerical scheme.\n",
"\n",
"More precisely, $u_{\\text{num},*} = (u_{\\text{num},*}^{(n)})_{n\\in\\{ 0 ,..., N\\}}$ will be a finite sequence of N + 1 (potentially complex) numbers such as $u_{\\text{num},*}^{(n)}$ should approximate the exact solution $u (n \\Delta t)$.\n",
"\n",
"In the following, Python codes will only need *numpy* and *matplotlib* packages."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Moreover we will use the following parameters in this part :"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"T = 10.0\n",
"N = 100\n",
"dt = T / N"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**1) Forwards Euler : ** We remind that Forwards Euler method (FE) is a first order method which computes the value of $u (t + \\Delta t)$ given $u(t)$ using the following Taylor expansion :\n",
"$$\n",
" u (t + \\Delta t) = u(t) + \\Delta t~f (u(t)) + \\mathcal{O}(\\Delta t ^2)\n",
"$$\n",
"Express $u_{\\text{num},\\text{FE}}^{(n + 1)}$ in terms of $u_{\\text{num},\\text{FE}}^{(n)}$. Deduce the expression of $u_{\\text{num},\\text{FE}}^{(n)}$ in terms of $u_{\\text{num},\\text{FE}}^{(0)} = u_0$.\n",
"\n",
"**2)** We define $q = \\alpha \\Delta t$. Program a Forwards Euler based solving method (using the recurrence relation) of the exponential system on $[0,T]$ and plot the numerical solution and the exact solution for several real values of $q$ and $u_0$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**3) Backwards Euler / Crank-Nicolson ** Same questions for Backwards Euler method and Crank-Nicolson method. We remind that Backwards Euler method relies on the following expansion :\n",
"$$\n",
" u (t + \\Delta t) = u(t) + \\Delta t~f(u(t + \\Delta t)) + \\mathcal{O}(\\Delta t ^2)\n",
"$$\n",
"and that Crank-Nicolson relies on the following expansion :\n",
"$$\n",
" u (t + \\Delta t) = u(t) + \\frac{\\Delta t}{2}\\left(f\\left(u\\left(t\\right)\\right) + f\\left(u\\left(t + \\Delta t\\right)\\right)\\right) + \\mathcal{O}(\\Delta t ^3)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**4)** Compare on this example the three previous methods by plotting the three numerical solutions next to the exact solution. We can use a value of $q$ which is included in the absolute stability region of every method."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**5)** Show that Crank-Nicolson scheme is indeed a second-order method."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Part 2 : The Harmonic Oscillator\n",
"\n",
"We now want to solve the following Cauchy problem $[0,T]$ :\n",
"$$\\begin{cases}\n",
" u^{\\prime\\prime} +\\omega^2 u =0\\\\\n",
" u(0) = u_0\\\\\n",
" u^\\prime(0) = u^\\prime_0\n",
"\\end{cases}\n",
"$$\n",
"where $u : [0,T] \\rightarrow \\mathbb{R}$, $\\omega\\in\\mathbb{R}$, $u_0 \\in \\mathbb{R}$ and $u^\\prime_0 \\in \\mathbb{R}$.\n",
"\n",
"**1)** Introducing $v = \\frac{u^\\prime}{\\omega}$, show that this sytem is equivalent to a complex exponential system such as defined in the previous part. We will express $\\alpha$ in terms of $\\omega$.\n",
"\n",
"**2)** Numerically solve this system for $\\omega = 1$ using Forwards Euler method, Backwards Euler method, Crank-Nicolson method."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**3)** Plot the previous numerical solutions in the $(u, v)$ plan. If we want to enforce energy conservation in this system, which numerical method should we use? Theoretically explain what we observe following the specified numerical method."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Part 3 : A chaotic dynamical system : the Lorenz system\n",
"\n",
"We are now interested in the Lorenz system, first studied by the mathematician and meteorologist Edward Lorenz to model the atmospheric convection (we will refer to the following [Wikipedia page](https://en.wikipedia.org/wiki/Lorenz_system)). The Lorenz system is defined as follows :\n",
"$$\\begin{cases}\n",
"x^\\prime = \\sigma(y - x)\\\\\n",
"y^\\prime = x (\\rho - z) - y\\\\\n",
"z^\\prime = xy - \\beta z\n",
"\\end{cases}\n",
"$$\n",
"where $x, y, z : \\mathbb{R} \\rightarrow \\mathbb{R}$ and $\\sigma, \\rho, \\beta \\in \\mathbb{R}$ are some numerical parameters. Lorenz first used the following values of the parameters :"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"rho = 28.0\n",
"sigma = 10.0\n",
"beta = 8/3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We want to solve the Lorenz system using Runge-Kutta multistep methods. We will call RK2, the following Runge-Kutta method :\n",
"$$\n",
" U_1 = u(t)\\\\\n",
" U_2 = u(t) + \\Delta t~f(U_1)\\\\\n",
" u(t + \\Delta t) = u(t) + \\frac{\\Delta t}{2}(f(U_1) + f(U_2)) + \\mathcal{O}(\\Delta t^3)\n",
"$$\n",
"\n",
"**1)** Show that RK2 is indeed a second-order method.\n",
"\n",
"**2)** Numerically solve the Lorenz system on the interval $[0,T]$ using RK2 method with the following initial values : $x (0) = y (0) = z (0) = 1$. We will use at least $T = 100$ and $N = 10000$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**3)** Plot the solution in the 3D-space $(x,y,z)$. We can use *[mpl_toolkits.mplot3d](https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html)* package for 3D plots."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We call this figure the *Lorenz attractor*.\n",
"The Lorenz system exhibits chaotic behaviors for some specific values of its parameters: this means that solutions are highly senstive to initial conditions.\n",
"\n",
"**4)** For the same values of the parameters, slightly change the initial conditions and observe this chaotic behavior."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}