% two columns generally contains more content than 1 col, despite the small gap % in the middle, because line breaks, and unfull lines cost less here. \documentclass[a4paper,twocolumn,enumboxcolor=cyan!50,10pt]{onepgnote} \usepackage{amsmath,amssymb} \usepackage{physics} \usepackage[margin=10pt]{geometry} \begin{document} \pagestyle{empty} \section{Basics} \subsection{Index Notation} \begin{enumerate} \item Whenever a quantity is summed over an index which appears exactly twice in each term in the sum, we leave out the summation sign. Such an index is called a \emph{dummy index}. \item An index appearing only once is called a \emph{free index}. \item No index may appear three times in a term. But one index can appear multiple times in an equation. A term is the basic unit in index notation. \end{enumerate} For example, \[ x_{ik}y_{jk} + a_{ik}b_{ik} = a_i + b_j \quad \hbox{means} \quad \sum_k x_{ik}y_{jk} + \sum_{k}a_{ik}b_{ik} = a_i + b_j \] \subsection{Definitions} The \emph{Kronecker delta} $\delta_{ij} := 1$ iff $i=j$ else $0$. We define the \emph{Levi-Civita symbol} $\varepsilon_{ijk}$ for $1 \le i,j,k \le 3$ to be \begin{enumerate} \item $0$ iff $(i,j,k)$ is not a permutation of $(1,2,3)$. \item $1$ iff $(i,j,k) \in \{(1,2,3),\allowbreak (2,3,1),\allowbreak (3,1,2)\}$. \item $-1$ iff $(i,j,k) \in \{(1,3,2),\allowbreak (2,1,3),\allowbreak (3,2,1)\}$. \end{enumerate} That is, $\varepsilon_{ijk}$ equals to the parity of the permutation $(1,2,3) \to (i,j,k)$. Similarly, in 2D, $\varepsilon_{ij}$ is $1$ whenever $(i,j) = (1,2)$, $-1$ when $(i,j) = (2,1)$, and $0$ otherwise. \begin{enumerate} \item The \emph{double dot product} is $\vb{A} : \vb{B} := \sum_{i}\sum_{j}A_{ij}B_{ij}$ \label{doubledot}. \item When one operand of the dot product is a matrix, it is interpreted as matrix multiplication. \item If we apply gradiant to a vector-valued function $\vb{f}(\vb{x})$, then we are putting the gradiant of each component of it into a column of the result matrix: $\nabla \vb{f} := [(\nabla f_1)^T, \dots, (\nabla f_n)^T]$ \label{gradvec}. \end{enumerate} \begin{enumerate} \item The \emph{divergence} of $\vb{P}$, $\nabla \cdot \vb{P} := \sum_{i}\pdv{P_i}{x_i}$, as if $\nabla = (\pdv{}{x_1},\allowbreak \dots,\allowbreak \pdv{}{x_n})$. It represents the source/sink of a $\vb{v}$ field. \mynote Distinguish $\nabla \cdot \vb{P}$ vs $\nabla \vb{P}$ \ref{gradvec}! \item The \emph{Laplace operator (Laplacian)} is defined as the divergence of the gradiant of function $f$, $\nabla^2 f := \nabla \cdot \nabla f$. Rarely, we may use $\Delta$ for it. \item The \emph{curl} of $\vb{P}$ is defined as $\nabla \times \vb{P} = (\pdv{v_3}{x_2} - \pdv{v_2}{x_3}, -(\pdv{v_3}{x_1} - \pdv{v_1}{x_3}), \pdv{v_2}{x_1} - \pdv{v_1}{x_2})$. It represents the vorticity of a $\vb{v}$ field. \end{enumerate} \subsection{Facts} \begin{enumerate} \item The area of a parallellogram equals $|\vb{a} \times \vb{b}|$. \item The volume of a parallelpiped equals $|\vb{a} \cdot (\vb{b} \times \vb{c})| = \det(\vb{a}, \vb{b}, \vb{c})$. \item $\vb{a} \times \vb{b} = \varepsilon_{ijk}\vb{e}_ia_jb_k$. \item $\varepsilon_{ijk}\varepsilon_{imn} = \delta_{jm}\delta_{kn} - \delta_{jn}\delta_{km}$. \end{enumerate} \section{Change of basis} For vectors and matrices. \subsection{For vectors} Let $U$ be a finite $n$-dimensional vector space over field $F$. Let $(\vb{u}_i)_{i=1}^n$ and $(\vb{v}_i)_{i=1}^n$ be two ordered bases of $U$. Let $\vb{x}$ be any vector in $U$. Then, there exists two unique scalar sequences $(a_i)_{i=1}^n$ and $(b_i)_{i=1}^n$ such that $\sum_{i=1}^n a_i\vb{u}_i = \sum_{i=1}^n b_i\vb{v}_i$. They are also vectors in their own, in the vector space $F^n$. We want to find the change of basis function $T: (a_i)_{i=1}^n \mapsto (b_i)_{i=1}^n$, which exists since $(\vb{v}_i)_{i=1}^n$ is a basis. $T$ may be broken down into $T = f \circ g$ where $g((a_i)_{i=1}^n) := a_i\vb{u}_i$, and $f(\vb{x} \in U) := (b_i)_{i=1}^n$ such that $b_i\vb{v}_i = \vb{x}$. It is clear that both $f$ and $g$ are linear, and thus $T$ is a linear function $F^n \to F^n$. We then seek its unique matrix representation under the standard basis $(\vb{e}_i)_{i=1}^n$ of $F^n$: $[T] = [T(\vb{e}_1)^T, \cdots, T(\vb{e}_n)^T]$. \mynote When $U = F^n$ and if we happen to use the standard basis of $F^n$, then $\vb{x}$ will appear exactly the same as $(a_i)_{i=1}^n$, which can be extremely confusing. We must constantly remind ourselves that $T$ operates on the coordinate vectors and not directly on $\vb{x}$, which is the same vector regardless of the basis used to represent it as coordinates. \subsubsection{Change back} Since $(\vb{u}_i)_{i=1}^n$ is a basis (it spans $U$ and each $(a_i)_{i=1}^n$ is unique), $T$ is invertible. Then, indeed, the function $T^{-1}$ will be the change-back funcstion. As another way, one may swap the places of $(\vb{u}_i)$ and $(\vb{v}_i)$ in the above discussion. \subsection{For matrices} Let $f$ be a linear function $U \to U$. Under any basis $B = (\vb{u}_i)_{i=1}^n$, $f$ has a unique matrix representation $[f]_{B}$. It is interesting to us how the matrix changes when we use another basis $B' = (\vb{v}_i)_{i=1}^n$ as the coordinate frame to get $[f]_{B'}$. Let $\vb{x} \in U$ be any vector, and denote $f(\vb{x}) =: y$. Let $(a_i)_{i=1}^n$ be the coordinates of $\vb{x}$ under basis $B$, $(b_i)_{i=1}^n$ be the coordinates of $\vb{y}$ under $B$, and $(a'_i)_{i=1}^n, (b'_i)_{i=1}^n$ be the coordinates for them under $B'$. Now, consider the change-of-basis function $T$ changing the coordinate frame under $B$ to $B'$, then, since $f$ and $\vb{x}$ don't change for the basis, we should have \[ [f]_{B}T^{-1}(a'_i)_{i=1}^n = T^{-1}(b'_i)_{i=1}^n \quad \to \quad T[f]_{B}T^{-1}(a'_i)_{i=1}^n = (b'_i)_{i=1}^n \] Thus $[f]_{B'} = T[f]_{B}T^{-1}$. \mynote If we are talking about orthognormal bases, then $T$ will be an orthognormal matrix, resulting in $T^{-1} = T^T$. \subsection{Invariants} Formally, we define an \emph{invariant} to be any function on such matrices such that $f(M) = f(TMT^{-1})$ for all applicable matrices $M, T$. Three frequent invariants are \begin{enumerate} \item $\trace M$ \item $\det M$ \item $ M_{11}M_{22} + M_{22}M_{33} + M_{11}M_{33} - M^2_{12} - M^2_{23} - M^2_{31}$. \end{enumerate} \section{Stress} Stress describes the forces present during the deformation of a material. It expresses the internal forces that neighbouring particles of a continuous material exert on each other. \subsection{Cauchy stress tensor} Cauchy observed that the stress vector across a surface will always be a linear function of the surface's normal vector. The matrix for this function (under a fixed basis) can be treated as a tensor and called the Cauchy stress tensor. For whatever reason, people chose to call the value of the Cauchy stress tensor function a \emph{traction vector}, and use stress to refer to the matrix (tensor) of the function. \mynote In practice we assume the function takes normalized normals. Under the standard basis, by the principle of conservation of angular momentum, the matrix can be shown to be symmetric \label{stresssym}. Thus, the traction vector can be computed by either $[\sigma]\vb{n}$ or $[\sigma]^T\vb{n}$. In the lecture the second way is used. \subsection{Normal \& shear tractions} For a traction vector, one can divide it into two components, \begin{enumerate} \item \emph{Normal traction}, which is parallel to the surface normal $\vb{n}$, \item \emph{Shear traction}, which is perpendicular to the normal $\vb{n}$. \end{enumerate} One can easily calculate the normal traction $\vb{t}_n = (\vb{t} \cdot \vb{n}) \vb{n}$. Then, the shear traction is $\vb{t}_s = \vb{t} - \vb{t}_n$. \mynote We may abuse notation sometimes to refer to the magnitude of the traction as traction. \mynote We also define the change in angle between two normalized vectors under a stress tensor $\sigma$ to be (by symmetry), $2\vb{v}_1[\sigma]\vb{v}_2^T = 2\vb{v}_2[\sigma]\vb{v}_1^T$. \subsection{Infinitesimal strain} The infinitesimal strain theory is an idealized theory in which one can arrive at simple formulae for deformation. Under this, the spatial and material coordinates are approximately the same, and we have the infinitesimal displacement/strain tensor $\epsilon = 1/2 ((\nabla_{\vb{x}} \vb{u})^T + \nabla_{\vb{x}} \vb{u})$ \label{inften1} and the infinitesimal rotation tensor $\omega = 1/2(\nabla_{\vb{x}} \vb{u} - (\nabla_{\vb{x}} \vb{u})^T)$, where $\vb{u}(\vb{x}, \dots)$ is the displacement vector field, and $\nabla_{\vb{x}}\vb{u} = [(\nabla_{\vb{x}} u_1)^T, \cdots, (\nabla_{\vb{x}} u_n)^T]$ \ref{gradvec}. About $\epsilon$, an original location vector $\vb{x}$'s deformation will be described by it in such a way that $\vb{x}' = \epsilon \vb{x}$. We have \begin{enumerate} \item The diagonal elements of $\epsilon$ represent fractional length changes. E.g., if $\vb{x} \parallel \vb{e}_1$, then $\epsilon_{11} = (|\vb{x}'| - |\vb{x}|)/|\vb{x}|)$. \item Thus, $\trace \epsilon = \nabla \cdot \vb{u}$ is the fractional change in volume. \item Off-diagonal elements represent changes in angle. This is because, the angle, for unit vectors, is $\arccos$ of $\vb{x} \cdot \vb{x}' = \vb{x}\cdot(\epsilon \vb{x}) = \vb{x}\epsilon\vb{x}^T$. \item As such, $2\epsilon_{ij}, i\ne j$ is the change in angle between $\vb{e}_i$ and $\vb{e}_j$ after the deformation. Also, given $\vb{p} \perp \vb{q}$, $2\vb{q}\epsilon\vb{p}^T$ is the change in angle between them. \end{enumerate} \section{Material vs spatial} Suppose we are interested in some physical property $\vb{P}$ of some material in space. \begin{enumerate} \item In material (Lagrangian) specification, the observer's eyes follows a particular particle, and the property is a function of the particle's initial location $\xi$ and the time $t$: $\vb{P}(\xi, t)$. \item In spatial (Eulerian) specification, the observer does not follow any particle but instead gives a global description of the space, resulting in function $\vb{P}(\vb{x}, t)$, giving the property for the particle at location $\vb{x}$ at time $t$. \end{enumerate} \subsection{Their link} Suppose we are given a spatial description $\vb{P}(\vb{x}, t)$, and would like to use this to follow a specific particle to give a material description to it. Then, $\vb{x}$, for the particle, is a function: $\vb{x}(\xi, t)$, and $\vb{P}$ becomes $\vb{P}(\vb{x}(\xi, t), t)$. Specifically, if we want to find $\pdv{\vb{P}}{t}$, then we need to use the chain rule to get \[ \pdv{\vb{P}}{\vb{x}}\pdv{\vb{x}}{t} + \pdv{\vb{P}}{t}\pdv{t}{t} = \pdv{\vb{P}}{\vb{x}}\vb{v}(\xi, t) + \pdv{\vb{P}}{t} = \sum_{i}\pdv{\vb{P}}{x_i}v_i(\xi, t) + \pdv{\vb{P}}{t} \] \section{Equations} \subsection{Terms} \begin{enumerate} \item \emph{Steady state} means everything is constant w.r.t.\ time $t$. \item \emph{No flow} means velocity $\vb{v} = \vb{0}$. \item \emph{No strain, stress} means the strain, stress tensors $\sigma, \epsilon = \vb{0}$. \end{enumerate} \subsection{of mass} For the conservation of mass, \begin{enumerate} \item In spatial description, we have $\pdv{\rho}{t} + \nabla \cdot (\rho\vb{v}) = 0$. \item In material description, it becomes $\frac{D\rho}{Dt} + \rho(\nabla \cdot \vb{v}) = 0$. \item If the material is incompressible, i.e., $\rho$ is constant, then the equation is reduced to $\nabla \cdot \vb{v} = 0$. \end{enumerate} \subsection{of momentum} \begin{enumerate} \item The conservation of angular momentum results in $\sigma$'s being symmetric \ref{stresssym}. \item As for that of linear momentum \label{conlinmom}, we have $ \rho\left(\pdv[2]{\vb{u}}{t} = \pdv{\vb{v}}{t}\right) = \vb{f} + \nabla\cdot\sigma = \sum_{i}\Big(f_i + \sum_j \pdv{\sigma_{ji}}{x_j}\Big). $ \end{enumerate} \subsection{of energy} The general form of conservation of energy in the lecture is $ \frac{D\rho C_pT}{Dt} = \nabla\cdot k(\nabla T) + A + \sigma:\vb{D}. $ , where $:$ is \ref{doubledot} and the terms from the left to the right are \begin{enumerate} \item change in temperature with time \item heat transfer by conduction (and radiation) \item heat production (including latent heat) \item heat generated by internal deformation. \end{enumerate} \subsection{Rheology} We have $\hbox{rheology} \cdot \hbox{deformation}\ (\epsilon) = \hbox{stress}\ (\sigma)$. \subsubsection{Elasticity} \begin{enumerate} \item Elasticity means a material's deformation under a force will be restored when the force is removed. Under perfect elasticity, Hook's law states that the distance of deformation is proportional (linear) to the force applied: $\sigma = C : \epsilon$ \ref{inften1}. \item Since it is linear, elasticity of a material is quantified by the \emph{elastic modulus}, defined as $\delta := \frac{\hbox{stress}}{\hbox{strain}}$. \item Young's modulus $E := \frac{\sigma_{11}}{\epsilon_{11}}$c, and Poisson's ratio $\nu := -\frac{\epsilon_{33}}{\epsilon_{11}}$. \item In homogeneous and isotropic materials, Lam\'e's constants $\lambda, \mu$ define Hooke's law in 3D $\sigma = 2\mu\epsilon + \lambda \trace (\epsilon) I_{3\times3}$. \item With the Bulk ($K$) and shear ($G$) modulus: $-p = K\theta$ (isotropic) $\sigma'_{ij} = \frac{\sigma_{kk}}{3} = K \epsilon_{kk} = 2G\epsilon'_{ij}$ (deviatoric), where $\sigma_{ij} + p\delta_{ij} =: \sigma'_{ij}$. Thus, $K = \lambda + 2\mu/3$, where the second RHS term is compressional and the third is shear. \item In Newtonian, general compressible fluids, $\sigma = (-p+\zeta\Delta)\vb{I} + 2\eta\vb{D}$, where $D$ is the strain rate, $\Delta = D_{kk} = \nabla \cdot \vb{v}$. We have the Navier-Stokes equation $ \nabla p + (\zeta + \eta) \nabla \Delta + \eta \nabla^{2} \vb{v} + \vb{f} = \rho \frac{D\vb{v}}{Dt} = \rho (\frac{\partial{\vb{v}}}{\partial{t}} + \vb{v} \cdot \nabla \vb{v}) $ If the fluid is incompressible, then $\Delta = \nabla \cdot \vb{v} = 0$ and it's simplified to $\sigma = -pI+2\eta D$, and we have the Navier-Stokes equation $ -\nabla p + \eta \nabla^{2} \vb{v} + \vb{f} = \rho (\frac{\partial{\vb{v}}}{\partial{t}} + \vb{v} \cdot \nabla \vb{v}) $. \end{enumerate} \subsubsection{Wave equation} Substituting Lam\'e's constants formula into the equation of conservation of linear momentum \ref{conlinmom}, we have $ \rho\pdv[2]{\vb{u}}{t} = \vb{f} + (\lambda + 2\mu)\nabla(\nabla\cdot\vb{u})-\mu\nabla\times\nabla\times\vb{u} $. \section{Interpolation} Let $(x_i, y_i)_{i=0}^{N}$ be $N+1$ data points. It can be shown that, provided $\forall i,j, x_i \ne x_j$, $\{(x_i^n)_{n=0}^N\}_{i=0}^{N}$ is linearly independent. Thus, the linear system $X\vb{a} = \vb{y}$, where $X$ has these vectors as rows and $\vb{y} = (y_i)_{i=0}^N$, has a unique solution, $\vb{a} = (a_n)_{n=0}^N$ which is the coefficient vector of the unique polynomial of degree $N$ passing through them. \subsection{Lagrange} Let $\mathcal{P}_N$ be the set of all polynomials of degree $N$ or less. It is a $(N+1)$-dimensional vector space. Lagrange found a basis $\{\ell_i(x)\}_{i=0}^N$, where $\ell_i(x) := \prod_{m \ne i}\frac{x-x_m}{x_i-x_m}$, and showed that $L(x) := \sum_{i=0}^Ny_i\ell_i(x)$ is the unique interpolating polynomial. \subsubsection{Remainder} Lagrange proved that, for any $f \in C^{N+1}[a,b]$, and datapoints $(x_i, y_i)_{i=0}^N$ that $f$ passes through, the unique interpolation polynomial $P_N(x)$ results in a remainder, $R(x) := f(x) - P_N(x)$, satisfying $\forall x \in [a,b],\ \exists c \in [a,b],\ R(x) = \Psi(x)f^{N+1}(c)/\allowbreak (N+1)!$, where $\Psi(x) := \prod_{i=0}^N(x-x_i)$. \mynote Thus, $\Psi(x) \cdot\allowbreak \max(f^{N+1}(c))/\allowbreak (N+1)!$ is an upper bound of it. \subsection{Variants} Observe that $\Psi(x)$ is much larger around $a$ or $b$ for evenly spaced datapoints. This, plus potentially large $f^{N+1}(c)$, gives very unstable results near the boundary. \begin{enumerate} \item Usually the choice of $x_i$ is the only thing we can control. Lanczos found that $\max_{x_i \in [-1,1]}\Psi(x)$ attains the minimum when $x_i$ are the roots of the Chebyshev polynomial $T_{N+1}(x)$, $x_i = \cos(\frac{2i-1}{2N}\pi)$. \item We may also interpolate $f$ piecewise to decrease the error, although it will make the interpolation function lose some smoothness. \end{enumerate} \section{Quadrature} \begin{enumerate} \item Midpoint $I_M := \sum_{i=0}^{n-1}\discretionary{}{}{} f \left ( \frac {x_{i+1}+x_i} {2} \right )\discretionary{}{}{} (x_{i+1}-x_i)$. \item Trapezoidal $\sum_{i=0}^{n-1}\discretionary{}{}{} \left(\frac{f(x_{i+1}) + f(x_{i})}{2}\right )\discretionary{}{}{} (x_{i+1}-x_i)$. \item Simpson's $I_S := \frac{2}{3}I_M + \frac{1}{3}I_T = \sum_i \frac{(x_{i+1}-x_i)}{6}(f(x_i)+4f(m)+f(x_{i+1}))$. \item Weddle's $I_W = I_{S_2} + \frac {\left (I_{S_2} - I_S \right )}{15}$, where $I_{S_2}$ is Simpson's applied on double amount of intervals. \item Composite trapezoidal $\frac{\Delta x}{2}[f(x_0) + 2f(x_1) + \dots + 2f(x_{n-1}) + f(x_n)]$. \item Composite Simpson's $\frac{\Delta x}{3}\left[ f \left ( x_0\right ) + 2\sum_{i=1}^{n/2 - 1} f\left(x_{2i}\right) + 4\sum_{i=1}^{n/2} f\left(x_{2i-1}\right) + f \left ( x_{n}\right ) \right]$. \end{enumerate} \mynote Composite Simpson actually uses 2 intervals as one, and the middle point is thus some $x_i$. \subsection{Error} All these rules can be regarded as doing piecewise polynomial interpolation on $f$ on evenly spaced datapoints, integrating $f$ minus the polynomial, summing over the intervals, and using the Lagrange remainder to find the error. We may find \begin{enumerate} \item Trapezoidal: $-\frac{1}{12} \Delta x^2 (b-a) \frac{1}{n} \sum_{i=0}^{n-1}\, f''\left(c_{x_i}\right)$ \item Midpoint: $\frac{1}{24} \Delta x^2 (b-a) \frac{1}{n} \sum_{i=0}^{n-1}\, f''\left(c_{x_i}\right)$. \item Simpson's: since the error $I-I_T \approx -2(I-I_M)$, we imagine there's a better approximation $I_S$ such that $I_S - I_T = -2(I-I_M)$, giving $I_S = \frac{2}{3}I_M + \frac{1}{3}I_T$. We have error $-\frac{\Delta x^4}{180} (b-a)f^{(4)}\left(c_x\right)$. \item Weddle's: it would be a fuss to derive the exact one but it should be proportional to $\Delta x^6$. \end{enumerate} \section{ODE} To approximately solve (satisfying convergence as $\Delta t \to 0$ and correct qualitative behaviour) $y'(t) = f(t, y(t))$, we have \begin{enumerate} \item By the definition of derivative or the Taylor series, we have $y_{n+1} \approx y_{n} + \Delta t y'(t_n)$. This the the forward Euler. \item By the definition of derivative, we may also say $y_{n+1} \approx y_n + \Delta t y'(t_{n+1})$, which is the backward Euler. \item We may also say $y'(t_n) \approx \frac{y_{n+1}-y_{n-1}}{2\Delta t}$ and as such $y_{n+1} \approx y_{n-1} + 2\Delta t y'(t_n)$, which is leapfrog. \item Recall the trapezoidal rule before, we may use it here to average the forward and backward Euler to obtain $y_{n+1} \approx y_{n} + \Delta t\frac{y'(t_n) + y'(t_{n+1})}{2}$. \end{enumerate} \subsection{Error} \begin{enumerate} \item The \emph{local error (LE)} is error committed at one step, assuming the previous step $y_{n}$ is exact. Thus it is $y_{n+1} - y_{n+1}'$. \item For example, using Taylor series, forward Euler has local error $(\frac{\Delta t^2}{2!}y''_n + \dots)$. \item The \emph{(local) truncation error (LTE)}, $\tau$, is defined by the local error scaled down: $\tau := LE / (\Delta t)$. \item A method is \emph{consistent} if $\lim_{\Delta t \to 0}\tau = 0$. \item The \emph{global error} is $E := \max_{t_0 \, \le \, t_n \, \le \, T}\| y_{n} - y(t_{n})\|$, only assuming the inital $y_0$ is exact. \item A method \emph{converges} iff $\lim_{\Delta t \to 0} E = 0$. \end{enumerate} \subsection{Stability} Stability in numerical methods of solving ODEs have different definitions, but in general we would want the numerial methods to exhibit the same important properties as the exact solution. We have \begin{enumerate} \item A numerical method is said to be \emph{A-stable}, if, when applied to the reference equation $y'=ky \land y(0) = 1$ for $k \in \mathbf{C}$, demonstrates that, provided $\real(k) < 0$, $\lim_{t \to \infty} \hbox{solution} \to 0$, the same property from the exact solution $y(t) = e^{kt}$. \item A numerical method is \emph{L-stable}, if it is A-stable, and that its stability function $\phi(x) \to 0$ as $z \to \infty$. \end{enumerate} \subsection{Adams-Bashforth} For general ODE of form $y'(t) = f(y(t), t)$, according to the fundamental theorem of calculus, $y_{n+1} - y_n = \int_{t_n}^{t_{n+1}}f(y(t),t) \dd{t}$. The A-B schemes uses combinations of $f_{i}$ for $k$ many $i$'s with $i \ne n+1$ to approximate the RHS integral to numerically solve the ODE. We have these AB schemes: \begin{enumerate} \item $k=0$-step: $y_{n+1} = y_n + \Delta t f_n$. \item $k=1$-step: $y_{n+1} = y_{n} +\frac{\Delta t}{2}\; (f_{n} + f_{n+1})$. \item $k=2$-step: $y_{n+2} = y_{n+1} + \frac{\Delta t}{12}\; (-f_{n} + 8 f_{n+1} + 5f_{n+2})$. \item $k=3$-step: $y_{n+3} = y_{n+2} + \frac{\Delta t}{24}\; (f_{n} - 5 f_{n+1} + 19 f_{n+2} + 9 f_{n+3})$. \item $k=4$-step: $y_{n+4} = y_{n+3} + \frac{\Delta t}{720}\; (-19 f_{n} + 106 f_{n+1} - 264 f_{n+2} + 646 f_{n+3} + 251 f_{n+4})$ \end{enumerate} \subsection{Runge-Kutta} Similar to the trapezoidal method, we have $y_{n+1} - y_n \approx \frac{1}{2}\allowbreak\Delta t\allowbreak(y'_n + y'_{n+1})$. However, here, $y'_{n+1} = f(y_{n+1},\allowbreak t_{n+1})$, where $y_{n+1}$ is not known. In Runge-Kutta 2-step (RK2), we first use forward Euler to get $y^* \approx y_{n+1}$, then use this to calculate the RHS, and finally gives an approximation to $y_{n+1}$ again. Perhaps surprisingly, RK2 is better in many ways than forward Euler. Specifically, it is L-stable. But the most common one is RK4, which has $y_{n+1} = y_n + \frac{h}{6}\cdot\allowbreak\left(k_1 + 2k_2 + 2k_3 + k_4 \right)$, where $h$ is the step size, and \begin{enumerate} \item $k_1 = \ f(t_n,\allowbreak y_n)$ \item $k_2 = \ f\!\left(t_n + \frac{h}{2},\allowbreak y_n + h \frac{k_1}{2}\right)$ \item $k_3 = \ f\!\left(t_n + \frac{h}{2},\allowbreak y_n + h \frac{k_2}{2}\right)$ \item $k_4 = \ f\!\left(t_n + h,\allowbreak y_n + h k_3\right)$ \end{enumerate}. \end{document}