{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "a6d9f65f",
            "cell_type": "markdown",
            "source": "<div class='alert alert-warning'>\n\nSciPy's interactive examples with Jupyterlite are experimental and may not always work as expected. Execution of cells containing imports may result in large downloads (up to 60MB of content for the first import from SciPy). Load times when importing from SciPy may take roughly 10-20 seconds. If you notice any problems, feel free to open an [issue](https://github.com/scipy/scipy/issues/new/choose).\n\n</div>",
            "metadata": {}
        },
        {
            "id": "a90dba59",
            "cell_type": "markdown",
            "source": "Basic exponential decay showing automatically chosen time points.\n",
            "metadata": {}
        },
        {
            "id": "39541963",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy.integrate import solve_ivp\ndef exponential_decay(t, y): return -0.5 * y\nsol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8])\nprint(sol.t)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[ 0.          0.11487653  1.26364188  3.06061781  4.81611105  6.57445806\n  8.33328988 10.        ]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "a76f368c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(sol.y)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[[2.         1.88836035 1.06327177 0.43319312 0.18017253 0.07483045\n  0.03107158 0.01350781]\n [4.         3.7767207  2.12654355 0.86638624 0.36034507 0.14966091\n  0.06214316 0.02701561]\n [8.         7.5534414  4.25308709 1.73277247 0.72069014 0.29932181\n  0.12428631 0.05403123]]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "8ca08668",
            "cell_type": "markdown",
            "source": "Specifying points where the solution is desired.\n",
            "metadata": {}
        },
        {
            "id": "7a225ba6",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8],\n                t_eval=[0, 1, 2, 4, 10])\nprint(sol.t)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[ 0  1  2  4 10]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "f3a7dfcd",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(sol.y)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[[2.         1.21305369 0.73534021 0.27066736 0.01350938]\n [4.         2.42610739 1.47068043 0.54133472 0.02701876]\n [8.         4.85221478 2.94136085 1.08266944 0.05403753]]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "a0a27fb1",
            "cell_type": "markdown",
            "source": "Cannon fired upward with terminal event upon impact. The ``terminal`` and\n``direction`` fields of an event are applied by monkey patching a function.\nHere ``y[0]`` is position and ``y[1]`` is velocity. The projectile starts\nat position 0 with velocity +10. Note that the integration never reaches\nt=100 because the event is terminal.\n",
            "metadata": {}
        },
        {
            "id": "77855622",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def upward_cannon(t, y): return [y[1], -0.5]\ndef hit_ground(t, y): return y[0]\nhit_ground.terminal = True\nhit_ground.direction = -1\nsol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground)\nprint(sol.t_events)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[array([40.])]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "8b977434",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(sol.t)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02\n 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "99f90bd4",
            "cell_type": "markdown",
            "source": "Use `dense_output` and `events` to find position, which is 100, at the apex\nof the cannonball's trajectory. Apex is not defined as terminal, so both\napex and hit_ground are found. There is no information at t=20, so the sol\nattribute is used to evaluate the solution. The sol attribute is returned\nby setting ``dense_output=True``. Alternatively, the `y_events` attribute\ncan be used to access the solution at the time of the event.\n",
            "metadata": {}
        },
        {
            "id": "18948572",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def apex(t, y): return y[1]\nsol = solve_ivp(upward_cannon, [0, 100], [0, 10],\n                events=(hit_ground, apex), dense_output=True)\nprint(sol.t_events)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[array([40.]), array([20.])]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "be65027b",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(sol.t)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02\n 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "a0054b4a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(sol.sol(sol.t_events[1][0]))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[100.   0.]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "5829cdd1",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(sol.y_events)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[array([[-5.68434189e-14, -1.00000000e+01]]),\n array([[1.00000000e+02, 1.77635684e-15]])]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "c1c106b5",
            "cell_type": "markdown",
            "source": "As an example of a system with additional parameters, we'll implement\nthe Lotka-Volterra equations [12].\n",
            "metadata": {}
        },
        {
            "id": "0bd5f86c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def lotkavolterra(t, z, a, b, c, d):\n    x, y = z\n    return [a*x - b*x*y, -c*y + d*x*y]\n",
            "outputs": []
        },
        {
            "id": "169fafd7",
            "cell_type": "markdown",
            "source": "We pass in the parameter values a=1.5, b=1, c=3 and d=1 with the `args`\nargument.\n",
            "metadata": {}
        },
        {
            "id": "79740e2a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "sol = solve_ivp(lotkavolterra, [0, 15], [10, 5], args=(1.5, 1, 3, 1),\n                dense_output=True)",
            "outputs": []
        },
        {
            "id": "fd6c68a7",
            "cell_type": "markdown",
            "source": "Compute a dense solution and plot it.\n",
            "metadata": {}
        },
        {
            "id": "24fe22ce",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "t = np.linspace(0, 15, 300)\nz = sol.sol(t)\nimport matplotlib.pyplot as plt\nplt.plot(t, z.T)\nplt.xlabel('t')\nplt.legend(['x', 'y'], shadow=True)\nplt.title('Lotka-Volterra System')\nplt.show()",
            "outputs": []
        },
        {
            "id": "f0d17002",
            "cell_type": "markdown",
            "source": "A couple examples of using solve_ivp to solve the differential\nequation ``y' = Ay`` with complex matrix ``A``.\n",
            "metadata": {}
        },
        {
            "id": "b9ab258b",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "A = np.array([[-0.25 + 0.14j, 0, 0.33 + 0.44j],\n              [0.25 + 0.58j, -0.2 + 0.14j, 0],\n              [0, 0.2 + 0.4j, -0.1 + 0.97j]])",
            "outputs": []
        },
        {
            "id": "941df9a3",
            "cell_type": "markdown",
            "source": "Solving an IVP with ``A`` from above and ``y`` as 3x1 vector:\n",
            "metadata": {}
        },
        {
            "id": "3d28e91a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def deriv_vec(t, y):\n    return A @ y\nresult = solve_ivp(deriv_vec, [0, 25],\n                   np.array([10 + 0j, 20 + 0j, 30 + 0j]),\n                   t_eval=np.linspace(0, 25, 101))\nprint(result.y[:, 0])",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[10.+0.j 20.+0.j 30.+0.j]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "51d01eb2",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(result.y[:, -1])",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[18.46291039+45.25653651j 10.01569306+36.23293216j\n -4.98662741+80.07360388j]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "efdddbf7",
            "cell_type": "markdown",
            "source": "Solving an IVP with ``A`` from above with ``y`` as 3x3 matrix :\n",
            "metadata": {}
        },
        {
            "id": "2274d18c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def deriv_mat(t, y):\n    return (A @ y.reshape(3, 3)).flatten()\ny0 = np.array([[2 + 0j, 3 + 0j, 4 + 0j],\n               [5 + 0j, 6 + 0j, 7 + 0j],\n               [9 + 0j, 34 + 0j, 78 + 0j]])",
            "outputs": []
        },
        {
            "id": "6c8711db",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "result = solve_ivp(deriv_mat, [0, 25], y0.flatten(),\n                   t_eval=np.linspace(0, 25, 101))\nprint(result.y[:, 0].reshape(3, 3))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[[ 2.+0.j  3.+0.j  4.+0.j]\n [ 5.+0.j  6.+0.j  7.+0.j]\n [ 9.+0.j 34.+0.j 78.+0.j]]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "7e1ba461",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(result.y[:, -1].reshape(3, 3))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[[  5.67451179 +12.07938445j  17.2888073  +31.03278837j\n    37.83405768 +63.25138759j]\n [  3.39949503 +11.82123994j  21.32530996 +44.88668871j\n    53.17531184+103.80400411j]\n [ -2.26105874 +22.19277664j -15.1255713  +70.19616341j\n   -38.34616845+153.29039931j]]"
                    },
                    "execution_count": null
                }
            ]
        }
    ]
}