Stablity_and_Root_Finding.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stablity and Root Finding\n",
    "__Author__ : Mohammad Rouintan , 400222042\n",
    "\n",
    "__Course__ : Undergraduate Numerical Analysis Course"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sympy import *"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1\n",
    "Prove by induction that first algorithm creates a sequence that second algorithm creates:"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First Algorithm:\n",
    "$$\n",
    "   f(n) = \n",
    "    \\begin{cases}\n",
    "      x_{0} = 1,\\; x_{1} = \\frac{1}{3}\\\\\n",
    "      x_{n + 1} = \\frac{13}{3}x_{n} - \\frac{4}{3}x_{n - 1}\\\\\n",
    "    \\end{cases}\n",
    "$$\n",
    "\n",
    "Second Algorithm:\n",
    "$$\n",
    "g(n) = x_{n} = (\\frac{1}{3})^{n}\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Basis__ :\n",
    "$$\n",
    "    \\begin{cases}\n",
    "      f(0) = x_{0} = 1\\\\\n",
    "      g(0) = x_{0} = (\\frac{1}{3})^{0} = 1\\\\\n",
    "    \\end{cases}\n",
    "    \\xrightarrow{} f(0) = g(0)\n",
    "\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Induction__ : We assume that $f(n) = g(n)$ for $n \\le k$. We must prove this equality for $n = k + 1$\n",
    "\n",
    "$$\n",
    "f(k + 1) = x_{k + 1} = \\frac{13}{3}x_{k} - \\frac{4}{3}x_{k - 1}\n",
    "$$\n",
    "$$\n",
    "= \\frac{13}{3}g(k) - \\frac{4}{3}g(k - 1)\n",
    "$$\n",
    "$$\n",
    "= \\frac{13}{3}(\\frac{1}{3})^{k} - \\frac{4}{3}(\\frac{1}{3})^{k - 1}\n",
    "$$\n",
    "$$\n",
    "= (\\frac{1}{3})^{k - 1}(\\frac{13}{3}(\\frac{1}{3}) - \\frac{4}{3})\n",
    "$$\n",
    "$$\n",
    "= (\\frac{1}{3})^{k - 1}(\\frac{13}{9} - \\frac{12}{9})\n",
    "$$\n",
    "$$\n",
    "= (\\frac{1}{3})^{k - 1}(\\frac{1}{9})\n",
    "$$\n",
    "$$\n",
    "= (\\frac{1}{3})^{k - 1}(\\frac{1}{3})^{2}\n",
    "$$\n",
    "$$\n",
    "= (\\frac{1}{3})^{k + 1} = g(k + 1)\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2\n",
    "Compute sequences for both algorithm for 15 terms with 7 float point precision."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def first_algorithm(n):\n",
    "    if n == 0:\n",
    "        return 1\n",
    "    elif n == 1:\n",
    "        return 1 / 3\n",
    "    \n",
    "    return (13 / 3) * first_algorithm(n - 1) - (4 / 3) * first_algorithm(n - 2)\n",
    "\n",
    "def second_algorithm(n):\n",
    "    return (1 / 3) ** n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First Algorithm\tSecond Algorithm\n",
      "1 1.0\n",
      "0.3333333333333333 0.3333333333333333\n",
      "0.11111111111111094 0.1111111111111111\n",
      "0.03703703703703626 0.03703703703703703\n",
      "0.012345679012342514 0.012345679012345677\n",
      "0.004115226337435884 0.004115226337448558\n",
      "0.0013717421124321456 0.0013717421124828527\n",
      "0.00045724737062478524 0.00045724737082761756\n",
      "0.00015241578946454185 0.0001524157902758725\n",
      "5.0805260179967644e-05 5.0805263425290837e-05\n",
      "1.6935074827137338e-05 1.693508780843028e-05\n",
      "5.644977344304949e-06 5.645029269476759e-06\n",
      "1.8814687224716613e-06 1.8816764231589195e-06\n",
      "6.263946716372672e-07 6.272254743863065e-07\n",
      "2.0575194713260943e-07 2.090751581287688e-07\n",
      "5.63988753916179e-08 6.969171937625627e-08\n"
     ]
    }
   ],
   "source": [
    "print('First Algorithm', 'Second Algorithm', sep='\\t')\n",
    "for i in range(16):\n",
    "    f_i = first_algorithm(i)\n",
    "    g_i = second_algorithm(i)\n",
    "    # print('{:.7f}'.format(f_i), '{:.7f}'.format(g_i), sep='\\t')\n",
    "    print(f_i, g_i)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3\n",
    "In this questions we want to compute $I_{n} = \\int_{0}^{1} x^{n}e^{x}\\, dx$. Write a recursive equation for computation with respect to part by part integration method. Compute all these integrals with your suggested recursive equation for $n = 0,1,...,15$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reqursive_equation(n):\n",
    "    if n == 0:\n",
    "        return np.exp(1) - 1\n",
    "        \n",
    "    return np.exp(1) - n * reqursive_equation(n - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.718282\n",
      "1.000000\n",
      "0.718282\n",
      "0.563436\n",
      "0.464536\n",
      "0.395600\n",
      "0.344685\n",
      "0.305490\n",
      "0.274362\n",
      "0.249028\n",
      "0.228002\n",
      "0.210265\n",
      "0.195100\n",
      "0.181983\n",
      "0.170519\n",
      "0.160496\n"
     ]
    }
   ],
   "source": [
    "for i in range(16):\n",
    "    ans = reqursive_equation(i)\n",
    "    print('{:.6f}'.format(ans))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4\n",
    "Solve following equation by Bisection, False Position, Newton and Improved Newton methods with\n",
    "threshold $|f(x_n)| < 10^{-16}$ with 20 precise floating point.(a and b are arbitrary) \n",
    "$$\n",
    "x + \\cos x = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def function(x):\n",
    "    return x + np.cos(x)\n",
    "\n",
    "def diffFunction(x):\n",
    "    return 1 - np.sin(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "a, b = -1.5, 0"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bisection Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def biSection(a, b):\n",
    "    if function(a) * function(b) >= 0:\n",
    "        print('There is not any root between a and b')\n",
    "        return\n",
    "\n",
    "    root = a\n",
    "    while np.abs(function(root)) >= 10**(-16):\n",
    "        root = (a + b) / 2\n",
    "        if function(root) == 0:\n",
    "            break\n",
    "\n",
    "        if function(root) * function(a) < 0:\n",
    "            b = root\n",
    "        else:\n",
    "            a = root\n",
    "\n",
    "    return root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Root = -0.73908513321516067229\n"
     ]
    }
   ],
   "source": [
    "print('Root =', '{:.20f}'.format(biSection(a, b)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### False Position Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def falsePosition(a, b):\n",
    "    if function(a) * function(b) >= 0:\n",
    "        print('There is not any root between a and b')\n",
    "        return\n",
    "\n",
    "    root = a\n",
    "    while np.abs(function(root)) >= 10**(-16):\n",
    "        root = (a * function(b) - b * function(a)) / (function(b) - function(a))\n",
    "        \n",
    "        if function(root) == 0:\n",
    "            break\n",
    "\n",
    "        if function(root) * function(a) < 0:\n",
    "            b = root\n",
    "        else:\n",
    "            a = root\n",
    "\n",
    "    return root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Root = -0.73908513321516067229\n"
     ]
    }
   ],
   "source": [
    "print('Root =', '{:.20f}'.format(falsePosition(a, b)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Newton Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newton(x0):\n",
    "    root = x0\n",
    "    while np.abs(function(root)) >= 10**(-16):\n",
    "        root = x0 - function(x0) / diffFunction(x0)\n",
    "        x0 = root\n",
    "\n",
    "    return root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Root = -0.73908513321516067229\n"
     ]
    }
   ],
   "source": [
    "print('Root =', '{:.20f}'.format(newton(-1)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Improved Newton Method (Secant Method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def secant(x0, x1):\n",
    "    root = x1\n",
    "    while np.abs(function(root)) >= 10**(-16):\n",
    "        root = x1 - function(x1) / ((function(x1) - function(x0)) / (x1 - x0))\n",
    "        x0 = x1\n",
    "        x1 = root\n",
    "        \n",
    "    return root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Root = -0.73908513321516067229\n"
     ]
    }
   ],
   "source": [
    "print('Root =', '{:.20f}'.format(secant(-1, -0.9)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 5\n",
    "Find roots of following equation:\n",
    "$$\n",
    "xe^{x} - 5 = 0\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We find roots of this equation with fixed point algorithm.\n",
    "$$\n",
    "x = ln(\\frac{5}{x})\\\\\n",
    "x_{0} = 2\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fixed_point(x):\n",
    "    if x == 0:\n",
    "        return 2\n",
    "    \n",
    "    return np.log(5 / fixed_point(x - 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "0.9162907318741551\n",
      "1.6968594842248554\n",
      "1.0806587320706331\n",
      "1.5318671201315661\n",
      "1.1829505810840526\n",
      "1.441426102541988\n",
      "1.243804939930723\n",
      "1.391262731110403\n",
      "1.2792261380252945\n",
      "1.3631825970022953\n",
      "1.2996158017062034\n",
      "1.3473692287953931\n",
      "1.311283940618623\n",
      "1.338431147870178\n",
      "1.3179397695471649\n",
      "1.333368175768879\n",
      "1.3217297084970905\n",
      "1.3304966484161829\n",
      "1.3238856200424909\n",
      "1.3288668483470925\n",
      "1.3251113270773656\n",
      "1.3279414360896864\n",
      "1.3258079616823177\n",
      "1.3274158564359095\n",
      "1.3262038252914221\n",
      "1.327117318085523\n",
      "1.3264287524872926\n",
      "1.3269477301838333\n",
      "1.326556547321885\n",
      "1.326851389782173\n",
      "1.326629152969372\n",
      "1.3267966588619082\n",
      "1.3266704025392237\n",
      "1.326765565825489\n",
      "1.326693837479721\n",
      "1.326747901503234\n",
      "1.326707151391614\n",
      "1.326737866147139\n",
      "1.3267147152958318\n",
      "1.3267321649020314\n",
      "1.3267190124953787\n",
      "1.326728925929449\n",
      "1.3267214538144496\n",
      "1.3267270858134925\n",
      "1.3267228407724005\n",
      "1.326726040411519\n",
      "1.326723628728263\n",
      "1.3267254465002631\n",
      "1.32672407638013\n",
      "1.326725109088777\n",
      "1.3267243306991603\n",
      "1.326724917399329\n",
      "1.3267244751823397\n",
      "1.3267248084971492\n",
      "1.3267245572658093\n",
      "1.3267247466278849\n",
      "1.3267246038988905\n",
      "1.3267247114788556\n",
      "1.3267246303919742\n",
      "1.3267246915100668\n",
      "1.3267246454431665\n",
      "1.3267246801654426\n",
      "1.326724653994017\n",
      "1.3267246737203597\n",
      "1.326724658851907\n",
      "1.3267246700587936\n",
      "1.326724661611761\n",
      "1.3267246679785925\n",
      "1.3267246631796832\n",
      "1.3267246667967934\n",
      "1.3267246640704475\n",
      "1.3267246661253924\n",
      "1.3267246645765065\n",
      "1.3267246657439575\n",
      "1.3267246648640076\n",
      "1.3267246655272575\n",
      "1.3267246650273423\n",
      "1.3267246654041465\n",
      "1.3267246651201356\n",
      "1.3267246653342046\n",
      "1.3267246651728533\n",
      "1.3267246652944695\n",
      "1.3267246652028029\n",
      "1.3267246652718954\n",
      "1.326724665219818\n",
      "1.3267246652590705\n",
      "1.3267246652294844\n",
      "1.3267246652517846\n",
      "1.3267246652349762\n",
      "1.3267246652476452\n",
      "1.3267246652380962\n",
      "1.3267246652452935\n",
      "1.3267246652398685\n",
      "1.3267246652439577\n",
      "1.3267246652408755\n",
      "1.3267246652431988\n",
      "1.3267246652414475\n",
      "1.3267246652427676\n",
      "1.3267246652417726\n",
      "1.3267246652425226\n",
      "1.326724665241957\n",
      "1.3267246652423834\n",
      "1.3267246652420621\n",
      "1.3267246652423044\n",
      "1.3267246652421216\n",
      "1.3267246652422595\n",
      "1.3267246652421556\n",
      "1.3267246652422338\n",
      "1.326724665242175\n",
      "1.3267246652422193\n",
      "1.3267246652421858\n",
      "1.326724665242211\n",
      "1.326724665242192\n",
      "1.3267246652422064\n",
      "1.3267246652421956\n",
      "1.3267246652422038\n",
      "1.3267246652421976\n",
      "1.3267246652422022\n",
      "1.3267246652421987\n",
      "1.3267246652422013\n",
      "1.3267246652421993\n",
      "1.326724665242201\n",
      "1.3267246652421998\n",
      "1.3267246652422005\n",
      "1.3267246652422\n",
      "1.3267246652422005\n",
      "1.3267246652422\n",
      "1.3267246652422005\n",
      "1.3267246652422\n"
     ]
    }
   ],
   "source": [
    "for i in range(130):\n",
    "    print(fixed_point(i))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion for this problem\n",
    "Write a conclusion and references which you've used in your homework"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}