{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem 3.1: SIS model with risk groups in Julia\n", "\n", "Author: Emma Accorsi @emmaaccorsi\n", "\n", "Date: 2018-10-02" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "using DifferentialEquations" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "sis_ode (generic function with 1 method)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function sis_ode(du,u,p,t)\n", " SH,IH,SL,IL = u\n", " betaHH,betaHL,betaLH,betaLL,gamma=p\n", " du[1]=-(betaHH*IH+betaHL*IL)*SH+gamma*IH\n", " du[2]=+(betaHH*IH+betaHL*IL)*SH-gamma*IH\n", " du[3]=-(betaLH*IH+betaLL*IL)*SL+gamma*IL\n", " du[4]=+(betaLH*IH+betaLL*IL)*SL-gamma*IL\n", " end" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.0, 15.0)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "parms =[10,0.1,0.1,1,1]\n", "init=[0.19999,0.00001,0.799,0.001]\n", "tspan=tspan = (0.0,15.0)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "retcode: Success\n", "Interpolation: 1st order linear\n", "t: 151-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", " 13.9\n", " 14.0\n", " 14.1\n", " 14.2\n", " 14.3\n", " 14.4\n", " 14.5\n", " 14.6\n", " 14.7\n", " 14.8\n", " 14.9\n", " 15.0\n", "u: 151-element Array{Array{Float64,1},1}:\n", " [0.19999, 1.0e-5, 0.799, 0.001] \n", " [0.199987, 1.31343e-5, 0.79902, 0.000980193] \n", " [0.199983, 1.65569e-5, 0.799039, 0.000960808]\n", " [0.19998, 2.0299e-5, 0.799058, 0.000941839] \n", " [0.199976, 2.43948e-5, 0.799077, 0.00092328] \n", " [0.199971, 2.88825e-5, 0.799095, 0.000905125]\n", " [0.199966, 3.38042e-5, 0.799113, 0.000887371]\n", " [0.199961, 3.92064e-5, 0.79913, 0.000870012] \n", " [0.199955, 4.51401e-5, 0.799147, 0.000853044]\n", " [0.199948, 5.16616e-5, 0.799164, 0.000836465]\n", " [0.199941, 5.88328e-5, 0.79918, 0.00082027] \n", " [0.199933, 6.67232e-5, 0.799196, 0.000804459]\n", " [0.199925, 7.54101e-5, 0.799211, 0.000789029]\n", " ⋮ \n", " [0.100111, 0.0998891, 0.775299, 0.024701] \n", " [0.100078, 0.0999224, 0.775082, 0.0249178] \n", " [0.100048, 0.0999524, 0.774871, 0.0251294] \n", " [0.100021, 0.0999794, 0.774664, 0.0253357] \n", " [0.0999962, 0.100004, 0.774463, 0.0255369] \n", " [0.0999737, 0.100026, 0.774267, 0.0257331] \n", " [0.0999526, 0.100047, 0.774076, 0.0259244] \n", " [0.0999326, 0.100067, 0.773889, 0.0261108] \n", " [0.0999141, 0.100086, 0.773708, 0.0262925] \n", " [0.0998972, 0.100103, 0.77353, 0.0264696] \n", " [0.0998817, 0.100118, 0.773358, 0.0266423] \n", " [0.0998676, 0.100132, 0.773189, 0.0268105] " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sis_prob = ODEProblem(sis_ode,init,tspan,parms)\n", "sis_sol = solve(sis_prob,saveat=0.1)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "using Plots" ] }, { "cell_type": "code", "execution_count": 6, "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", "0\n", "\n", "\n", "5\n", "\n", "\n", "10\n", "\n", "\n", "15\n", "\n", "\n", "0.0\n", "\n", "\n", "0.2\n", "\n", "\n", "0.4\n", "\n", "\n", "0.6\n", "\n", "\n", "0.8\n", "\n", "\n", "Time (Years)\n", "\n", "\n", "Proportion of Population\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "u1(t)\n", "\n", "\n", "\n", "u2(t)\n", "\n", "\n", "\n", "u3(t)\n", "\n", "\n", "\n", "u4(t)\n", "\n", "\n" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(sis_sol,xlabel=\"Time (Years)\",ylabel=\"Proportion of Population\")" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.3", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.3" } }, "nbformat": 4, "nbformat_minor": 2 }