{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# Exercises 8: Numerical Solution of ODE's 1\n",
" \n",
"**Alvar WallstrÃ¶m Lindell**\n",
"\n",
"Date: **Mar 10, 2021**\n",
"Submission deadline: **Mar 29, 2021**\n",
"\n",
"If you want to have a nicer theme for your jupyter notebook,\n",
"download the [cascade stylesheet file tma4125.css](https://www.math.ntnu.no/emner/TMA4125/2020v/part_II/notebooks/tma4125.css)\n",
"and execute the next cell:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" try:\n",
" with open(\"tma4125.css\", \"r\") as f:\n",
" styles = f.read()\n",
" return HTML(styles)\n",
" except FileNotFoundError:\n",
" pass #Do nothing\n",
"\n",
"# Comment out next line and execute this cell to restore the default notebook style \n",
"css_styling()"
]
},
{
"cell_type": "code",
"execution_count": 182,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib\n",
"matplotlib.rcParams.update({'font.size': 12})\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this exercise set you will be analyzing and implementing the following three explicit Runge-Kutta methods:\n",
"\n",
"Explicit Euler\n",
"$$\n",
"\\begin{array}{c|c}\n",
"0 & 0 \\\\ \\hline\n",
"& 1\n",
"\\end{array}\n",
"$$\n",
"Heun's Method\n",
"$$\n",
"\\begin{array}{c|cc}\n",
"0 & 0 & 0 \\\\ \n",
"1 & 1 & 0 \\\\ \\hline\n",
"& \\tfrac{1}{2} & \\tfrac{1}{2}\n",
"\\end{array}\n",
"$$\n",
"Classic 4-stage Runge-Kutta (RK4)\n",
"$$\n",
"\\begin{array}{c|cccc}\n",
"0 & 0 & 0 & 0 & 0 \\\\ \n",
"\\tfrac{1}{2} & \\tfrac{1}{2} & 0 & 0 & 0 \\\\ \n",
"\\tfrac{1}{2} & 0 &\\tfrac{1}{2} & 0 & 0 \\\\\n",
"1 & 0 & 0 & 1 & 0\\\\ \\hline\n",
"& \\tfrac{1}{6} & \\tfrac{1}{3} & \\tfrac{1}{3} & \\tfrac{1}{6}\n",
"\\end{array}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### **Exercise 1** \n",
"#### Convergence orders\n",
"Calculate analytically the convergence order's of the Explicit Euler method and Heun's method. If you want, you can also check that RK4 has order 4."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### **Exercise 2** \n",
"#### Implementing and testing the methods\n",
"In this exercise we will numerically solve the ODE\n",
"$$\n",
"y'(t) = f(y), \\quad y(0) = y_0\n",
"$$\n",
"in the interval $t \\in [0,T]$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**a)** Implement the three RK-methods, for example in the following way:\n",
"\n",
"1. For each method, implement a function which takes one step with a given step size. The function should take as arguments:\n",
" * The current solution value $y_n \\approx y(t_n)$\n",
" * The time point $t_n$\n",
" * The step size $h$\n",
" * The right-hand side function $f$\n",
" \n",
" The function should return the new solution value $y_{n+1} \\approx y(t_{n+1}) = y(t_n + h).$\n",
" \n",
"2. Implement a function which repeatedly applies a RK-method in the interval $[0,T]$. It should take as arguments:\n",
" * The initial value y_0\n",
" * The step size h\n",
" * The final time T\n",
" * The Runge-Kutta method\n",
" \n",
" The function should return two arrays:\n",
" * One array TT containing all the time-points \n",
"$0 = t_0,t_1,...,t_N = T$\n",
" \n",
" * One array YY containing all the function values \n",
"$y_0,y_1,...,y_N$\n",
"\n",
"Test the methods on the ODE\n",
"$$\n",
"y'(t) = -y(t), \\quad y(0) = 1, \\quad t \\in [0,1].\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**b)**\n",
"We will now numerically investigate the RK-methods. We can do this since we know what the exact solution to the ODE above is. We assume that the error $e = |y(T) - y_N|$ when using step size $h$ is approximately\n",
"$$\n",
"e \\approx Ch^p\n",
"$$\n",
"for some $C>0$ and $p.$ Note that $p$ is what we call the convergence order. We assume that $p$ and $C$ is the same when using different step sizes $h$. Let $e_1$ and $e_2$ be the errors when using step sizes $h_1$ and $h_2$. Then we have\n",
"\n",
"$$\n",
"\\frac{e_1}{e_2} \\approx \\frac{h_1^p}{h_2^p} = \\left(\\frac{h_1}{h_2}\\right)^p.\n",
"$$\n",
"Taking logarithms on both sides we get\n",
"$$\n",
"\\log(e_1/e_2) \\approx p \\log(h_1/h_2)\n",
"$$\n",
"or\n",
"$$\n",
"p \\approx \\frac{\\log(e_1/e_2)}{\\log(h_1/h_2)}.\n",
"$$\n",
"The value on the right-hand side of this equation is what we call the Experimental Order of Convergence, or EOC. We will now try to estimate the order of convergence using EOC-values.\n",
"\n",
"Do the following for each method:\n",
"1. For j=0,...,5, find the numerical solution $y_N$ of the ODE at $t = 10$ using step size $h_j = 2^{-j}$.\n",
"2. Calculate the error $e_j = |y(10) - y_N|$.\n",
"3. Calculate the EOC for neighbouring step sizes, that is using equation for $p$ above above with $e_j$ and $e_{j+1}$ for $j=0,...,4$. This should give you $5$ different approximations. \n",
"\n",
"Draw a conclusion about the order of convergence $p$ for each method. Does it agree with the result in exercise 1?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**c)**\n",
"We will finally test the three methods on the ODE\n",
"\n",
"$$\n",
"y'(t) = -2ty(t), \\quad y(0) = 1, \\quad t \\in [0,0.4].\n",
"$$\n",
"\n",
"This has exact solution $e^{-t^2}$. Find the approximate value of $y(0.4)$ using \n",
" * The Explicit Euler method with $h = 0.1$\n",
" * Heun's method with $h = 0.2$\n",
" * The RK4 method with $h = 0.4$ \n",
" \n",
"The step sizes are chosen such that each method needs to perform 4 function evaluations. How do the errors compare?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### **Exercise 3**\n",
"#### SIR Model\n",
"\n",
"The SIR model is a system of first oder ODE's which model the dynamics of disease. There are $S(t)$ susceptible/healthy individuals, $I(t)$ infected individuals and $R(t)$ recovered individuals. Each susceptible person has a risk of becoming infected, a risk which is proportional to the number of infected people $I(t)$, with a proportinality constant $\\beta>0.$ Each infected person also has a chance of recovering, with a recovery constant $\\gamma>0$. This leads to the coupled system of first-order ODE's\n",
"\n",
"$$\n",
"S'(t) = -\\beta S(t) I(t),\\\\\n",
"I'(t) = \\beta S(t) I(t) - \\gamma I(t),\\\\\n",
"R'(t) = \\gamma I(t).\n",
"$$\n",
"\n",
"We can rewrite this in vector form as\n",
"\n",
"$$\n",
"\\mathbf{u}'(t) = \\mathbf{f}(\\mathbf{u}(t))\n",
"$$\n",
"\n",
"where we have defined\n",
"\n",
"$$\n",
"\\mathbf{u}(t) = \n",
"\\begin{pmatrix}\n",
"S(t)\\\\\n",
"I(t)\\\\\n",
"R(t)\n",
"\\end{pmatrix}, \\quad \n",
"\\mathbf{f}(\\mathbf{u}(t)) = \n",
"\\begin{pmatrix}\n",
"-\\beta S(t) I(t)\\\\\n",
"\\beta S(t) I(t) - \\gamma I(t)\\\\\n",
"\\gamma I(t)\n",
"\\end{pmatrix}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**a)** Show that the system is conservative, that is that the total number of individuals $S(t)+I(t)+R(t)$ is constant."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**b)** Numerically solve the system \n",
"\n",
"$$\n",
"\\mathbf{u}'(t) = \\mathbf{f}(\\mathbf{u}(t)),\\quad \\mathbf{u}(0) = \\mathbf{u}_0,\\quad t\\in [0,T]\n",
"$$\n",
"\n",
"with either of the RK-methods. You may pick end-time $T$ and step size $h$, as well as initial number of individuals. A suitable inital conditions could e.g. be $\\mathbf{u}_0 = (50,10,0)^T$. Plot the solution as a function of time. Also plot the total number of individuals. Is the total number conserved? To check this, you might calculate the maximum total and the minimum total over the interval."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### **Exercise 4**\n",
"#### A second-order system (Not obligatory)\n",
"\n",
"For this final exercise, consider three weights, connected through springs, hanging from the ceiling. The weights have masses $m_0$, $m_1$ and $m_2$, the springs have spring constants $k_0$, $k_1$ and $k_2$, and natural lengths $l_0$, $l_1$ and $l_2$. The weights are hanging at a distance $z_0$, $z_1$ and $z_2$ from the ceiling."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The total kinetic energy is \n",
"$$T = \\sum_{i=0}^2 \\tfrac{1}{2}m_i\\dot{z_i}^2,$$\n",
"the potential energy arizing from the gravitaional forces is\n",
"$$\n",
"V_G = -\\sum_{i=0}^2 m_i g z_i\n",
"$$\n",
"and the potential energy arizing from the spring deformation is\n",
"$$\n",
"V_S = \\sum_{i=0}^2 \\tfrac{1}{2}k_i(z_i-z_{i-1}-l_i)^2.\n",
"$$\n",
"Here we have defined $z_{-1} = 0$ and the dot denoting derivative with respect to time. The constant $g$ is the gravitational acceleration, which in general is approximately $9.81$. The total energy is then $H = T+V_G+V_S$. By either balancing forces or (for those who know it) using Lagrangian mechanics we can derive the equations of motion:\n",
"\n",
"$$\n",
"\\ddot{z_0} =g - \\tfrac{k_0}{m_0}(z_0-l_0) + \\tfrac{k_1}{m_0}(z_1 - z_0 - l_1),\\\\\n",
"\\ddot{z_1} =g - \\tfrac{k_1}{m_1}(z_1 - z_0 - l_1) + \\tfrac{k_2}{m_1}(z_2-z_1-l_2),\\\\\n",
"\\ddot{z_2} =g - \\tfrac{k_2}{m_2}(z_2-z_1-l_2).\n",
"$$\n",
"\n",
"This is a system of three coupled second-order ODE's. We want to rewrite this as a system of six first-order ODE's. We do this by introducing the variables $w_i = \\dot{z_i}$. The system can then be rewritten as \n",
"\n",
"$$\n",
"\\dot{w_0} =g - \\tfrac{k_0}{m_0}(z_0-l_0) + \\tfrac{k_1}{m_0}(z_1 - z_0 - l_1),\\\\\n",
"\\dot{w_1} =g - \\tfrac{k_1}{m_1}(z_1 - z_0 - l_1) + \\tfrac{k_2}{m_1}(z_2-z_1-l_2),\\\\\n",
"\\dot{w_2} =g - \\tfrac{k_2}{m_2}(z_2-z_1-l_2),\\\\\n",
"\\dot{z_0} =w_0,\\\\\n",
"\\dot{z_1} =w_1, \\\\\n",
"\\dot{z_2} =w_2,\n",
"$$\n",
"\n",
"or, in vector form \n",
"\n",
"$$\n",
"\\dot{\\mathbf{u}} = \\mathbf{F}(\\mathbf{u})\n",
"$$ \n",
"\n",
"where \n",
"\n",
"$$\n",
"\\mathbf{u} = \n",
"\\begin{pmatrix}\n",
"u_0\\\\\n",
"u_1\\\\\n",
"u_2\\\\\n",
"u_3\\\\\n",
"u_4\\\\\n",
"u_5\n",
"\\end{pmatrix} = \n",
"\\begin{pmatrix}\n",
"w_0\\\\\n",
"w_1\\\\\n",
"w_2\\\\\n",
"z_0\\\\\n",
"z_1\\\\\n",
"z_2\n",
"\\end{pmatrix}\n",
"$$\n",
"\n",
"and\n",
"\n",
"$$\n",
"\\mathbf{F}(\\mathbf{u}) = \n",
"\\begin{pmatrix}\n",
"g - \\tfrac{k_0}{m_0}(u_3-l_0) + \\tfrac{k_1}{m_0}(u_4 - u_3 - l_1)\\\\\n",
"g - \\tfrac{k_1}{m_1}(u_4 - u_3 - l_1) + \\tfrac{k_2}{m_1}(u_5-u_4-l_2)\\\\\n",
"g - \\tfrac{k_2}{m_2}(u_5-u_4-l_2)\\\\\n",
"u_0\\\\\n",
"u_1\\\\\n",
"u_2\n",
"\\end{pmatrix}.\n",
"$$\n",
"\n",
"We can now solve this system using our RK-methods. In order for the system to have a unique solution we do not only need initial positions $z_0(0), z_1(0)$ and $z_2(0)$, but also initial velocities $w_0(0), w_1(0)$ and $w_2(0)$.\n",
"\n",
"\n",
"**a)** Try to solve the system with some initial conditions and parameters $m_0,m_1,m_2$, $k_0,k_1,k_2$ and $l_0,l_1,l_2$. Plot the positions $z_0(t), z_1(t)$ and $z_2(t)$ versus time $t$. Play around with different methods and step sizes $h$. What happens when you chose $h$ too large?\n",
"\n",
"**b)** Plot the total energy $H$ as a function of time. In particular, calculate the maximum difference in the energy. How does this change when you use different methods and different step sizes $h$?"
]
}
],
"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.7.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}