{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# a 4-stage Runge-Kutta method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The purpose of this notebook is to illustrate a general function to apply a 4-stage Runge-Kutta method to solve an initial value problem."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The system of first-order ordinary differential equations is defined by the right hand side vector. This right hand side is defined by a function of two arguments: one number (the independent variable) and one vector (the dependent variables). "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"using Printf"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"(process:35296): GLib-GIO-WARNING **: 18:35:16.930: 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:35296): GLib-GIO-WARNING **: 18:35:16.977: 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,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"rk4"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" rk4(F::Function,\n",
" h::Float64,endT::Float64,\n",
" init::Array{Float64})\n",
"\n",
"applies a 4-stage Runge-Kutta method with step size h,\n",
"for the time interval [0, endT] with initial values in init.\n",
"\n",
"The function F has two input arguments:\n",
"1. a number, the value of the independent variable (time),\n",
"2. a vector with the values of the dependent variables.\n",
"\n",
"Returns an array of size length(init) + 1 with arrays all\n",
"of the same length: the number of time steps.\n",
"\"\"\"\n",
"function rk4(F::Function,\n",
" h::Float64,endT::Float64,\n",
" init::Array{Float64})\n",
" dim = length(init)\n",
" k1 = zeros(dim)\n",
" k2 = zeros(dim)\n",
" k3 = zeros(dim)\n",
" k4 = zeros(dim)\n",
" n = Int(ceil(endT/h))\n",
" result = [zeros(n+1) for k=1:dim+1]\n",
" for k=1:dim\n",
" result[k+1][1] = init[k]\n",
" end\n",
" y = init\n",
" for i=1:n\n",
" t = i*h\n",
" k1 = F(t, y)\n",
" k2 = F(t + h/2, y + (h/2)*k1)\n",
" k3 = F(t + h/2, y + (h/2)*k2)\n",
" k4 = F(t + h, y + h*k3)\n",
" y = y + (h/6)*(k1 + 2*k2 + 2*k3 + k4)\n",
" result[1][i+1] = t\n",
" for k=1:dim\n",
" result[k+1][i+1] = y[k]\n",
" end\n",
" end\n",
" return result\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us test the function ``rk4`` on the pendulum problem.\n",
"Written as a system of first-order differential equations, the pendulum problem appears in standard form as\n",
"\n",
"$$\n",
" \\begin{eqnarray*}\n",
" y'_1 & = & y_2 \\\\\n",
" y'_2 & = & - \\frac{g}{\\ell} \\sin(y_1)\n",
" \\end{eqnarray*}\n",
"$$\n",
"\n",
"For the length $\\ell$ of the string we take $\\ell = 1$\n",
"and the value for the gravitational constant $g = 9.807$.\n",
"\n",
"The first variable $y_1(t)$ represents an angle. When zero, then the pendulum is at rest.\n",
"The second variable $y_2(t)$ is the velocity of the pendulum."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The right hand side of the system of first-order differential equations is defined by the ``f`` below."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"f (generic function with 1 method)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(t, y) = [y[2], -9.807*sin(y[1])]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us test ``f`` on the initial conditions: $y_1(0) = \\pi/4$ (initial angle of deviation) and zero velocity $y_2(0) = 0$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Vector{Float64}:\n",
" 0.0\n",
" -6.934596203096471"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"initc = [pi/4, 0]\n",
"f(0, initc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have the initial value problem correctly defined, we select a time interval ``T`` and a step size ``h``. Then we can run the function."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Vector{Float64}}:\n",
" [0.0, 0.2617993877991494, 0.5235987755982988, 0.7853981633974483, 1.0471975511965976, 1.308996938995747, 1.5707963267948966, 1.832595714594046, 2.0943951023931953, 2.356194490192345 … 3.926990816987241, 4.1887902047863905, 4.45058959258554, 4.71238898038469, 4.974188368183839, 5.235987755982988, 5.497787143782138, 5.759586531581287, 6.021385919380436, 6.283185307179586]\n",
" [0.7853981633974483, 0.5577022297935539, -0.001390003368037518, -0.5576223291179553, -0.7803391428962481, -0.551077450418934, 0.005972007651994571, 0.5573614896924796, 0.7752715503816151, 0.5436232946861215 … 0.5591568455428995, 0.7650154106355238, 0.5262556049004672, -0.026570786902836208, -0.5611198316981506, -0.7597682634544577, -0.5163484524790499, 0.03559487675366424, 0.5637258740137575, 0.7543960410076648]\n",
" [0.0, -1.6597334564648372, -2.3896421246361914, -1.6464723381575372, 0.011519456902007175, 1.6593276148371194, 2.3749551309457595, 1.626165190602221, -0.026504002586837494, -1.6614356540635091 … 1.5777909648432975, -0.0664686625045996, -1.6726710103948312, -2.3299770462820866, -1.5497731488884736, 0.09124910132818576, 1.6815144032265286, 2.314353655234497, 1.5192212239336786, -0.11909428196752514]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"T = 2*pi\n",
"h = T/24\n",
"r = rk4(f, h, T, [pi/4, 0.0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us print the output better."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1 0.00 7.853982e-01 0.000000e+00\n",
" 2 0.26 5.577022e-01 -1.659733e+00\n",
" 3 0.52 -1.390003e-03 -2.389642e+00\n",
" 4 0.79 -5.576223e-01 -1.646472e+00\n",
" 5 1.05 -7.803391e-01 1.151946e-02\n",
" 6 1.31 -5.510775e-01 1.659328e+00\n",
" 7 1.57 5.972008e-03 2.374955e+00\n",
" 8 1.83 5.573615e-01 1.626165e+00\n",
" 9 2.09 7.752716e-01 -2.650400e-02\n",
" 10 2.36 5.436233e-01 -1.661436e+00\n",
" 11 2.62 -1.172142e-02 -2.360173e+00\n",
" 12 2.88 -5.578880e-01 -1.603262e+00\n",
" 13 3.14 -7.701727e-01 4.485365e-02\n",
" 14 3.40 -5.353478e-01 1.665929e+00\n",
" 15 3.67 1.860030e-02 2.345213e+00\n",
" 16 3.93 5.591568e-01 1.577791e+00\n",
" 17 4.19 7.650154e-01 -6.646866e-02\n",
" 18 4.45 5.262556e-01 -1.672671e+00\n",
" 19 4.71 -2.657079e-02 -2.329977e+00\n",
" 20 4.97 -5.611198e-01 -1.549773e+00\n",
" 21 5.24 -7.597683e-01 9.124910e-02\n",
" 22 5.50 -5.163485e-01 1.681514e+00\n",
" 23 5.76 3.559488e-02 2.314354e+00\n",
" 24 6.02 5.637259e-01 1.519221e+00\n",
" 25 6.28 7.543960e-01 -1.190943e-01\n"
]
}
],
"source": [
"t, y1, y2 = r\n",
"for i=1:length(t)\n",
" stri = @sprintf(\"%3d\", i)\n",
" strt = @sprintf(\"%.2f\", t[i])\n",
" stry1 = @sprintf(\"%.6e\", y1[i])\n",
" stry2 = @sprintf(\"%.6e\", y2[i])\n",
" println(\"$stri $strt $stry1 $stry2\")\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A plot summarizes the computations best."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(t,y1,legend=:bottomright)\n",
"plot!(t,y2)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# savefig(\"figrk4pend.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The jagged nature of the plot shows that we have used a minimum amount of steps to reach a plausible correct plot of perpetual motion."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The 3-Body Problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A more interesting problem is the motion of three bodies. For simplicity we restrict to the plane."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We consider three bodies with respective masses $m_1$, $m_2$, $m_3$\n",
"in the plane with positions\n",
"$(x_1(t),y_1(t))$, $(x_2(t),y_2(t))$, $(x_3(t),y_3(t))$\n",
"evolving over time $t$, governed by a system of second order\n",
"differential equations, shown below for the movement of the first body:\n",
"\\begin{eqnarray*}\n",
" \\frac{d^2 x_1(t)}{dt^2} & \\!\\!=\\!\\! &\n",
" - \\frac{m_2 (x_1(t) - x_2(t))}\n",
" {\\left( (x_1(t) - x_2(t))^2 + (y_1(t) - y_2(t))^2 \\right)^{3/2}}\n",
" - \\frac{m_3 (x_1(t) - x_3(t))}\n",
" {\\left( (x_1(t) - x_3(t))^2 + (y_1(t) - y_3(t))^2 \\right)^{3/2}} \\\\\n",
" \\frac{d^2 y_1(t)}{dt^2} & \\!\\!=\\!\\! &\n",
" - \\frac{m_2 (y_1(t) - y_2(t))}\n",
" {\\left( (x_1(t) - x_2(t))^2 + (y_1(t) - y_2(t))^2 \\right)^{3/2}}\n",
" - \\frac{m_3 (y_1(t) - y_3(t))}\n",
" {\\left( (x_1(t) - x_3(t))^2 + (y_1(t) - y_3(t))^2 \\right)^{3/2}} \\\\\n",
"\\end{eqnarray*}\n",
"Four additional equations determine the positions for the second\n",
"and third body."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To turn this into a system of first order differential\n",
"equations we introduce new variables\n",
"$u_i$, $v_i$ for the velocities of $x_i$, $y_i$ so we have\n",
"\\begin{eqnarray*}\n",
" \\frac{d x_1(t)}{dt} & \\!\\!=\\!\\! & u_1(t) \\\\\n",
" \\frac{d u_1(t)}{dt^2} & \\!\\!=\\!\\! &\n",
" - \\frac{m_2 (x_1(t) - x_2(t))}\n",
" {\\left( (x_1(t) - x_2(t))^2 + (y_1(t) - y_2(t))^2 \\right)^{3/2}}\n",
" - \\frac{m_3 (x_1(t) - x_3(t))}\n",
" {\\left( (x_1(t) - x_3(t))^2 + (y_1(t) - y_3(t))^2 \\right)^{3/2}} \\\\\n",
" \\frac{d y_1(t)}{dt} & \\!\\!=\\!\\! & v_1(t) \\\\\n",
" \\frac{d v_1(t)}{dt} & \\!\\!=\\!\\! &\n",
" - \\frac{m_2 (y_1(t) - y_2(t))}\n",
" {\\left( (x_1(t) - x_2(t))^2 + (y_1(t) - y_2(t))^2 \\right)^{3/2}}\n",
" - \\frac{m_3 (y_1(t) - y_3(t))}\n",
" {\\left( (x_1(t) - x_3(t))^2 + (y_1(t) - y_3(t))^2 \\right)^{3/2}} \\\\\n",
"\\end{eqnarray*}\n",
"Applying this rewriting leads to a system of 12 first order differential\n",
"equations in the positions and velocities of the three bodies.\n",
"Defining the masses, initial positions and velocities leads to an\n",
"initial value problem."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"F"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" Right hand size for the 3-body problem with unit masses.\n",
"\"\"\"\n",
"function F(t, y)\n",
" x1, u1, y1, v1, x2, u2, y2, v2, x3, u3, y3, v3 = y\n",
" r = [u1, -(x1-x2)/((x1-x2)^2 + (y1-y2)^2)^(3/2) - (x1-x3)/((x1-x3)^2 + (y1-y3)^2)^(3/2),\n",
" v1, -(y1-y2)/((x1-x2)^2 + (y1-y2)^2)^(3/2) - (y1-y3)/((x1-x3)^2 + (y1-y3)^2)^(3/2),\n",
" u2, -(x2-x1)/((x2-x1)^2 + (y2-y1)^2)^(3/2) - (x2-x3)/((x2-x3)^2 + (y2-y3)^2)^(3/2),\n",
" v2, -(y2-y1)/((x2-x1)^2 + (y2-y1)^2)^(3/2) - (y2-y3)/((x2-x3)^2 + (y2-y3)^2)^(3/2),\n",
" u3, -(x3-x1)/((x3-x1)^2 + (y3-y1)^2)^(3/2) - (x3-x2)/((x3-x2)^2 + (y3-y2)^2)^(3/2),\n",
" v3, -(y3-y1)/((x3-x1)^2 + (y3-y1)^2)^(3/2) - (y3-y2)/((x3-x2)^2 + (y3-y2)^2)^(3/2)]\n",
" return r\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us test the function on the initial conditions."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"12-element Vector{Float64}:\n",
" 0.466203685\n",
" -1.2125054397049002\n",
" 0.43236573\n",
" 0.3038594099200009\n",
" 0.466203685\n",
" 1.2125054397049002\n",
" 0.43236573\n",
" -0.3038594099200009\n",
" -0.93240737\n",
" 0.0\n",
" -0.86473146\n",
" 0.0"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# the initial positions\n",
"ip1 = [0.97000436, -0.24308753]\n",
"ip2 = [-ip1[1], -ip1[2]]\n",
"ip3 = [0, 0]\n",
"# the initial velocities\n",
"iv3 = [-0.93240737, -0.86473146]\n",
"iv2 = [-iv3[1]/2, -iv3[2]/2]\n",
"iv1 = iv2\n",
"# The initial conditions\n",
"init = [ip1[1], iv1[1], ip1[2], iv1[2],\n",
" ip2[1], iv2[1], ip2[2], iv2[2],\n",
" ip3[1], iv3[1], ip3[2], iv3[2]]\n",
"y = F(0.0, init)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"T = 2.1\n",
"h = 0.01\n",
"r = rk4(F, h, T, init)\n",
"t, x1, u1, y1, v1, x2, u2, y2, v2, x3, u3, y3, v3 = r;"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(x1,y1,label=\"1\")\n",
"plot!(x2,y2,label=\"2\")\n",
"plot!(x3,y3,label=\"3\")"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# savefig(\"figrk4celestial.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"
}
},
"nbformat": 4,
"nbformat_minor": 4
}