{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "66d6dea5",
   "metadata": {},
   "source": [
    "# Interval Arithmetic"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c3f12d0",
   "metadata": {},
   "source": [
    "We illustrate interval arithmetic, available in ``mpmath``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3549addd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "import mpmath"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2045e4ae",
   "metadata": {},
   "source": [
    "### 1. a motivating example"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51c970de",
   "metadata": {},
   "source": [
    "The example is taken from the paper of Stefano Taschini:\n",
    "*Interval Arithmetic: Python Implementation and Applications*\n",
    "Proceedings of the 7th Python in Science Conference (SciPy 2008).\n",
    "\n",
    "While we get the final correct approximation when we use 36\n",
    "decimal places using mpmath of SymPy, the choice of 36 is not obvious."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c151715f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def exact_versus_numeric(A, B):\n",
    "    \"\"\"\n",
    "    Compares the exact to the numeric evaluation\n",
    "    of a expression at (A, B).\n",
    "    \"\"\"\n",
    "    F = lambda x, y: (333.75 - x**2)*y**6 \\\n",
    "        + x**2*(11*x**2*y**2 - 121*y**4 - 2) \\\n",
    "        + 5.5*y**8 + x/(2*y)\n",
    "    Z = F(float(A), float(B))\n",
    "    # to check the answer with SymPy :\n",
    "    X, Y = sp.var('x, y')\n",
    "    G = (sp.Rational(33375)/100 - X**2)*Y**6 \\\n",
    "        + X**2*(11*X**2*Y**2 - 121*Y**4 - 2) \\\n",
    "        + sp.Rational(55)/10*Y**8 \\\n",
    "        + sp.Rational(1)*X/(2*Y)\n",
    "    print('evaluating', G, 'at', (A, B))\n",
    "    E = sp.Subs(G, (X, Y), (A, B)).doit()\n",
    "    E15 = E.evalf(15)\n",
    "    print('numerical value :', Z)\n",
    "    print('exact value :', E, '~', E15)\n",
    "    print('error :', abs(E15 - Z))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e41d4db4",
   "metadata": {},
   "source": [
    "We compare the numerical evaluation to the exact evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f319b546",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "evaluating x**2*(11*x**2*y**2 - 121*y**4 - 2) + x/(2*y) + 11*y**8/2 + y**6*(1335/4 - x**2) at (77617, 33096)\n",
      "numerical value : 1.1726039400531787\n",
      "exact value : -54767/66192 ~ -0.827396059946821\n",
      "error : 2.00000000000000\n"
     ]
    }
   ],
   "source": [
    "(A, B) = (77617, 33096)\n",
    "exact_versus_numeric(A, B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "46d12f0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def variable_precision(A, B):\n",
    "    \"\"\"\n",
    "    The expression is evaluated with variable precision\n",
    "    at (A, B).\n",
    "    \"\"\"\n",
    "    # let us double the precision to 30\n",
    "    mpmath.mp.dps = 30\n",
    "    MP_F = lambda x, y: (mpmath.mpf('333.75') \\\n",
    "         - x**2)*y**6 \\\n",
    "         + x**2*(11*x**2*y**2 - 121*y**4 - 2) \\\n",
    "         + mpmath.mpf('5.5')*y**8 + x/(2*y)\n",
    "    MP_A30 = mpmath.mpf(str(A))\n",
    "    MP_B30 = mpmath.mpf(str(B))\n",
    "    Z30 = MP_F(MP_A30, MP_B30)\n",
    "    print('using 30 digits :', Z30)\n",
    "    # adjusting the precision to 35\n",
    "    mpmath.mp.dps = 35\n",
    "    MP_A35 = mpmath.mpf(str(A))\n",
    "    MP_B35 = mpmath.mpf(str(B))\n",
    "    Z35 = MP_F(MP_A35, MP_B35)\n",
    "    print('using 35 digits :', Z35)\n",
    "    # adjusting the precision to 36\n",
    "    mpmath.mp.dps = 36\n",
    "    MP_A36 = mpmath.mpf(str(A))\n",
    "    MP_B36 = mpmath.mpf(str(B))\n",
    "    Z36 = MP_F(MP_A36, MP_B36)\n",
    "    print('using 36 digits :', Z36)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3e0ef7c0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "using 30 digits : 1.17260394005317863185883490452\n",
      "using 35 digits : 1.1726039400531786318588349045201837\n",
      "using 36 digits : -0.827396059946821368141165095479816292\n"
     ]
    }
   ],
   "source": [
    "variable_precision(A, B)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84ce9aa4",
   "metadata": {},
   "source": [
    "We see a catastrophic change in the numbers, going from 35 to 36 digits."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de5e9263",
   "metadata": {},
   "source": [
    "### 2. illustration of using interval arithmetic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "72b2288e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpmath import iv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4fac1b45",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.0, 3.0] has type <class 'mpmath.ctx_iv.ivmpf'>\n",
      "[3.0, 3.0] / [3.0, 4.0] = [0.75, 1.0]\n",
      "Observe strings!\n",
      "[0.10000000000000000555, 0.10000000000000000555]\n",
      "[0.099999999999999991673, 0.10000000000000000555]\n",
      "some properties of [0.099999999999999991673, 0.10000000000000000555]\n",
      "middle : [0.10000000000000000555, 0.10000000000000000555]\n",
      "width : [1.3877787807814456755e-17, 1.3877787807814456755e-17]\n",
      "left bound : [0.099999999999999991673, 0.099999999999999991673]\n",
      "right bound : [0.10000000000000000555, 0.10000000000000000555]\n",
      "internal representation of [0.099999999999999991673, 0.10000000000000000555]\n",
      "{'_mpi_': ((0, 7205759403792793, -56, 53), (0, 3602879701896397, -55, 52))}\n",
      "fraction = 7205759403792793\n",
      "exponent = -56\n",
      "0.09999999999999999\n"
     ]
    }
   ],
   "source": [
    "iv.dps = 15\n",
    "X = iv.mpf(3)\n",
    "print(X, 'has type', type(X))\n",
    "Y = iv.mpf([3, 4])\n",
    "Z = X/Y\n",
    "print(X, '/', Y, '=', Z)\n",
    "# 0.1 /= '0.1'\n",
    "A = iv.mpf(0.1)\n",
    "B = iv.mpf('0.1')\n",
    "print('Observe strings!')\n",
    "print(A)\n",
    "print(B)\n",
    "print('some properties of', B)\n",
    "print('middle :', B.mid)\n",
    "print('width :', B.delta)\n",
    "print('left bound :', B.a)\n",
    "print('right bound :', B.b)\n",
    "print('internal representation of', B)\n",
    "print(B.__dict__)\n",
    "FRACTION = B.__dict__['_mpi_'][0][1]\n",
    "EXPONENT = B.__dict__['_mpi_'][0][2]\n",
    "print('fraction =', FRACTION)\n",
    "print('exponent =', EXPONENT)\n",
    "LB = FRACTION*2.0**EXPONENT\n",
    "print(LB)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11176cad",
   "metadata": {},
   "source": [
    "We can evaluate functions at intervals, as illustrated below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bfe976f0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "e : [2.7182818284590450908, 2.7182818284590455349]\n",
      "log(e) : [0.99999999999999988898, 1.000000000000000222]\n",
      "sin(e) : [0.41078129050290840274, 0.41078129050290884683]\n",
      "cos(e) : [-0.91173391478696530488, -0.91173391478696497181]\n",
      "pi : [3.141592653589793116, 3.1415926535897935601]\n",
      "sin(pi) : [-3.216245299353273201e-16, 1.2246467991473532072e-16]\n",
      "cos(pi) : [-1.0, -0.99999999999999988898]\n"
     ]
    }
   ],
   "source": [
    "iv.dps = 15\n",
    "E = iv.exp(1)\n",
    "print('e :', E)\n",
    "print('log(e) :', iv.log(E))\n",
    "print('sin(e) :', iv.sin(E))\n",
    "print('cos(e) :', iv.cos(E))\n",
    "P = iv.pi\n",
    "print('pi :', P)\n",
    "print('sin(pi) :', iv.sin(P))\n",
    "print('cos(pi) :', iv.cos(P))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e32bf7ef",
   "metadata": {},
   "source": [
    "### 3. the motivating example solved"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b5dfa34",
   "metadata": {},
   "source": [
    "We show how to obtain a correct numerical approximation for this example\n",
    "using the interval arithmetic in the mpmath of SymPy.\n",
    "The width of the interval gives us an upper bound for the error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a048cd28",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "evaluating x**2*(11*x**2*y**2 - 121*y**4 - 2) + x/(2*y) + 11*y**8/2 + y**6*(1335/4 - x**2) at (77617, 33096)\n",
      "exact value : -54767/66192 ~ -0.827396059946821\n",
      "using 35 decimal places ...\n",
      "[-6.827396059946821368141165095479816292382, 1.172603940053178631858834904520183709123]\n",
      "using 36 decimal places ...\n",
      "[-0.82739605994682136814116509547981629200549, -0.82739605994682136814116509547981629181741]\n"
     ]
    }
   ],
   "source": [
    "import sympy as sp\n",
    "X, Y = sp.var('x,y')\n",
    "G = (sp.Rational(33375)/100 - X**2)*Y**6 \\\n",
    "  + X**2*(11*X**2*Y**2 - 121*Y**4 - 2) \\\n",
    "  + sp.Rational(55)/10*Y**8 \\\n",
    "  + sp.Rational(1)*X/(2*Y)\n",
    "(A, B) = (77617, 33096)\n",
    "print('evaluating', G, 'at', (A, B))\n",
    "E = sp.Subs(G, (X, Y), (A, B)).doit()\n",
    "E15 = E.evalf(15)\n",
    "print('exact value :', E, '~', E15)\n",
    "# now we use interval arithmetic\n",
    "from mpmath import iv\n",
    "print('using 35 decimal places ...')\n",
    "iv.dps = 35\n",
    "IV_F = lambda x, y: (iv.mpf('333.75') \\\n",
    "- x**2)*y**6 \\\n",
    "+ x**2*(iv.mpf('11')*x**2*y**2 \\\n",
    "- iv.mpf('121')*y**4 - iv.mpf('2')) \\\n",
    "+ iv.mpf('5.5')*y**8 + x/(iv.mpf('2')*y)\n",
    "IV_A = iv.mpf(str(A))\n",
    "IV_B = iv.mpf(str(B))\n",
    "IV_Z = IV_F(IV_A, IV_B)\n",
    "print(IV_Z)\n",
    "print('using 36 decimal places ...')\n",
    "iv.dps = 36\n",
    "IV_F = lambda x, y: (iv.mpf('333.75') \\\n",
    "- x**2)*y**6 \\\n",
    "+ x**2*(iv.mpf('11')*x**2*y**2 \\\n",
    "- iv.mpf('121')*y**4 - iv.mpf('2')) \\\n",
    "+ iv.mpf('5.5')*y**8 + x/(iv.mpf('2')*y)\n",
    "IV_A = iv.mpf(str(A))\n",
    "IV_B = iv.mpf(str(B))\n",
    "IV_Z = IV_F(IV_A, IV_B)\n",
    "print(IV_Z)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef087e0c",
   "metadata": {},
   "source": [
    "### 4. Newton's method with interval arithmetic"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dcbcd2e",
   "metadata": {},
   "source": [
    "First we show a naive use of interval arithmetic in Newton's method.\n",
    "with the mpmath package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7256c132",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "naive interval Newton :\n",
      "[1.3107142857142854986, 1.514285714285714457]\n",
      "[1.1989198131568696848, 1.6218713507201252266]\n",
      "[0.93598868491625553112, 1.8564955830691347582]\n",
      "[0.16323583369742133975, 2.4568902244900931997]\n",
      "[-12.200194069435248423, 8.5013781837177671008]\n"
     ]
    }
   ],
   "source": [
    "from mpmath import iv\n",
    "iv.dps = 15\n",
    "print('naive interval Newton :')\n",
    "X = iv.mpf(['1.4', '1.5'])\n",
    "for i in range(5):\n",
    "    X = X - (X**2 - 2)/(2*X)\n",
    "    print(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bc2bb14",
   "metadata": {},
   "source": [
    "We see the interval getting wider, so the naive version of Newton's method does not work."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "457e77fd",
   "metadata": {},
   "source": [
    "Now we use interval arithmetic in the proper manner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d5ad52c4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "proper interval Newton :\n",
      "[1.4146551724137930162, 1.4146551724137932382]\n",
      "[1.4142136313013509152, 1.4142136313013513593]\n",
      "[1.4142135623730964777, 1.4142135623730969218]\n",
      "[1.4142135623730949234, 1.4142135623730951455]\n",
      "[1.4142135623730949234, 1.4142135623730951455]\n"
     ]
    }
   ],
   "source": [
    "print('proper interval Newton :')\n",
    "X = iv.mpf(['1.4', '1.5'])\n",
    "for i in range(5):\n",
    "    X = X.mid\n",
    "    X = X - (X**2 - 2)/(2*X)\n",
    "    print(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44df5442",
   "metadata": {},
   "source": [
    "So we see that Newton's method converges well."
   ]
  }
 ],
 "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}