Lecture-14-Ordinary-Differential-Equations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 14:  Solutions to Ordinary Differential Equations and Viscoelasticity\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Background\n",
    "----\n",
    "\n",
    "What are differential equations?  Grant Sanderson from 3Blue1Brown gives us this conceptual description on the origin of differential equations, \"Differential equations arise whenever it is easier to describe change than absolute amounts.\"  The analysis of mechanics of moving bodies, vibrating strings, and heat flow were some of the areas where differential equations were first used.  In this lecture, viscoelasticity will be the topic used to demonstrate the solution to differential equations.  The invention of Calculus and use of the deriviative and integral naturally lead to differential equations.\n",
    "\n",
    "For us, an equation that relates an independent variable (such as time) to some quantities and the derivatives of those quantities with respect to the independent variable is an ordinary differental equation (ODE).  Using x as the independent variable, a general ordinary differential equation looks like this (see Wikipedia page on ODEs):\n",
    "$$\n",
    "y^{(n)}=\\sum _{i=0}^{n-1}a_{i}(x)y^{(i)}+r(x)\n",
    "$$\n",
    "where $y'$ are the derivatives in the independent variable $x$ and the $a$ could be functions of $x$.\n",
    "\n",
    "In physics, you have probably solved this differential equation:\n",
    "$$\n",
    "m{\\frac {\\mathrm {d} ^{2}x(t)}{\\mathrm {d} t^{2}}}=F(x(t))\n",
    "$$\n",
    "where t is now the independent variable.  This is a differential equation restatement of $F= ma$ where it is easier to relate the acceleration and velocity as a differential quantity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## What skills will I learn?\n",
    "\n",
    "* You will identify differential equations and practice solving them using Sympy and numerical methods.\n",
    "* Solutions to standard and more general viscoelastic problems will be developed.\n",
    "\n",
    "## What steps should I take?\n",
    "\n",
    "1. Review the solutions to ODEs from physics and electronics.\n",
    "1. Practice using sympy and python tools to solve these equations in an organized way.\n",
    "1. Review viscoelasticity concepts.\n",
    "1. Apply the techniquest of ODE solutions to viscoelastic problems.  Review the Maxwell model and solve the Kelvin model for yourself.\n",
    "1. Examine the solutions to two important ODEs that we will use when solving Fick's law.  \n",
    "\n",
    "## Reading and Reference\n",
    "\n",
    "* 3Blue1Brown's playlist on differential equations:  https://www.youtube.com/playlist?list=PLZHQObOWTQDNPOjrT6KVlfJuKtYTftqH6\n",
    "* https://en.wikipedia.org/wiki/Ordinary_differential_equation\n",
    "* https://en.wikipedia.org/wiki/Differential_equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1:  Review of Physics and Electronics ODEs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example using DC Circuits\n",
    "\n",
    "A simple DC circuit includes a switch (initially open), an inductor (L=2) and a resistor (R=6) in series.  The DC voltage is 100V.  What is the current in the circuit after the switch is closed?  We know that:\n",
    "\n",
    "$$ E_L + E_R = 100V $$\n",
    "\n",
    "and that \n",
    "\n",
    "$$ E_L = L \\frac{dI(t)}{dt} $$\n",
    "\n",
    "and \n",
    "\n",
    "$$ E_R = RI(t) $$  \n",
    "\n",
    "(You will see a strong correlation between these problems and viscoelastic problems later in this lecture.)  Typically, one would solve this problem by separation of variables.  In `sympy` there is a function called `dsolve()` that can help us with this.  Let us import `sympy` into the current namespace in a way that helps us keep track of what functions are coming from `sympy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import sympy as sp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# can also run quietly using:\n",
    "#sp.init_session(quiet=True)\n",
    "\n",
    "# set up some common symbols and report back to the user.\n",
    "sp.init_session()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create some additional symbols that are a bit easier to remember in this problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These will be functions of time, thus the classification.\n",
    "current, voltage = sp.symbols('current voltage', cls=Function)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assign to the variable `equationToSolve` the equation we are going to pass to `dsolve()`.  At this point I have to say that Mathematica is really way more advanced than `sympy`.  If you find yourself really doing some heavy computing you may want to default to numerical methods if you are desperate to stay with Python or consider switching to Mathematica to do symbolic calculations.  I'm really caught between two worlds on this one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This creates a relation between two objects.\n",
    "?Eq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "equationToSolve = Eq(2*current(t).diff(t) + 9*current(t), 100)\n",
    "equationToSolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we solve the equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solutionToEquation = dsolve(equationToSolve, current(t))\n",
    "solutionToEquation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, when we use `dsolve()` the integration constants are not automatically created for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('C1 C2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We substitute the initial condition and then solve for C1.  Afterwards we substitute C1 back into the general solution and plot the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "particularSolution = solutionToEquation.subs([(current(t),0),(t,0)])\n",
    "particularSolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "solutionSet = solveset(particularSolution,C1)\n",
    "solutionList = [a for a in solutionSet]\n",
    "solutionSet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`solveset` returns a set object.  When the set is a `FiniteSet` we can convert to a list and de-reference the list as I've done here.  Not sure this is good practice, in general.  So - I would only do this in an interactive session when I can inspect the return value.  I would not write library code like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solutionToEquation.subs(C1,-100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(solutionToEquation.subs(C1,-100).rhs,(t,0,2));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time with no substitutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L, R, t, V0 = symbols('L R t V0', positive=True)\n",
    "A, V = symbols('A V', cls=Function)\n",
    "\n",
    "rLCircuit = A(t).diff(t) + (R/L)*A(t) - V(t)\n",
    "rLCircuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = dsolve(rLCircuit, A(t))\n",
    "solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "particularSolution = dsolve(rLCircuit, A(t)).subs([(A(t),0),(t,0),(V(t),V0)]).doit()\n",
    "particularSolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "const1 = -L*V0/R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simplify(solution.subs([(C1,const1),(V(t),V0)]).doit())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Second Example First Order ODE: Skydiving"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example is taking from Arfken's book.  Example problem 7.2.1.  The model describes a falling skydiver under a parachute.  We are trying to find the velocity profile as a function of time as well as the terminal velocity.  The equation of motion is:\n",
    "\n",
    "$$ m \\dot{v} = m g - b v^2 $$\n",
    "\n",
    "where v is $v(t)$, a function of time.  Later in the problem we will provide values for $m$, $g$, and $b$ to find a particular solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#from sympy import *\n",
    "#init_session()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mass, g, b, t = symbols('mass g b t', real=True)\n",
    "v = symbols('v', cls=Function)\n",
    "\n",
    "paraChute = mass*v(t).diff(t)-mass*g+b*(v(t))**2\n",
    "paraChute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = dsolve(paraChute,v(t))\n",
    "solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "particularSolution = solution.subs([(v(t),0),(t,0)])\n",
    "particularSolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C1 = symbols('C1')\n",
    "\n",
    "const = solve(particularSolution, C1)\n",
    "const"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to think a bit about our identities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1/tanh(t - I*pi/2).expand(complex=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1/tanh(t - I*pi/2).expand(complex=True).simplify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_solution = (solution.subs(C1,const[1]).rhs).expand().simplify()\n",
    "final_solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can substitute and plot the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%matplotlib notebook\n",
    "\n",
    "plot(final_solution.subs([(g,9.8),(mass,70000),(b,700000)]),(t,0,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy and SciPy give us the ability to do this work numerically, too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%matplotlib notebook\n",
    "\n",
    "from scipy.integrate import odeint  # for integrate.odeint\n",
    "from pylab import *  # for plotting commands\n",
    "\n",
    "def deriv(velocity,time):\n",
    "    mass = 70000.\n",
    "    drag = 700000.\n",
    "    gravity = 9.8\n",
    "    dvdt = gravity-(drag/mass)*velocity**2\n",
    "    return array(dvdt)\n",
    "\n",
    "time = linspace(0.0,1.0,1000)\n",
    "yinit = array([0.0])  # initial values\n",
    "y = odeint(deriv,yinit,time)\n",
    "\n",
    "figure()\n",
    "plot(time, y[:,0])\n",
    "xlabel('t')\n",
    "ylabel('y')\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DIY:  Nonseparable Exact ODE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve the following differential equation and plot the family of solutions for different constants:\n",
    "\n",
    "$$\n",
    "\\frac{d f(x)}{dx} + \\frac{1 + f(x)}{x} = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your dsolve code goes here!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# it is good practice to put your solution back into the differential equation to check the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your plotting code goes here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RL Circuit with Forcing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve Example 7.2.5 from Arfken.  Assume that $I(t=0) = 0$.  Find the solution first, then allow $V(t) = V0$ (a constant).  While we won't do it here, you could let $V(t)$ be a function like $V0 \\sin(\\omega t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%matplotlib notebook\n",
    "#from sympy import *\n",
    "#init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L, R, t = symbols('L R t')\n",
    "A, V = symbols('A V', cls=Function)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rLCircuit = L*A(t).diff(t) + R*A(t) - V(t)\n",
    "rLCircuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.symbols('V0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rLCircuit_constantV = rLCircuit.subs(V(t),V0)\n",
    "rLCircuit_constantV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = dsolve(rLCircuit_constantV, A(t))\n",
    "solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution_with_ics = solution.subs([(A(t),0),(t,0)])\n",
    "solution_with_ics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.symbols('C1')\n",
    "const = sp.solve(solution_with_ics.rhs,C1)\n",
    "const"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "particularSolution = (solution.rhs).subs(C1, ((L/R)*sp.log(-V0))).simplify()\n",
    "particularSolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.plot(particularSolution.subs([(R,20),(V0,10),(L,2)]),(t,0,10));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Postscript:  Using `dsolve` in an Organized Way"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `dsolve` documentation can be found [here](https://docs.sympy.org/latest/modules/solvers/ode.html) and covers the types of equations that can be solved currently.  In a simple case the function `dsolve` is called/used by supplying the differental equation of interest and the function for which the solution is desired.  The results are then used in subsequent calculations to determine the constants in the problem.  We will use the examples below to reinforce this simple view of differential equation solving.  \n",
    "\n",
    "My recommendation for developing your solutions is as follows:\n",
    "\n",
    "* express the equation/model you are interested in solving in the most general format\n",
    "* substitute any known functional relationships that can eliminate unknown terms or derivatives\n",
    "* call `dsolve` to find the general solution for your model\n",
    "* substitute initial conditions into the results to determine the constants\n",
    "* substitute the constants into the general solution to get your particular solution\n",
    "* functionalize the results and visualize the results in a plot (if appropriate)\n",
    "* use `interact` to permit exploration of free parameters and continue to develop an intuititve understanding of your model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_session()\n",
    "sp.init_printing()\n",
    "f = Function('f')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "modelToSolve = Eq(Derivative(f(x), x, x) + 9*f(x),0)\n",
    "modelToSolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "generalSolution = dsolve(modelToSolve, f(x))\n",
    "generalSolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are other ways to specify the derivative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsolve(sp.diff(f(x),x,2) + 9*f(x), f(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Part 2:  Review of Viscoelasticity\n",
    "----\n",
    "\n",
    "The spring obeys Hooke's law: \n",
    "\n",
    "$$\\sigma=M\\epsilon $$\n",
    "\n",
    "where $\\sigma$ is the stress, $M$ is the modulus of elasticity in the linear regime, and $\\epsilon$ is the elastic strain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The dashpot obeys: \n",
    "\n",
    "$$\\sigma = \\eta \\frac{d \\epsilon}{dt}$$\n",
    "\n",
    "where $\\sigma$, and $\\epsilon$ are defined as above and $\\eta$ is related to the viscous response."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Some materials exhibit time dependent behavior under loading.  In polymers this behavior originates from chain elasticity due to coiling and uncoiling and non-reversible motion of the chains relative to one another.  The elastic and viscous behaviors are typically modeled using the spring and [dashpot](https://en.wikipedia.org/wiki/Dashpot).\n",
    "\n",
    "In this lecture we will discuss the Maxwell model and the Kelvin model of viscoelastic materials.  Combinations of springs and dashpots can provide insight into materials behavior.  Combining many such Maxwell and Kelvin models can give insight into the spectrum of dynamic mechanical behavior seen in polymer systems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Maxwell Model\n",
    "\n",
    "![](./images/Maxwell_diagram.png)\n",
    "\n",
    "### Kelvin Model\n",
    "\n",
    "![](./images/Kelvin_Voigt_diagram.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To give you some concept for what we will discuss, this section covers a short derivation for the spring and dashpot in series.  Our first assumption is that for a spring and a dashpot in series that the stress is equally carried by each member in the load train and that the strain is a sum of strains in each element.  This is stated as:  \n",
    "\n",
    "$$ \\epsilon_{total} = \\epsilon_s + \\epsilon_d $$\n",
    "\n",
    "Taking the time derivative for each term produces:  \n",
    "\n",
    "$$ \\frac{d\\epsilon_{\\mathrm{total}}}{dt} = \\frac{d\\epsilon_s}{dt} + \\frac{d\\epsilon_d}{dt} $$\n",
    "\n",
    "We assume that the spring response is instantaneous so that we can write:  \n",
    "\n",
    "$$\\frac{d\\sigma}{dt} = M\\frac{d\\epsilon}{dt}  $$\n",
    "\n",
    "Using the condition that the stresses are the same reduces the number of variables by one.  Using the condition that the strains add up to the total strain results in:\n",
    "\n",
    "$$ \\frac{d \\epsilon_{total}(t)}{dt} = \\frac{1}{M} \\frac{d \\sigma(t)}{dt} + \\frac{\\sigma(t)}{\\eta}$$\n",
    "\n",
    "Don't forget that both stress and strain are functions of time. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models of Viscoelastic Behavior:  Maxwell Model - Constant Stress\n",
    "\n",
    "The physical experiment is one where we apply a constant stress to the Maxwell model.  The model responds by setting the stress in every element as constant and the strain is the sum of the individual strains.  The thought experiment is as follows:\n",
    "\n",
    "* A stress is applied to the series (chain) of elements at rest\n",
    "* The spring responds by Hooke's law instantaneously\n",
    "* The dashpot is unresponsive at $t=0$\n",
    "* From $t>0$ the spring's strain is fixed but the dashpot begins to strain\n",
    "* There is no limit on how far the dashpot can extend\n",
    "* A good exercise is to sketch a diagram of these steps\n",
    "\n",
    "![](./images/Maxwell_diagram.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by setting up our environment.  We import `interact` and `fixed` to help us plot and visualize the results of our computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "from ipywidgets import interact, fixed\n",
    "import sympy as sp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.init_session()\n",
    "sp.var('omega sigma epsilon eta M sigma0 C1 C2')\n",
    "sp.var('sigma epsilon', cls=Function)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Maxwell model is defined as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxwellModel = Eq(epsilon(t).diff(t), sigma(t).diff(t)/M + sigma(t)/eta)\n",
    "maxwellModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "First we will work out the solution for a constant applied stress, $\\sigma_0$.  We start with the general solution and make the substitution that $\\sigma(t)$ is constant by changing the symbol from a function of `t` to $\\sigma_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "generalSolution = dsolve(maxwellModel,epsilon(t))\n",
    "generalSolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "constantStressSolution = dsolve(maxwellModel,epsilon(t)).subs(sigma(t),sigma0).doit()\n",
    "constantStressSolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know that at $t=0$ in this system, Hooke's law defines the strain at $t=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "constantStressSolution.subs([(epsilon(t),sigma0/M),(t,0)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solutionToPlot = constantStressSolution.subs(C1,sigma0/M).rhs\n",
    "solutionToPlot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot the solution with appropriate substitutions for the two parameters and one constant.  We extract the RHS and use `subs` and `ipywidgets` to interactively plot the solution.  Alternatively you could `lambdafy` the solution and use `numpy` to develop the interactive plot (you might be more satisfied with the results)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotSolution(eta0,M0):\n",
    "    plot(solutionToPlot.subs([(eta,eta0),(M,M0),(sigma0,100)]),\n",
    "         (t,0,5),\n",
    "         #ylim=(0,1000),\n",
    "         xlabel='time',\n",
    "         ylabel=r'$\\epsilon(t)$'\n",
    "        )\n",
    "    \n",
    "interact(plotSolution,eta0=(1,100,1),M0=(1,100,1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The major takeawy here is that strain rises linearly with time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models of Viscoelastic Behavior:  Maxwell Model - Constant Strain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the strain is constant the time rate of change is zero and we can make ths substitution right from the start."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "from ipywidgets import interact, fixed\n",
    "import sympy as sp\n",
    "sp.init_session()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.var('t, M, eta, epsilon0, sigma0, C1, C2')\n",
    "sp.var('epsilon, sigma', cls=Function)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Maxwell model is defined as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxwellModel = sp.Eq(epsilon(t).diff(t), sigma(t).diff(t)/M + sigma(t)/eta)\n",
    "maxwellModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "constantStrainModel = maxwellModel.subs(diff(epsilon(t),t),0)\n",
    "constantStrainModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solutionToConstantStrainModel1 = dsolve(constantStrainModel.rhs,sigma(t))\n",
    "solutionToConstantStrainModel1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsolve(maxwellModel,sigma(t)).subs(epsilon(t),0).doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "solutionToConstantStrainModel2 = dsolve(maxwellModel,sigma(t))\n",
    "solution = solutionToConstantStrainModel2.subs(epsilon(t),0).doit()\n",
    "solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution.subs([(t,0),(sigma(0),sigma0)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we make the final substitutions and get our final equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution.subs([(C1,sigma0)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does this say?  We apply an initial stress at time zero.  The material strains to whatever is predicted by Hooke's law.  And then, without changing the strain again through the experiment, the stress drops exponentially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solutionToPlot = solution.subs([(C1,sigma0)]).rhs\n",
    "solutionToPlot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Functionalize and plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import *\n",
    "\n",
    "def plotSolution(M0, eta0):\n",
    "    plot(solutionToPlot.subs([(M,M0),(eta,eta0),(sigma0,100)]),\n",
    "         (t,0,5),\n",
    "         ylim=(0,100),\n",
    "         xlabel='time',\n",
    "         ylabel=r'$\\sigma(t)$'\n",
    "        )\n",
    "    \n",
    "interact(plotSolution,M0=(1,100,1),eta0=(1,100,1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DIY:  Kelvin Model - Constant Stress\n",
    "\n",
    "![](./images/Kelvin_Voigt_diagram.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following similar logic as above, except this time the strains are all taken to be the same and the stresses are additive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "from ipywidgets import interact, fixed\n",
    "import sympy as sp\n",
    "sp.init_session()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.var('t, M, eta, epsilon0, sigma0, C1, C2')\n",
    "sp.var('epsilon, sigma', cls=Function)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kelvinModel = Eq(sigma(t),M*epsilon(t)+eta*epsilon(t).diff(t))\n",
    "kelvinModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this problem we will let stress be a constant for all time:\n",
    "\n",
    "$$ \\sigma(t) = \\sigma_0 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Work on your solution here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Preparation for Solving Fick's Law\n",
    "\n",
    "In this section I pose two important differential equation forms that result from the decomposition of Fick's Law.  As stated above, the steps to follow are:\n",
    "\n",
    "* express the equation/model you are interested in solving in the most general format\n",
    "* substitute any known functional relationships that can eliminate unknown terms or derivatives\n",
    "* call `dsolve` to find the general solution for your model\n",
    "* substitute initial conditions into the results to determine the constants\n",
    "* substitute the constants into the general solution to get your particular solution\n",
    "* functionalize the results and visualize the results in a plot (if appropriate)\n",
    "* use `interact` to permit exploration of free parameters and continue to develop an intuititve understanding of your model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "I'll help you get started by writing down the first model equation to solve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Restart the kernel here.\n",
    "\n",
    "%matplotlib notebook\n",
    "import sympy as sp\n",
    "import numpy as np\n",
    "from ipywidgets import interact, fixed\n",
    "sp.init_session()\n",
    "sp.init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "t = sp.symbols('t')\n",
    "x = sp.symbols('x', cls=Function)\n",
    "k, alpha = sp.symbols('k alpha', positive=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Important ODE One\n",
    "\n",
    "This is the first ODE of interest:\n",
    "\n",
    "$$\n",
    "\\frac{d x(t)}{dt } \\frac{1}{\\alpha x(t)} + k = 0\n",
    "$$\n",
    "\n",
    "Solve this ODE with boundary conditions: $x(t) = 1$ at $t = 0$ and $x(t) = 0$ as $t \\rightarrow \\infty$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "firstModelToSolve = x(t).diff(t)/(x(t)*alpha) + k\n",
    "functionWeAreLookingFor = x(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "firstModelToSolve, functionWeAreLookingFor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "hintList = sp.classify_ode(firstModelToSolve)\n",
    "hintList"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solutionList = [sp.dsolve(firstModelToSolve, functionWeAreLookingFor, hint=hint) for hint in hintList]\n",
    "solutionList"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solutionToFirstModel = sp.dsolve(firstModelToSolve, functionWeAreLookingFor, hint=hintList[2])\n",
    "solutionToFirstModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "sp.var('C1'), solutionToFirstModel.subs([(C1,1)])\n",
    "sp.plot(sp.exp(-2*t),(t,0,10));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Important ODE Two\n",
    "\n",
    "This is the other ODE of interest:\n",
    "\n",
    "$$\n",
    "\\frac{d^2 x(t)}{dt^2} \\frac{1}{x(t)} + k = 0\n",
    "$$\n",
    "\n",
    "Solve this ODE with boundary conditions:\n",
    "\n",
    "$$\n",
    "x(t) = 0\n",
    "$$\n",
    "\n",
    "at $t = 0$ and $t = 2\\pi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import sympy as sp\n",
    "import numpy as np\n",
    "from ipywidgets import interact, fixed\n",
    "sp.init_session()\n",
    "sp.init_printing()\n",
    "\n",
    "t = sp.symbols('t')\n",
    "x = sp.symbols('x', cls=Function)\n",
    "k, alpha = sp.symbols('k alpha', positive=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "secondModelToSolve = x(t).diff(t,2)/x(t) + k**2\n",
    "sp.var('C1 C2')\n",
    "solutionToSecondModel = sp.dsolve(secondModelToSolve).subs(C2,0)\n",
    "solutionToSecondModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Homework\n",
    "----\n",
    "\n",
    "Siolve the Maxwell model for a time varying input stress.  Solve in `sympy` and plot with `numpy` (you'll get better plotting performance.  You will be graded on your solution as well as your visualization and written response.  Choose a forcing function such as:\n",
    "\n",
    "$$\\sigma(t) = \\sigma_0 \\sin(\\omega t)$$\n",
    "\n",
    "for this problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}