{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## SIR model in Julia using DifferentialEquations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Author*: Simon Frost @sdwfrost\n", "*Editor*: Chris Rackauckas @ChrisRackauckas\n", "\n", "*Date*: 2018-07-12\n", "\n", "### Using the Macro DSL" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "using DifferentialEquations" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "┌ Warning: Symbolic calculations could not initiate. Likely there's a function which is not differentiable by SymEngine.\n", "└ @ ParameterizedFunctions /home/simon/.julia/packages/ParameterizedFunctions/ozDxQ/src/ode_def_opts.jl:244\n" ] }, { "data": { "text/plain": [ "(::SIRModel{getfield(Main, Symbol(\"##3#4\")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any}) (generic function with 2 methods)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sir_ode = @ode_def SIRModel begin\n", " dS = -β*S*I\n", " dI = β*S*I-γ*I\n", " dR = γ*I\n", " end β γ" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.0, 200.0)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "parms = [0.1,0.05]\n", "init = [0.99,0.01,0.0]\n", "tspan = (0.0,200.0)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[36mODEProblem\u001b[0m with uType \u001b[36mArray{Float64,1}\u001b[0m and tType \u001b[36mFloat64\u001b[0m. In-place: \u001b[36mtrue\u001b[0m\n", "timespan: (0.0, 200.0)\n", "u0: [0.99, 0.01, 0.0]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sir_prob = ODEProblem(sir_ode,init,tspan,parms)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "sir_sol = solve(sir_prob,saveat = 0.1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualisation" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "using Plots" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "0\n", "\n", "\n", "50\n", "\n", "\n", "100\n", "\n", "\n", "150\n", "\n", "\n", "200\n", "\n", "\n", "0.00\n", "\n", "\n", "0.25\n", "\n", "\n", "0.50\n", "\n", "\n", "0.75\n", "\n", "\n", "1.00\n", "\n", "\n", "Time\n", "\n", "\n", "Number\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "S(t)\n", "\n", "\n", "\n", "I(t)\n", "\n", "\n", "\n", "R(t)\n", "\n", "\n" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(sir_sol,xlabel=\"Time\",ylabel=\"Number\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using the Function Interface" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "retcode: Success\n", "Interpolation: 1st order linear\n", "t: 2001-element Array{Float64,1}:\n", " 0.0\n", " 0.1\n", " 0.2\n", " 0.3\n", " 0.4\n", " 0.5\n", " 0.6\n", " 0.7\n", " 0.8\n", " 0.9\n", " 1.0\n", " 1.1\n", " 1.2\n", " ⋮ \n", " 198.9\n", " 199.0\n", " 199.1\n", " 199.2\n", " 199.3\n", " 199.4\n", " 199.5\n", " 199.6\n", " 199.7\n", " 199.8\n", " 199.9\n", " 200.0\n", "u: 2001-element Array{Array{Float64,1},1}:\n", " [0.99, 0.01, 0.0] \n", " [0.989901, 0.0100491, 5.01227e-5] \n", " [0.989801, 0.0100985, 0.000100492]\n", " [0.989701, 0.010148, 0.000151108] \n", " [0.9896, 0.0101979, 0.000201972] \n", " [0.989499, 0.0102479, 0.000253087]\n", " [0.989397, 0.0102982, 0.000304452]\n", " [0.989295, 0.0103487, 0.000356069]\n", " [0.989193, 0.0103995, 0.000407939]\n", " [0.989089, 0.0104504, 0.000460064]\n", " [0.988986, 0.0105017, 0.000512444]\n", " [0.988882, 0.0105531, 0.000565081]\n", " [0.988777, 0.0106049, 0.000617976]\n", " ⋮ \n", " [0.21018, 0.0149422, 0.774878] \n", " [0.210148, 0.014899, 0.774953] \n", " [0.210117, 0.0148558, 0.775027] \n", " [0.210086, 0.0148128, 0.775101] \n", " [0.210055, 0.01477, 0.775175] \n", " [0.210024, 0.0147272, 0.775249] \n", " [0.209993, 0.0146845, 0.775322] \n", " [0.209962, 0.014642, 0.775396] \n", " [0.209932, 0.0145996, 0.775469] \n", " [0.209901, 0.0145573, 0.775542] \n", " [0.209871, 0.0145151, 0.775614] \n", " [0.20984, 0.0144731, 0.775687] " ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function sir_ode2(du,u,p,t)\n", " S,I,R = u\n", " b,g = p\n", " du[1] = -b*S*I\n", " du[2] = b*S*I-g*I\n", " du[3] = g*I\n", "end\n", "parms = [0.1,0.05]\n", "init = [0.99,0.01,0.0]\n", "tspan = (0.0,200.0)\n", "sir_prob2 = ODEProblem(sir_ode2,init,tspan,parms)\n", "sir_sol = solve(sir_prob2,saveat = 0.1)" ] } ], "metadata": { "@webio": { "lastCommId": "1d6f38bbda8746869f917ca72588a0e4", "lastKernelId": "46125098-ed1c-451b-91a7-126686e11ad3" }, "kernelspec": { "display_name": "Julia 1.0.1", "language": "julia", "name": "julia-1.0" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.0.1" } }, "nbformat": 4, "nbformat_minor": 2 }