{
 "cells": [
  {
   "cell_type": "raw",
   "id": "21753ddf",
   "metadata": {
    "vscode": {
     "languageId": "raw"
    }
   },
   "source": [
    "---\n",
    "title: \"Solution of Initial value problems\"\n",
    "format:\n",
    "  revealjs:\n",
    "    slide-level: 3\n",
    "    embed-resources: true\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fbc26bc",
   "metadata": {},
   "source": [
    "### What are *Initial Value Problems*?\n",
    "\n",
    "Considering an independent variable $t$, usually **time**,\n",
    "\n",
    "Given a system of differential equations that relate a number of functions $y(t)$ to their first derivatives\n",
    "\n",
    "$dy / dt = f(t, y)$\n",
    "\n",
    "and a set of initial values for $y(t)$ at time $t_0$\n",
    "\n",
    "$y(t_0) = y_0$\n",
    "\n",
    "we want to find the solutions $y(t)$.\n",
    "\n",
    "A numerical solution is an approximation of at several time points (ususally to make a plot of the solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eeeaebbb",
   "metadata": {},
   "source": [
    "### The importance of SODE\n",
    "\n",
    "- They appear in physics (because a = dv /dt = F/m and v = dx /dt)\n",
    "- They appear in chemistry because of rate laws of reactions and mass conservation\n",
    "- They appear everywhere in science, where time derivatives are known, but time courses need to be computed"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5d67be8",
   "metadata": {},
   "source": [
    "### Warm up example: mass bound to spring (no friction)\n",
    "\n",
    "Mass $m$ is bound to spring (bound to wall) contact with the ground is frictionless\n",
    "\n",
    "Spring has constant $k$, such that deviation $x$ from rest, produces force $F = - k x$\n",
    "\n",
    "Differencial equations:\n",
    "\n",
    "$ d x / dt = v$\n",
    "\n",
    "$ d v / dt = a = F / m = -k x / m$\n",
    "\n",
    "and let's pull the mass to position $x_0$ at time 0, although $v_0 = 0$\n",
    "\n",
    "Let's start with $m = 1$, $k = 1$ and $x_0 = 1$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d2420b7",
   "metadata": {},
   "source": [
    "### `solve_ivp()` function\n",
    "\n",
    "The `solve_ivp()` function finds numerical solutions of *Initial Value Problems*.\n",
    "\n",
    "It has the \"signature\"\n",
    "\n",
    "```\n",
    "solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)\n",
    "```\n",
    "\n",
    "Consult the [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) if needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43467d57",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the function\n",
    "from scipy.integrate import solve_ivp\n",
    "# and import numpy\n",
    "import numpy as np\n",
    "\n",
    "m = 1\n",
    "k = 1\n",
    "x0 = 1\n",
    "\n",
    "def spring_sode(t, y):\n",
    "    # position x is y[0] and velocity v is y[1]\n",
    "    dx = y[1]\n",
    "    dv = - k * y[0] / m\n",
    "    return np.array([dx, dv])\n",
    "\n",
    "sol = solve_ivp(spring_sode, (0, 10), (x0, 0))\n",
    "\n",
    "print(sol.t)\n",
    "print(sol.y[0])\n",
    "print(sol.y[1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ae6e3e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# don't ask...\n",
    "%config InlineBackend.figure_formats = ['svg']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dca0a194",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Obviously much better if we plot\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "f, ax = plt.subplots(1,1)\n",
    "\n",
    "ax.plot(sol.t, sol.y[0], label='$x$')\n",
    "ax.plot(sol.t, sol.y[1], label='$v$')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a25b7f6",
   "metadata": {},
   "source": [
    "Obviously ugly, especially because the time points are too few\n",
    "\n",
    "The times that we want the solution can be specified to `solve_ipv()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1045bb87",
   "metadata": {},
   "outputs": [],
   "source": [
    "times = np.linspace(0, 10, 200)\n",
    "sol = solve_ivp(spring_sode, (0, 10), (x0, 0), t_eval=times)\n",
    "\n",
    "# print(sol.t)\n",
    "# print(sol.y[0])\n",
    "# print(sol.y[1])\n",
    "f, ax = plt.subplots(1,1)\n",
    "\n",
    "ax.plot(sol.t, sol.y[0], label='$x$')\n",
    "ax.plot(sol.t, sol.y[1], label='$v$')\n",
    "ax.set(xlim=(0,10), xlabel=\"$t$\")\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09ab6110",
   "metadata": {},
   "source": [
    "### Mass bound to spring with friction\n",
    "\n",
    "What happens if there is friction, so that F has on more term:\n",
    "\n",
    "$F = -k x - \\mu v $\n",
    "\n",
    "Obtain the solution for $\\mu = 0.5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3dc7318c",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 1\n",
    "k = 1\n",
    "x0 = 1\n",
    "\n",
    "mu = 0.5\n",
    "\n",
    "def spring_sode_friction(t, y):\n",
    "    # position x is y[0] and velocity v is y[1]\n",
    "    dx = y[1]\n",
    "    dv = (- k * y[0] -mu * y[1] ) / m\n",
    "    return np.array([dx, dv])\n",
    "\n",
    "times = np.linspace(0, 20, 200)\n",
    "final_t = times[-1]\n",
    "sol = solve_ivp(spring_sode_friction, (0, final_t), (x0, 0), t_eval=times)\n",
    "\n",
    "f, ax = plt.subplots(1,1)\n",
    "\n",
    "ax.plot(sol.t, sol.y[0], label='$x$')\n",
    "ax.plot(sol.t, sol.y[1], label='$v$')\n",
    "ax.set(xlim=(0, final_t), xlabel=\"$t$\")\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dc00b3b",
   "metadata": {},
   "source": [
    "Exercise 1.\n",
    "\n",
    "For the frictionless case explore k = 0.2, 1, 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "090cca3d",
   "metadata": {},
   "source": [
    "Exercise 2.\n",
    "\n",
    "For the case with friction, explore k = 1 and $\\mu = 0.1, 0.5, 5$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a469c8f4",
   "metadata": {},
   "source": [
    "Exercise 3.\n",
    "\n",
    "In a closed, well-stirred vessel,\n",
    "\n",
    "$A \\rightarrow B$, first-order reaction with $k = 0.1$\n",
    "\n",
    "$B \\rightarrow A$, first-order reaction with $k = 0.2$\n",
    "\n",
    "$B \\rightarrow C$, first-order reaction with $k = 0.5$\n",
    "\n",
    "$C \\rightarrow B$, first-order reaction with $k = 0.1$\n",
    "\n",
    "At $t_0$ there is only $A_0 = 2$\n",
    "\n",
    "Compute the three chemical species time dependence $A(t)$ $B(t)$ $C(t)$ up to time 100."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65973d45",
   "metadata": {},
   "source": [
    "Exercise 4.\n",
    "\n",
    "What happens in the last system if there is a continuous inflow of A, with zero order\n",
    "\n",
    "$\\rightarrow A$, zero-order inflow with $k = 0.1$\n",
    "\n",
    "Compute the three chemical species time dependence $A(t)$ $B(t)$ $C(t)$ up to time 100."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PMN",
   "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.12.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}