Miscellaneous updates and fixups for docs

docs/src/examples/beeler_reuter.md

@@ -8,19 +8,23 @@ Note that this tutorial does not use the [higher order IMEX methods built into D
 
 There are hundreds of ionic models that describe cardiac electrical activity in various degrees of detail. Most are based on the classic [Hodgkin-Huxley model](https://en.wikipedia.org/wiki/Hodgkin%E2%80%93Huxley_model) and define the time-evolution of different state variables in the form of nonlinear first-order ODEs. The state vector for these models includes the transmembrane potential, gating variables, and ionic concentrations. The coupling between cells is through the transmembrane potential only and is described as a reaction-diffusion equation, which is a parabolic PDE,
 
-$$\partial V / \partial t = \nabla (D  \nabla V) - \frac {I_\text{ion}} {C_m},$$
+```math
+\frac{\partial V}{\partial t} = \nabla (D  \nabla V) - \frac {I_\text{ion}} {C_m},
+```
 
-where $V$ is the transmembrane potential, $D$ is a diffusion tensor, $I_\text{ion}$ is the sum of the transmembrane currents and is calculated from the ODEs, and $C_m$ is the membrane capacitance, usually assumed to be constant. Here, we model a uniform and isotropic medium. Therefore, the model can be simplified to,
+where ``V`` is the transmembrane potential, ``D`` is a diffusion tensor, ``I_\text{ion}`` is the sum of the transmembrane currents and is calculated from the ODEs, and ``C_m`` is the membrane capacitance, usually assumed to be constant. Here, we model a uniform and isotropic medium. Therefore, the model can be simplified to,
 
-$$\partial V / \partial t = D \Delta{V} - \frac {I_\text{ion}} {C_m},$$
+```math
+\frac{\partial V}{\partial t} = D \nabla^2 V - \frac {I_\text{ion}} {C_m},
+```
 
-where $D$ is now a scalar. By nature, these models have to deal with different time scales and are therefore classified as *stiff*. Commonly, they are solved using the explicit Euler method, typically with a closed form for the integration of the gating variables (the Rush-Larsen method, see below). We can also solve these problems using implicit or semi-implicit PDE solvers (e.g., the [Crank-Nicholson method](https://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method) combined with an iterative solver). Higher order explicit methods such as Runge-Kutta and linear multistep methods cannot overcome the stiffness and are not particularly helpful.
+where ``D`` is now a scalar. By nature, these models have to deal with different time scales and are therefore classified as *stiff*. Commonly, they are solved using the explicit Euler method, typically with a closed form for the integration of the gating variables (the Rush-Larsen method, see below). We can also solve these problems using implicit or semi-implicit PDE solvers (e.g., the [Crank-Nicholson method](https://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method) combined with an iterative solver). Higher order explicit methods such as Runge-Kutta and linear multistep methods cannot overcome the stiffness and are not particularly helpful.
 
 In this tutorial, we first develop a CPU-only IMEX solver and then show how to move the explicit part to a GPU.
 
 ### The Beeler-Reuter Model
 
-We have chosen the [Beeler-Reuter ventricular ionic model](https://pmc.ncbi.nlm.nih.gov/articles/PMC1283659/) as our example. It is a classic model first described in 1977 and is used as a base for many other ionic models. It has eight state variables, which makes it complicated enough to be interesting without obscuring the main points of the exercise. The eight state variables are: the transmembrane potential ($V$), sodium-channel activation and inactivation gates ($m$ and $h$, similar to the Hodgkin-Huxley model), with an additional slow inactivation gate ($j$), calcium-channel activation and deactivations gates ($d$ and $f$), a time-dependent inward-rectifying potassium current gate ($x_1$), and intracellular calcium concentration ($c$). There are four currents: a sodium current ($i_{Na}$), a calcium current ($i_{Ca}$), and two potassium currents, one time-dependent ($i_{x_1}$) and one background time-independent ($i_{K_1}$).
+We have chosen the [Beeler-Reuter ventricular ionic model](https://pmc.ncbi.nlm.nih.gov/articles/PMC1283659/) as our example. It is a classic model first described in 1977 and is used as a base for many other ionic models. It has eight state variables, which makes it complicated enough to be interesting without obscuring the main points of the exercise. The eight state variables are: the transmembrane potential (``V``), sodium-channel activation and inactivation gates (``m`` and ``h``, similar to the Hodgkin-Huxley model), with an additional slow inactivation gate (``j``), calcium-channel activation and deactivations gates (``d`` and ``f``), a time-dependent inward-rectifying potassium current gate (``x_1``), and intracellular calcium concentration (``c``). There are four currents: a sodium current (``i_\text{Na}``), a calcium current (``i_\text{Ca}``), and two potassium currents, one time-dependent (``i_{x_1}``) and one background time-independent (``i_{K_1}``).
 
 ## CPU-Only Beeler-Reuter Solver
 
@@ -137,27 +141,39 @@ We use an explicit solver for all the state variables except for the transmembra
 
 The [Rush-Larsen](https://ieeexplore.ieee.org/document/4122859/) method replaces the explicit Euler integration for the gating variables with direct integration. The starting point is the general ODE for the gating variables in Hodgkin-Huxley style ODEs,
 
-$$\frac{dg}{dt} = \alpha(V) (1 - g) - \beta(V) g$$
+```math
+\frac{dg}{dt} = (1 - g) \alpha(V) - g \beta(V)
+```
 
-where $g$ is a generic gating variable, ranging from 0 to 1, and $\alpha$ and $\beta$ are reaction rates. This equation can be written as,
+where ``g`` is a generic gating variable, ranging from 0 to 1, and ``α`` and ``β`` are reaction rates. This equation can be written as,
 
-$$\frac{dg}{dt} = (g_{\infty} - g) / \tau_g,$$
+```math
+\frac{dg}{dt} = \frac{g_{\infty} - g}{\tau_g},
+```
 
-where $g_\infty$ and $\tau_g$ are
+where ``g_∞`` and ``τ_g`` are
 
-$$g_{\infty} = \frac{\alpha}{(\alpha + \beta)},$$
+```math
+g_{\infty} = \frac{\alpha}{\alpha + \beta},
+```
 
 and,
 
-$$\tau_g = \frac{1}{(\alpha + \beta)}.$$
+```math
+\tau_g = \frac{1}{\alpha + \beta}.
+```
 
-Assuming that $g_\infty$ and $\tau_g$ are constant for the duration of a single time step ($\Delta{t}$), which is a reasonable assumption for most cardiac models, we can integrate directly to have,
+Assuming that ``g_∞`` and ``τ_g`` are constant for the duration of a single time step (``Δt``), which is a reasonable assumption for most cardiac models, we can integrate directly to have,
 
-$$g(t + \Delta{t}) = g_{\infty} - \left(g_{\infty} - g(\Delta{t})\right)\,e^{-\Delta{t}/\tau_g}.$$
+```math
+g(t + Δt) = g_∞ - \left(g_∞ - g(Δt)\right)\,e^{-Δt/τ_g}.
+```
 
-This is the Rush-Larsen technique. Note that as $\Delta{t} \rightarrow 0$, this equation morphs into the explicit Euler formula,
+This is the Rush-Larsen technique. Note that as ``Δt → 0``, this equation morphs into the explicit Euler formula,
 
-$$g(t + \Delta{t}) = g(t) + \Delta{t}\frac{dg}{dt}.$$
+```math
+g(t + Δt) = g(t) + Δt \frac{dg}{dt}.
+```
 
 `rush_larsen` is a helper function that use the Rush-Larsen method to integrate the gating variables.
 
@@ -169,7 +185,7 @@ $$g(t + \Delta{t}) = g(t) + \Delta{t}\frac{dg}{dt}.$$
 end
 ```
 
-The gating variables are updated as below. The details of how to calculate $\alpha$ and $\beta$ are based on the Beeler-Reuter model and not of direct interest to this tutorial.
+The gating variables are updated as below. The details of how to calculate ``α`` and ``β`` are based on the Beeler-Reuter model and not of direct interest to this tutorial.
 
 ```julia
 function update_M_cpu(g, v, Δt)
@@ -225,8 +241,9 @@ end
 
 ### Implicit Solver
 
-Now, it is time to define the derivative function as an associated function of **BeelerReuterCpu**. We plan to use the CVODE_BDF solver as our implicit portion. Similar to other iterative methods, it calls the deriv function with the same $t$ multiple times. For example, these are consecutive $t$s from a representative run:
+Now, it is time to define the derivative function as an associated function of **BeelerReuterCpu**. We plan to use the `CVODE_BDF` solver as our implicit portion. Similar to other iterative methods, it calls the deriv function with the same ``t`` multiple times. For example, these are consecutive ``t``s from a representative run:
 
+```
 0.86830
 0.86830
 0.85485
@@ -240,8 +257,9 @@ Now, it is time to define the derivative function as an associated function of *
 0.87233
 0.88598
 ...
+```
 
-Here, every time step is called three times. We distinguish between two types of calls to the deriv function. When $t$ changes, the gating variables are updated by calling `update_gates_cpu`:
+Here, every time step is called three times. We distinguish between two types of calls to the deriv function. When ``t`` changes, the gating variables are updated by calling `update_gates_cpu`:
 
 ```julia
 function update_gates_cpu(u, XI, M, H, J, D, F, C, Δt)
@@ -265,7 +283,7 @@ function update_gates_cpu(u, XI, M, H, J, D, F, C, Δt)
 end
 ```
 
-On the other hand, du is updated at each time step, since it is independent of $\Delta{t}$.
+On the other hand, du is updated at each time step, since it is independent of ``Δt``.
 
 ```julia
 # iK1 is the inward-rectifying potassium current
@@ -367,7 +385,7 @@ deriv_cpu = BeelerReuterCpu(u0, 1.0);
 prob = DE.ODEProblem(deriv_cpu, u0, (0.0, 50.0));
 ```
 
-For stiff reaction-diffusion equations, CVODE_BDF from Sundial library is an excellent solver.
+For stiff reaction-diffusion equations, `CVODE_BDF` from Sundial library is an excellent solver.
 
 ```julia
 @time sol = DE.solve(prob, Sundials.CVODE_BDF(linear_solver = :GMRES), saveat = 100.0);
@@ -405,7 +423,7 @@ Some libraries, such as [ArrayFire](https://github.com/arrayfire/arrayfire), hid
 
 The key to fast CUDA programs is to minimize CPU/GPU memory transfers and global memory accesses. The implicit solver is currently CPU only, but it only requires access to the transmembrane potential. The rest of state variables reside on the GPU memory.
 
-We modify ``BeelerReuterCpu`` into ``BeelerReuterGpu`` by defining the state variables as *CuArray*s instead of standard Julia *Array*s. The name of each variable defined on the GPU is prefixed by *d_* for clarity. Note that $\Delta{v}$ is a temporary storage for the Laplacian and stays on the CPU side.
+We modify `BeelerReuterCpu` into `BeelerReuterGpu` by defining the state variables as *CuArray*s instead of standard Julia *Array*s. The name of each variable defined on the GPU is prefixed by *d_* for clarity. Note that ``Δv`` is a temporary storage for the Laplacian and stays on the CPU side.
 
 ```julia
 import CUDA
@@ -556,7 +574,7 @@ A CUDA program does not directly deal with GPCs and SMs. The logical view of a C
 
 Each thread can find its logical coordinate by using few pre-defined indexing variables (*threadIdx*, *blockIdx*, *blockDim* and *gridDim*) in C/C++ and the corresponding functions (e.g., `threadIdx()`) in Julia. Their variables and functions are defined automatically for each thread and may return a different value depending on the calling thread. The return value of these functions is a 1-, 2-, or 3-dimensional structure whose elements can be accessed as `.x`, `.y`, and `.z` (for a 1-dimensional case, `.x` reports the actual index and `.y` and `.z` simply return 1). For example, if we deploy a kernel in 128 blocks and with 256 threads per block, each thread will see
 
-```
+```c
     gridDim.x = 128;
     blockDim=256;
 ```
@@ -565,13 +583,13 @@ while `blockIdx.x` ranges from 0 to 127 in C/C++ and 1 to 128 in Julia. Similarl
 
 A C/C++ thread can calculate its index as
 
-```
+```c
     int idx = blockDim.x * blockIdx.x + threadIdx.x;
 ```
 
 In Julia, we have to take into account base 1. Therefore, we use the following formula
 
-```
+```julia
     idx = (blockIdx().x-UInt32(1)) * blockDim().x + threadIdx().x
 ```
 

docs/src/examples/classical_physics.md

@@ -6,9 +6,11 @@ If you're getting some cold feet to jump in to DiffEq land, here are some handcr
 
 #### Radioactive Decay of Carbon-14
 
-$$f(t,u) = \frac{du}{dt}$$
+```math
+\frac{du}{dt} = f(t,u) = - λ u
+```
 
-The Radioactive decay problem is the first order linear ODE problem of an exponential with a negative coefficient, which represents the half-life of the process in question. Should the coefficient be positive, this would represent a population growth equation.
+The Radioactive decay problem is the first order linear ODE problem of an exponential with a negative coefficient, which represents the half-life of the process in question. Should the coefficient be positive, this would represent a population growth equation. ``λ`` is known as the decay rate, and can be related to the half-life as ``λ = \ln(2)/t_{1/2}``.
 
 ```@example physics
 import OrdinaryDiffEq as ODE, Plots
@@ -41,21 +43,27 @@ Plots.plot!(sol.t, t -> 2^(-t / t½), lw = 3, ls = :dash, label = "Analytical So
 
 Another classical example is the harmonic oscillator, given by:
 
-$$\ddot{x} + \omega^2 x = 0$$
+```math
+\ddot{x} + \omega^2 x = 0
+```
 
 with the known analytical solution
 
-$$\begin{align*}
+```math
+\begin{align*}
 x(t) &= A\cos(\omega t - \phi) \\
 v(t) &= -A\omega\sin(\omega t - \phi),
-\end{align*}$$
+\end{align*}
+```
 
 where
 
-$$A = \sqrt{c_1 + c_2} \qquad\text{and}\qquad \tan \phi = \frac{c_2}{c_1}$$
+```math
+A = \sqrt{c_1 + c_2} \qquad\text{and}\qquad \tan \phi = \frac{c_2}{c_1}
+```
 
-with $c_1, c_2$ constants determined by the initial conditions such that
-$c_1$ is the initial position and $\omega c_2$ is the initial velocity.
+with ``c_1``, ``c_2`` constants determined by the initial conditions such that
+``c_1`` is the initial position and ``\omega c_2`` is the initial velocity.
 
 Instead of transforming this to a system of ODEs to solve with `ODEProblem`,
 we can use `SecondOrderODEProblem` as follows.
@@ -98,22 +106,28 @@ Thus, if we want the first series to be `x`, we have to flip the order with `var
 
 #### Simple Pendulum
 
-We will start by solving the pendulum problem. In the physics class, we often solve this problem by small angle approximation, i.e. $ sin(\theta) \approx \theta$, because otherwise, we get an elliptic integral which doesn't have an analytic solution. The linearized form is
+We will start by solving the pendulum problem. In the physics class, we often solve this problem by small angle approximation, i.e. ``\sin(θ) \approx θ``, because otherwise, we get an elliptic integral which doesn't have an analytic solution. The linearized form is
 
-$$\ddot{\theta} + \frac{g}{L}{\theta} = 0$$
+```math
+\ddot{θ} + \frac{g}{L} θ = 0
+```
 
 But we have numerical ODE solvers! Why not solve the *real* pendulum?
 
-$$\ddot{\theta} + \frac{g}{L}{\sin(\theta)} = 0$$
+```math
+\ddot{θ} + \frac{g}{L} \sin(θ) = 0
+```
 
 Notice that now we have a second order ODE.
 In order to use the same method as above, we need to transform it into a system
-of first order ODEs by employing the notation $d\theta = \dot{\theta}$.
+of first order ODEs by employing the notation ``ω(t) = \dot{θ}``.
 
-$$\begin{align*}
-&\dot{\theta} = d{\theta} \\
-&\dot{d\theta} = - \frac{g}{L}{\sin(\theta)}
-\end{align*}$$
+```math
+\begin{align*}
+\dot{θ} &= ω \\
+\dot{ω} &= - \frac{g}{L} \sin(θ)
+\end{align*}
+```
 
 ```@example physics
 # Simple Pendulum Problem
@@ -129,9 +143,8 @@ tspan = (0.0, 6.3)
 
 #Define the problem
 function simplependulum(du, u, p, t)
-    θ = u[1]
-    dθ = u[2]
-    du[1] = dθ
+    θ, ω = u
+    du[1] = ω
     du[2] = -(g / L) * sin(θ)
 end
 
@@ -167,16 +180,17 @@ Plots.plot(p, xlims = (-9, 9))
 A more complicated example is given by the double pendulum. The equations governing
 its motion are given by the following (taken from this [Stack Overflow question](https://mathematica.stackexchange.com/questions/40122/help-to-plot-poincar%C3%A9-section-for-double-pendulum))
 
-$$\frac{d}{dt}
-\begin{pmatrix}
-\alpha \\ l_\alpha \\ \beta \\ l_\beta
-\end{pmatrix}=
+```math
+\frac{d}{dt}
+\begin{pmatrix} \alpha \\ l_\alpha \\ \beta \\ l_\beta \end{pmatrix}
+=
 \begin{pmatrix}
 2\frac{l_\alpha - (1+\cos\beta)l_\beta}{3-\cos 2\beta} \\
 -2\sin\alpha - \sin(\alpha + \beta) \\
 2\frac{-(1+\cos\beta)l_\alpha + (3+2\cos\beta)l_\beta}{3-\cos2\beta}\\
 -\sin(\alpha+\beta) - 2\sin(\beta)\frac{(l_\alpha-l_\beta)l_\beta}{3-\cos2\beta} + 2\sin(2\beta)\frac{l_\alpha^2-2(1+\cos\beta)l_\alpha l_\beta + (3+2\cos\beta)l_\beta^2}{(3-\cos2\beta)^2}
-\end{pmatrix}$$
+\end{pmatrix}
+```
 
 ```@example physics
 #Double Pendulum Problem
@@ -197,7 +211,7 @@ function polar2cart(sol; dt = 0.02, l1 = L₁, l2 = L₂, vars = (2, 4))
     x1 = l1 * sin.(p1)
     y1 = l1 * -cos.(p1)
     (u, (x1 + l2 * sin.(p2),
-        y1 - l2 * cos.(p2)))
+         y1 - l2 * cos.(p2)))
 end
 
 #Define the Problem
@@ -231,7 +245,7 @@ Instead of looking at the full phase space, we can look at Poincaré sections,
 which are sections through a higher-dimensional phase space diagram.
 This helps to understand the dynamics of interactions and is wonderfully pretty.
 
-The Poincaré section in this is given by the collection of $(β,l_β)$ when $α=0$ and $\frac{dα}{dt}>0$.
+The Poincaré section in this is given by the collection of ``(β,l_β)`` when ``α=0`` and ``\frac{dα}{dt}>0``.
 
 ```@example physics
 #Constants and setup
@@ -241,10 +255,7 @@ tspan2 = (0.0, 500.0)
 
 #Define the problem
 function double_pendulum_hamiltonian(udot, u, p, t)
-    α = u[1]
-    lα = u[2]
-    β = u[3]
-    lβ = u[4]
+    α, lα, β, lβ = u
     udot .= [2(lα - (1 + cos(β))lβ) / (3 - cos(2β)),
         -2sin(α) - sin(α + β),
         2(-(1 + cos(β))lα + (3 + 2cos(β))lβ) / (3 - cos(2β)),
@@ -284,22 +295,30 @@ Plots.plot(p, xlabel = "\\beta", ylabel = "l_\\beta", ylims = (0, 0.03))
 
 The Hénon-Heiles potential occurs when non-linear motion of a star around a galactic center, with the motion restricted to a plane.
 
-$$\begin{align}
+```math
+\begin{align*}
 \frac{d^2x}{dt^2}&=-\frac{\partial V}{\partial x}\\
 \frac{d^2y}{dt^2}&=-\frac{\partial V}{\partial y}
-\end{align}$$
+\end{align*}
+```
 
 where
 
-$$V(x,y)={\frac {1}{2}}(x^{2}+y^{2})+\lambda \left(x^{2}y-{\frac {y^{3}}{3}}\right).$$
+```math
+V(x,y) = \frac {1}{2} (x^2+y^2) + λ \left(x^2 y - \frac{y^3}{3}\right).
+```
 
-We pick $\lambda=1$ in this case, so
+We pick ``λ=1`` in this case, so
 
-$$V(x,y) = \frac{1}{2}(x^2+y^2+2x^2y-\frac{2}{3}y^3).$$
+```math
+V(x,y) = \frac{1}{2} \left(x^2 + y^2 + 2 x^2 y - \frac{2}{3} y^3\right).
+```
 
 Then the total energy of the system can be expressed by
 
-$$E = T+V = V(x,y)+\frac{1}{2}(\dot{x}^2+\dot{y}^2).$$
+```math
+E = T+V = V(x,y) + \frac{1}{2} \left(\dot{x}^2+\dot{y}^2\right).
+```
 
 The total energy should conserve as this system evolves.
 
@@ -317,10 +336,7 @@ E(x, y, dx, dy) = V(x, y) + 1 // 2 * (dx^2 + dy^2);
 
 #Define the function
 function Hénon_Heiles(du, u, p, t)
-    x = u[1]
-    y = u[2]
-    dx = u[3]
-    dy = u[4]
+    x, y, dx, dy = u
     du[1] = dx
     du[2] = dy
     du[3] = -x - 2x * y