{ "cells": [ { "cell_type": "markdown", "id": "0276fde7-633a-4c42-87f5-522a9469c574", "metadata": {}, "source": [ "# Computing definite integrals\n", "\n", "In this short tutorial, we will show how it is possible to use heyoka.py to compute definite integrals.\n", "\n", "Consider the integral\n", "\n", "$$\n", "\\int_A^B f\\left( t \\right)\\,dt, \\tag{1}\n", "$$\n", "\n", "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", "\n", "$$\n", "\\frac{dx}{dt} = f\\left( t \\right), \\tag{2}\n", "$$\n", "\n", "where $x$ is a dummy state variable. Indeed, the integration of (2) by quadrature between $t=A$ and $t=B$ yields:\n", "\n", "$$\n", "x\\left(B\\right) - x\\left(A\\right) = \\int_A^B f\\left(t\\right) \\, dt.\n", "$$\n", "\n", "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", "\n", "* we set $t=A$ as initial time coordinate,\n", "* we set $x\\left( A \\right) = 0$ as initial condition,\n", "\n", "then the definite integral (1) is the value of $x$ at $t=B$.\n", "\n", "## Examples\n", "\n", "Let us start easy:" ] }, { "cell_type": "code", "execution_count": 1, "id": "dfd57957-4b76-41e4-97bc-70e0e8dbf3f5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.0" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import heyoka as hy, numpy as np\n", "x, = hy.make_vars(\"x\")\n", "\n", "# Integrate sin(t) between 0 and pi.\n", "ta = hy.taylor_adaptive([(x, hy.sin(hy.time))], [0])\n", "ta.propagate_until(np.pi)\n", "\n", "# Print the result.\n", "ta.state[0]" ] }, { "cell_type": "markdown", "id": "b7904907-2d5e-494c-bff8-074333b5c828", "metadata": {}, "source": [ "As expected, $\\int_0^\\pi \\sin t\\, dt = 2$. Let's try to change the integration range:" ] }, { "cell_type": "code", "execution_count": 2, "id": "8919b3d6-4fb5-4975-97ef-0f4bfae4b6c8", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.9564491424152821" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Reset the state.\n", "ta.state[0] = 0\n", "\n", "# New integration limits: from 1 to 2.\n", "ta.time = 1\n", "ta.propagate_until(2)\n", "\n", "ta.state[0]" ] }, { "cell_type": "markdown", "id": "bca411e5-b972-4bce-988a-8be03e7648e1", "metadata": {}, "source": [ "Indeed, $\\int_1^2 \\sin t\\, dt = \\cos 1 - \\cos 2 = 0.9564491424152821\\ldots$.\n", "\n", "Let us try with a more complicated function:\n", "\n", "$$\n", "\\int_\\sqrt{2}^\\sqrt{3} \\frac{\\sin \\left( \\cos t \\right)\\cdot \\operatorname{erf}{t}}{\\log\\left(\\sqrt{t}\\right)}\\,dt.\n", "$$" ] }, { "cell_type": "code", "execution_count": 3, "id": "48d63fa5-0c64-48df-9dca-fd577ec37c35", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.012382281847117892" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ta = hy.taylor_adaptive([(x, hy.sin(hy.cos(hy.time))*hy.erf(hy.time)/hy.log(hy.sqrt(hy.time)))],\n", " [0], time = np.sqrt(2))\n", "ta.propagate_until(np.sqrt(3))\n", "ta.state[0]" ] }, { "cell_type": "markdown", "id": "23fceebd-ecb2-4b11-bffd-046835ae2163", "metadata": {}, "source": [ "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", "\n", "Note that, since heyoka.py supports integration backwards in time, flipping around the integration limits also works as expected:" ] }, { "cell_type": "code", "execution_count": 4, "id": "7a643902-1a07-4046-8266-0776d5f31b04", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.012382281847117878" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ta.state[0] = 0\n", "ta.time = np.sqrt(3)\n", "ta.propagate_until(np.sqrt(2))\n", "ta.state[0]" ] }, { "cell_type": "markdown", "id": "383ad99c-282d-4c08-b9ea-19156ac2c0ff", "metadata": {}, "source": [ "Let us also perform the integration in extended precision:" ] }, { "cell_type": "code", "execution_count": 5, "id": "34cf8f4b-3636-412f-86fe-f6c2549f8772", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.012382281847117883605" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ta = hy.taylor_adaptive([(x, hy.sin(hy.cos(hy.time))*hy.erf(hy.time)/hy.log(hy.sqrt(hy.time)))],\n", " [np.longdouble(0)], time = np.sqrt(np.longdouble(3)), fp_type=\"long double\")\n", "ta.propagate_until(np.sqrt(np.longdouble(2)))\n", "ta.state[0]" ] }, { "cell_type": "markdown", "id": "9dea8d4b-091d-421b-ac86-d92d4736ad46", "metadata": {}, "source": [ "## Limitations and caveats\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 6, "id": "8f1c1684-03c4-431d-9f21-00989502b337", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "nan" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Cannot compute the area of a semi-circle!\n", "ta = hy.taylor_adaptive([(x, hy.sqrt(1 - hy.time**2))],\n", " [0], time = -1)\n", "ta.propagate_until(1)\n", "ta.state[0]" ] } ], "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.10" } }, "nbformat": 4, "nbformat_minor": 5 }