newton_iter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before you turn this problem in, make sure everything runs as expected. First, **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart) and then **run all cells** (in the menubar, select Cell$\\rightarrow$Run All).\n",
    "\n",
    "Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\", as well as your github nickname below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NICKNAME = \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "faffb78b7d48e753358e7f2b080e58e1",
     "grade": false,
     "grade_id": "cell-d238cbae9896563e",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "87e98d582358e0470c1593a461300f5f",
     "grade": false,
     "grade_id": "cell-e24f93d1f5c3d8e8",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# Часть I. Метод Ньютона нахождения корней уравнения."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "461b52f42752997f98438e039802ad00",
     "grade": false,
     "grade_id": "cell-dcdac7e559daf63a",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Напишите алгоритм, который выполняет итерации Ньютона для заданной функции $f(x)$ с известной производной $f'(x)$. Ваша функция должна находить корни $f(x)$ с заданной точностью $\\epsilon$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "e97b768e6a3dc7e5496abd0004b542aa",
     "grade": false,
     "grade_id": "cell-e47469fedb8cc24a",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def newton_iteration(f, fder, x0, eps=1e-5, maxiter=100):\n",
    "    \"\"\"Нахождение корней $f(x) = 0$ через итерации Ньютона.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    f : callable\n",
    "        Функция, корни которой мы хотим найти.\n",
    "    fder : callable\n",
    "        Производная `f`.\n",
    "    x0 : float\n",
    "        Начальное приближение итераций Ньютона.\n",
    "    eps : float\n",
    "        Заданная точность. \n",
    "        Алгоритм прекращает работу когда расстояние между последовательными приближениями меньше `eps`.\n",
    "        По умолчанию 1e-5.\n",
    "    maxiter : int\n",
    "        Максимальное число итераций (по умолчанию 100).\n",
    "        Алгоритм прекращается, когда число итераций достигает `maxiter`.\n",
    "        Этот параметр нужен лишь для предотвращения бесконечного зацикливания.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    x : float\n",
    "        Найденное приближение к корню.\n",
    "    niter : int\n",
    "        Количество итераций.\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "62b80a064f6eb5c08d15643180c8d5df",
     "grade": true,
     "grade_id": "cell-0c23dfbb8447f5aa",
     "locked": true,
     "points": 2,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "from numpy.testing import assert_allclose\n",
    "\n",
    "xx, nit = newton_iteration(lambda x: x**2 - 1, lambda x: 2.*x, x0=4)\n",
    "assert_allclose(xx, 1.0, atol=1e-5)\n",
    "assert nit < 10\n",
    "\n",
    "xx, nit = newton_iteration(lambda x: x**2 - 1, lambda x: 2.*x, x0=-4)\n",
    "assert_allclose(xx, -1.0, atol=1e-5)\n",
    "assert nit < 10\n",
    "\n",
    "\n",
    "from math import log, exp\n",
    "\n",
    "xx, nit = newton_iteration(lambda x: exp(x)-2, lambda x: exp(x), x0=4, eps=1e-10)\n",
    "assert_allclose(xx, log(2), atol=1e-10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "7a6d66171d8f71ad126a9061605b3023",
     "grade": false,
     "grade_id": "cell-5df8b5b63ea89558",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Используя функцию `newton_iteration` для функции $f(x) = x^2 - 1$, постройте зависимость логарифма ошибки нахождения корня методом Ньютона от числа итераций. Какова сходимость метода (линейная, квадратичная, иная)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "8b9c677729824c82366781be99506a15",
     "grade": true,
     "grade_id": "cell-d429579b202c13d0",
     "locked": false,
     "points": 2,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "c4d8c23c1ec39e503217971bf3f2a5e8",
     "grade": false,
     "grade_id": "cell-f588bf8c5619986c",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# Часть II. Модифицированный метод Ньютона для кратных корней.\n",
    "\n",
    "Рассмотрим функцию, имеющую кратный корень (например, $f(x) = (x^2 - 1)^2$). Реализуйте модифицированную итерацию Ньютона, \n",
    "\n",
    "$$\n",
    "x_{n+1} = x_{n} - m \\frac{f(x_n)}{f'(x_n)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "585d03c9ede970d8b84d11324180b5ab",
     "grade": false,
     "grade_id": "cell-2f641d195eae18f0",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def mod_newton(f, fder, x0, m, eps=1e-5, maxiter=100):\n",
    "    \"\"\"Нахождение корней $f(x) = 0$ через итерации Ньютона.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    f : callable\n",
    "        Функция, корни которой мы хотим найти.\n",
    "    fder : callable\n",
    "        Производная `f`.\n",
    "    x0 : float\n",
    "        Начальное приближение итераций Ньютона.\n",
    "    eps : float\n",
    "        Заданная точность. \n",
    "        Алгоритм прекращает работу когда расстояние между последовательными приближениями меньше `eps`.\n",
    "        По умолчанию 1e-5.\n",
    "    maxiter : int\n",
    "        Максимальное число итераций (по умолчанию 100).\n",
    "        Алгоритм прекращается, когда число итераций достигает `maxiter`.\n",
    "        Этот параметр нужен лишь для преcдотвращения бесконечного зацикливания.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    x : float\n",
    "        Найденное приближение к корню.\n",
    "    niter : int\n",
    "        Количество итераций.\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "    return x, it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "f90f9328b8bb1dcef23530b126b54a79",
     "grade": true,
     "grade_id": "cell-d2c8d42721e39f15",
     "locked": true,
     "points": 2,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "for m in [1, 2, 3, 4, 5, 6]:\n",
    "    xx, nit = mod_newton(lambda x: (x**2 - 1)**4,\n",
    "                         lambda x: 4*(x**2 - 1)**3 * 2 * x,\n",
    "                         x0=2, m=m, maxiter=10000, eps=1e-9)\n",
    "    assert_allclose(xx, 1.0, atol=1e-8)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "e23163e7ecd9cd6453bcc7f992765a8c",
     "grade": false,
     "grade_id": "cell-6cfe913974c653c8",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Используйте модифицированную итерацию Ньютона для функции $f(x) = (x^2 - 1)^2$ со значениями  $m= 1, 2, 3, 4, 5, 6$. Сколько итераций требуется для достижения сходимости с фиксированным значением $\\epsilon$? Согласуются ли ваши наблюдения с предсказанием о квадратичной сходимости метода при значении $m$ равном кратности корня, и линейной сходимости для остальных значений $m$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "455ad6f64a12c655dbda83b423e7f70e",
     "grade": true,
     "grade_id": "cell-8e82e12fb5659e3b",
     "locked": false,
     "points": 2,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "b9ad26f1248014359154e8daf5387813",
     "grade": false,
     "grade_id": "cell-abdbc2a612344c2f",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# Part III. Newton's fractal.\n",
    "\n",
    "(4 дополнительных балла)\n",
    "\n",
    "Consider the equation\n",
    "$$ x^3 = 1 $$\n",
    "\n",
    "It has three solutions in the complex plane, $x_k = \\exp(i\\, 2\\pi k/ 3)$, $k = 0, 1, 2$.\n",
    "\n",
    "The Newton's iterations converge to one of these solutions, depending on the starting point in the complex plane (to converge to a complex-valued solution, the iteration needs a complex-valued starting point).\n",
    "\n",
    "Plot the **{basins of attraction}** of these roots on the complex plane of $x$ (i.e., on the plane $\\mathrm{Re}x$ -- $\\mathrm{Im}x$). To this end, make a series of calculations, varying the initial conditions on a grid of points. Color the grid in three colors, according to the root, to which iterations converged.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "822cf69d84ca6e84f8f0923be0fa47ec",
     "grade": true,
     "grade_id": "cell-dca407cdfc9774b4",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "025263dde967693b6468d3df7b9d36ab",
     "grade": false,
     "grade_id": "cell-afcddb4086fe5494",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# Часть IV. Метод Ньютона для системы нелинейных уравнений."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "e4f5f4089635cba3d6d371acc61386ac",
     "grade": false,
     "grade_id": "cell-43a22a68039359f3",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Напишите алгоритм, который выполняет итерации Ньютона для заданной системы нелинейных уравнений $\\vec{f}(x) = 0$ с известной матрицей Якоби $J(x)$. Ваша функция должна находить корни вектор-функции $\\vec{f}(x)$ с заданной точностью $\\epsilon$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "5c1866fe37d7d5cda77d61daa4c02c59",
     "grade": false,
     "grade_id": "cell-4580b8353412b2d3",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def newton_system(F, jcb, X0, eps=1e-5, maxiter=100):\n",
    "    \"\"\"Нахождение корней $f(x) = 0$ через итерации Ньютона.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    F : callable\n",
    "        Вектор-функция системы, которую мы хотим решить.\n",
    "    jcb : callable\n",
    "        Якобиан `f`.\n",
    "    X0 : array-like of floats, shape (n,)\n",
    "        Начальное приближение итераций Ньютона.\n",
    "    eps : float\n",
    "        Заданная точность. \n",
    "        Алгоритм прекращает работу когда расстояние между последовательными приближениями меньше `eps`.\n",
    "        По умолчанию 1e-5.\n",
    "    maxiter : int\n",
    "        Максимальное число итераций (по умолчанию 100).\n",
    "        Алгоритм прекращается, когда число итераций достигает `maxiter`.\n",
    "        Этот параметр нужен лишь для предотвращения бесконечного зацикливания.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    X : array-like of floats, shape (n,)\n",
    "        Найденное приближение к корню.\n",
    "    niter : int\n",
    "        Количество итераций.\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "480e4de969c51ac16e69d602e708c959",
     "grade": true,
     "grade_id": "cell-32a522217549fd0c",
     "locked": true,
     "points": 3,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "from numpy.testing import assert_allclose\n",
    "\n",
    "def func(X):\n",
    "    return np.array([X[0]**2 - 2*X[1]**4 + 1,\n",
    "                     X[0] - X[1]**3 + 1])\n",
    "\n",
    "def jac(X):\n",
    "    return np.array([[2.*X[0], -8.*X[1]**3],\n",
    "                     [1., -3.*X[1]**2]])\n",
    "\n",
    "xx, nit = newton_system(func, jac, np.array([3., 2.]))\n",
    "assert_allclose(xx, np.array([4.02879487, 1.71325223]), atol=1e-5)\n",
    "assert nit < 10\n",
    "\n",
    "xx, nit = newton_system(func, jac, np.asarray([1., 2.]))\n",
    "assert_allclose(xx, np.asarray([-0.35119943,  0.86570596]), atol=1e-5)\n",
    "assert nit < 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "4111148c7af39bf662bc113a2afeba56",
     "grade": false,
     "grade_id": "cell-d6571d426b977d3d",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Напишите функцию, возвращающую матрицу якобиана функции трех переменных `func2`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "26139d771a38d202a923a0ac9f1613b9",
     "grade": false,
     "grade_id": "cell-d091fffbc05825d7",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def func2(X):\n",
    "    x, y, z = X\n",
    "    return np.array([x**2 + y**2 + z**2 - 25,\n",
    "                     x*y + y*z + z*x - 5,\n",
    "                     x + y - 3])\n",
    "\n",
    "def jac2(X):\n",
    "    r\"\"\"Return the Jacobian of `func2(X)`.\n",
    "    \n",
    "    \n",
    "    The Jacobian matrix is defined as \n",
    "    \n",
    "    $$\n",
    "    J_{ij} = \\partial f_i / \\partial x_j\n",
    "    $$\n",
    "    \n",
    "    so that the first row contains the derivatives of $f_0$\n",
    "    with respect to the first, second etc coordinates; the second\n",
    "    row contains the derivatives of $f_1$ with respect to\n",
    "    the first, second etc coordinates; and so on.\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "16ef7b8539b67984802da101b1bdb64b",
     "grade": true,
     "grade_id": "cell-3a1ae706ab35a030",
     "locked": true,
     "points": 2,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "X = [1, 2, 3]\n",
    "assert jac2(X).shape == (3, 3)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "957a02fd59297e4f10bbca791b124183",
     "grade": true,
     "grade_id": "cell-7f2820844f3069af",
     "locked": true,
     "points": 2,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# more testing in this cell, keep it intact\n"
   ]
  }
 ],
 "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.8.13"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": true,
   "autocomplete": false,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}