{
"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 04- Linear systems and least squares "
]
},
{
"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": "5afb622e-c357-4a77-9356-e75ad0400743",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-44ae0a875bc8bfaa",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy as sp\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib\n",
"import timeit\n",
"import sys"
]
},
{
"cell_type": "markdown",
"id": "f5bf47d4-6137-42c1-9715-7c9a0fae632f",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-02fc1f88d5beaba2",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Convergence of bisection method\n",
"\n",
"Let us consider the function represented in the figure which appears when you run the cell below."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a4e832c8-dfb3-454e-b7e5-dfa1db4cad6f",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-830c563a25d5bbad",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"x_lin = np.linspace(-2, 2.5, 100)\n",
"f = lambda x: 1/3 - x/3 - x**2/3 + x**3/3\n",
"plt.plot(x_lin, f(x_lin))\n",
"plt.scatter([-1, 1], [0, 0], s=50, color=\"black\")\n",
"plt.text(-1, 0, r\"$\\alpha$\", fontsize=20, horizontalalignment=\"right\", verticalalignment=\"bottom\")\n",
"plt.text(1, 0, r\"$\\beta$\", fontsize=20, horizontalalignment=\"left\", verticalalignment=\"top\")\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "a8b916f4-d970-4482-8db5-f3f249bb3f56",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-840787d820d0fc38",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 1:** If we use the bisection method to find one of the roots of the function using $[-2, 2]$ as the initial interval and $10^{-5}$ as the tolerance, to which root does the method converge and in how many iterations? [Choose one answer]\n",
"\n",
"1. The method converges to $\\alpha$ in at most $18$ iterations.\n",
"2. The method converges to $\\beta$ in at most $18$ iterations.\n",
"3. The method converges to $\\alpha$ in at most $10$ iterations. \n",
"4. The method converges to $\\beta$ in at most $10$ iterations. \n",
"5. It is not possible to tell to which root the method will converge, but we will need at most $18$ iterations.\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"Answer 1 is the correct choice.\n",
"\n",
"The method cannot converge to $\\beta$, as the function doesn't change signs around $\\beta$. However, the function changes signs around $\\alpha$, therefore the method converges to $\\alpha$.\n",
"\n",
"To compute the number of iterations needed to guarantee that the bisecion method finds a root closer to the actual root than $tol$, we can use the formula\n",
"$$\n",
"k > \\log_2 \\left(\\frac{b-a}{tol}\\right)-1 \n",
"$$\n",
"to see that we need at most $18$ iterations. \n",
"\n",
"This formular comes from the fact, that the error of the method in the $k$-th iteration is bounded by *half* of the sub-interval size, i.e., the maximum distance of the sub-interval midpoint to the sub-interval bonudaries. Since the sub-interval length in iteration $k$ is\n",
"$$\n",
"\\left(\\frac{1}{2}\\right)^{k}(b-a),\n",
"$$\n",
"we need to ensure that $k$ satisfies \n",
"$$\n",
"\\frac{1}{2} \\left(\\frac{1}{2}\\right)^{k}(b-a) \\leq tol \\implies k > \\log_2 \\left(\\frac{b-a}{tol}\\right)-1.\n",
"$$\n",
"\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "5c382083-4412-48d0-9537-5b1166329c43",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-b106366fbd018302",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Bisection method for a problem in statics"
]
},
{
"cell_type": "markdown",
"id": "ed77abe1-0eb9-430e-b629-064500999c81",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-5085771c0d7d0f54",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"A certain planar mechanical system consists of four rigid rods of length $a_1$, $a_2$, $a_3$, and $a_4$ linked together. The configuration in which this rods are in with respect to each other is uniquely determined by the angles between the first and second rod, denoted with $\\theta$, and the first and fourth rod, denoted with $\\omega$. We provide you with a function which visualizes such a configuration. Run the below two cells to see the result."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "acafdd2e-cb50-4ee4-9929-c33c9b9d7efe",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-c936918b2618e06a",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"def plot_configuration(theta, omega, a1, a2, a3, a4):\n",
" fig, ax = plt.subplots()\n",
" ax.axis(\"equal\")\n",
" ax.axis(\"off\")\n",
" ax.plot([- a1 / 3, 4 * a1 / 3], [0, 0], linewidth=4, linestyle=\"--\", color=\"black\")\n",
" ax.plot([0, np.cos(theta) * a2, a1 + np.cos(omega) * a3, a1, 0], [0, np.sin(theta) * a2, np.sin(omega) * a3, 0, 0], linewidth=7, marker=\"o\", markersize=15, color=\"black\")\n",
" ax.text(a1 / 2, 0, \"$a_1$\", fontsize=20, horizontalalignment=\"center\", verticalalignment=\"center\", bbox=dict(facecolor='white', edgecolor='white', boxstyle='round,pad=0.1'))\n",
" ax.text(np.cos(theta) * a2 / 2, np.sin(theta) * a2 / 2, \"$a_2$\", fontsize=20, horizontalalignment=\"center\", verticalalignment=\"center\", bbox=dict(facecolor='white', edgecolor='white', boxstyle='round,pad=0.1'))\n",
" ax.text(a1 + np.cos(omega) * a3 / 2, np.sin(omega) * a3 / 2, \"$a_3$\", fontsize=20, horizontalalignment=\"center\", verticalalignment=\"center\", bbox=dict(facecolor='white', edgecolor='white', boxstyle='round,pad=0.1'))\n",
" ax.text((np.cos(theta) * a2 + a1 + np.cos(omega) * a3) / 2, (np.sin(theta) * a2 + np.sin(omega) * a3) / 2, \"$a_4$\", fontsize=20, horizontalalignment=\"center\", verticalalignment=\"center\", bbox=dict(facecolor='white', edgecolor='white', boxstyle='round,pad=0.1'))\n",
" ax.add_patch(matplotlib.patches.Arc((a1, 0), a1 / 3, a1 / 3, theta1=0, theta2=omega * 180 / np.pi, linewidth=4))\n",
" ax.text(6 * a1 / 5, 0, r\"$\\omega$\", fontsize=20, horizontalalignment=\"left\", verticalalignment=\"bottom\")\n",
" ax.add_patch(matplotlib.patches.Arc((0, 0), a1 / 3, a1 / 3, theta1=0, theta2=theta * 180 / np.pi, linewidth=4))\n",
" ax.text(a1 / 5, 0, r\"$\\theta$\", fontsize=20, horizontalalignment=\"left\", verticalalignment=\"bottom\")\n",
" plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "951ab36e-d6a9-4e0c-b2d5-83127d0eeccf",
"metadata": {},
"outputs": [],
"source": [
"theta = np.pi / 4\n",
"omega = np.pi / 2\n",
"a1 = 1\n",
"a2 = 1\n",
"a3 = 1 / np.sqrt(2)\n",
"a4 = 1 - 1 / np.sqrt(2)\n",
"plot_configuration(theta, omega, a1, a2, a3, a4)"
]
},
{
"cell_type": "markdown",
"id": "4aeed61d-7544-4b9f-9c32-254ed8f1f179",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-91918afab664d118",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Given the lengths of the rods $a_1$, $a_2$, $a_3$, $a_4$, and the angle $\\omega$, we would like to find out what value the angle $\\theta$ takes. In the lecture you have seen that every valid configuration satisfies the non-linear equation: \n",
"\n",
"$$\n",
"\\frac{a_1}{a_2} \\cos(\\omega) - \\frac{a_1}{a_4} \\cos(\\theta) - \\cos(\\omega - \\theta) + \\frac{a_1^2 + a_2^2 - a_3^2 + a_4^2}{2 a_2 a_4} = 0\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "3c2bd5fd-6a62-438f-a999-d3957029b045",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-72cca06cba72ce38",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 2:** Complete the function `constraint`, which returns the value of the left-hand side of the above equation for given $\\theta$, $\\omega$, $a_1$, $a_2$, $a_3$, and $a_4$.\n",
"\n",
"*Hint*: You can use `np.cos` to compute the cosine of a number.\n",
"
\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "411518a8-2e74-4de4-a5e4-c71dfc5c5802",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-a1549b554aabd0df",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def constraint(theta, omega, a1, a2, a3, a4):\n",
" \"\"\" constraint function \"\"\"\n",
"\n",
" ### BEGIN SOLUTION\n",
" return a1/a2 * np.cos(omega) - a1/a4 * np.cos(theta) - np.cos(omega - theta) + (a1**2 + a2**2 - a3**2 + a4**2) / (2*a2*a4)\n",
" ### END SOLUTION"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3bc2a82c-d863-4d63-884b-61611eea61c5",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-3b0338d318d428aa",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"assert not (constraint(np.pi / 2, np.pi / 2, 1, 1, 1, 1) is None), f\"the function 'constraint' returned nothing, make sure to 'return' your result\"; assert (C := constraint(np.pi / 2, np.pi / 2, 1, 1, 1, 1)) == 0.0, f\"'constraint(np.pi / 2, np.pi / 2, 1, 1, 1, 1)' should return 0 but got {C}\"; assert (C := constraint(np.pi / 2, np.pi / 2, 1, 1, 1, 1)) == 0.0, f\"'constraint(np.pi / 4, np.pi / 4, 1, 1, 1, 1)' should return 0 but got {C}\"; assert (C := constraint(np.pi / 2, 0, 2, 1, 1, 1)) == 4.5, f\"'constraint(np.pi / 2, 0, 1, 1, 1, 1)' should return 4.5 but got {C}\"; print(\"Great job! Your implementation passes our checks.\")"
]
},
{
"cell_type": "markdown",
"id": "320c53e3-b820-4502-b2b3-fb6b25bfc90a",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-d66bbc1091619ca0",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"As you have seen during the lecture, this non-linear equation cannot be solved analytically. Hence, we will need to solve this nonlinear equation analytically."
]
},
{
"cell_type": "markdown",
"id": "c35c51d9-790d-40a1-8d3a-62bd6ee65cd7",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-be8dabbf96f6a8e4",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 3:** In the below cell, complete the `bisection` function which takes as input a function $f:\\mathbb{R} \\to \\mathbb{R}$, the startpoint `a` and endpoint `b` of the search-interval $[a, b]$, and a tolerance `tol`, and uses the bisection method (Algorithm 1.2 in the lecture notes) to output the location of a zero `alpha` and the number of iterations `niter`.\n",
"\n",
"*Hint*: Additionally to Algorithm 1.2, our implementation also checks in every iteration if any of $f(a)$, $f(b)$, or $f(x^{(k)})$ are zero, in which case the root has already been found!\n",
"
\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fdc7b3c1-5aeb-4901-bc5c-7f7cce81ed5a",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-8c4f56d698407194",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def bisection(f, a, b, tol):\n",
" \"\"\" bisection method with stopping criterion \"\"\"\n",
"\n",
" alpha = a # approximate root\n",
" k_min = int(np.ceil(np.log2((b - a) / tol) - 1)) # compute number of iterations needed\n",
" x_k = (a + b) / 2\n",
" for k in range(k_min):\n",
" if f(x_k) == 0:\n",
" return x_k, k\n",
" elif f(a) == 0:\n",
" return a, k\n",
" elif f(b) == 0:\n",
" return b, k\n",
" \n",
" ### BEGIN SOLUTION\n",
" if f(x_k) * f(a) < 0:\n",
" b = x_k\n",
" else:\n",
" a = x_k\n",
" x_k = (a + b) / 2\n",
" ### END SOLUTION\n",
"\n",
" alpha = x_k\n",
" niter = k + 1\n",
" \n",
" return alpha, niter"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "78ee3826-3fae-49c8-8421-fe364003177a",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-ff847361f31cadb6",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"assert not (bisection(lambda x: x, -np.pi, 1, 0.1) is None), f\"the function 'bisection' returned nothing, make sure to 'return' your result\"; assert isinstance(bisection(lambda x: x, -np.pi, 1, 1e-10), tuple), f\"'bisection' should return 2 items, but only got 1\"; assert np.isclose(B := bisection(lambda x: x, -1e-10, 1, 1e-10)[0], 0), f\"'bisection(f, -1e-10, 1, 1e-10)' for 'f(x) = x' should return 'alpha = 0', but got {B}\"; assert np.isclose(B := bisection(lambda x: x, -1, 1e-10, 1e-10)[0], 0), f\"'bisection(f, -1, 1e-10, 1e-10)' for 'f(x) = x' should return 'alpha = 0', but got {B}\"; assert (B := bisection(lambda x: x, -1e-10, 1, 0.5)[1]) == 1, f\"'bisection(f, -1e-10, 1, 0.5)' for 'f(x) = x' should return 'niter = 1', but got {B}\"; print(\"Great job! It seems like your function works correctly.\")"
]
},
{
"cell_type": "markdown",
"id": "1b851792-bafe-4b73-87ca-e20ef975b5ab",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-2de6bc98f5eeca46",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Now, we would like to see how our method works in practice. We take a simple configuration where all rods are of length $1$ and $\\omega = \\pi / 4$. We will use the bisection method to look for a valid $\\theta \\in [0, \\pi]$ and plot the result below."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2c66a80f-5f47-4895-a23e-fa2ffbc08d69",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-0fa083ff54b68cc5",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"a1 = 1\n",
"a2 = 1\n",
"a3 = 1\n",
"a4 = 1\n",
"omega = np.pi / 4\n",
"\n",
"target = lambda theta: constraint(theta, omega, a1, a2, a3, a4)\n",
"theta, _ = bisection(target, 0, np.pi, 1e-10)\n",
"print(f\"found valid configuration with angle theta = {theta}\")\n",
"plot_configuration(theta, omega, a1, a2, a3, a4)"
]
},
{
"cell_type": "markdown",
"id": "decbb182-246a-4b38-ac9f-a578f14474e7",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-2dc462c53842052a",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Of course $\\theta = 0$ is a valid configuration, where the first and second, and the third and fourth rod overlap. However, if we suppose the rods cannot overlap, then we will need to exclude $\\theta = 0$ from our search interval. We will do this by slightly shrinking the interval to $[10^{-10}, \\pi]$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dcfc24be-c1f7-48cc-9ac6-8e2bd6670465",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-ba0f56fcdce721e5",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"a1 = 1\n",
"a2 = 1\n",
"a3 = 1\n",
"a4 = 1\n",
"omega = np.pi / 4\n",
"\n",
"target = lambda theta: constraint(theta, omega, a1, a2, a3, a4)\n",
"theta, _ = bisection(target, 1e-10, np.pi, 1e-10)\n",
"print(f\"found valid configuration with angle theta = {theta}\")\n",
"plot_configuration(theta, omega, a1, a2, a3, a4)"
]
},
{
"cell_type": "markdown",
"id": "aa8d4904-9632-42bf-9141-83d46a47b8dc",
"metadata": {},
"source": [
"And now a bit more interesting configuration. The angle $\\omega$ is now larger than $\\pi / 2$, while all rods are of length $1$, except for the first one which is half that lenght. Feel free to play around with the parameters a bit with this example. Be aware, that many combinations of parameters will be invalid, in which case often $\\theta = \\pi$ is returned."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "afc1adb8-79a9-462a-9d08-7af387fd3484",
"metadata": {},
"outputs": [],
"source": [
"a1 = 0.5\n",
"a2 = 1\n",
"a3 = 1\n",
"a4 = 1\n",
"omega = 0.75 * np.pi\n",
"\n",
"target = lambda theta: constraint(theta, omega, a1, a2, a3, a4)\n",
"theta, _ = bisection(target, 0, np.pi, 1e-10)\n",
"print(f\"found valid configuration with angle theta = {theta}\")\n",
"plot_configuration(theta, omega, a1, a2, a3, a4)"
]
},
{
"cell_type": "markdown",
"id": "d978e131-0ac6-45b8-9cd4-d0167075bd81",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-59a71f058e303216",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Only few configurations of parametrs are valid. For example, when the length of one rod is larger than the sum of the lengths of the other rods. Then we can see that the `constraint` function will never be zero. We can already see this by plotting it as a function of $\\theta$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8ed4021c-f3a6-4da5-b5ec-c801f957bc0b",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-988f67240ece3380",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"a1 = 1\n",
"a2 = 1\n",
"a3 = 1\n",
"a4 = 4\n",
"omega = np.pi / 4\n",
"x_lin = np.linspace(0, np.pi)\n",
"target = lambda theta: constraint(theta, omega, a1, a2, a3, a4)\n",
"plt.plot(x_lin, target(x_lin))\n",
"plt.xlabel(r\"$\\theta$\")\n",
"plt.ylabel(\"constraint function\")\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "8b8588c4-ab8b-4bc1-b9a4-5f689649dd81",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-abfcc407541f339f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Newton's method for systems of equations"
]
},
{
"cell_type": "markdown",
"id": "d786821e-8fa6-4a4a-8f1e-5625d07617ed",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-cae03da525c3bcaa",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 4:** Complete the following implementation of Newton's method for systems of equations (Algorithm 1.7 in the lecture notes).\n",
"\n",
"*Hint:* To solve an equation of the form $Ax = b$ for $x$, where $A \\in \\mathbb{R}^{n \\times n}$, $b \\in \\mathbb{R}^{n}$, and $x \\in \\mathbb{R}^{n}$, use `x = np.linalg.solve(A, b)`. To compute the norm $\\lVert x \\rVert$ of a vector $x$, use `np.linalg.norm(x)`.\n",
"
\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fd0af4a0-ae3f-4e22-a064-7907f90ce0c4",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-a226ef1dfd4e3222",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def newtonsys(F, JF, 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",
" ### BEGIN SOLUTION\n",
" dx = np.linalg.solve(JF(x[-1]), -F(x[-1]))\n",
" x.append(x[-1] + dx)\n",
" r.append(np.linalg.norm(x[-1] - x[-2]))\n",
" k = k + 1\n",
" ### END SOLUTION\n",
" return x, r, k"
]
},
{
"cell_type": "markdown",
"id": "401cb305-5c62-4821-80cc-300e1c737e52",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-ffb52f9122775c6a",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Let us consider the nonlinear system $\\mathbf{F}({\\bf x})=\\mathbf{0}$ with\n",
"\n",
"$$\n",
"\\mathbf{F}({\\bf x})\n",
"%F(x_1,x_2)\n",
"=\\left[\n",
"\\begin{array}{lcl}\n",
"e^{x_1^2+x_2^2}-\\alpha \\\\\n",
"e^{x_1^2-x_2^2}-1\n",
"\\end{array}\n",
"\\right],\n",
"$$\n",
"where the parameter $\\alpha$ takes the values $1$ or $e$."
]
},
{
"cell_type": "markdown",
"id": "ed372459-120d-45f9-9818-3dd0b56ddbb6",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-b8fdc0909564f89f",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 5 (Theoretical):** Compute the Jacobian matrix $J_\\mathbf{F}$ associated to the nonlinear system and the first Newton iteration in the case $\\alpha=e$ with ${\\bf x}^{(0)}=(1,1)^\\top$ as initial point.\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"We have\n",
"$$\n",
"J_\\mathbf{F}({\\bf x})=\n",
"\\left[\n",
"\\begin{array}{lc}\n",
"2 x_1 e^{x_1^2 + x_2^2} & 2 x_2 e^{x_1^2 + x_2^2} \\\\\n",
"2 x_1 e^{x_1^2 - x_2^2} & - 2 x_2 e^{x_1^2 - x_2^2}\n",
"\\end{array}\n",
"\\right].\n",
"$$\n",
"\n",
"The first Newton iteration ${\\bf x}^{(1)}$ is given by ${\\bf x}^{(1)}={\\bf x}^{(0)}+\\delta{\\bf x}$, where $\\delta{\\bf x}$ is the solution of\n",
"\n",
"$$\n",
"J_{\\bf F}({\\bf x}^{(0)})\\delta{\\bf x}=-\\mathbf{F}({\\bf x}^{(0)}).\n",
"$$\n",
"\n",
"For ${\\bf x}^{(0)}=(1,1)^\\top$ we have\n",
"\n",
"$${\\bf F}({\\bf x}^{(0)})=\\left[\n",
"\\begin{array}{c}\n",
"e(e-1) \\\\\n",
"0\n",
"\\end{array}\n",
"\\right]\n",
"\\quad \\text{and} \\quad\n",
"J_{\\bf F}({\\bf x}^{(0)})=\\left[\n",
"\\begin{array}{lc}\n",
"2e^2 & 2e^2 \\\\\n",
"2 & - 2\n",
"\\end{array}\n",
"\\right].\n",
"$$\n",
"\n",
"This means we have to solve the system\n",
"\n",
"$$\n",
"\\left[\n",
"\\begin{array}{lc}\n",
"2e^2 & 2e^2 \\\\\n",
"2 & - 2\n",
"\\end{array}\n",
"\\right]\\left[\n",
"\\begin{array}{c}\n",
"\\delta x_1 \\\\\n",
"\\delta x_2\n",
"\\end{array}\n",
"\\right]=\\left[\n",
"\\begin{array}{c}\n",
"e(1-e) \\\\\n",
"0\n",
"\\end{array}\n",
"\\right].\n",
"$$\n",
"\n",
"The second equation gives us $\\delta x_1=\\delta x_2$ and, replacing this result in the first one, we get $\\delta x_1=\\frac{1-e}{4e}$, namely $\\delta{\\bf x}=\\frac{1-e}{4e}{\\bf x}^{(0)}$. Finally, we find ${\\bf x}^{(1)}=(\\frac{1-e}{4e}+1,\\frac{1-e}{4e}+1)^\\top\\approx (0.84197,0.84197)^\\top$.\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "212d764a-b436-4935-8f6c-f81cf9ab8b44",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-196f6af8f39ef6f7",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
"\n",
"**Exercise 6:** Fix a tolerance of `tol = 1e-8` and a maximal number of iterations `nmax=1000`. Apply the vector Newton method for both cases $\\alpha=1$ and $\\alpha=e$, using as initial value ${\\bf x}_1^{(0)}=(1/10,1/10)^\\top$ and ${\\bf x}_e^{(0)}=(1,1)^\\top$, respectively. How do the result and the number of iterations change? What is the order of convergence of the method in each case?\n",
"\n",
"*Hint:* Try to work as much with NumPy arrays as possible, since you can do mathematical operations with them (`+`, `-`, ...), which are not possible between lists. To convert a Python list `l = [1, 2, 3]` to a NumPy array, use `np.array(l)`.\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "103268a9-14c3-4c2e-832e-8c883b7b5f60",
"metadata": {},
"source": [
"\n",
"\n",
"**Note:** Choosing the initial point can strongly influence if, and to which zero, Newton's method converges. In general, the closer to the zero the initial point is chosen, the better the method will converge to that zero. For one-dimentional functions with one input, we can often visually determine a good starting point. Once we move to two or more inputs, we may need to use a few iterations of another method (bisection, ...) to find a suitable stating point.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4c259442-eb57-4301-bfe9-eb327f21dad3",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-0c27b83ac14a5ea4",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def F(x, alpha):\n",
" ### BEGIN SOLUTION\n",
" return np.array([np.exp(x[0] ** 2 + x[1] ** 2) - alpha, np.exp(x[0] ** 2 - x[1] ** 2) - 1])\n",
" ### END SOLUTION\n",
"\n",
"def JF(x):\n",
" ### BEGIN SOLUTION\n",
" return np.array([[2 * x[0] * np.exp(x[0] ** 2 + x[1] ** 2),\n",
" 2 * x[1] * np.exp(x[0] ** 2 + x[1] ** 2)],\n",
" [2 * x[0] * np.exp(x[0] ** 2 - x[1] ** 2),\n",
" -2 * x[1] * np.exp(x[0] ** 2 - x[1] ** 2)]])\n",
" ### END SOLUTION\n",
"\n",
"# the function F when α = 1\n",
"def F_1(x):\n",
" return F(x, alpha=1)\n",
"\n",
"# the function F when α = e\n",
"def F_e(x):\n",
" return F(x, alpha=np.e)\n",
"\n",
"### BEGIN SOLUTION\n",
"x0_1 = [1/10, 1/10]\n",
"x0_e = [1, 1]\n",
"\n",
"tol = 1e-8\n",
"nmax = 1000\n",
"\n",
"x_1, r_1, niter_1 = newtonsys(F_1, JF, x0_1, tol, nmax)\n",
"x_e, r_e, niter_e = newtonsys(F_e, JF, x0_e, tol, nmax)\n",
"\n",
"print(f\"converged to x_1 = {x_1[-1]} after {niter_1} iterations\")\n",
"print(f\"converged to x_e = {x_e[-1]} after {niter_e} iterations\")\n",
"\n",
"plt.semilogy(range(niter_1), np.array(r_1[1:]) / np.array(r_1[:-1]), marker=\"o\", label=r\"$p=1$\")\n",
"plt.semilogy(range(niter_1), np.array(r_1[1:]) / np.array(r_1[:-1]) ** 2, marker=\"o\", label=r\"$p=2$\")\n",
"plt.ylabel(r\"increment ratio $\\|| x^{(k + 1)} - x^{(k)} \\|| / \\|| x^{(k)} - x^{(k - 1)} \\||^{p}$\")\n",
"plt.xlabel(r\"iteration $k$\")\n",
"plt.legend()\n",
"plt.grid()\n",
"plt.show()\n",
"plt.semilogy(range(niter_e), np.array(r_e[1:]) / np.array(r_e[:-1]), marker=\"o\", label=r\"$p=1$\")\n",
"plt.semilogy(range(niter_e), np.array(r_e[1:]) / np.array(r_e[:-1]) ** 2, marker=\"o\", label=r\"$p=2$\")\n",
"plt.ylabel(r\"increment ratio $\\|| x^{(k + 1)} - x^{(k)} \\|| / \\|| x^{(k)} - x^{(k - 1)} \\||^{p}$\")\n",
"plt.xlabel(r\"iteration $k$\")\n",
"plt.legend()\n",
"plt.grid()\n",
"plt.show()\n",
"\n",
"# We see that the Newton method converges with second order in the case α = e\n",
"# whereas in the case α = 1, it only converges with first order. This is due\n",
"# to the fact that det(Jf([0, 0])) = 0.\n",
"\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "2443aa45-0fca-4113-995c-70e9adc40d89",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-a26fae9678ff662c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"## The end\n",
"\n",
"Congratulations! You have made it to the end of the fourth 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
}