{
"cells": [
{
"cell_type": "markdown",
"id": "56889fcd-f08a-4475-9928-af817748f5c8",
"metadata": {
"editable": true,
"nbgrader": {
"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": {
"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": {
"editable": true,
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"editable": true,
"nbgrader": {
"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": {
"editable": true,
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" return x - h * f((i + 1) * h, x) - u[i]\n",
" ### END SOLUTION\n",
" # Jacobian of the function\n",
" def Jg(x):\n",
" ### BEGIN SOLUTION\n",
" return np.eye(n) - h * Jf((i + 1) * h, x)\n",
" ### END SOLUTION\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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"editable": true,
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" return np.array([u[1], -u[0]])\n",
" ### END SOLUTION\n",
"\n",
"# Jacobian of f\n",
"def Jf(t, u):\n",
" # Returns the value of the right-hand side of the ODE\n",
" ### BEGIN SOLUTION\n",
" return np.array([[0, 1], [-1, 0]])\n",
" ### END SOLUTION\n",
"\n",
"\n",
"# Parameters\n",
"tau = 20\n",
"u_0 = np.array([1, 0])\n",
"m = 200\n",
"\n",
"### BEGIN SOLUTION\n",
"# Solve ODE with the implicit Euler method \n",
"u_forward_euler = forward_euler(f, u_0, tau, m)\n",
"u_implicit_euler = implicit_euler(f, Jf, u_0, tau, m)\n",
"\n",
"# Visualize the result\n",
"plot_solution(t, u_forward_euler)\n",
"plot_solution(t, u_implicit_euler)\n",
"\n",
"# Both methods diverge from the exact solution:\n",
"# The forward Euler method spirals outwards, while the Backward Euler\n",
"# spirals inward from the unit circle, on which the exact solution lies.\n",
"\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "30dea496-5a71-4534-96d3-a3b2f70ca919",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" return x - h / 2 * (f(i * h, u[i]) + f((i + 1) * h, x)) - u[i]\n",
" ### END SOLUTION\n",
" def Jg(x):\n",
" ### BEGIN SOLUTION\n",
" return np.eye(n) - h / 2 * (Jf(i * h, u[i]) + Jf((i + 1) * h, x))\n",
" ### END SOLUTION\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": {
"nbgrader": {
"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": {
"nbgrader": {
"grade": false,
"grade_id": "cell-536bae331b3b7a4a",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"u_crank_nicolson = crank_nicolson(f, Jf, u_0, tau, m)\n",
"plot_solution(t, u_crank_nicolson)\n",
"\n",
"# The approximate solution with the Crank-Nicolson method approximately\n",
"# coincides with the exact solution, as it stays on the unit circle.\n",
"\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "35759787-aec8-43d2-bee0-38370b9ba304",
"metadata": {
"nbgrader": {
"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": {
"editable": true,
"nbgrader": {
"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",
"### BEGIN SOLUTION \n",
"h_list = tau / m_list\n",
"error_implicit_euler_list = []\n",
"error_crank_nicolson_list = []\n",
"for m in m_list:\n",
" u_implicit_euler = implicit_euler(f, Jf, u_0, tau, m)\n",
" u_crank_nicolson = crank_nicolson(f, Jf, u_0, tau, m)\n",
" t = np.linspace(0, tau, m + 1)\n",
" u_exact = u(t, u_0)\n",
" error_implicit_euler_list.append(np.linalg.norm(u_exact - u_implicit_euler))\n",
" error_crank_nicolson_list.append(np.linalg.norm(u_exact - u_crank_nicolson))\n",
"\n",
"plt.loglog(h_list, error_implicit_euler_list, label=r\"Implicit Euler\")\n",
"plt.loglog(h_list, error_crank_nicolson_list, label=r\"Crank-Nicolson\")\n",
"plt.loglog(h_list, h_list, color=\"black\", linestyle=\"--\", label=r\"$h$\")\n",
"plt.loglog(h_list, h_list ** 2, color=\"black\", linestyle=\"dotted\", label=r\"$h^2$\")\n",
"plt.xlabel(r\"$h$\")\n",
"plt.grid()\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"\n",
"# The implicit Euler method is roughly parallel to the line (h, h), hence of order 1.\n",
"# The Crank-Nicolson error is roughly parallel to the line (h, h^2), hence of order 2.\n",
"\n",
"### END SOLUTION "
]
},
{
"cell_type": "markdown",
"id": "c879e37f-6292-449e-811d-1c919a53500a",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"Using the relation $w(t)=v'(t)$, we have $w'(t) = v''(t)$ and therefore the second order ordinary differential equation can be written as\n",
"\n",
"$$\n",
"LC w'(t) + RC w(t) + v(t) = F,\n",
"$$\n",
"\n",
"which leads to the system of first order equations\n",
"\n",
"$$\n",
" \\begin{cases}\n",
" w'(t) = \\frac{1}{LC}\\,F - \\frac{R}{L}\\,w(t) - \\frac{1}{LC}\\,v(t),\\\\[8pt]\n",
" v'(t) = w(t)\n",
" \\end{cases}\n",
"$$\n",
"\n",
"We can rewrite it in matrix form\n",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
" w(t) \\\\ v(t)\n",
"\\end{pmatrix}'\n",
"=\n",
"\\begin{pmatrix}\n",
" % -\\frac{R}{L} & -\\frac{1}{LC} \\\\\n",
" -R/L & -1/(LC) \\\\\n",
" 1 & 0 \n",
"\\end{pmatrix}\n",
"\\begin{pmatrix}\n",
" w(t) \\\\ v(t)\n",
"\\end{pmatrix}\n",
"+\n",
"\\begin{pmatrix}\n",
" F/(LC) \\\\ 0\n",
"\\end{pmatrix}\n",
"$$\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "03ed4af7-39f1-4989-8426-62af2f283d1c",
"metadata": {
"editable": true,
"nbgrader": {
"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",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"For the forward Euler method has the following stability criterion:\n",
"\n",
"$$\n",
"|1 + h \\lambda_i(A) | < 1\n",
"$$\n",
"where $\\lambda_i(A)$, $i=1,2$ are the eigenvalues of the matrix $A$. We notice here that this\n",
"criterion applies only if the eigenvalues of the matrix are real and negative.\n",
"With the given values, we have\n",
"\n",
"$$\n",
"A=\n",
"\\left(\\begin{array}{cc}\n",
" -10 & -10 \\\\\n",
" 1 & 0 \n",
" \\end{array} \\right),\n",
"$$\n",
"\n",
"whose eigenvalues are $\\lambda_1(A) = -8.873$ and $\\lambda_2(A) = -1.127$, as can be computed with the following code:\n",
"\n",
"```python\n",
"L = 0.01\n",
"C = 10\n",
"R = 0.1\n",
"F = 1.0\n",
"A = np.array([\n",
" [-R / L, -1 / (L * C)],\n",
" [1, 0]\n",
"])\n",
"\n",
"eigenvalues = np.linalg.eigvals(A)\n",
"print(\"Eigenvalues of A:\", eigenvalues)\n",
"```\n",
"\n",
"The stability criterion is therefore satisfied when $h > 0$ is chosen such that \n",
"\n",
"$$\n",
"|1 - 8.873 h | < 1 \\iff h < 2 / 8.873 = 0.225.\n",
"$$\n",
"\n",
"This can be computed with the following code snippet:\n",
"\n",
"```python\n",
"lambda_max = max(abs(eigenvalues))\n",
"print(\"Maximum magnitude of eigenvalues:\", lambda_max)\n",
"\n",
"dt_max = 2 / lambda_max\n",
"print(\"Maximum stable time step Δt:\", dt_max)\n",
"```\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "40263770-63fa-44ce-9454-80e5fee12382",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" A = np.array([[-R/L, -1/(L*C)], [1, 0]])\n",
" b = np.array([F/(L*C), 0])\n",
" return A @ u + b\n",
" ### END SOLUTION\n",
"\n",
"tau = 10\n",
"u_0 = np.array([0, 1])\n",
"\n",
"m_list = [43, 46, 500]\n",
"\n",
"### BEGIN SOLUTION\n",
"for m in m_list:\n",
" u_forward_euler = forward_euler(f, u_0, tau, m)\n",
" t = np.linspace(0, tau, m + 1)\n",
" plt.plot(t, np.asarray(u_forward_euler)[:, 1], label=r\"$m = {:d}$\".format(m))\n",
"\n",
"plt.ylabel(r\"Forward Euler approximation $v(t)$\")\n",
"plt.xlabel(r\"time $t$\")\n",
"plt.legend()\n",
"plt.grid()\n",
"plt.show()\n",
"\n",
"# We notice that the solution computed by the forward Euler method\n",
"# oscillates and does not approach 0 if we take m = 43 for which h = 0.233\n",
"# does not satisfy the stability condition. However, with m = 46, namely \n",
"# h = 0.217 < 0.225, we get a numerical solution which approaches 0\n",
"# (oscillating). Finally, if we take m = 500, namely h = 0.02 which is way\n",
"# smaller than the minimum necessary value, we find a good approximation.\n",
" \n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "740b14c0-69a3-4849-b561-1c103d492a4e",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" A = np.array([[-R/L, -1/(L*C)], [1, 0]])\n",
" return A\n",
" ### END SOLUTION\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",
"\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 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
}