{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "9847dee6",
            "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": "d980ee1f",
            "cell_type": "markdown",
            "source": "In the first example, we solve Bratu's problem\n```\n\ny'' + k * exp(y) = 0\ny(0) = y(1) = 0\n\n```\nfor k = 1.\n\nWe rewrite the equation as a first-order system and implement its\nright-hand side evaluation\n```\n\ny1' = y2\ny2' = -exp(y1)\n\n```",
            "metadata": {}
        },
        {
            "id": "a2454bc3",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\ndef fun(x, y):\n    return np.vstack((y[1], -np.exp(y[0])))",
            "outputs": []
        },
        {
            "id": "7a171dbd",
            "cell_type": "markdown",
            "source": "Implement evaluation of the boundary condition residuals:\n",
            "metadata": {}
        },
        {
            "id": "d19a4f03",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def bc(ya, yb):\n    return np.array([ya[0], yb[0]])",
            "outputs": []
        },
        {
            "id": "923ed3f4",
            "cell_type": "markdown",
            "source": "Define the initial mesh with 5 nodes:\n",
            "metadata": {}
        },
        {
            "id": "201dc453",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "x = np.linspace(0, 1, 5)",
            "outputs": []
        },
        {
            "id": "0b9a0dca",
            "cell_type": "markdown",
            "source": "This problem is known to have two solutions. To obtain both of them, we\nuse two different initial guesses for y. We denote them by subscripts\na and b.\n",
            "metadata": {}
        },
        {
            "id": "c0d75434",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "y_a = np.zeros((2, x.size))\ny_b = np.zeros((2, x.size))\ny_b[0] = 3",
            "outputs": []
        },
        {
            "id": "b05fdb8a",
            "cell_type": "markdown",
            "source": "Now we are ready to run the solver.\n",
            "metadata": {}
        },
        {
            "id": "ce4e61eb",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.integrate import solve_bvp\nres_a = solve_bvp(fun, bc, x, y_a)\nres_b = solve_bvp(fun, bc, x, y_b)",
            "outputs": []
        },
        {
            "id": "91d3d224",
            "cell_type": "markdown",
            "source": "Let's plot the two found solutions. We take an advantage of having the\nsolution in a spline form to produce a smooth plot.\n",
            "metadata": {}
        },
        {
            "id": "f6da42fa",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "x_plot = np.linspace(0, 1, 100)\ny_plot_a = res_a.sol(x_plot)[0]\ny_plot_b = res_b.sol(x_plot)[0]\nimport matplotlib.pyplot as plt\nplt.plot(x_plot, y_plot_a, label='y_a')\nplt.plot(x_plot, y_plot_b, label='y_b')\nplt.legend()\nplt.xlabel(\"x\")\nplt.ylabel(\"y\")\nplt.show()",
            "outputs": []
        },
        {
            "id": "651ab942",
            "cell_type": "markdown",
            "source": "We see that the two solutions have similar shape, but differ in scale\nsignificantly.\n\nIn the second example, we solve a simple Sturm-Liouville problem\n```\n\ny'' + k**2 * y = 0\ny(0) = y(1) = 0\n\n```\nIt is known that a non-trivial solution y = A * sin(k * x) is possible for\nk = pi * n, where n is an integer. To establish the normalization constant\nA = 1 we add a boundary condition\n```\n\ny'(0) = k\n\n```\nAgain, we rewrite our equation as a first-order system and implement its\nright-hand side evaluation\n```\n\ny1' = y2\ny2' = -k**2 * y1\n\n```",
            "metadata": {}
        },
        {
            "id": "ff250a15",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def fun(x, y, p):\n    k = p[0]\n    return np.vstack((y[1], -k**2 * y[0]))",
            "outputs": []
        },
        {
            "id": "0cd6c87c",
            "cell_type": "markdown",
            "source": "Note that parameters p are passed as a vector (with one element in our\ncase).\n\nImplement the boundary conditions:\n",
            "metadata": {}
        },
        {
            "id": "e780a796",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def bc(ya, yb, p):\n    k = p[0]\n    return np.array([ya[0], yb[0], ya[1] - k])",
            "outputs": []
        },
        {
            "id": "c4c34047",
            "cell_type": "markdown",
            "source": "Set up the initial mesh and guess for y. We aim to find the solution for\nk = 2 * pi, to achieve that we set values of y to approximately follow\nsin(2 * pi * x):\n",
            "metadata": {}
        },
        {
            "id": "7d6f02f5",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "x = np.linspace(0, 1, 5)\ny = np.zeros((2, x.size))\ny[0, 1] = 1\ny[0, 3] = -1",
            "outputs": []
        },
        {
            "id": "7951e31f",
            "cell_type": "markdown",
            "source": "Run the solver with 6 as an initial guess for k.\n",
            "metadata": {}
        },
        {
            "id": "d1f278d2",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "sol = solve_bvp(fun, bc, x, y, p=[6])",
            "outputs": []
        },
        {
            "id": "5e530142",
            "cell_type": "markdown",
            "source": "We see that the found k is approximately correct:\n",
            "metadata": {}
        },
        {
            "id": "42911ef6",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "sol.p[0]",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "6.28329460046"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "a49c3557",
            "cell_type": "markdown",
            "source": "And, finally, plot the solution to see the anticipated sinusoid:\n",
            "metadata": {}
        },
        {
            "id": "e3954bc0",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "x_plot = np.linspace(0, 1, 100)\ny_plot = sol.sol(x_plot)[0]\nplt.plot(x_plot, y_plot)\nplt.xlabel(\"x\")\nplt.ylabel(\"y\")\nplt.show()",
            "outputs": []
        }
    ]
}