1{
2 "cells": [
3  {
4   "cell_type": "markdown",
5   "id": "0276fde7-633a-4c42-87f5-522a9469c574",
6   "metadata": {},
7   "source": [
8    "# Computing definite integrals\n",
9    "\n",
10    "In this short tutorial, we will show how it is possible to use heyoka.py to compute definite integrals.\n",
11    "\n",
12    "Consider the integral\n",
13    "\n",
14    "$$\n",
15    "\\int_A^B f\\left( t \\right)\\,dt, \\tag{1}\n",
16    "$$\n",
17    "\n",
18    "where $f\\left( t \\right)$ is a differentiable function of $t$ and $A, B\\in\\mathbb{R}$. This integral can be seen as as the solution of the time-dependent differential equation\n",
19    "\n",
20    "$$\n",
21    "\\frac{dx}{dt} = f\\left( t \\right), \\tag{2}\n",
22    "$$\n",
23    "\n",
24    "where $x$ is a dummy state variable. Indeed, the integration of (2) by quadrature between $t=A$ and $t=B$ yields:\n",
25    "\n",
26    "$$\n",
27    "x\\left(B\\right) - x\\left(A\\right) = \\int_A^B f\\left(t\\right) \\, dt.\n",
28    "$$\n",
29    "\n",
30    "Note that we are always free to choose $x\\left( A \\right) = 0$, because the dynamics of $x$ in (2) does not depend on the value of $x$ itself. Thus, provided that we set up a numerical integration of (2) in which\n",
31    "\n",
32    "* we set $t=A$ as initial time coordinate,\n",
33    "* we set $x\\left( A \\right) = 0$ as initial condition,\n",
34    "\n",
35    "then the definite integral (1) is the value of $x$ at $t=B$.\n",
36    "\n",
37    "## Examples\n",
38    "\n",
39    "Let us start easy:"
40   ]
41  },
42  {
43   "cell_type": "code",
44   "execution_count": 1,
45   "id": "dfd57957-4b76-41e4-97bc-70e0e8dbf3f5",
46   "metadata": {},
47   "outputs": [
48    {
49     "data": {
50      "text/plain": [
51       "2.0"
52      ]
53     },
54     "execution_count": 1,
55     "metadata": {},
56     "output_type": "execute_result"
57    }
58   ],
59   "source": [
60    "import heyoka as hy, numpy as np\n",
61    "x, = hy.make_vars(\"x\")\n",
62    "\n",
63    "# Integrate sin(t) between 0 and pi.\n",
64    "ta = hy.taylor_adaptive([(x, hy.sin(hy.time))], [0])\n",
65    "ta.propagate_until(np.pi)\n",
66    "\n",
67    "# Print the result.\n",
68    "ta.state[0]"
69   ]
70  },
71  {
72   "cell_type": "markdown",
73   "id": "b7904907-2d5e-494c-bff8-074333b5c828",
74   "metadata": {},
75   "source": [
76    "As expected, $\\int_0^\\pi \\sin t\\, dt = 2$. Let's try to change the integration range:"
77   ]
78  },
79  {
80   "cell_type": "code",
81   "execution_count": 2,
82   "id": "8919b3d6-4fb5-4975-97ef-0f4bfae4b6c8",
83   "metadata": {},
84   "outputs": [
85    {
86     "data": {
87      "text/plain": [
88       "0.9564491424152821"
89      ]
90     },
91     "execution_count": 2,
92     "metadata": {},
93     "output_type": "execute_result"
94    }
95   ],
96   "source": [
97    "# Reset the state.\n",
98    "ta.state[0] = 0\n",
99    "\n",
100    "# New integration limits: from 1 to 2.\n",
101    "ta.time = 1\n",
102    "ta.propagate_until(2)\n",
103    "\n",
104    "ta.state[0]"
105   ]
106  },
107  {
108   "cell_type": "markdown",
109   "id": "bca411e5-b972-4bce-988a-8be03e7648e1",
110   "metadata": {},
111   "source": [
112    "Indeed, $\\int_1^2 \\sin t\\, dt = \\cos 1 - \\cos 2 = 0.9564491424152821\\ldots$.\n",
113    "\n",
114    "Let us try with a more complicated function:\n",
115    "\n",
116    "$$\n",
117    "\\int_\\sqrt{2}^\\sqrt{3} \\frac{\\sin \\left( \\cos t \\right)\\cdot \\operatorname{erf}{t}}{\\log\\left(\\sqrt{t}\\right)}\\,dt.\n",
118    "$$"
119   ]
120  },
121  {
122   "cell_type": "code",
123   "execution_count": 3,
124   "id": "48d63fa5-0c64-48df-9dca-fd577ec37c35",
125   "metadata": {},
126   "outputs": [
127    {
128     "data": {
129      "text/plain": [
130       "0.012382281847117892"
131      ]
132     },
133     "execution_count": 3,
134     "metadata": {},
135     "output_type": "execute_result"
136    }
137   ],
138   "source": [
139    "ta = hy.taylor_adaptive([(x, hy.sin(hy.cos(hy.time))*hy.erf(hy.time)/hy.log(hy.sqrt(hy.time)))],\n",
140    "                        [0], time = np.sqrt(2))\n",
141    "ta.propagate_until(np.sqrt(3))\n",
142    "ta.state[0]"
143   ]
144  },
145  {
146   "cell_type": "markdown",
147   "id": "23fceebd-ecb2-4b11-bffd-046835ae2163",
148   "metadata": {},
149   "source": [
150    "The result matches the value produced by [Wolfram Alpha](https://www.wolframalpha.com/input/?i=N%5BIntegrate%5BSin%5BCos%5Bx%5D%5D*Erf%5Bx%5D%2FLog%5BSqrt%5Bx%5D%5D%2C+%7Bx%2C+Sqrt%5B2%5D%2C+Sqrt%5B3%5D%7D%5D%2C+16%5D).\n",
151    "\n",
152    "Note that, since heyoka.py supports integration backwards in time, flipping around the integration limits also works as expected:"
153   ]
154  },
155  {
156   "cell_type": "code",
157   "execution_count": 4,
158   "id": "7a643902-1a07-4046-8266-0776d5f31b04",
159   "metadata": {},
160   "outputs": [
161    {
162     "data": {
163      "text/plain": [
164       "-0.012382281847117878"
165      ]
166     },
167     "execution_count": 4,
168     "metadata": {},
169     "output_type": "execute_result"
170    }
171   ],
172   "source": [
173    "ta.state[0] = 0\n",
174    "ta.time = np.sqrt(3)\n",
175    "ta.propagate_until(np.sqrt(2))\n",
176    "ta.state[0]"
177   ]
178  },
179  {
180   "cell_type": "markdown",
181   "id": "383ad99c-282d-4c08-b9ea-19156ac2c0ff",
182   "metadata": {},
183   "source": [
184    "Let us also perform the integration in extended precision:"
185   ]
186  },
187  {
188   "cell_type": "code",
189   "execution_count": 5,
190   "id": "34cf8f4b-3636-412f-86fe-f6c2549f8772",
191   "metadata": {},
192   "outputs": [
193    {
194     "data": {
195      "text/plain": [
196       "-0.012382281847117883605"
197      ]
198     },
199     "execution_count": 5,
200     "metadata": {},
201     "output_type": "execute_result"
202    }
203   ],
204   "source": [
205    "ta = hy.taylor_adaptive([(x, hy.sin(hy.cos(hy.time))*hy.erf(hy.time)/hy.log(hy.sqrt(hy.time)))],\n",
206    "                        [np.longdouble(0)], time = np.sqrt(np.longdouble(3)), fp_type=\"long double\")\n",
207    "ta.propagate_until(np.sqrt(np.longdouble(2)))\n",
208    "ta.state[0]"
209   ]
210  },
211  {
212   "cell_type": "markdown",
213   "id": "9dea8d4b-091d-421b-ac86-d92d4736ad46",
214   "metadata": {},
215   "source": [
216    "## Limitations and caveats\n",
217    "\n",
218    "This method for the computation of definite integrals inherits all the peculiarities and caveats of heyoka.py. For instance, the computation will fail if the derivative of $f\\left( t \\right)$ becomes infinite within the integration interval:"
219   ]
220  },
221  {
222   "cell_type": "code",
223   "execution_count": 6,
224   "id": "8f1c1684-03c4-431d-9f21-00989502b337",
225   "metadata": {},
226   "outputs": [
227    {
228     "data": {
229      "text/plain": [
230       "nan"
231      ]
232     },
233     "execution_count": 6,
234     "metadata": {},
235     "output_type": "execute_result"
236    }
237   ],
238   "source": [
239    "# Cannot compute the area of a semi-circle!\n",
240    "ta = hy.taylor_adaptive([(x, hy.sqrt(1 - hy.time**2))],\n",
241    "                        [0], time = -1)\n",
242    "ta.propagate_until(1)\n",
243    "ta.state[0]"
244   ]
245  }
246 ],
247 "metadata": {
248  "kernelspec": {
249   "display_name": "Python 3 (ipykernel)",
250   "language": "python",
251   "name": "python3"
252  },
253  "language_info": {
254   "codemirror_mode": {
255    "name": "ipython",
256    "version": 3
257   },
258   "file_extension": ".py",
259   "mimetype": "text/x-python",
260   "name": "python",
261   "nbconvert_exporter": "python",
262   "pygments_lexer": "ipython3",
263   "version": "3.8.10"
264  }
265 },
266 "nbformat": 4,
267 "nbformat_minor": 5
268}
269