{
"cells": [
{
"cell_type": "markdown",
"id": "d448df07",
"metadata": {},
"source": [
"# Solving the pendulum equation"
]
},
{
"cell_type": "markdown",
"id": "749f620e",
"metadata": {},
"source": [
"We experiment with Euler's method and the modified Euler method\n",
"on the pendulum problem.\n",
"For the gravitational constant, we take 9.807."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "84e34877",
"metadata": {},
"outputs": [],
"source": [
"using Printf"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "0ba538f8",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"(process:36004): GLib-GIO-WARNING **: 16:53:10.737: Unexpectedly, UWP app `Clipchamp.Clipchamp_2.5.1.0_neutral__yxz26nhyzhsrt' (AUMId `Clipchamp.Clipchamp_yxz26nhyzhsrt!App') supports 41 extensions but has no verbs\n",
"\n",
"(process:36004): GLib-GIO-WARNING **: 16:53:10.800: Unexpectedly, UWP app `Evernote.Evernote_10.47.7.0_x64__q4d96b2w5wcc2' (AUMId `Evernote.Evernote_q4d96b2w5wcc2!Evernote') supports 1 extensions but has no verbs\n"
]
}
],
"source": [
"using Plots"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "3047fdbb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"eulerpend"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" eulerpend(n::Int64,verbose::Bool=true)\n",
"\n",
"applies Euler's method with n*p evaluations\n",
"on the interval [0,p*2*pi] on the pendulum problem.\n",
"\"\"\"\n",
"function eulerpend(p::Int,n::Int,verbose::Bool=true)\n",
" h = 2*p*pi/n\n",
" if verbose\n",
" println(\" i t y1 y2 \")\n",
" end\n",
" (rt, ry1, ry2) = (zeros(n+1), zeros(n+1), zeros(n+1))\n",
" y1 = pi/4\n",
" y2 = 0\n",
" (rt[1], ry1[1], ry2[1]) = (0, y1, y2)\n",
" for i=1:n\n",
" t = i*h\n",
" (y1, y2) = (y1 + h*y2, y2 + h*(-9.807*sin(y1)))\n",
" (rt[i+1], ry1[i+1], ry2[i+1]) = (t, y1, y2)\n",
" if verbose\n",
" stri = @sprintf(\"%3d\", i)\n",
" strt = @sprintf(\"%5.2f\", t)\n",
" stry1 = @sprintf(\"%9.6e\", y1)\n",
" stry2 = @sprintf(\"%9.6e\", y2)\n",
" println(\"$stri $strt $stry1 $stry2\")\n",
" end\n",
" end\n",
" return (rt, ry1, ry2)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "455f1d66",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"eulermodpend"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" eulermodpend(n::Int64,verbose::Bool=true)\n",
"\n",
"applies the modified Euler method with n*p evaluations\n",
"on the interval [0,p*2*pi] on the pendulum problem.\n",
"\"\"\"\n",
"function eulermodpend(p::Int,n::Int,verbose::Bool=true)\n",
" h = 2*p*pi/n\n",
" if verbose\n",
" println(\" i t y1 y2 \")\n",
" end\n",
" (rt, ry1, ry2) = (zeros(n+1), zeros(n+1), zeros(n+1))\n",
" y1 = pi/4\n",
" y2 = 0\n",
" (rt[1], ry1[1], ry2[1]) = (0, y1, y2)\n",
" for i=1:n\n",
" t = i*h\n",
" (y1a, y2a) = (y1 + h*y2, y2 + h*(-9.807*sin(y1)))\n",
" (y1b, y2b) = (y1, y2)\n",
" y1 = y1 + (h/2)*(y2a + y2b)\n",
" y2 = y2 + (h/2)*(-9.807*sin(y1a)-9.807*sin(y1b))\n",
" (rt[i+1], ry1[i+1], ry2[i+1]) = (t, y1, y2)\n",
" if verbose\n",
" stri = @sprintf(\"%3d\", i)\n",
" strt = @sprintf(\"%.2f\", t)\n",
" stry1 = @sprintf(\"%.6e\", y1)\n",
" stry2 = @sprintf(\"%.6e\", y2)\n",
" println(\"$stri $strt $stry1 $stry2\")\n",
" end\n",
" end\n",
" return (rt, ry1, ry2)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "83339ebc",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(3, 100)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(p, n) = (3, 100)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "dd939b31",
"metadata": {},
"outputs": [],
"source": [
"t, y1, y2 = eulerpend(p,p*n,false);"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "fd2bc797",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(t,y1,label=\"position\")\n",
"plot!(t,y2,label=\"velocity\")"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "ffe41db4",
"metadata": {},
"outputs": [],
"source": [
"# savefig(\"figeulerpend.png\")"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "e1cc5f61",
"metadata": {},
"outputs": [],
"source": [
"mt, my1, my2 = eulermodpend(p,p*n,false);"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "59267b59",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(mt,my1,label=\"position\")\n",
"plot!(mt,my2,label=\"velocity\")"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "f6ca1078",
"metadata": {},
"outputs": [],
"source": [
"# savefig(\"figeulermodpend.png\")"
]
},
{
"cell_type": "markdown",
"id": "a3264538",
"metadata": {},
"source": [
"# A Predator-Prey Simulation"
]
},
{
"cell_type": "markdown",
"id": "451a6c0d",
"metadata": {},
"source": [
"A simple model of the evolution of the population sizes of a prey and a predator species is defined by a system of two nonlinear differential equations:"
]
},
{
"cell_type": "markdown",
"id": "708cf707",
"metadata": {},
"source": [
" $$\n",
" \\left\\{\n",
" \\begin{array}{lcr}\n",
" {\\displaystyle \\frac{d r(t)}{d t}}\n",
" & = & a~\\! r(t) + b~\\! r(t) f(t), \\quad a > 0, b < 0,\n",
" \\\\\n",
" {\\displaystyle \\frac{d f(t)}{d t}}\n",
" & = & c~\\! f(t) + d~\\! r(t) f(t), \\quad c < 0, d > 0.\n",
" \\end{array}\n",
" \\right.\n",
" $$"
]
},
{
"cell_type": "markdown",
"id": "179c367e",
"metadata": {},
"source": [
"We can apply the modified Euler method in our computer simulation. The constants $a$, $b$, $c$, $d$ are defined by ``A``, ``B``, ``C``, ``D`` below:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "a957a80c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.25"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = 1.1 # growth rate of the prey\n",
"B = -0.5 # correction to prey from encounter with predator\n",
"C = -0.75 # decay rate of the predator\n",
"D = 0.25 # correction to predator from encounter with prey"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "f3988edb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"modeulerpredprey"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" modeulerpredprey(T::Float64,n::Int,\n",
" y1t0::Float64,y2t0::Float64,\n",
" verbose::Bool=true)\n",
"\n",
"applies the modified Euler method to the predator-prey model\n",
"with time span [0, T], n steps, and initial values (y1t0, y2t0),\n",
"respectively for y1 and y2.\n",
"\"\"\"\n",
"function modeulerpredprey(T::Float64,n::Int,\n",
" y1t0::Float64,y2t0::Float64,\n",
" verbose::Bool=true)\n",
" h = T/n\n",
" if verbose\n",
" println(\" i t y1 y2 \")\n",
" end\n",
" (rt, ry1, ry2) = (zeros(n+1), zeros(n+1), zeros(n+1))\n",
" y1 = y1t0\n",
" y2 = y2t0\n",
" (rt[1], ry1[1], ry2[1]) = (0, y1, y2)\n",
" for i=1:n\n",
" t = i*h\n",
" F1p = A*y1 + B*y1*y2\n",
" F2p = C*y2 + D*y1*y2\n",
" (y1p, y2p) = (y1 + h*F1p, y2 + h*F2p)\n",
" F1c = A*y1p + B*y1p*y2p\n",
" F2c = C*y2p + D*y1p*y2p\n",
" y1 = y1 + (h/2)*(F1p + F1c)\n",
" y2 = y2 + (h/2)*(F2p + F2c)\n",
" (rt[i+1], ry1[i+1], ry2[i+1]) = (t, y1, y2)\n",
" if verbose\n",
" stri = @sprintf(\"%3d\", i)\n",
" strt = @sprintf(\"%.2f\", t)\n",
" stry1 = @sprintf(\"%.6e\", y1)\n",
" stry2 = @sprintf(\"%.6e\", y2)\n",
" println(\"$stri $strt $stry1 $stry2\")\n",
" end\n",
" end\n",
" return (rt, ry1, ry2)\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "91eadfeb",
"metadata": {},
"source": [
"To run the simulation, we need to define a time interval and the number of steps. Also, we need two initial values for the population sizes. The setup of the simulation is defined below:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "fd8ae8c1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(3.0, 2.0)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(endtime, nbpts) = (20.0, 100)\n",
"(y1valt0, y2valt0) = (3.0, 2.0)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "53ddc896",
"metadata": {},
"outputs": [],
"source": [
"t, y1, y2 = modeulerpredprey(endtime,nbpts,y1valt0,y2valt0,false);"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "48aea8c0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1 0.00 3.000000e+00 2.000000e+00 0.00e+00\n",
" 11 2.00 3.341181e+00 2.256135e+00 3.41e-01\n",
" 21 4.00 2.812097e+00 2.375841e+00 1.88e-01\n",
" 31 6.00 2.766754e+00 2.055365e+00 2.33e-01\n",
" 41 8.00 3.299205e+00 2.090896e+00 2.99e-01\n",
" 51 10.00 3.079936e+00 2.408659e+00 7.99e-02\n",
" 61 12.00 2.669408e+00 2.195260e+00 3.31e-01\n",
" 71 14.00 3.069206e+00 2.002071e+00 6.92e-02\n",
" 81 16.00 3.314841e+00 2.299640e+00 3.15e-01\n",
" 91 18.00 2.758366e+00 2.348322e+00 2.42e-01\n"
]
}
],
"source": [
"freq = 10\n",
"idx = 1\n",
"step = Int(nbpts/freq)\n",
"for k=1:freq\n",
" stri = @sprintf(\"%5d\", idx)\n",
" strt = @sprintf(\"%5.2f\", t[idx])\n",
" stry1 = @sprintf(\"%13.6e\", y1[idx])\n",
" stry2 = @sprintf(\"%13.6e\", y2[idx])\n",
" serr = @sprintf(\"%.2e\", abs(y1[idx] - y1[1]))\n",
" println(\"$stri $strt $stry1 $stry2 $serr\")\n",
" idx = idx + step\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "5c610841",
"metadata": {},
"source": [
"The numbers look plausable, but only the plots will show the expected periodic behavior."
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "ae925e3e",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(t, y1, color=\"green\", label=\"prey\")\n",
"plot!(t, y2, color=\"red\", label=\"predator\")"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "ec4725d0",
"metadata": {},
"outputs": [],
"source": [
"savefig(\"figpredpreysim.png\")"
]
},
{
"cell_type": "markdown",
"id": "2db50a8b",
"metadata": {},
"source": [
"Let us plot the phase portrait."
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "84fd5982",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(y1,y2,color=\"blue\", xlabel=\"prey\", ylabel=\"predator\", label=\"\")"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "3801e624",
"metadata": {},
"outputs": [],
"source": [
"# savefig(\"figpredpreyphase.png\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.8.0",
"language": "julia",
"name": "julia-1.8"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.0"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": false
}
},
"nbformat": 4,
"nbformat_minor": 5
}