diffeq-uncertainties.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to solve the differential equation `du/dt = f(t, u)`, where `f(t, u) = α * u`, in the interval `t ∈ [0, 1]`.  We know that the solution to this equation is `u(t) = u_0 exp(α * t)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using DifferentialEquations, Measurements # Load the packages\n",
    "α = (1.01 ± 0.01) # Define the α parameter as a number with uncertainty\n",
    "f(t, u) = α * u # Define the f function\n",
    "u0 = 1/2 ± 0 # The initial value of the Cauchy problem\n",
    "tspan = (0.0 ± 0, 1.0 ± 0) # The extrema of the domain in which we want to solve the equation\n",
    "prob = ODEProblem(f, u0, tspan) # Define the ODEProblem\n",
    "sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8); # Solve the problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is the analytic solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "u = u0 .* exp.(α .* sol.t);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the solution found by `DifferentialEquations.jl` and the analytic solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "At time 0.0:\n",
      "0.5 ± 0.0\n",
      "0.5 ± 0.0\n",
      "\n",
      "At time 0.012407826196308189:\n",
      "0.5063053789114713 ± 3.769289486273972e-5\n",
      "0.5063053789114712 ± 3.7692894862736955e-5\n",
      "\n",
      "At time 0.042500886481028156:\n",
      "0.52193026853015 ± 6.72494218852177e-5\n",
      "0.521930268530135 ± 6.724942188505842e-5\n",
      "\n",
      "At time 0.0817803932926208:\n",
      "0.5430526603660298 ± 8.15081766960003e-5\n",
      "0.5430526603659422 ± 8.1508176695602e-5\n",
      "\n",
      "At time 0.12887323477424767:\n",
      "0.5695064255096831 ± 0.00010898724594719012\n",
      "0.5695064255093786 ± 0.00010898724594568129\n",
      "\n",
      "At time 0.18409719135533997:\n",
      "0.6021738925732998 ± 0.00013105341682102388\n",
      "0.602173892572426 ± 0.00013105341681784647\n",
      "\n",
      "At time 0.24627365251891972:\n",
      "0.6412019663929366 ± 0.00015886106381869243\n",
      "0.64120196639087 ± 0.00015886106381198824\n",
      "\n",
      "At time 0.31479180297608794:\n",
      "0.6871467088572084 ± 0.00018569231191653082\n",
      "0.6871467088529182 ± 0.00018569231190505639\n",
      "\n",
      "At time 0.38859490981815525:\n",
      "0.7403247618946622 ± 0.0002155136694726928\n",
      "0.740324761886693 ± 0.00021551366945420482\n",
      "\n",
      "At time 0.46685971668181403:\n",
      "0.8012206780339247 ± 0.0002462164149989751\n",
      "0.8012206780203366 ± 0.0002462164149718729\n",
      "\n",
      "At time 0.5487136259449821:\n",
      "0.8702746579779278 ± 0.00027921963537337163\n",
      "0.8702746579563262 ± 0.0002792196353356099\n",
      "\n",
      "At time 0.6334320104824218:\n",
      "0.9480188907838683 ± 0.00031402163556843257\n",
      "0.9480188907514014 ± 0.0003140216355183908\n",
      "\n",
      "At time 0.72035921761696:\n",
      "1.0350146997549017 ± 0.00035123419392624525\n",
      "1.035014699708293 ± 0.00035123419386214993\n",
      "\n",
      "At time 0.8089540967699169:\n",
      "1.1318986650134468 ± 0.0003908984317072107\n",
      "1.1318986649489977 ± 0.0003908984316274656\n",
      "\n",
      "At time 0.8987609607927629:\n",
      "1.2393677739920392 ± 0.0004334066410309232\n",
      "1.2393677739056368 ± 0.00043340664093383026\n",
      "\n",
      "At time 0.9894116178526793:\n",
      "1.3581976365396038 ± 0.0004789992744179338\n",
      "1.3581976364267052 ± 0.00047899927430177354\n",
      "\n",
      "At time 1.0:\n",
      "1.3728005076225709 ± 0.013728005076303318\n",
      "1.3728005075084582 ± 0.013728005075084582\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for i in eachindex(sol)\n",
    "    println(\"At time \", Measurements.value(sol.t[i]), \":\\n\", sol[i], \"\\n\", u[i], \"\\n\")\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}