{ "cells": [ { "cell_type": "markdown", "id": "56889fcd-f08a-4475-9928-af817748f5c8", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "a70d6ddc354b5213a64848beb11c2028", "grade": false, "grade_id": "cell-5e2b58ff19f0997c", "locked": true, "schema_version": 3, "solution": false, "task": false }, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## Numerical Analysis - Fall semester 2025\n", "# Serie 12 - Numerical solution of ODEs" ] }, { "cell_type": "markdown", "id": "5d0c899b-58da-49e8-ace0-828c8d5b60ca", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "10edd7f8d5f42ef44d458651f4e462e8", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "575f9cc0d9b7147c8556ab42de3acc33", "grade": false, "grade_id": "cell-6bbe893666203fcf", "locked": true, "schema_version": 3, "solution": false, "task": false }, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import matplotlib as mpl" ] }, { "cell_type": "markdown", "id": "29297903-3767-4724-a9fd-26446df83c6b", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "95c52e93d8981d6303bdbc394245a3fb", "grade": false, "grade_id": "cell-1ee29b227cd9bc29", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "### Implicit methods for simple ODE" ] }, { "cell_type": "markdown", "id": "18718143-81d6-4dd8-a680-ba60ae8c288c", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "5fdfc9287632817678b637f0ab653e68", "grade": false, "grade_id": "cell-4fd86941baed31c7", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "Given is the linear system of ordinary differential equations\n", "$$\n", "\\begin{cases}\n", " u'(t) = \\begin{pmatrix} u_x'(t) \\\\ u_y'(t) \\end{pmatrix} = \\begin{pmatrix} u_y(t) \\\\ - u_x(t) \\end{pmatrix}, & t \\in (0, 20], \\\\\n", " {u}(0)=\\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix}.\n", "\\end{cases}\n", "$$\n", "\n", "The exact solution of this system is\n", "\n", "$$\n", "{u}(t) = \\begin{pmatrix} \\cos(t) \\\\ - \\sin(t) \\end{pmatrix}\n", "$$" ] }, { "cell_type": "markdown", "id": "5eede764-1f2e-4205-b8f7-babc55a57972", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "3df34c50c12745c9affa1596deda7d3f", "grade": false, "grade_id": "cell-2fecdbf7092cce30", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "We provide you with the function `plot_solution`, which visualizes the approximate solutions ${u}_i$ at the time steps $t_i$ for $i=0, 1, 2, \\dots, m$. It takes as input the time steps $t_0, t_1, t_2, \\dots, t_m$ in a NumPy array `t` of size $m+1$, and the corresponding solutions ${u}_0, {u}_1, {u}_2, \\dots, {u}_m$ collected in a list." ] }, { "cell_type": "code", "execution_count": null, "id": "6b60ce76-92a1-4768-8969-a2cf3e21a8c5", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "40e69d817914d9c7e238e4eab6624950", "grade": false, "grade_id": "cell-226919bf0ecfd5cd", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "outputs": [], "source": [ "def plot_solution(t, u):\n", " fig, ax = plt.subplots()\n", "\n", " u = np.asarray(u)\n", " \n", " # Determine x- and y-values as well as midpoints\n", " x = u[:, 0]\n", " y = u[:, 1]\n", " x_midpts = np.hstack((x[0], 0.5 * (x[1:] + x[:-1]), x[-1]))\n", " y_midpts = np.hstack((y[0], 0.5 * (y[1:] + y[:-1]), y[-1]))\n", "\n", " # Define line segments\n", " coord_start = np.column_stack((x_midpts[:-1], y_midpts[:-1]))[:, np.newaxis, :]\n", " coord_mid = np.column_stack((x, y))[:, np.newaxis, :]\n", " coord_end = np.column_stack((x_midpts[1:], y_midpts[1:]))[:, np.newaxis, :]\n", " segments = np.concatenate((coord_start, coord_mid, coord_end), axis=1)\n", " linewidths = np.linspace(0, 6, len(t))[::-1]\n", " lc = mpl.collections.LineCollection(segments, linewidths=linewidths, cmap=\"plasma\", capstyle=\"butt\")\n", " lc.set_array(t)\n", "\n", " # Plot the line segments\n", " lines = ax.add_collection(lc)\n", " fig.colorbar(lines, label=r\"parameter value $t$\")\n", " ax.scatter(x[0], y[0], color=mpl.colormaps[\"plasma\"](0), s=70)\n", " ax.annotate(r\"$u_0$\", (x[0], y[0]), xytext=(-17, 5), textcoords='offset points')\n", " ax.margins(0.2, 0.2)\n", " ax.add_patch(mpl.patches.Circle((0, 0), radius=1, fill=False, edgecolor=\"black\", linestyle=\"--\", linewidth=2))\n", " plt.xlabel(r\"$u_x(t)$\")\n", " plt.ylabel(r\"$u_y(t)$\")\n", " plt.grid()\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "6e69a566-ad49-4dba-a666-1f9090e659ec", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "df1e44e558209e278c71c741b9876939", "grade": false, "grade_id": "cell-2167b1fa3c4822c3", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "Below, we visualize the exact solution for $m = 200$ time steps in $[0, 20]$, which we can obtain from the function `u`." ] }, { "cell_type": "code", "execution_count": null, "id": "ba6bee86-c014-4221-b8e4-e1f8efa2c9b1", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "073250209fd3bb0ef227877ee8b3b423", "grade": false, "grade_id": "cell-38af7845654bf536", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "outputs": [], "source": [ "# Function which returns the exact solutions at the values t\n", "def u(t, u_0):\n", " return np.array([np.cos(t) * u_0[0] - np.sin(t) * u_0[1], - np.sin(t) * u_0[0] - np.cos(t) * u_0[1]]).T\n", "\n", "# Parameters\n", "tau = 20\n", "u_0 = np.array([1, 0])\n", "m = 200\n", "t = np.linspace(0, tau, m + 1)\n", "\n", "# Plots the exact solution\n", "u_exact = u(t, u_0)\n", "plot_solution(t, u_exact)" ] }, { "cell_type": "markdown", "id": "bdccf1ce-87e4-4369-88b7-d4c6abdeca8b", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "99831f428c6b9e4ad8280c4d5011fc5a", "grade": false, "grade_id": "cell-747105b80de40b4f", "locked": true, "schema_version": 3, "solution": false, "task": false }, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "We now want to apply different methods for approximating this ODE. To this extent, we provide you with the function `newton_method` which you have implemented in `serie04.ipynb`. Given a function $g : \\mathbb{R}^n \\to \\mathbb{R}^n$ and its Jacobian $Jg : \\mathbb{R}^n \\to \\mathbb{R}^{n \\times n}$, it iteratively determines a point $x$ for which $g(x) = 0$. This function will be useful in the implementation of the implicit ODE solvers." ] }, { "cell_type": "code", "execution_count": null, "id": "555bca84-36ef-4617-8c9b-e4dfe5e60c6d", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "5a13859c1a643da8ad1010eecc1ba65f", "grade": false, "grade_id": "cell-51a05e07f7705eac", "locked": true, "schema_version": 3, "solution": false, "task": false }, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "def newton_method(g, Jg, x0, tol, nmax):\n", " x = [] # list of iterates\n", " x.append(x0)\n", " r = [] # list of increments\n", " r.append(tol + 1)\n", " k = 0 # iteration counter\n", " while r[-1] > tol and k < nmax:\n", " dx = np.linalg.solve(Jg(x[-1]), -g(x[-1]))\n", " x.append(x[-1] + dx)\n", " r.append(np.linalg.norm(x[-1] - x[-2]))\n", " k = k + 1\n", " return x[-1]" ] }, { "cell_type": "markdown", "id": "2411d565-068b-42a5-aba7-7718cfd214b4", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "e82b449df986ebe292f919254c877196", "grade": false, "grade_id": "cell-762ceb3167a23709", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 1:** Complete the function `implicit_euler` which implements the implicit Euler method to compute the approximations $u_i \\approx u(t_i)$ for $i = 0, 1, 2, \\dots, m$. That is, identify a function $g : \\mathbb{R}^n \\to \\mathbb{R}^n$ along with its Jacobian $Jg : \\mathbb{R}^n \\to \\mathbb{R}^{n \\times n}$ whose zero gives you the next step of the implicit Euler method.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "6a487b49-91f0-4529-b002-8c76aed71550", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "20764b6ce8ed54ea76581ac55ec8f94d", "grade": false, "grade_id": "cell-163468dc3546c755", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "def implicit_euler(f, Jf, u_0, tau, m):\n", " n = 1 if isinstance(u_0, (int, float, complex)) else len(u_0)\n", " h = tau / m\n", " u = []\n", " u.append(u_0)\n", " for i in range(m):\n", " # Function whose zero is the next step of the implicit Euler method\n", " def g(x):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", " # Jacobian of the function\n", " def Jg(x):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", " # Method which finds the zero of the function g \n", " u.append(newton_method(g, Jg, u[i], 1e-8, 100))\n", " return u" ] }, { "cell_type": "markdown", "id": "960f1696-831d-4d54-9061-e0e722f83b3c", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "14bc50ae770d19717f6d519d4ac190cf", "grade": false, "grade_id": "cell-110ec4ed8354e501", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "We further provide you with the implementation of the forward Euler method you have completed last week in `serie11.ipynb`." ] }, { "cell_type": "code", "execution_count": null, "id": "7a0539d8-859b-4b3e-a60e-489b4fce8d61", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "2c1cb70059b35b8da377462b18467be4", "grade": false, "grade_id": "cell-30eaef1d562f8d9a", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "outputs": [], "source": [ "def forward_euler(f, u_0, tau, m):\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", " return u" ] }, { "cell_type": "markdown", "id": "b0fa6bda-60a2-43c1-b8b7-e3fab06d315d", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "4f38f333980931587ed749535d2338fa", "grade": false, "grade_id": "cell-46a79c7c49f2da6c", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 2:** Plot the approximate solutions of the system of ordinary differential equations from the implicit Euler method for $m=200$ time steps. How well do the methods approximate the exact solution?\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "a85d6724-ff8b-4a2f-8553-f09e9ce82368", "metadata": { "deletable": false, "editable": true, "nbgrader": { "cell_type": "code", "checksum": "d568cffa984f8b1237c24066132a3171", "grade": false, "grade_id": "cell-38cdcd5fcc268677", "locked": false, "schema_version": 3, "solution": true, "task": false }, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "# Right-hand side of the ODE\n", "def f(t, u):\n", " # Returns the value of the right-hand side of the ODE\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", "\n", "# Jacobian of f\n", "def Jf(t, u):\n", " # Returns the value of the right-hand side of the ODE\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", "\n", "\n", "# Parameters\n", "tau = 20\n", "u_0 = np.array([1, 0])\n", "m = 200\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "30dea496-5a71-4534-96d3-a3b2f70ca919", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "4f9b8081bbee0e6a0752707c14c56f9a", "grade": false, "grade_id": "cell-8cc76854f83abbce", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 3:** Complete the function `crank_nicolson` which implements the Crank-Nicolson method to compute the approximations $u_i \\approx u(t_i)$ for $i = 0, 1, 2, \\dots, m$. That is, identify a function $g : \\mathbb{R}^n \\to \\mathbb{R}^n$ along with its Jacobian $Jg : \\mathbb{R}^{n} \\to \\mathbb{R}^{n \\times n}$ whose zero gives you the next step of the Crank-Nicolson method\n", "
\n", "" ] }, { "cell_type": "code", "execution_count": null, "id": "cb457469-01a9-40cf-8808-d460283c8f09", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "6e30375a1bf739b1dc2ce3cb4ba7f1c7", "grade": false, "grade_id": "cell-659c6747be1d07e7", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "def crank_nicolson(f, Jf, u_0, tau, m):\n", " n = 1 if isinstance(u_0, (int, float, complex)) else len(u_0)\n", " h = tau / m\n", " u = []\n", " u.append(u_0)\n", " for i in range(m):\n", " # Function whose zero is the next step of the Crank-Nicolson method\n", " def g(x):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", " def Jg(x):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", " # Method which finds the zero of the function g \n", " u.append(newton_method(g, Jg, u[i], 1e-8, 100))\n", " \n", " return u" ] }, { "cell_type": "markdown", "id": "5ad9a9d4-640e-445d-8aec-ae572cd3ca5e", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "013f36e0134797e18d1db769f5b7b4e4", "grade": false, "grade_id": "cell-5f2d71cea5de2d3a", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 4:** As for the forward and implicit Euler method, plot the approximate solutions of the system of ordinary differential equations from the Crank-Nicolson method for $m=200$ time steps. How well do the methods approximate the exact solution?\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "6887381f-6e82-4fcd-acf7-5c748b5dc943", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "9c00e9eb66af6821d9da06aa49fc3bd2", "grade": false, "grade_id": "cell-536bae331b3b7a4a", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "35759787-aec8-43d2-bee0-38370b9ba304", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "c94c8604b8d20be3db7164b9df9a2ce6", "grade": false, "grade_id": "cell-739597547e26d8fd", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 5:** Visually determine the order of convergence of the implicit Euler and Crank-Nicolson methods by plotting the Euclidean error of the approximation $u_m$ from the exact value $u(\\tau)$ for $m = [100, 200, 400, 800, 1600]$.\n", "\n", "*Hint:* To compute the Euclidean norm of a vector `x`, use `np.linalg.norm(x)`. \n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "f5fbe281-e191-4867-874c-943bc30256de", "metadata": { "deletable": false, "editable": true, "nbgrader": { "cell_type": "code", "checksum": "7df49310f550e8a587a0ca67a64a8123", "grade": false, "grade_id": "cell-2505b9528e502b9f", "locked": false, "schema_version": 3, "solution": true, "task": false }, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "# Parameters\n", "tau = 20\n", "u_0 = np.array([1, 0])\n", "\n", "m_list = np.array([100, 200, 400, 800, 1600])\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "c879e37f-6292-449e-811d-1c919a53500a", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "abcdb24e537620448e1b0b73196e3621", "grade": false, "grade_id": "cell-e4bdb77687ae74a9", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "### Electrical circuit" ] }, { "cell_type": "markdown", "id": "04aae1c2-1351-4b42-8c88-c425e5ace787", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "4c25c75d7af47b335b998c1f8649bdd1", "grade": false, "grade_id": "cell-9013e0e3b764c3f8", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "We want to study an electrical circuit, which you can visualize by running the cell below." ] }, { "cell_type": "code", "execution_count": null, "id": "6ebc8831-1614-4fcb-ad81-0ee6c290e9a4", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "3ec85c6c0deff6bf80405d15df1f8721", "grade": false, "grade_id": "cell-e8223d2610a407c2", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "outputs": [], "source": [ "fig, ax = plt.subplots(); ax.axis(\"equal\"); ax.axis(\"off\"); ax.plot([0.4, 1, 1], [-1, -1, -0.1], linewidth=4, color=\"black\"); ax.plot([1, 1, 0.4], [0.1, 1, 1], linewidth=4, color=\"black\"); ax.plot([-0.4, -1, -1], [1, 1, 0.1], linewidth=4, color=\"black\"); ax.plot([-1, -1, -0.4], [-0.1, -1, -1], linewidth=4, color=\"black\"); ax.plot([0.8, 1.2], [-0.1, -0.1], linewidth=4, color=\"black\"); ax.plot([0.8, 1.2], [0.1, 0.1], linewidth=4, color=\"black\"); ax.text(0.65, 0.25, \"$C$\", fontsize=30, horizontalalignment=\"center\", verticalalignment=\"center\"); ax.plot([-0.85, -1.15], [-0.1, -0.1], linewidth=4, color=\"black\"); ax.plot([-0.7, -1.3], [0.1, 0.1], linewidth=4, color=\"black\"); ax.text(-0.7, 0.3, \"$+$\", fontsize=30, horizontalalignment=\"center\", verticalalignment=\"center\"); ax.text(-0.7, -0.2, \"$-$\", fontsize=30, horizontalalignment=\"center\", verticalalignment=\"center\"); ax.add_patch(mpl.patches.Arc((-0.3, 1), 0.2, 0.4, theta1=0, theta2=180, linewidth=4)); ax.add_patch(mpl.patches.Arc((-0.1, 1), 0.2, 0.4, theta1=0, theta2=180, linewidth=4)); ax.add_patch(mpl.patches.Arc((0.1, 1), 0.2, 0.4, theta1=0, theta2=180, linewidth=4)); ax.add_patch(mpl.patches.Arc((0.3, 1), 0.2, 0.4, theta1=0, theta2=180, linewidth=4)); ax.text(0, 0.75, \"$L$\", fontsize=30, horizontalalignment=\"center\", verticalalignment=\"center\"); ax.plot([-0.4, -0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35, 0.4], [-1, -0.85, -1.15, -0.85, -1.15, -0.85, -1.15, -0.85, -1.15, -1], linewidth=4, color=\"black\"); ax.text(0, -0.65, \"$R$\", fontsize=30, horizontalalignment=\"center\", verticalalignment=\"center\"); plt.show()" ] }, { "cell_type": "markdown", "id": "ad4efae2-851f-4a49-abe5-c9c99c8871a6", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "4caad121f4e253b52b6c80eea727d372", "grade": false, "grade_id": "cell-1d2e8df8d105092f", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "The circuit satisfies the following equation for the voltage $v$ around the capacitor:\n", "\n", "$$\n", "LC v''(t) + RC v'(t) + v(t) = F\n", "$$\n", "\n", "where $L$ is the inductance of the coil, $C$ is the capacity of the capacitor, $R$ the resistance\n", "and $F$ the motor force of the generator, which we assume to be constant in time. To apply the methods from the course, we will need to transform this second order equation into a system of first order equations." ] }, { "cell_type": "markdown", "id": "9bc60503-7e43-4321-834b-9b2178f2d11b", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "8ae4a21d3cb5b2204807721b8c75b5fc", "grade": false, "grade_id": "cell-0309da2cf2703165", "locked": true, "points": 0, "schema_version": 3, "solution": false, "task": true } }, "source": [ "
\n", " \n", "**Exercise 6 (Theoretical):** Introduce a new unknown $w$ such as $w(t)=v'(t)$. Use this relation to rewrite the equation above as a linear system of first order equations\n", "\n", "$$\n", "u'(t) = A {u}(t) + {b} = {f}(t, {u}(t))\n", "$$\n", "where ${u}(t) = ( w(t), v(t) )^\\top$ and $A$ is a $2 \\times 2$ matrix.\n", "
\n" ] }, { "cell_type": "markdown", "id": "03ed4af7-39f1-4989-8426-62af2f283d1c", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "1236dbbedb5bb22e93a672e2695408b8", "grade": false, "grade_id": "cell-ad0ae62c96f52c36", "locked": true, "points": 0, "schema_version": 3, "solution": false, "task": true }, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "
\n", " \n", "**Exercise 7 (Theoretical):** We want to appoximate the solution of this system using the forward Euler method. Give the general stability criterion for the forward Euler method applied to linear systems. Then, compute the condition on $h$ for the forward Euler method to be stable for the system at hand, for the values $L=0.01$, $C=10$ and $R=0.1$.\n", "\n", "*Hint:* Compute the eigenvalues $\\lambda_i(A)$ of the matrix $A$ found in the previous exercise. You can use Python to help you with this. The forward Euler method is said to be stable when \n", "$$\n", "|1 + h \\lambda_i(A) | < 1.\n", "$$\n", "
\n" ] }, { "cell_type": "markdown", "id": "40263770-63fa-44ce-9454-80e5fee12382", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "f53af6ad2ee9d392c7c97a3c80ef0def", "grade": false, "grade_id": "cell-8c9605039bcf22b8", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 8:** Consider the initial condition ${u}_0 = (0,1)^\\top$ and $F=0$. Using the forward Euler method (function `forward_euler` which we provide you above), compute an approximation of the solution of the system on the interval $[0,\\tau]$ with $\\tau=10$ and $m=43, 46,$ and $500$ sub-intervals, which corresponds to a time-step of $h = 0.233, 0.217,$ and $0.02$, respectively. Visualise the obtained approximations for $v(t)$ and comment on the obtained results.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "d0799479-d28e-426a-843e-3cbfc750176f", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "db90bb6982b91fe5fbf9baa378d702ba", "grade": false, "grade_id": "cell-c37258905c66ed30", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "def f(t, u, R=0.1, L=0.01, C=10, F=0):\n", " # Returns the value of the right-hand side of the ODE\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", "\n", "tau = 10\n", "u_0 = np.array([0, 1])\n", "\n", "m_list = [43, 46, 500]\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "740b14c0-69a3-4849-b561-1c103d492a4e", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "551d3713749ddeb0c0fa15d9727906c7", "grade": false, "grade_id": "cell-defd00e2e566bc91", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 9:** Repeat the previous exercise using the implicit Euler method (function `implicit_euler` which you have implemented above). Comment on the obtained results.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "81fd1ffd-751f-46f1-8c7d-29340f0338bd", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "e6eecd21f38dbffafc12a4994d9ffcb7", "grade": false, "grade_id": "cell-7df38cdbb93e9961", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "def Jf(t, u, R=0.1, L=0.01, C=10, F=0):\n", " # Returns the value of the right-hand side of the ODE\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", "\n", "tau = 10\n", "u_0 = np.array([0, 1])\n", "\n", "m_list = [43, 46, 500]\n", "for m in m_list:\n", " u_implicit_euler = implicit_euler(f, Jf, u_0, tau, m)\n", " t = np.linspace(0, tau, m + 1)\n", " plt.plot(t, np.asarray(u_implicit_euler)[:, 1], label=r\"$m = {:d}$\".format(m))\n", "\n", "plt.ylabel(r\"Implicit Euler approximation $v(t)$\")\n", "plt.xlabel(r\"time $t$\")\n", "plt.legend()\n", "plt.grid()\n", "plt.show()\n", "\n", "# Contrary to what was observed with the forward Euler method, we note that the\n", "# backward Euler method gives stable results for the three considered time-steps.\n" ] }, { "cell_type": "markdown", "id": "3065e471-f22d-42dd-a7dd-7d6156139e36", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "5d8b085a14bab8e8ce777995d1d753ec", "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 twelfth 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 }