{ "cells": [ { "cell_type": "markdown", "id": "56889fcd-f08a-4475-9928-af817748f5c8", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "38095bee14f44cc9567ae387caffb751", "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": { "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": "5afb622e-c357-4a77-9356-e75ad0400743", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "896fc9aab44afbe84f4354fdd6801bb3", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "20871369e34903ca238e60207d2c19c9", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "b4a36921354fd0ced2972e010f4f6573", "grade": false, "grade_id": "cell-830c563a25d5bbad", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "outputs": [], "source": [ "x_lin = np.linspace(-2, 2, 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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "18f06ae84a22d4ad40b3f776fc383c1b", "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" ] }, { "cell_type": "markdown", "id": "5c382083-4412-48d0-9537-5b1166329c43", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "bc1b4351642790530b8f9f0109e093cb", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "0d570b193f5c62330bcf15ea61531b95", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "0cbbc2d7c02e3eac709c34cac4895324", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "ed5fa77a83fe5a057932d62f6e809e5f", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "9b65ea2ad5edf343085eb3c975c6fa12", "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": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "2fc261105c65d532865b9d9d069c0b84", "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", " # YOUR CODE HERE\n", " raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "id": "3bc2a82c-d863-4d63-884b-61611eea61c5", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "68793232d9d3dcc7957e7a82546a0a90", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "7a5be7bb60d2ae710033221813a99b34", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "b6a0ef39d35f2d5b51d9f5942205c9fe", "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": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "a39d22a76ccf5225818e285729f088df", "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", " # YOUR CODE HERE\n", " raise NotImplementedError()\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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "850ba6738546b263a179a661a872997e", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "1b6c857de9e88812caa6a21db4e601fd", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "c77d8ed80cac498361aac6f0ee53044a", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "723d70c91cd8338af7c7a989f46227ed", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "bb4671d2187d80d5617c11206509076c", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "5cd1397a963f380fb0daeb4a01e5e060", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "797d1a0d7da6046e575d5211ad1be239", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "fe1a4f0fd6b629b402f27114f66eef5a", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "6415e3f1133a8df80580244f114e8447", "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": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "2dee3d375b0e4c0d9b8427710d817e6e", "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", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", " return x, r, k" ] }, { "cell_type": "markdown", "id": "401cb305-5c62-4821-80cc-300e1c737e52", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "b36739fbb70ff06a2589e857b0c5656d", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "aa65ecf787814435ea24c09c9c878e0d", "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" ] }, { "cell_type": "markdown", "id": "212d764a-b436-4935-8f6c-f81cf9ab8b44", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "07f4c91299adcd913573c731be1cf84d", "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": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "4ab2ddb48d777c684f1de03057182b2c", "grade": false, "grade_id": "cell-0c27b83ac14a5ea4", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "def F(x, alpha):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", "\n", "def JF(x):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\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", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "2443aa45-0fca-4113-995c-70e9adc40d89", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "e1c1420917f8f2c9b00a9f315707fe44", "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 }