Calc2-Lab4-ODEs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculus 2, Lab 4: Differential Equations\n",
    "\n",
    "Today we will be using the Numpy and SciPy Python libraries in order to numerical approximate solutions to _Ordinary Differential Equations_ using the Euler method. The session will build on the Jupyter notebooks which you have been using in class this week."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slope fields\n",
    "\n",
    "The following code was introduced in class to demonstrate how to use Python to generate a slope field for a differential equation of the form $y'=f(x,y)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dfield(f):\n",
    "    for j in x:\n",
    "        for k in y:\n",
    "            slope = f(k,j)\n",
    "            domain = np.linspace(j-0.05,j+0.05,2)\n",
    "            def fun(y1,x1):\n",
    "                z = slope*(domain-x1)+y1\n",
    "                return z\n",
    "            plt.plot(domain,fun(k,j),solid_capstyle='projecting',solid_joinstyle='bevel')\n",
    "\n",
    "    plt.title(\"Slope field for y'=f(x,y)\")\n",
    "    plt.grid(True)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 1\n",
    "\n",
    "Produce slope fields for the following differential equations:\n",
    "\n",
    "(a) $y' = y\\sin(x)$\n",
    "\n",
    "(b) $y' = x^2(1-y)$\n",
    "\n",
    "For each equation, you will need to:\n",
    "- define the function `f(y,x)` (remember that we have to write y,x and not x,y!)\n",
    "- define x and y using the `np.linspace` command\n",
    "- run the code above by putting `dfield(f)` as the last line in your code cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exact solutions\n",
    "\n",
    "Recall that for a separable differential equation of the form $y'=f(x)g(y)$, we can solve by rearranging this as $\\dfrac{1}{g(y)}\\dfrac{dy}{dx}=f(x)$ and then integrating both sides with respect to $x$.\n",
    "\n",
    "### Problem 2\n",
    "\n",
    "Solve the initial value problems below (the equations are the same as in problem 1). \n",
    "\n",
    "(a) $y'=y\\sin(x)$, with $y(0)=\\pi/2$.\n",
    "\n",
    "(b) $y'=x^2(1-y)$, with $y(0)=3$.\n",
    "\n",
    "You can either solve them by hand, or use SymPy. The SymPy code is as follows: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sy\n",
    "sy.init_printing()\n",
    "x = sy.symbols('x')\n",
    "y = sy.Function('y')\n",
    "\n",
    "yprime = sy.diff(y(x),x)\n",
    "eq = sy.Eq(yprime,x-y(x))\n",
    "eq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What this does:\n",
    "- defines `x` as a symbolic variable\n",
    "- defines `y` (symbolically) as a function of `x`\n",
    "- defines `yprime` as the derivative of `y` with respect to `x`\n",
    "- defines a differential equation \"`eq`\" of the form $y'=x-y$. (You should replace `x-y` with an approprirate function.)\n",
    "- displays the differential equation\n",
    "\n",
    "You can now ask the computer to solve the equation using the command `sy.dsolve(eq)`. If the output looks messy, you can try wrapping this command in `sy.simplify()`. To specify that you want to solve for $y$, you can use `sy.dsolve(eq,y(x))`, and to enter an intial condition $y(0)=1$, try `sy.dsolve(eq,y,ics={y(0):1})` as an optional argument (after the `y(x)` in `sy.dsolve(eq,y)`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Euler's Method\n",
    "\n",
    "- Essentially a numerical counterpart to slope fields.\n",
    "- Recall the linear approximation $$g(x_0+\\Delta x) \\approx g(x_0) + g'(x_0)\\Delta x.$$\n",
    "- Given $y'=f(x,y)$, look for solution $y=g(x)$ with initial condition $g(x_0)=y_0$.\n",
    "- Notice $g'(x_0) = f(x_0,y_0)$.\n",
    "- Proceed in \"steps\", with fixed **step size** $\\Delta x$ as follows:\n",
    "\n",
    "1. $(x_0,y_0)$ is given.\n",
    "2. $x_1 = x_0+\\Delta x$ and $y_1 = y_0+f(x_0,y_0)\\Delta x$\n",
    "3. $x_2 = x_0+2\\Delta x$ and $y_2 = y_1+f(x_1,y_1)\\Delta x$\n",
    "4. $x_3 = x_0+3\\Delta x$ and $y_3 = y_2+f(x_2,y_2)\\Delta x$\n",
    "\n",
    "etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code for Euler's Method\n",
    "\n",
    "The following is borrowed from [Patrick Walls at UBC](https://www.math.ubc.ca/~pwalls/math-python/differential-equations/first-order/). We define a function called odeEuler that implements the steps outlined above.\n",
    "\n",
    "Input is an ODE of the form $y'=f(y,x)$. Syntax is `odeEuler(f,y0,x)`, where:\n",
    "- `f` is our function\n",
    "- `y0` is the initial value for $y$, corresponding to the value in `x` at index 0.\n",
    "- `x` is a NumPy array (the values to be used for the steps)\n",
    "\n",
    "Before you can use it, you will need to define the function `f(y,x)`, and the values for `x` (for example, `x=np.linspace(0,2,100)` will use 100 samples for `x` between 0 and 2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def odeEuler(f,y0,x):\n",
    "    y = np.zeros(len(x))\n",
    "    y[0] = y0\n",
    "    for n in range(0,len(x)-1):\n",
    "        y[n+1] = y[n] + f(y[n],x[n])*(x[n+1] - x[n])\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output is a NumPy array containing the $y$ values correpsonding to the $x$ values in `x`.<br><br>\n",
    "\n",
    "If you want to plot these points, you can do the following. (You can copy this code into the empty cells below, and then edit as needed.)\n",
    "\n",
    "```\n",
    "def f(y,x):\n",
    "    return 2*x-3*y #you need to replace 2*x-3*y with your function here\n",
    "x = np.linspace(a,b,n) #replace a, b, and n with numbers\n",
    "y0 = 0 #replace 0 with the initial value for y here\n",
    "y = odeEuler(f,y0,x)\n",
    "plt.plot(x,y)\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 3\n",
    "\n",
    "Apply Euler's method to the two initial value problems from Problem 2. For each problem, you will need to define the function `f(y,x)` and the initial value `y0` before running the `odeEuler` and plot commands. For each differential equation, use $x$ values from 0 to 2 with 50 steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking the solution\n",
    "\n",
    "How well does our algorithm work? We can plot our values against the exact ones.<br>\n",
    "\n",
    "For each initial value problem done in problems 2 and 3: first define the function `f(y,x)` used by the `odeEuler` function, and then define a function `g(x)` that gives the exact solution. You can use your answers from Problem 2, but you will need to enter them as NumPy functions rather than SymPy functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You will need to change the definitions of f(y,x), y0 and g(x) below to the ones for part (a)\n",
    "# You may also need to define x as an np.linspace\n",
    "def f(y,x):\n",
    "    return 2*x-3*y\n",
    "def g(x):\n",
    "    return (2/3)*x-2/9+(2/9)*np.exp(-3*x)\n",
    "y0 = 0\n",
    "x = np.linspace(0,2,50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_exact = g(x)\n",
    "y_euler = odeEuler(f,y0,x)\n",
    "plt.plot(x,y_euler,'b.-',x,y_exact,'r-')\n",
    "plt.legend(['Euler','Exact'])\n",
    "plt.axis([0,2,0,2])\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If your graphs don't fit well in the window, you can adjust the $x$ and $y$ values shown by changing the numbers in the `plt.axis` line.\n",
    "\n",
    "Next, repeat the above for part (b)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter f(y,x) and g(x) for part (b) here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# copy,paste, and re-run the plotting code above here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional Challenge\n",
    "\n",
    "If you have time, see if you can produce a plot that shows the slope field, the exact solution, and the Euler's method approximation, all in the same plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}