{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2e9478f5-3f70-4811-975f-efcbc2a6fbaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "using IntervalArithmetic;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a4b63c69-6c31-47c2-88a2-0713a91d875f",
   "metadata": {},
   "outputs": [],
   "source": [
    "setdisplay(:full);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "e04106da-5cf8-442d-b3d9-e6197761e985",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.884058535717245e28"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's try Rump's example of yesterday lecture -- please note that indices start from 1 and that the power is ^ and not **\n",
    "function f(x)\n",
    "    (333.75 - x[1]^2) * x[2]^6 + x[1]^2 * (11*x[1]^2*x[2]^2 - 121*x[2]^4 - 2) + 5.5*x[2]^8 + (x[1]/(2*x[2]))\n",
    "end\n",
    "# Floating-point evaluation\n",
    "f([77617,33096])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "521b21c6-4b11-4bd9-9c88-428305928408",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{Float64}(-3.541774862152234e21, 3.5417748621522344e21, com, false)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Interval evaluation, using Float64\n",
    "f([interval(77617),interval(33096)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "33c9dfa5-7a86-4be2-8c5a-d6e01bc49278",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{BigFloat}(-0.8273960599468213681411650954798162919990331157843848199178148416727096930142628, -0.8273960599468213681411650954798162919990331157843848199178148416727096930142455, com, false)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# And now, let's ask for \"BigFloat\" (and I do not know how big they are)\n",
    "f([interval(BigFloat,77617),interval(BigFloat,33096)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "13a6909b-5967-43cb-9511-50f5ea973ead",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Pfactorized (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now let's compare several evaluation schemes: first, the basis is (cheating) using the factored form of a polynomial\n",
    "function Pfactorized(x)\n",
    "    (x-1)*(x-2)*(x-3)*(x-4)*(x-5)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "eef595f3-7605-4895-9330-e697273a3aaf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Pfifth (generic function with 1 method)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# and its expanded form, followed by its successive derivatives\n",
    "function P(x)\n",
    "       x^5-15*x^4+85*x^3-225*x^2+274*x-120\n",
    "       end\n",
    "function Pprime(x)\n",
    "    5*x^4 - 60*x^3 + 255*x^2 - 450*x + 274\n",
    "end\n",
    "function Psecond(x)\n",
    "    20*x^3 - 180*x^2 + 510*x - 450\n",
    "end\n",
    "function Pthird(x)\n",
    "    60*x^2 - 360*x + 510\n",
    "end\n",
    "function Pfourth(x)\n",
    "    120*x - 360\n",
    "end\n",
    "function Pfifth(x)\n",
    "    120\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "76e35710-9e9e-418d-8b04-1b00b7768448",
   "metadata": {},
   "outputs": [],
   "source": [
    "# let us define some intervals: x=[0,6] and it is subdivided into 6 parts, y=[0,1], z=[-5,5]\n",
    "x=interval(0,6); x1=interval(0,1);x2=interval(1,2);x3=interval(2,3);x4=interval(3,4);x5=interval(4,5);x6=interval(5,6);y=interval(0,1);z=interval(-5,5);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "189e2857-1a6f-431d-afe3-d221ded0f543",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{Float64}(-1200.0, 1200.0, com, false)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let us compare the accuracy of the evaluation using the different forms: first, the factored one\n",
    "Pfactorized(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "663bcd0e-3a37-4839-bd30-9a40014db09a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{Float64}(-27660.0, 27660.0, com, false)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# now the naive evaluation of the whole interval\n",
    "P(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "806fc5b7-ff47-412b-9a3d-d25bbd89a5fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{Float64}(-12540.0, 12660.0, trv, false)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the union of the naive evaluation on each subinterval\n",
    "union_interval(P(x0P(x1),P(x2),P(x3),P(x4),P(x5),P(x6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "5f2d80b0-c916-4ecc-ac5c-0ca57e1c4a60",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Taylor1 (generic function with 1 method)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function Taylor1(x,f,fprime)\n",
    "    f(mid(x))+fprime(x)*(x-mid(x))\n",
    "    end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "7e1f8caf-c895-45b3-a5e5-7331e3c055ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{Float64}(-47802.0, 47802.0, com, false)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Taylor1(x,P,Pprime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "c442916e-8d64-4322-9a75-af39db7f5f70",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Taylor2 (generic function with 1 method)"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function Taylor2(x,f,fprime,fsecond)\n",
    "    f(mid(x))+fprime(mid(x))*(x-mid(x))+fsecond(x)/2*(x-mid(x))^2\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "3c92e49c-2a63-4a66-a342-35bda43e9360",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{Float64}(-31197.0, 31197.0, com, false)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Taylor2(x,P,Pprime,Psecond)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "72875c14-33ab-45a2-bf53-2bfd4396b1c7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Taylor3 (generic function with 1 method)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function Taylor3(x,f,fprime,fsecond,fthird)\n",
    "    f(mid(x))+fprime(mid(x))*(x-mid(x))+fsecond(mid(x))/2*(x-mid(x))^2+fthird(x)/6*(x-mid(x))^3\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "ba79897e-a27f-4ce5-9f03-010a66166263",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{Float64}(-12027.0, 12027.0, com, false)"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Taylor3(x,P,Pprime,Psecond,Pthird)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "8670570e-4847-4d47-8a61-fd1dec64dba7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Taylor4 (generic function with 1 method)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function Taylor4(x,f,fprime,fsecond,fthird,ffourth)\n",
    "    f(mid(x))+fprime(mid(x))*(x-mid(x))+fsecond(mid(x))/2*(x-mid(x))^2+fthird(mid(x))/6*(x-mid(x))^3+ffourth(x)/24*(x-mid(x))^4\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "8a63a05c-f60c-44b6-94ef-ff7acdebff61",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{Float64}(-1362.0, 1362.0, com, false)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Taylor4(x,P,Pprime,Psecond,Pthird,Pfourth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "e3d127a5-f3b9-4895-88e8-419db4fff2e7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Taylor5 (generic function with 1 method)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function Taylor5(x,f,fprime,fsecond,fthird,ffourth,ffifth)\n",
    "    f(mid(x))+fprime(mid(x))*(x-mid(x))+fsecond(mid(x))/2*(x-mid(x))^2+fthird(mid(x))/6*(x-mid(x))^3+ffourth(mid(x))/24*(x-mid(x))^4+ffifth(x)/120*(x-mid(x))^5\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "cd00cf35-c371-4291-b052-0def7d52ce4f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{Float64}(-390.0, 390.0, com, false)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Taylor5(x,P,Pprime,Psecond,Pthird,Pfourth,Pfifth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "50b4a929-af7a-47c7-b4ad-eef7980d5f37",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "hthird (generic function with 1 method)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# To show how dependent it is on the expression of the function, let's start with another function with complicated derivatives\n",
    "\n",
    "function h(x)\n",
    "    exp(x)*sin(x)/sqrt(x^2+1)\n",
    "end\n",
    "function hprime(x)\n",
    "    exp(x)*sin(x)/sqrt(x^2 + 1) + exp(x)*cos(x)/sqrt(x^2 + 1) - exp(x)*sin(x)*x/sqrt(x^2 + 1)^3\n",
    "end\n",
    "function hsecond(x)\n",
    "    2*exp(x)*cos(x)/sqrt(x^2 + 1) - 2*exp(x)*sin(x)*x/sqrt(x^2 + 1)^3 - 2*exp(x)*cos(x)*x/sqrt(x^2 + 1)^3 + 3*exp(x)*sin(x)*x^2/sqrt(x^2 + 1)^5 - exp(x)*sin(x)/sqrt(x^2 + 1)^3\n",
    "end\n",
    "function hthird(x)\n",
    "    2*exp(x)*cos(x)/sqrt(x^2 + 1) - 2*exp(x)*sin(x)/sqrt(x^2 + 1) - 6*exp(x)*cos(x)*x/sqrt(x^2 + 1)^3 + 9*exp(x)*sin(x)*x^2/sqrt(x^2 + 1)^5 - 3*exp(x)*sin(x)/sqrt(x^2 + 1)^3 + 9*exp(x)*cos(x)*x^2/sqrt(x^2 + 1)^5 - 3*exp(x)*cos(x)/sqrt(x^2 + 1)^3 - 15*exp(x)*sin(x)*x^3/sqrt(x^2 + 1)^7 + 9*exp(x)*sin(x)*x/sqrt(x^2 + 1)^5\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "81b55921-6da8-4394-a2fb-884aaa1f202d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{Float64}(-148.41315910257663, 148.41315910257663, com, false)"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "470ea918-1c71-4dbb-9b94-ba17b6b1dcba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{Float64}(-46.93236175050935, 14.202619362158785, trv, false)"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "union_interval(h(interval(-5,-3)),h(interval(-3,-1)),h(interval(-1,1)),h(interval(1,3)),h(interval(3,5)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "b788399d-86a5-435c-9c4e-0d40e95708bd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interval{Float64}(-35.995478306560166, 8.167787037086395, trv, false)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "union_interval(h(interval(-5,-4)),(interval(-4,-3)),h(interval(-3,-2)),h(interval(-2,-1)),h(interval(-1,0)),h(interval(0,1)),h(interval(1,2)),h(interval(2,3)),h(interval(3,4)),h(interval(4,5)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab1bfef6-ce30-464d-ab36-c387d17d97c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now it is your turn: compare with Taylor evaluations of order 1, 2 and 3 avec the whole interval z"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.12",
   "language": "julia",
   "name": "julia-1.12"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}