{
"cells": [
{
"cell_type": "markdown",
"id": "56889fcd-f08a-4475-9928-af817748f5c8",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-5e2b58ff19f0997c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## Numerical Analysis - Fall semester 2025\n",
"# Serie 11 - Romberg integration and ODEs"
]
},
{
"cell_type": "markdown",
"id": "5d0c899b-58da-49e8-ace0-828c8d5b60ca",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-e56f5ef08590e711",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"First, we will need to import some of the usual packages. You will have to run this cell every time you restart your notebook."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7d640d3e-d688-47ad-95b6-96eed7152684",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-6bbe893666203fcf",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import scipy as sp"
]
},
{
"cell_type": "markdown",
"id": "f63dc15a-a707-43e2-9722-59d5ba8eb67f",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-cf5c15060376e2e7",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Romberg integration method"
]
},
{
"cell_type": "markdown",
"id": "27fd5981-14a6-408f-9f16-00a774a2be5f",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-95f85d5601bd8908",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Given $a, b \\in \\mathbb{R}$, $a < b$, $n \\in \\mathbb{N} \\setminus \\{0, 1\\}$, and $f \\in \\mathrm{C}^6([a, b], \\mathbb{R})$, the Euler--Maclaurin expansion ensures the existence of $c_1, c_2 \\in \\mathbb{R}$ such that\n",
"$$\n",
"\\int_a^b f = T(n) + c_1 h^2 + c_2 h^4 + \\mathrm{O}(h^6),\n",
"$$\n",
"where $h := (b-a)/n$ and $T(n) := h\\left((f(a)+f(b))/2+\\sum_{i=1}^{n-1} f(a+hi)\\right)$ is the composite trapezium formula. "
]
},
{
"cell_type": "markdown",
"id": "4ceed02c-e5db-4c1d-9b9f-9be98c53aeb6",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-6378ec916260fd6f",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 1 (Theoretical):** \n",
"Propose a formula $T_2(n)$ relying only on $T(n)$, $T(2n)$, and $T(4n)$ such that\n",
"$$\n",
"\\int_a^b f = T_2(n) + \\mathrm{O}(h^6).\n",
"$$\n",
"How many evaluations of $f$ does $T_2(n)$ require?\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"Define $T_1(n) := (4T(2n)-T(n))/3$. By proceeding as in the lecture, we obtain\n",
"$$\n",
"\\int_a^b f = T_1(n) + \\frac{3}{4}c_2 h^4 + \\mathrm{O}(h^6).\n",
"$$\n",
"Thus,\n",
"$$\n",
"\\int_a^b f = T_1(2n) + \\frac{3}{64}c_2 h^4 + \\mathrm{O}(h^6).\n",
"$$\n",
"Therefore,\n",
"$$\n",
"\\int_a^b f = T_2(n) + \\mathrm{O}(h^6)\n",
"$$\n",
"with\n",
"$$\n",
"T_2(n) := \\frac{16T_1(2n)-T_1(n)}{15}.\n",
"$$\n",
"Computing $T(4n)$ requires evaluating $f$ at $a+\\frac{h}{4}i$ for all $i \\in \\{0, \\dots, 4n\\}$, which represents $4n+1$ evaluations of $f$.\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "c8c1936d-5795-41cf-845d-3fb1a490f02a",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-d1011fe74ddbe047",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Theoretical tools for ODEs"
]
},
{
"cell_type": "markdown",
"id": "a7651f50-0778-46d9-a531-c96c24ce8038",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-3c49e73133deb1ac",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 2 (Theoretical):** \n",
"Consider the autonomous vector field associated with the logistic growth model, namely $f : \\mathbb{R} \\to \\mathbb{R} : x \\mapsto x(1-x)$. For every $u_0 \\in [0, \\infty)$, compute the maximal integral curve $u$ of $f$ such that $u(0) = u_0$.\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"The function $f$ is continuously differentiable, hence locally Lipschitz continuous. Thus, the initial-value problem admits a unique maximal solution $u$. If $u_0 \\in \\{0, 1\\}$, i.e., $u_0$ is an equilibrium point of $f$, then the maximal solution is the constant function $\\mathbb{R} \\to \\mathbb{R} : t \\mapsto u_0$. Otherwise, by global uniqueness, for every $t$ in the interval of existence, $u(t) \\in (1, \\infty)$ if $u_0 \\in (1, \\infty)$ and $u(t) \\in (0, 1)$ if $u_0 \\in (0, 1)$. Thus, on the interval of existence,\n",
"\\begin{equation*}\n",
"\\frac{u'}{u(1-u)} = 1\n",
"\\end{equation*}\n",
"and, for every $t$ in the interval of existence,\n",
"\\begin{equation*}\n",
"\\int_0^t \\frac{u'}{u(1-u)} = t\n",
"\\end{equation*}\n",
"or, equivalently,\n",
"\\begin{equation*}\n",
"\\int_{u_0}^{u(t)} \\frac{1}{x(1-x)} \\mathrm{d}x = t.\n",
"\\end{equation*}\n",
"Since, for all $x \\in \\mathbb{R} \\setminus \\{0, 1\\}$,\n",
"\\begin{equation*}\n",
"\\frac{1}{x(1-x)} = \\frac{1}{2} \\left(\\frac{1}{x}+\\frac{1}{1-x}\\right),\n",
"\\end{equation*}\n",
"it holds that\n",
"\\begin{align*}\n",
"\\int_{u_0}^{u(t)} \\frac{1}{x(1-x)} \\mathrm{d}x\n",
"&= \\left\\{\\begin{array}{ll}\n",
"\\left[\\ln x - \\ln(1-x)\\right]_{x=u_0}^{x=u(t)} & \\text{if } u_0 \\in (0, 1),\\\\[2mm]\n",
"\\left[\\ln x - \\ln(x-1)\\right]_{x=u_0}^{x=u(t)} & \\text{if } u_0 \\in (1, \\infty).\n",
"\\end{array}\\right.\\\\\n",
"&= \\left\\{\\begin{array}{ll}\n",
"\\left[\\ln\\left(\\dfrac{x}{1-x}\\right)\\right]_{x=u_0}^{x=u(t)} & \\text{if } u_0 \\in (0, 1),\\\\[4mm]\n",
"\\left[\\ln\\left(\\dfrac{x}{x-1}\\right)\\right]_{x=u_0}^{x=u(t)} & \\text{if } u_0 \\in (1, \\infty).\n",
"\\end{array}\\right.\\\\\n",
"&= \\ln\\left(\\frac{u(t)(u_0-1)}{(u(t)-1)u_0}\\right).\n",
"\\end{align*}\n",
"Hence,\n",
"\\begin{equation*}\n",
"\\ln\\left(\\frac{u(t)(u_0-1)}{(u(t)-1)u_0}\\right) = t\n",
"\\end{equation*}\n",
"or, equivalently,\n",
"\\begin{equation*}\n",
"u(t) = \\frac{u_0}{u_0+(1-u_0)\\exp(-t)}.\n",
"\\end{equation*}\n",
"In conclusion, the maximal solution is:\n",
"- $u : \\mathbb{R} \\to \\mathbb{R} : t \\mapsto u_0$ if $u_0 \\in \\{0, 1\\}$;\n",
"- $u : \\mathbb{R} \\to \\mathbb{R} : t \\mapsto \\dfrac{u_0}{u_0+(1-u_0)\\exp(-t)}$ if $u_0 \\in (0, 1)$;\n",
"- $u : (\\ln(1-\\frac{1}{u_0}), \\infty) \\to \\mathbb{R} : t \\mapsto \\dfrac{u_0}{u_0+(1-u_0)\\exp(-t)}$ if $u_0 \\in (1, \\infty)$.\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "4e019a05-0667-4229-8d66-77f1f60f3cda",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-b6258d4954827230",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 3 (Theoretical):** \n",
"Compute the solution to the IVP\n",
"\\begin{align*}\n",
"\\left\\{\\begin{array}{l}\n",
"u'(t) = Au(t) \\text{ for all } t \\in \\mathbb{R}\\\\\n",
"u(0) = u_0\n",
"\\end{array}\\right.&&\n",
"\\text{where}&&\n",
"A \\coloneq \\begin{bmatrix} 1 & -1\\\\ 1 & 1 \\end{bmatrix}&&\n",
"\\text{and}&&\n",
"u_0 \\in \\mathbb{R}^2.\n",
"\\end{align*}\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"The eigendecomposition\n",
"\\begin{equation*}\n",
"A = \\frac{1}{2} \\begin{bmatrix} 1 & \\imath\\\\ \\imath & 1 \\end{bmatrix} \\mathrm{diag}(1-\\imath, 1+\\imath) \\begin{bmatrix} 1 & -\\imath\\\\ -\\imath & 1 \\end{bmatrix}\n",
"\\end{equation*}\n",
"gives, after a few calculations,\n",
"\\begin{equation*}\n",
"u(t) = \\exp(t) \\begin{bmatrix} \\cos t & -\\sin t\\\\ \\sin t & \\cos t \\end{bmatrix} u_0.\n",
"\\end{equation*}\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "2a23c5ca-7434-4632-b3a5-1d37c1b24a69",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-4e32ab71461f3c03",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 4 (Theoretical):** \n",
"Consider the autonomous vector field associated with a damped pendulum, namely $f : \\mathbb{R}^2 \\to \\mathbb{R}^2 : (x_1, x_2) \\mapsto (x_2, - \\sin x_1 - x_2)$. What are the equilibrium points of $f$? Study their stability using Lyapunov's first method.\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"The set of equilibrium points is $\\pi \\mathbb{Z} \\times \\{0\\}$. The Jacobian matrix of $f$ at $(x_1, x_2) \\in \\mathbb{R}^2$ is\n",
"\\begin{equation*}\n",
"J(x_1, x_2) = \\begin{bmatrix} 0 & 1\\\\ -\\cos x_1 & -1 \\end{bmatrix}.\n",
"\\end{equation*}\n",
"For all $k \\in \\mathbb{Z}$,\n",
"\\begin{equation*}\n",
"J(k\\pi, 0) = \\begin{bmatrix} 0 & 1\\\\ (-1)^{k+1} & -1 \\end{bmatrix}.\n",
"\\end{equation*}\n",
"For every $k \\in \\mathbb{Z}$, the eigenvalues of $J(k\\pi, 0)$ are $(-1\\pm\\sqrt{5})/2$ if $k$ is odd and $(-1\\pm\\imath\\sqrt{3})/2$ if $k$ is even. In conclusion, for every $k \\in \\mathbb{Z}$, the equilibrium point $(k\\pi, 0)$ is asymptotically stable if $k$ is even (the pendulum is at its minimum height) and unstable if $k$ is odd (the pendulum is at its maximum height).\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "7f504b66-6e8b-454d-8a9b-7810e53e49b0",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-2f76ac45431af2b3",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Explicit methods for population model"
]
},
{
"cell_type": "markdown",
"id": "40813ae6-1d05-46ac-93e5-6a9cccdd71d9",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-c92451a8da9db056",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Let $f : \\mathbb{R} \\times \\mathbb{R}^n \\to \\mathbb{R}^n$ and $(t_0, u_0) \\in \\operatorname{dom} f$. Let us consider an ordinary differential equation of the form\n",
"\n",
"$$\n",
"\\begin{cases}\n",
" u'(t) = f(t, u(t)), & t \\in (0, \\tau], \\\\\n",
" u(0)=u_0.\n",
"\\end{cases}\n",
"$$\n",
"\n",
"We want to approximate $u(t_i)$ at the discrete time-steps $t_i = i h$, $i = 0, 1, 2, \\ldots,m$, where $h = \\tau/m$ is the step size and $m$ the number of time-steps. "
]
},
{
"cell_type": "markdown",
"id": "4b008597-b7bf-4761-a266-5c6b696a90ea",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-c21721b11e77180b",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 5:** Complete the function `forward_euler` which implements the forward Euler method to compute the approximations $u_i \\approx u(t_i)$ for $i = 0, 1, 2, \\dots, m$.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "654958ce-11b9-4b26-a87c-a92d60668008",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-a41bad7f587b2cb7",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def forward_euler(f, u_0, tau, m):\n",
" ### BEGIN SOLUTION\n",
" h = tau / m\n",
" u = []\n",
" u.append(u_0)\n",
" for i in range(m):\n",
" u.append(u[i] + h * f(i * h, u[i]))\n",
" ### END SOLUTION\n",
" return u"
]
},
{
"cell_type": "markdown",
"id": "79c04a9d-0932-4de2-9508-8c64be7f13b0",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-a9a51259c1344ec6",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 6:** Complete the function `heun` which implements the Heun method to compute the approximations $u_i \\approx u(t_i)$ for $i = 0, 1, 2, \\dots, m$.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c2c7cb2a-fd72-4ca2-9888-88b4da58951c",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-3e22f0ffd534b765",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def heun(f, u_0, tau, m):\n",
" ### BEGIN SOLUTION\n",
" h = tau / m\n",
" u = []\n",
" u.append(u_0)\n",
" for i in range(m):\n",
" u.append(u[i] + h / 2 * f(i * h, u[i]) + h / 2 * f((i + 1) * h, u[i] + h * f(i * h, u[i])))\n",
" ### END SOLUTION\n",
" return u"
]
},
{
"cell_type": "markdown",
"id": "96ca18ab-baf9-4424-9584-7b5b246ef134",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-546d0e879a54168d",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Let us consider a population of $u$ individuals in an environment\n",
"where at most $u_{\\max}=1000$ individuals can coexist. We suppose that initially\n",
"the number of individuals is $u_0 = 100$ and that the growth factor equals\n",
"a constant $C=2/15$. The considered model for the evolution of the\n",
"population is the following:\n",
"\n",
"$$\n",
"\\left\\{\\begin{array}{ll}\n",
" u'(t)=Cu(t)\\left(1-\\frac{u(t)}{u_{\\max}}\\right), & t \\in (0,100], \\\\\n",
" u(0)=u_0. & \n",
"\\end{array}\n",
"\\right.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "9783ac31-f06b-4a8f-b84b-4327fbf300d7",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-be34d926e36a1dbd",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 7:** Use the two functions you have implemented above to compute the approximate solution $u_i$, $i=0,1,...,m$ for $m=20$ time-steps and plot on the same graph the obtained numerical solutions in terms of time. Compare the obtained approximations for both forward Euler and Heun with the exact solution given by\n",
" $$\n",
" u(t) = \\frac{u_0u_\\max}{u_0+(u_\\max-u_0)\\exp(-Ct)}.\n",
" $$\n",
" \n",
" Which method gives a better approximation?\n",
"
\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d85d7412-63ff-4bb5-920d-f6e3fcb8aea5",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-ee50b5b4981d323b",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"\n",
"# RHS of ODE\n",
"def f(t, u, u_max=1000, C=2/15):\n",
" return C * u * (1 - u / u_max)\n",
"\n",
"# Exact solution\n",
"def u(t, u0, u_max=1000, C=2/15):\n",
" return 1 / ((1 - np.exp(- C * t)) / u_max + np.exp(- C * t) / u0)\n",
"\n",
"# Parameters\n",
"tau = 100\n",
"u_0 = 100\n",
"m = 20\n",
"\n",
"# Solve ODEs with the two methods\n",
"u_forward_euler = forward_euler(f, u_0, tau, m)\n",
"u_heun = heun(f, u_0, tau, m)\n",
"\n",
"# Determine exact solution\n",
"t = np.linspace(0, tau, m + 1)\n",
"u_exact = u(t, u_0)\n",
"\n",
"# Visualize the results\n",
"plt.plot(t, u_forward_euler, label=\"forward Euler\")\n",
"plt.plot(t, u_heun, label=\"Heun\")\n",
"plt.plot(t, u_exact, color=\"black\", linestyle=\"--\", label=\"exact\")\n",
"plt.ylabel(r\"$u(t)$\")\n",
"plt.xlabel(r\"$t$\")\n",
"plt.grid()\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"# Comparing the obtained approximations with the exact solution, we notice that\n",
"# Heun gives a better approximation of the true solution than forward Euler. \n",
"# This is expected, because Heun's method is more accurate than forward Euler\n",
"# as it uses a corrector step to refine the solution by averaging slopes and\n",
"# the forward Euler method tends to under- or overshoot, especially for larger\n",
"# time steps, due to its simple one-step approximation.\n",
"\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "e6b8708e-df6d-46d0-a55a-632402c287a8",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-e07241bfa6f26b16",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 8:** Repeat the last exercise using $m=2000$ time-steps. Is the obtained approximation with $m=2000$ better than the one obtained with $m=20$?\n",
"
\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "62cca3be-75ce-4b98-b278-a460c05e8394",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-ecc3d186889f7cfa",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"\n",
"# Parameters\n",
"tau = 100\n",
"u_0 = 100\n",
"m = 2000\n",
"\n",
"# Solve ODEs with the two methods\n",
"u_forward_euler = forward_euler(f, u_0, tau, m)\n",
"u_heun = heun(f, u_0, tau, m)\n",
"\n",
"# Determine exact solution\n",
"t = np.linspace(0, tau, m + 1)\n",
"u_exact = u(t, u_0)\n",
"\n",
"# Visualize the results\n",
"plt.plot(t, u_forward_euler, label=\"forward Euler\")\n",
"plt.plot(t, u_heun, label=\"Heun\")\n",
"plt.plot(t, u_exact, color=\"black\", linestyle=\"--\", label=\"exact\")\n",
"plt.ylabel(r\"$u(t)$\")\n",
"plt.xlabel(r\"$t$\")\n",
"plt.grid()\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"# The lines for the forward Euler, Heun's method, and the exact solution are almost\n",
"# indistinguishable from each other. This suggests that the numerical methods with\n",
"# a large number of time steps (m = 2000) closely approximate the exact solution.\n",
"# The fact that the numerical solutions align so well with the exact solution\n",
"# indicates that the step size (h = tau/m) is sufficiently small to reduce numerical\n",
"# error significantly.\n",
"\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "3065e471-f22d-42dd-a7dd-7d6156139e36",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-5b39b5b957674ceb",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"## The end\n",
"\n",
"Great! You have reached the end of the eleventh exercise notebook."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.12.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}