{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# ODEs 3\n",
"* System of equations\n",
"* Higher order derivatives\n",
"* Decoupled ODEs\n",
"* Adaptive time steps"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## ODE system\n",
"* So far we have considered $dy/dt = f(y,t)$, with one equation in one variable.\n",
"* For a system of ODEs, we have $d\\vec{y}/{dt} = \\vec{f}(\\vec{y},t).$\n",
" * For example, for 2 ODEs in 2 variables:\n",
" $$ \\frac{dx}{dt} = g(x,z,t), $$\n",
" $$ \\frac{dz}{dt} = h(x,z,t).$$\n",
" * We can write this as \n",
" $$ \\frac{d\\vec{y}}{dt} = \\vec{f}(\\vec{y},t),$$\n",
" \n",
" where $x=y[0]$, $z=y[1]$, $g=f[0]$ and $h=f[1]$.\n",
" * Note that each rate depends (in general) on all the variables.\n",
" * Note, below, we'll leave off the vector arrow symbols for simplicity.\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Explicit solution is a straightforward extension of the one equation case. Python's array functionality even allows nearly identical codes for systems of equations and for one equation.\n",
"* Implicit solutions require solution of a linear system of equations at each step for linear ODE systems. For nonlinear ODE systems, a nonlinear system must be solved at each step.\n",
" * Linear:\n",
" $$f(y) = Ay + b,$$\n",
" $$ y_{k+1} = y_k + \\Delta t(Ay_{k+1}+b),$$\n",
" $$(I-\\Delta tA)y_{k+1} = (y_k + \\Delta tb).$$\n",
" This last equation has the form $By_{k+1}=c$, which is a linear system solved for $y_{k+1}$ at each step.\n",
" * Note, $A$ and $b$ can depend on time, which is not explicitly shown.\n",
" * Nonlinear:\n",
" $$y_{k+1} = y_k + \\Delta t f(y_{k+1},t).$$\n",
" * Rearrange and solve the following nonlinear system for the $y_{k+1}$ vector at each step:\n",
" $$ F(y_{k+1}) = y_{k+1}-y_k - \\Delta tf(y_{k+1},t) = 0.$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Example\n",
"The following reactions are given:\n",
"$$A + B \\rightarrow C,$$\n",
"$$A + C \\rightarrow D.$$\n",
"* Reactions have rate constants $k_1 = 1$ and $k_2 = 2$.\n",
"* Let $A_0=B_0=1$ and $C_0=D_0=0$.\n",
"* Solve to 5 seconds.\n",
"* Species concentrations are given by the following rate equations:\n",
"\\begin{align*}\n",
"\\frac{dA}{dt} &= -k_1AB - k_2AC, \\\\\n",
"\\frac{dB}{dt} &= -k_1AB, \\\\\n",
"\\frac{dC}{dt} &= k_1AB - k_2AC, \\\\\n",
"\\frac{dD}{dt} &= k_2AC.\n",
"\\end{align*}\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Solve this system using the fourth order RK method:\n",
"\\begin{align*}\n",
"y_{k+1} &= y_k + h\\left(\\frac{1}{6}S_1 + \\frac{2}{6}S_2 + \\frac{2}{6}S_3 + \\frac{1}{6}S_4\\right),\\\\\n",
"& S_1 = f(y_k), \\\\\n",
"& S_2 = f(y_k+\\frac{h}{2}S_1), \\\\\n",
"& S_3 = f(y_k+\\frac{h}{2}S_2), \\\\\n",
"& S_4 = f(y_k+hS_3). \n",
"\\end{align*}"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"import Plots as plt"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"function odeRK4(f, y0, t)\n",
" ns = length(t)-1\n",
" y = zeros(length(t), length(y0))\n",
" y[1,:] = y0\n",
" \n",
" for k in 1:ns\n",
" h = t[k+1] - t[k]\n",
" S1 = f(y[k,:], t[k])\n",
" S2 = f(y[k,:]+0.5*h*S1, t[k]+0.5*h)\n",
" S3 = f(y[k,:]+0.5*h*S2, t[k]+0.5*h)\n",
" S4 = f(y[k,:]+ h*S3, t[k]+ h)\n",
" y[k+1,:] = y[k,:] + h/6*(S1 + 2*S2 + 2*S3 + S4)\n",
" end\n",
" \n",
" return y\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"function rhsf(ABCD, t)\n",
" A = ABCD[1]\n",
" B = ABCD[2]\n",
" C = ABCD[3]\n",
" D = ABCD[4]\n",
" \n",
" k1 = 1\n",
" k2 = 2\n",
" \n",
" dAdt = -k1*A*B - k2*A*C\n",
" dBdt = -k1*A*B\n",
" dCdt = k1*A*B - k2*A*C\n",
" dDdt = k2*A*C\n",
" \n",
" return [dAdt, dBdt, dCdt, dDdt]\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"ABCD_initial = [1,1,0,0]\n",
"tend = 5\n",
"t = LinRange(0,tend,100)\n",
"\n",
"ABCD = odeRK4(rhsf, ABCD_initial, t);"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plt.resetfontsizes(); plt.scalefontsizes(1.5)\n",
"plt.plot(t,ABCD, lw=2, label=[\"A\" \"B\" \"C\" \"D\"])\n",
"plt.plot!(foreground_color_legend=nothing)\n",
"plt.plot!(xlabel=\"t\", ylabel=\"Concentration\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Higher order derivatives\n",
"* A single higher order ODE results in a system of first order ODEs.\n",
"$$y^{\\prime\\prime\\prime}=\\frac{d^3y}{dt^3} = f(y,y^{\\prime},y^{\\prime\\prime},t).$$\n",
"\n",
"Question: can we solve this using the techniques we have used for first order ODE's? If so, how?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Can you convert this one third-order equation to a system of three first-order equations?\n",
"\n",
"$$y^{\\prime\\prime\\prime}=\\frac{d^3y}{dt^3} = f(y,y^{\\prime},y^{\\prime\\prime},t).$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* let $x=y^{\\prime\\prime}$, then $y^{\\prime\\prime\\prime}=x^{\\prime}$.\n",
"* let $z=y^{\\prime}$, then $y^{\\prime\\prime}=z^{\\prime}$.\n",
"* Then the resulting system of ODEs is:\n",
"\n",
"\\begin{align}\n",
"x^{\\prime} &= \\frac{dx}{dt} = f(y,z,x,t), \\\\\n",
"y^{\\prime} &= \\frac{dy}{dt} = z, \\\\\n",
"z^{\\prime} &= \\frac{dz}{dt} = x.\n",
"\\end{align}\n",
"\n",
"* This is a system of three first order ODEs in three variables. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Decoupled equations\n",
"\n",
"* Consider a linear system of ODEs.\n",
"$$y^{\\prime} = Ay+b.$$\n",
"\n",
"* In general, each rate equation depends on all the variables. \n",
"\n",
"* In our linear algebra review, we discussed how to decouple a system of linear equations.\n",
" * How did we do that?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* To decouple the system, change to an eigenvector basis so that each component equation depends only on its own variable (in the new basis) \n",
"* Let $V$ be a matrix whose columns are the eigenvectors of matrix $A$.\n",
"* Let $\\Lambda$ be a diagonal matrix whose diagonal elements are the eigenvalues of $A$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Then,\n",
"\n",
"$$AV = V\\Lambda,$$\n",
"$$A = V\\Lambda V^{-1}.$$\n",
"\n",
"* Insert this into the ODE:\n",
"$$y^{\\prime} = V\\Lambda V^{-1}y + b.$$\n",
"* Multiply through by $V^{-1}$:\n",
"$$(V^{-1}y^{\\prime}) = \\Lambda(V^{-1}y) + (V^{-1}b).$$\n",
"* Now, let $\\hat{y}=V^{-1}y$, and $\\hat{b}=V^{-1}b$:\n",
"$$\\hat{y}^{\\prime} = \\Lambda \\hat{y} + \\hat{b}.$$\n",
"* Because $\\Lambda$ is diagonal, this system is decoupled. That is component $i$ is given by\n",
"$$\\hat{y}^{\\prime}_i = \\lambda_i\\hat{y}_i + \\hat{b}_i.$$\n",
"* This equation has a simple analytic solution.\n",
" * When solved, all the $\\hat{y}_i(t)$ are known. \n",
" * Then $y(t)$ are given by\n",
" $$y(t) = V\\hat{y}(t).$$\n",
" \n",
"This analysis can be useful for solving ODEs analytically, but also for analyzing (and modifying) stability properties of ODEs."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## ODE Step Size\n",
"* See [Numerical Recipes](http://s3.amazonaws.com/nrbook.com/book_C210.html) section 16.2 for more detailed explanations.\n",
"* When solving an ODE, we needed a step size $\\Delta t$ or $h$.\n",
"* How should we select the step size?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Consider integrating $dy/dt = \\tanh(t)$\n",
"* What part of the solution will dictate the step size?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"* Do we have to use this step size everywhere?\n",
"* If not, how can we make the computer choose the stepsize in an \"intelligent\" way?\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Suppose we know the error $\\Delta$ for a given step size $\\Delta t=h_1$.\n",
"* We'll show how to get $\\Delta$ below.\n",
"\n",
"Suppose we set some desired error that we are okay with on a given step, like $|\\Delta|\\le\\epsilon = atol + |y|rtol.$ \n",
"* A large $y\\rightarrow$ rtol controls; \n",
"* A small $y\\rightarrow$ atol controls.\n",
"\n",
"Given a known error for a known step, how can we change our step to get the desired error? Assume we are using the RK4 method.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Adjust $h$ to get the desired error.\n",
"* For a globally $4^{th}$ order method $\\rightarrow$ $\\Delta = \\mathcal{O}(h^5)$ $\\rightarrow$ $\\Delta\\sim h^5$.\n",
"* Then \n",
"$$\\frac{\\Delta_2}{\\Delta_1} = \\frac{\\epsilon}{\\Delta_1} = \\left(\\frac{h_2}{h_1}\\right)^5,$$\n",
"\n",
"$$\\rightarrow h_2 = h_1\\left(\\frac{\\epsilon}{\\Delta_1}\\right)^{1/5}.$$\n",
"\n",
"* So, guess an initial $h_1$.\n",
" * If the error $\\Delta_1$ is too big, then redo the step using a smaller $h$ as computed using the above equation.\n",
" * If the error is too small, take the step, but do the *next* step with a larger $h$ computed using the equation above."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### How to compute $\\Delta$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"* See N.R.\n",
"* Two approaches:\n",
" 1. Step doubling\n",
" 2. Felberg.\n",
" \n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"#### Step doubling\n",
"\n",
"\n",
"* Consider two grids where we can either take two size $h$ steps, or one size $2h$ step:\n",
"``` \n",
"(A) <----h---->|<----h---->\n",
"(B) <---------2h---------->\n",
"```\n",
"* Now, let $\\Delta$ = $y_B-y_A$.\n",
" * Recall, if the error of the method (per step) is $\\mathcal{O}(h^5)$, then \n",
" \n",
" \\begin{align}\n",
" y_{exact} &= y_A + 2\\mathcal{O}(h^5), \\\\\n",
" y_{exact} &= y_B + \\mathcal{O}((2h)^5) = y_B + 32\\mathcal{O}(h^5).\n",
" \\end{align}\n",
" * The error in $y_{B}$ is 16 times larger than the error in $y_{A}$, so we consider $y_A$ to be exact (compared to $y_B$), and evaluate $\\Delta = y_B-y_A$. \n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* ***Cost*** \n",
" * Each RK step requires 4 function evaluations.\n",
" * 3 total steps (two for grid (A) and one for grid (B)) $\\rightarrow$ 12 function evaluations. Actually 11, since the grids share a starting point.\n",
" * We compare 11 required using step doubling to 8 required without step doubling.\n",
" * 37.5%=(11-8)/8 is the cost increase.\n",
" * This cost increase pays for itself in terms of allowing (ideally) at least 37.5% fewer overall steps due to the adaptive stepsize control.\n",
" * For a given error, without adaptive stepsize, the *whole* time domain would have to use the most stringent step size.\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"#### [Felberg](https://onlinelibrary.wiley.com/doi/abs/10.1002/zamm.19660460102)\n",
"\n",
"* You can write an RK method where one linear combination of slopes results in an $\\mathcal{O}(h^6)$ method, and another combination of the same slopes gives an $\\mathcal{O}(h^5)$ method.\n",
"\n",
"\\begin{align}\n",
"y_{k+1} &= y_k + h(a_1S_1 + a_2S_2 + a_3S_3 + a_4S_4 + a_5S_5 + a_6S_6 ) + \\mathcal{O}(h^6),\\\\\n",
"\\hat{y}_{k+1} &= y_k + h(b_1S_1 + b_2S_2 + b_3S_3 + b_4S_4 + b_5S_5 + b_6S_6 ) + \\mathcal{O}(h^5).\n",
"\\end{align}\n",
"\n",
"* Then let $\\Delta=y_{k+1}-\\hat{y}_{k+1}$.\n",
"\n",
"[See Hoffman](https://ignite.byu.edu/cbe541/lectures/felberg.pdf)"
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Julia 1.8.3",
"language": "julia",
"name": "julia-1.8"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.3"
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}