Appendix A

Formula Handbook

Rendered from uploaded appendix source with KaTeX. No raw LaTeX displayed.

Chapter 1: Floating-Point Arithmetic

Floating-point number

x=\pm m\beta^e.

Here x is the represented number, m is the significand or mantissa, \beta\ge2 is the base, and e\in\Z is the exponent. In a normalized system, m is scaled so that its leading digit is nonzero.

Unit roundoff and machine epsilon

u=\frac12\beta^{1-p}, \qquad \epsmach=\beta^{1-p}, \qquad u=\frac{\epsmach}{2}.

Here p is the number of base-\beta digits in the significand, u is the unit roundoff for rounding to nearest, and \epsmach is the distance from 1 to the next larger floating-point number.

Standard rounding model

\fl(x)=x(1+\delta), \qquad |\delta|\le u.

Here \fl(x) is the floating-point rounding of x, \delta is the relative rounding error, and u is the unit roundoff. The model assumes x is in the normalized range and no overflow or underflow occurs.

Floating-point arithmetic model

\fl(x\circ y)=(x\circ y)(1+\delta), \qquad |\delta|\le u, \qquad \circ\in\{+,-,\times,/\}.

Here x and y are floating-point inputs, \circ is one arithmetic operation, and \delta is the local relative rounding error.

Accumulated product of roundoff factors

\prod_{j=1}^{n}(1+\delta_j)=1+\theta_n, \qquad |\theta_n|\le\gamma_n, \qquad \gamma_n=\frac{nu}{1-nu}.

Here \delta_j is the rounding error in the j-th operation, n is the number of operations being grouped, and the bound is meaningful when nu<1.

Relative condition number of a scalar function

\kappa_f(x)= \left| \frac{x f'(x)}{f(x)} \right|.

Here f is a scalar function, x is the input, f'(x) is the derivative, and \kappa_f(x) measures amplification of small relative perturbations in x.

Matrix condition number

\kappa(A)=\norm{A}\norm{A^{-1}}.

Here A\in\R^{n\times n} is nonsingular, \norm{\cdot} is a chosen matrix norm, and \kappa(A) measures sensitivity of the linear system Ax=b.

Stable quadratic formula

For

ax^2+bx+c=0,

compute

q=-\frac12\left(b+\sgn(b)\sqrt{b^2-4ac}\right), \qquad x_1=\frac{q}{a}, \qquad x_2=\frac{c}{q}.

Here a,b,c are the polynomial coefficients, x_1,x_2 are the roots, and q is an auxiliary quantity chosen to avoid catastrophic cancellation.

Chapter 2: Nonlinear Equations

Bisection method

c_k=\frac{a_k+b_k}{2}.

Here [a_k,b_k] is the current bracketing interval, c_k is its midpoint, and the method requires

f(a_k)f(b_k)<0.

The error bound after k bisections is

|x_\ast-c_k|\le \frac{b_0-a_0}{2^{k+1}},

where x_\ast is a root in the original interval [a_0,b_0].

Fixed-point iteration

x_{k+1}=g(x_k).

Here g is the iteration function, x_k is the current approximation, and a fixed point satisfies x_\ast=g(x_\ast). If

|g'(x)|\le L<1

on an invariant interval, then

|x_k-x_\ast|\le L^k|x_0-x_\ast|.

Newton method for one equation

x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}.

Here f(x)=0 is the nonlinear equation, x_k is the current approximation, f'(x_k)\ne0, and x_{k+1} is the next approximation. For a simple root,

e_{k+1} = \frac{f''(x_\ast)}{2f'(x_\ast)}e_k^2+O(e_k^3), \qquad e_k=x_k-x_\ast.

Modified Newton method for a multiple root

x_{k+1}=x_k-m\frac{f(x_k)}{f'(x_k)}.

Here m is the known multiplicity of the root x_\ast. This correction restores quadratic convergence when the root multiplicity is known.

Secant method

x_{k+1} = x_k- f(x_k) \frac{x_k-x_{k-1}}{f(x_k)-f(x_{k-1})}.

Here x_{k-1} and x_k are two previous approximations. The derivative f'(x_k) is replaced by the divided difference

\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}.

Aitken acceleration

\widehat x_k= x_k- \frac{(x_{k+1}-x_k)^2} {x_{k+2}-2x_{k+1}+x_k}.

Here x_k,x_{k+1},x_{k+2} are three consecutive terms of a linearly convergent sequence, and \widehat x_k is the accelerated approximation.

Newton method for nonlinear systems

J_F(x_k)s_k=-F(x_k), \qquad x_{k+1}=x_k+s_k.

Here F:\R^n\to\R^n, J_F(x_k) is the Jacobian matrix with entries

(J_F)_{ij}(x_k)=\frac{\partial F_i}{\partial x_j}(x_k),

s_k is the Newton correction, and x_{k+1} is the new iterate.

Damped Newton method

x_{k+1}=x_k+\alpha_k s_k, \qquad J_F(x_k)s_k=-F(x_k), \qquad 0<\alpha_k\le1.

Here \alpha_k is a damping parameter, often chosen by line search to reduce a merit function such as

\Phi(x)=\frac12\norm{F(x)}_2^2.

Broyden update

B_{k+1} = B_k+ \frac{(y_k-B_ks_k)s_k^T}{s_k^Ts_k}.

Here B_k approximates J_F(x_k),

s_k=x_{k+1}-x_k, \qquad y_k=F(x_{k+1})-F(x_k).

The update enforces the secant equation B_{k+1}s_k=y_k.

Chapter 3: Interpolation and Approximation

Data notation for interpolation

x_j=\text{interpolation node}, \qquad f_j=f(x_j), \qquad j=0,\ldots,n.

Here x_j are distinct points and f_j are the sampled data values. Every interpolation formula below assumes x_i\ne x_j for i\ne j.

Lagrange cardinal polynomial

\ell_j(x)= \prod_{\substack{0\le m\le n\\m\ne j}} \frac{x-x_m}{x_j-x_m}.

Here j is the index of the cardinal function, x_m are all interpolation nodes, and

\ell_j(x_i)=\delta_{ij}.

Lagrange interpolating polynomial

p_n(x)=\sum_{j=0}^{n}f_j\ell_j(x).

Here p_n is the unique polynomial of degree at most n satisfying

p_n(x_j)=f_j=f(x_j), \qquad j=0,\ldots,n.

Interpolation error formula

f(x)-p_n(x) = \frac{f^{(n+1)}(\xi_x)}{(n+1)!} \prod_{j=0}^{n}(x-x_j).

Here \xi_x is a point between the smallest and largest of

x,x_0,\ldots,x_n,

and the formula assumes f\in C^{n+1} on the interpolation interval.

Newton divided difference

f[x_i,\ldots,x_{i+k}] = \frac{ f[x_{i+1},\ldots,x_{i+k}] - f[x_i,\ldots,x_{i+k-1}] } {x_{i+k}-x_i}.

Here f[x_i]=f_i=f(x_i), and f[x_i,\ldots,x_{i+k}] is the divided difference of order k.

Newton interpolation formula

p_n(x)= f[x_0]+ \sum_{k=1}^{n} f[x_0,\ldots,x_k]\prod_{j=0}^{k-1}(x-x_j).

Here p_n is the same interpolation polynomial as the Lagrange polynomial, written in nested Newton form.

Barycentric interpolation weights

w_j^{\rm bar}= \frac{1}{\displaystyle\prod_{\substack{0\le m\le n\\m\ne j}}(x_j-x_m)}.

Here w_j^{\rm bar} is the barycentric weight associated with node x_j. The superscript ``bar'' is used here to distinguish barycentric weights from quadrature weights.

Barycentric interpolation formula

p_n(x)= \frac{\displaystyle\sum_{j=0}^{n}\frac{w_j^{\rm bar}}{x-x_j}f_j} {\displaystyle\sum_{j=0}^{n}\frac{w_j^{\rm bar}}{x-x_j}}, \qquad x\ne x_j.

Here f_j=f(x_j), w_j^{\rm bar} are barycentric weights, and if x=x_j one sets p_n(x_j)=f_j directly.

Chebyshev--Lobatto nodes on \([-1,1]\)

x_j=\cos\left(\frac{j\pi}{n}\right), \qquad j=0,\ldots,n.

Here n is the polynomial degree and n+1 is the number of nodes. These nodes include both endpoints -1 and 1.

Affine mapping of nodes to \([a,b]\)

t_j=\frac{a+b}{2}+\frac{b-a}{2}x_j.

Here x_j\in[-1,1] is a reference node, t_j\in[a,b] is the physical node, and a<b are the interval endpoints.

Cubic spline second-derivative unknowns

M_j=S''(x_j), \qquad j=0,\ldots,n.

Here S is the cubic spline interpolant. For a natural spline,

M_0=0, \qquad M_n=0.

Bernstein basis and Bernstein approximation

B_{k,n}(x)= {n\choose k}x^k(1-x)^{n-k}, \qquad k=0,\ldots,n,

\mathcal B_n f(x)= \sum_{k=0}^{n} f\left(\frac{k}{n}\right)B_{k,n}(x).

Here B_{k,n} is the Bernstein basis polynomial and \mathcal B_n f is the Bernstein approximation on [0,1].

Chapter 4: Direct Linear Systems

Linear system and residual

Ax=b, \qquad r=b-A\widetilde x.

Here A\in\R^{n\times n}, b\in\R^n, x is the exact solution, \widetilde x is a computed solution, and r is the residual.

LU factorization without pivoting

A=LU.

Here L is lower triangular, usually with diagonal entries 1, and U is upper triangular. Solve

Ly=b, \qquad Ux=y.

LU factorization with partial pivoting

PA=LU.

Here P is a permutation matrix, L is lower triangular, and U is upper triangular. Solve

Ly=Pb, \qquad Ux=y.

Gaussian elimination multiplier

\ell_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}.

Here a_{ik}^{(k)} is the current entry in row i, column k, after the first k-1 elimination steps, and a_{kk}^{(k)} is the pivot.

Cholesky factorization

A=LL^T, \qquad A=A^T\succ0.

Here L is lower triangular with positive diagonal entries. The diagonal and off-diagonal entries are

\ell_{ii}= \left(a_{ii}-\sum_{k=1}^{i-1}\ell_{ik}^2\right)^{1/2},

\ell_{ji}= \frac{1}{\ell_{ii}} \left(a_{ji}-\sum_{k=1}^{i-1}\ell_{jk}\ell_{ik}\right), \qquad j>i.

Schur complement

For a block matrix

\begin{pmatrix} A&B\\ C&D \end{pmatrix},

the Schur complement of A is

S=D-CA^{-1}B.

Here A must be nonsingular, and B,C,D are conforming blocks.

Householder reflector

H=I-2vv^T, \qquad \norm{v}_2=1.

Here v is the Householder direction and H is orthogonal:

H^TH=I.

QR factorization

A=QR, \qquad Q^TQ=I.

Here Q has orthonormal columns and R is upper triangular. For square nonsingular A, solve

Rx=Q^Tb.

Singular value decomposition

A=U\Sigma V^T.

Here U and V are orthogonal, and

\Sigma=\diag(\sigma_1,\ldots,\sigma_r), \qquad \sigma_1\ge\cdots\ge\sigma_r>0

contains the singular values.

Sherman--Morrison formula

(A+uv^T)^{-1} = A^{-1}- \frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}.

Here A is nonsingular, u,v\in\R^n, and the formula requires

1+v^TA^{-1}u\ne0.

Woodbury formula

(A+UCV)^{-1} = A^{-1}-A^{-1}U(C^{-1}+VA^{-1}U)^{-1}VA^{-1}.

Here A and C are nonsingular, and U,V are rectangular matrices with compatible dimensions.

Iterative refinement

r^{(k)}=b-Ax^{(k)}, \qquad A d^{(k)}=r^{(k)}, \qquad x^{(k+1)}=x^{(k)}+d^{(k)}.

Here x^{(k)} is the current solution, r^{(k)} is its residual, and d^{(k)} is the correction.

Chapter 5: Iterative Linear Systems

Stationary iteration from a splitting

A=M-N, \qquad x^{(k+1)}=M^{-1}Nx^{(k)}+M^{-1}b.

Here M must be easy to invert, N=M-A, and the iteration matrix is

G=M^{-1}N.

The method converges for every initial guess if

\rho(G)<1.

Jacobi method

x_i^{(k+1)} = \frac1{a_{ii}} \left( b_i-\sum_{\substack{j=1\\j\ne i}}^{n}a_{ij}x_j^{(k)} \right), \qquad i=1,\ldots,n.

Here A=(a_{ij}), b=(b_i), and x_i^{(k)} is the i-th component at iteration k.

Gauss--Seidel method

x_i^{(k+1)} = \frac1{a_{ii}} \left( b_i-\sum_{j<i}a_{ij}x_j^{(k+1)} -\sum_{j>i}a_{ij}x_j^{(k)} \right).

Here newly updated components x_j^{(k+1)} are used immediately for j<i.

SOR method

x_i^{(k+1)} = (1-\omega)x_i^{(k)} + \frac{\omega}{a_{ii}} \left( b_i-\sum_{j<i}a_{ij}x_j^{(k+1)} -\sum_{j>i}a_{ij}x_j^{(k)} \right).

Here \omega is the relaxation parameter. \omega=1 gives Gauss--Seidel.

Richardson iteration

x^{(k+1)}=x^{(k)}+\alpha(b-Ax^{(k)}).

Here \alpha is a scalar step parameter. For A=A^T\succ0, a standard optimal choice is

\alpha_{\rm opt} = \frac{2}{\lambda_{\min}(A)+\lambda_{\max}(A)}.

Krylov subspace

\mathcal K_k(A,r_0) = \spann\{r_0,Ar_0,A^2r_0,\ldots,A^{k-1}r_0\}.

Here r_0=b-Ax_0 is the initial residual and k is the Krylov dimension.

Conjugate gradient method

\alpha_k=\frac{r_k^Tr_k}{p_k^TAp_k}, \qquad x_{k+1}=x_k+\alpha_kp_k,

r_{k+1}=r_k-\alpha_kAp_k, \qquad \beta_k=\frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k},

p_{k+1}=r_{k+1}+\beta_kp_k.

Here A=A^T\succ0, r_k=b-Ax_k, and p_k is the search direction.

CG error estimate

\norm{e_k}_A \le 2 \left( \frac{\sqrt{\kappa(A)}-1}{\sqrt{\kappa(A)}+1} \right)^k \norm{e_0}_A.

Here e_k=x_\ast-x_k,

\norm{e}_A=(e^TAe)^{1/2},

and \kappa(A)=\lambda_{\max}(A)/\lambda_{\min}(A) for A=A^T\succ0.

Arnoldi relation for GMRES

AV_k=V_{k+1}\overline H_k.

Here V_k=[v_1,\ldots,v_k] has orthonormal columns, \overline H_k\in\R^{(k+1)\times k} is upper Hessenberg, and the columns of V_k span \mathcal K_k(A,r_0).

GMRES least-squares problem

x_k=x_0+V_ky_k, \qquad y_k=\arg\min_y \norm{\beta e_1-\overline H_ky}_2.

Here \beta=\norm{r_0}_2, e_1=(1,0,\ldots,0)^T, and r_0=b-Ax_0.

Preconditioning

Left preconditioning:

M^{-1}Ax=M^{-1}b.

Right preconditioning:

AM^{-1}y=b, \qquad x=M^{-1}y.

Here M is a preconditioner chosen so that M^{-1}A or AM^{-1} is better conditioned than A, and applying M^{-1} is inexpensive.

Chapter 6: Eigenvalue Problems

Eigenvalue problem

Ax=\lambda x, \qquad x\ne0.

Here A\in\C^{n\times n}, \lambda\in\C is an eigenvalue, and x is a corresponding eigenvector.

Rayleigh quotient

\rho(x)=\frac{x^TAx}{x^Tx} \quad \text{or} \quad \rho(x)=\frac{x^\ast Ax}{x^\ast x} \text{ in the complex case.}

Here x\ne0. If Ax=\lambda x, then \rho(x)=\lambda.

Power method

y_k=Ax_k, \qquad x_{k+1}=\frac{y_k}{\norm{y_k}}.

Here x_k is the current normalized vector. The method targets the eigenvalue of largest modulus when it is simple and dominant.

Inverse iteration

(A-\sigma I)y_k=x_k, \qquad x_{k+1}=\frac{y_k}{\norm{y_k}}.

Here \sigma is a shift, I is the identity matrix, and the method targets the eigenvalue closest to \sigma.

Rayleigh quotient iteration

\mu_k=\rho(x_k), \qquad (A-\mu_k I)y_k=x_k, \qquad x_{k+1}=\frac{y_k}{\norm{y_k}}.

Here the shift is updated at each iteration by the Rayleigh quotient.

QR iteration

A_k=Q_kR_k, \qquad A_{k+1}=R_kQ_k.

Here Q_k is orthogonal/unitary, R_k is upper triangular, and A_{k+1} is orthogonally similar to A_k.

Shifted QR iteration

A_k-\mu_k I=Q_kR_k, \qquad A_{k+1}=R_kQ_k+\mu_kI.

Here \mu_k is the shift, often chosen to accelerate convergence.

Gershgorin disks

\lambda\in \bigcup_{i=1}^{n} \left\{ z\in\C: |z-a_{ii}|\le R_i \right\}, \qquad R_i=\sum_{\substack{j=1\\j\ne i}}^{n}|a_{ij}|.

Here A=(a_{ij}), and every eigenvalue of A lies in at least one disk.

Eigenpair residual

r=Ax-\mu x.

Here \mu is an approximate eigenvalue, x is an approximate eigenvector, and r measures the defect in the eigenvalue equation.

Hermitian residual bound

\dist(\mu,\sigma(A))\le \norm{r}_2.

Here A=A^\ast, \sigma(A) is the spectrum of A, r=Ax-\mu x, and \norm{x}_2=1.

Bauer--Fike theorem

\min_{\lambda\in\sigma(A)} |\widetilde\lambda-\lambda| \le \kappa(V)\norm{E}, \qquad A=V\Lambda V^{-1}.

Here \widetilde\lambda is an eigenvalue of A+E, E is the perturbation, and \kappa(V)=\norm{V}\norm{V^{-1}}.

Chapter 7: Least Squares

Linear least-squares problem

\min_{x\in\R^n}\norm{Ax-b}_2.

Here A\in\R^{m\times n}, usually m\ge n, b\in\R^m, and x is the unknown coefficient vector.

Residual and orthogonality

r=b-Ax, \qquad A^Tr=0.

Here r is the least-squares residual. Orthogonality means the residual is perpendicular to the column space of A.

Normal equations

A^TAx=A^Tb.

Here A^TA\in\R^{n\times n}. If A has full column rank, then A^TA is symmetric positive definite.

QR least-squares solution

A=QR, \qquad Rx=Q^Tb.

Here Q\in\R^{m\times n} has orthonormal columns, R\in\R^{n\times n} is upper triangular, and x is found by back substitution.

SVD least-squares solution

A=U\Sigma V^T, \qquad x=A^\dagger b=V\Sigma^\dagger U^Tb.

Here \Sigma^\dagger is obtained by replacing each nonzero singular value \sigma_i by 1/\sigma_i.

Projection matrix onto \((A)\)

P_A=A(A^TA)^{-1}A^T.

Here A has full column rank and P_A b is the orthogonal projection of b onto the column space of A.

Weighted least squares

\min_x (Ax-b)^TW(Ax-b).

Here W=W^T\succ0 is a weight matrix. The normal equations are

A^TWAx=A^TWb.

Tikhonov regularization

\min_x \left( \norm{Ax-b}_2^2+\lambda\norm{x}_2^2 \right).

Here \lambda>0 is the regularization parameter. The system is

(A^TA+\lambda I)x=A^Tb.

Truncated SVD solution

x_k= \sum_{i=1}^{k} \frac{u_i^Tb}{\sigma_i}v_i.

Here u_i and v_i are left and right singular vectors, \sigma_i are singular values, and k is the truncation index.

Nonlinear least squares and Gauss--Newton

\min_x \frac12\norm{r(x)}_2^2, \qquad J(x_k)^TJ(x_k)s_k=-J(x_k)^Tr(x_k).

Here r(x)\in\R^m is the nonlinear residual vector and

J_{ij}(x)=\frac{\partial r_i}{\partial x_j}(x)

is its Jacobian.

Chapter 8: Numerical Differentiation

Forward difference

f'(x)\approx\frac{f(x+h)-f(x)}{h}.

Here h>0 is the step size. The truncation error is O(h) for smooth f.

Backward difference

f'(x)\approx\frac{f(x)-f(x-h)}{h}.

Here h>0. This formula is first order.

Central difference

f'(x)\approx\frac{f(x+h)-f(x-h)}{2h}.

Here h>0. This formula is second order:

f'(x)=\frac{f(x+h)-f(x-h)}{2h}+O(h^2).

Second derivative central difference

f''(x)\approx \frac{f(x-h)-2f(x)+f(x+h)}{h^2}.

Here h>0. The truncation error is O(h^2).

Fourth-order first derivative

f'(x)\approx \frac{-f(x+2h)+8f(x+h)-8f(x-h)+f(x-2h)}{12h}.

Here h>0. The formula uses the five points

x-2h,\quad x-h,\quad x,\quad x+h,\quad x+2h.

Fourth-order second derivative

f''(x)\approx \frac{-f(x+2h)+16f(x+h)-30f(x)+16f(x-h)-f(x-2h)} {12h^2}.

Here h>0. This is a five-point fourth-order approximation.

Richardson extrapolation

D(h)=D+C h^p+O(h^{p+1}),

then

D_{\rm rich}= \frac{2^pD(h/2)-D(h)}{2^p-1}.

Here D is the exact derivative quantity, D(h) is an approximation with step h, and p is the leading order of the error.

Roundoff-truncation balance

Forward difference:

E(h)\approx C_1h+C_2\frac{u}{h}, \qquad h_{\rm opt}=O(\sqrt u).

Central difference:

E(h)\approx C_1h^2+C_2\frac{u}{h}, \qquad h_{\rm opt}=O(u^{1/3}).

Here u is the unit roundoff.

Complex-step derivative

f'(x)\approx \frac{\operatorname{Im}(f(x+ih))}{h}.

Here i^2=-1, h>0, and f must be extendable to complex arguments. This formula avoids subtractive cancellation.

Finite-difference weights by moment matching

For nodes x_j=x_0+c_jh, choose weights a_j such that

\sum_j a_j c_j^m = \begin{cases} 0, & m\ne q,\\ q!, & m=q, \end{cases}

for m=0,\ldots,r. Then

f^{(q)}(x_0)\approx \frac{1}{h^q}\sum_j a_j f(x_0+c_jh).

Here q is the derivative order and r controls the accuracy.

Chebyshev differentiation matrix

For

x_j=\cos\left(\frac{j\pi}{N}\right), \qquad j=0,\ldots,N,

define c_0=c_N=2 and c_j=1 otherwise. For i\ne j,

D_{ij} = \frac{c_i}{c_j} \frac{(-1)^{i+j}}{x_i-x_j}.

The diagonal is

D_{ii}=-\sum_{j\ne i}D_{ij}.

Chapter 9: Numerical Integration

Quadrature notation

I(f)=\int_a^b f(x)\,\dd x, \qquad Q(f)=\sum_{j=1}^{q}w_j^{\rm quad} f(x_j).

Here x_j are quadrature nodes, w_j^{\rm quad} are quadrature weights, and q is the number of quadrature points.

Midpoint rule

\int_a^b f(x)\,\dd x \approx (b-a)f\left(\frac{a+b}{2}\right).

Here a,b are interval endpoints and the midpoint is m=(a+b)/2.

Trapezoidal rule

\int_a^b f(x)\,\dd x \approx \frac{b-a}{2}[f(a)+f(b)].

Here the nodes are x_1=a, x_2=b, and both weights are

w_1^{\rm quad}=w_2^{\rm quad}=\frac{b-a}{2}.

Composite trapezoidal rule

Let

x_j=a+jh, \qquad h=\frac{b-a}{n}, \qquad f_j=f(x_j).

Then

\int_a^b f(x)\,\dd x \approx h\left[ \frac12f_0+\sum_{j=1}^{n-1}f_j+\frac12f_n \right].

Here n is the number of subintervals and n+1 is the number of nodes.

Simpson rule

Let

m=\frac{a+b}{2}.

Then

\int_a^b f(x)\,\dd x \approx \frac{b-a}{6}[f(a)+4f(m)+f(b)].

The nodes are a,m,b, and the weights are

\frac{b-a}{6},\qquad \frac{4(b-a)}{6},\qquad \frac{b-a}{6}.

Composite Simpson rule

Let n be even,

x_j=a+jh, \qquad h=\frac{b-a}{n}, \qquad f_j=f(x_j).

Then

\int_a^b f(x)\,\dd x \approx \frac{h}{3} \left[ f_0+4\sum_{\substack{1\le j\le n-1\\j\ {\rm odd}}}f_j +2\sum_{\substack{2\le j\le n-2\\j\ {\rm even}}}f_j +f_n \right].

Gauss--Legendre rule on \([a,b]\)

Given reference nodes \xi_j\in[-1,1] and reference weights \omega_j^{GL}, set

x_j=\frac{b-a}{2}\xi_j+\frac{a+b}{2}, \qquad w_j^{\rm quad}=\frac{b-a}{2}\omega_j^{GL}.

Then

\int_a^b f(x)\,\dd x \approx \sum_{j=1}^{q}w_j^{\rm quad} f(x_j).

Here q is the number of Gauss points.

Two-point Gauss--Legendre on \([a,b]\)

Reference nodes and weights:

\xi_1=-\frac1{\sqrt3}, \qquad \xi_2=\frac1{\sqrt3}, \qquad \omega_1^{GL}=\omega_2^{GL}=1.

Physical nodes:

x_j=\frac{b-a}{2}\xi_j+\frac{a+b}{2}.

Physical weights:

w_1^{\rm quad}=w_2^{\rm quad}=\frac{b-a}{2}.

Three-point Gauss--Legendre on \([a,b]\)

Reference nodes:

\xi_1=-\sqrt{\frac35}, \qquad \xi_2=0, \qquad \xi_3=\sqrt{\frac35}.

Reference weights:

\omega_1^{GL}=\omega_3^{GL}=\frac59, \qquad \omega_2^{GL}=\frac89.

Map to [a,b] by

x_j=\frac{b-a}{2}\xi_j+\frac{a+b}{2}, \qquad w_j^{\rm quad}=\frac{b-a}{2}\omega_j^{GL}.

Romberg extrapolation

Let R_{k,0}=T(h_k), where T(h_k) is the composite trapezoidal rule with step h_k. Then

R_{k,j}=R_{k,j-1} + \frac{R_{k,j-1}-R_{k-1,j-1}}{4^j-1}.

Here j is the extrapolation level.

Adaptive Simpson error estimate

Let S(a,b) be Simpson's rule on [a,b], and let m=(a+b)/2. Then

E\approx \frac{S(a,m)+S(m,b)-S(a,b)}{15}.

Here E estimates the local integration error.

Monte Carlo integration

For X_j\sim U(a,b),

\int_a^b f(x)\,\dd x \approx \frac{b-a}{N}\sum_{j=1}^{N}f(X_j).

Here N is the number of random samples. The statistical error is typically O(N^{-1/2}).

Chapter 10: Initial-Value Problems for ODEs

IVP and step notation

y'(t)=f(t,y(t)), \qquad y(t_0)=y_0.

Time grid:

t_n=t_0+nh, \qquad y_n\approx y(t_n), \qquad f_n=f(t_n,y_n).

Here h>0 is the step size.

Exact integral form

y(t_{n+1}) = y(t_n)+ \int_{t_n}^{t_{n+1}}f(s,y(s))\,\dd s.

Here t_{n+1}=t_n+h. Every one-step method approximates this integral.

Explicit Euler

y_{n+1}=y_n+h f_n.

Here f_n=f(t_n,y_n). This is first order globally.

Backward Euler

y_{n+1}=y_n+h f(t_{n+1},y_{n+1}).

Here y_{n+1} is unknown on both sides. Define the nonlinear residual

G(Y)=Y-y_n-hf(t_{n+1},Y),

and solve G(Y)=0.

Trapezoidal method

y_{n+1} = y_n+\frac h2 \left[ f(t_n,y_n)+f(t_{n+1},y_{n+1}) \right].

This is implicit and second order.

Explicit midpoint RK2

k_1=f(t_n,y_n),

k_2=f\left(t_n+\frac h2,\ y_n+\frac h2k_1\right),

y_{n+1}=y_n+h k_2.

Here k_1,k_2 are stage slopes.

Heun method

k_1=f(t_n,y_n), \qquad k_2=f(t_n+h,y_n+hk_1),

y_{n+1}=y_n+\frac h2(k_1+k_2).

This is the explicit trapezoidal method.

Classical RK4

\begin{aligned} k_1&=f(t_n,y_n),\\ k_2&=f\left(t_n+\frac h2,y_n+\frac h2k_1\right),\\ k_3&=f\left(t_n+\frac h2,y_n+\frac h2k_2\right),\\ k_4&=f(t_n+h,y_n+hk_3),\\ y_{n+1}&=y_n+\frac h6(k_1+2k_2+2k_3+k_4). \end{aligned}

Here k_1,\ldots,k_4 are stage slopes and h is the step size.

Stability test equation

y'=\lambda y, \qquad z=h\lambda.

A method gives

y_{n+1}=R(z)y_n.

Here R(z) is the stability function. The absolute stability region is

\mathcal S=\{z\in\C:\ |R(z)|\le1\}.

Common stability functions

Explicit Euler:

R(z)=1+z.

Backward Euler:

R(z)=\frac1{1-z}.

Trapezoidal:

R(z)=\frac{1+z/2}{1-z/2}.

RK4:

R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}.

AB2 method

y_{n+1}=y_n+\frac h2(3f_n-f_{n-1}).

Here f_n=f(t_n,y_n) and f_{n-1}=f(t_{n-1},y_{n-1}).

BDF2 method

\frac{3y_{n+1}-4y_n+y_{n-1}}{2h} = f(t_{n+1},y_{n+1}).

This is implicit. The unknown is y_{n+1}.

Embedded adaptive step update

h_{\rm new}=h\,s\,E^{-1/(p+1)}.

Here h is the current step size, s\in(0,1) is a safety factor, E is the normalized local error estimate, and p is the order of the lower method in the embedded pair.

Chapter 11: Boundary-Value Problems for ODEs

Second-order BVP

y''=f(x,y,y'), \qquad y(a)=\alpha, \qquad y(b)=\beta.

Here a,b are endpoints and \alpha,\beta are boundary data.

Shooting residual

F(s)=y(b;s)-\beta.

Here s=y'(a) is the guessed initial slope, y(b;s) is the terminal value obtained by solving the IVP, and the correct slope satisfies F(s)=0.

Newton update for shooting

s_{k+1}=s_k-\frac{F(s_k)}{F'(s_k)}.

Here F'(s_k) may be computed by finite differences or by the variational equation.

Centered finite difference for \(y''\)

With

x_i=a+ih, \qquad Y_i\approx y(x_i),

y''(x_i)\approx \frac{Y_{i-1}-2Y_i+Y_{i+1}}{h^2}.

Here h=(b-a)/N.

Centered finite difference for \(y'\)

y'(x_i)\approx \frac{Y_{i+1}-Y_{i-1}}{2h}.

Here Y_i\approx y(x_i) and h is the mesh spacing.

Weak form for a linear second-order problem

For

-(p y')'+q y=g, \qquad y(a)=y(b)=0,

the weak form is

\int_a^b p y'v'\,\dd x+ \int_a^b q yv\,\dd x = \int_a^b gv\,\dd x.

Here v is a test function satisfying v(a)=v(b)=0.

Galerkin approximation

y_h(x)=\sum_{j=1}^{N}Y_j\phi_j(x).

Here \phi_j are basis functions, Y_j are unknown coefficients, and N is the number of degrees of freedom.

Stiffness and load entries

K_{ij}= \int_a^b p\phi_j'\phi_i'\,\dd x+ \int_a^b q\phi_j\phi_i\,\dd x,

F_i=\int_a^b g\phi_i\,\dd x.

Here K=(K_{ij}) is the stiffness matrix and F=(F_i) is the load vector.

Element quadrature

For an element [x_e,x_{e+1}],

\int_{x_e}^{x_{e+1}}G(x)\,\dd x \approx \frac{h_e}{2} \sum_{\ell=1}^{q}w_\ell^{\rm ref} G\left(x_{e,c}+\frac{h_e}{2}\xi_\ell\right).

Here h_e=x_{e+1}-x_e, x_{e,c}=(x_e+x_{e+1})/2, \xi_\ell\in[-1,1] are reference Gauss nodes, and w_\ell^{\rm ref} are reference Gauss weights.

Nonlinear Galerkin Newton system

If R(Y)=0 is the nonlinear residual, Newton's method is

J(Y_k)s_k=-R(Y_k), \qquad Y_{k+1}=Y_k+s_k.

Here

J_{ij}(Y_k)=\frac{\partial R_i}{\partial Y_j}(Y_k)

is the Jacobian matrix.

Chapter 12: Partial Differential Equations

Grid notation in one space dimension

x_i=a+i\Delta x, \qquad t^n=n\Delta t, \qquad U_i^n\approx u(x_i,t^n).

Here \Delta x is the spatial mesh size and \Delta t is the time step.

Explicit heat equation scheme

For

u_t=\kappa u_{xx},

U_i^{n+1} = U_i^n+\mu(U_{i-1}^n-2U_i^n+U_{i+1}^n), \qquad \mu=\frac{\kappa\Delta t}{\Delta x^2}.

Here \kappa>0 is diffusivity. Stability in one dimension requires

\mu\le\frac12.

Crank--Nicolson heat scheme

\frac{U_i^{n+1}-U_i^n}{\Delta t} = \frac{\kappa}{2} \left[ \frac{U_{i-1}^{n+1}-2U_i^{n+1}+U_{i+1}^{n+1}}{\Delta x^2} + \frac{U_{i-1}^{n}-2U_i^{n}+U_{i+1}^{n}}{\Delta x^2} \right].

Here all U_i^{n+1} are solved simultaneously.

Centered wave equation scheme

For

u_{tt}=c^2u_{xx},

U_i^{n+1} = 2U_i^n-U_i^{n-1} + \lambda^2(U_{i-1}^n-2U_i^n+U_{i+1}^n),

\lambda=\frac{c\Delta t}{\Delta x}.

Here c is wave speed. Stability requires \lambda\le1 in one dimension.

Upwind advection for \(a>0\)

For

u_t+a u_x=0,

U_i^{n+1}=U_i^n-\nu(U_i^n-U_{i-1}^n), \qquad \nu=\frac{a\Delta t}{\Delta x}.

Here a>0 is the advection velocity. Stability requires 0\le\nu\le1.

Five-point Laplacian

For a square grid with spacing h,

-\Delta_hU_{i,j} = \frac{ 4U_{i,j}-U_{i-1,j}-U_{i+1,j}-U_{i,j-1}-U_{i,j+1} }{h^2}.

Here U_{i,j}\approx u(x_i,y_j).

Poisson weak form

For

-\Delta u=f \quad\text{in }\Omega, \qquad u=0\quad\text{on }\partial\Omega,

the weak form is

\int_\Omega\nabla u\cdot\nabla v\,\dd x = \int_\Omega f v\,\dd x.

Here v\in H_0^1(\Omega) is a test function.

Finite-volume update

\overline U_j^{\,n+1} = \overline U_j^{\,n} - \frac{\Delta t}{\Delta x} \left( F_{j+1/2}^{\,n}-F_{j-1/2}^{\,n} \right).

Here \overline U_j^n is the cell average in cell j, and F_{j\pm1/2}^n are numerical fluxes across cell interfaces.

Method of lines

\frac{\dd U(t)}{\dd t}=F(t,U(t)).

Here U(t) is the vector of spatial degrees of freedom obtained after spatial discretization, and an ODE solver is then applied in time.

Chapter 13: Spectral and Pseudospectral Methods

Fourier series notation

u(x)=\sum_{k=-K}^{K}\widehat u_k e^{ikx}.

Here k is the Fourier wavenumber, K is the maximum retained wavenumber, \widehat u_k is the Fourier coefficient, and e^{ikx} is the Fourier basis.

Fourier spectral differentiation

\widehat{u_x}_k=ik\widehat u_k, \qquad \widehat{u_{xx}}_k=-k^2\widehat u_k.

Here \widehat{u_x}_k and \widehat{u_{xx}}_k are Fourier coefficients of the first and second derivatives.

Chebyshev polynomial

T_k(x)=\cos(k\arccos x), \qquad -1\le x\le1.

Here T_k is the Chebyshev polynomial of degree k.

Chebyshev--Lobatto nodes

x_j=\cos\left(\frac{j\pi}{N}\right), \qquad j=0,\ldots,N.

Here N is the polynomial degree and N+1 is the number of nodes.

Chebyshev differentiation matrix

Let

c_0=c_N=2, \qquad c_j=1,\quad 1\le j\le N-1.

For i\ne j,

D_{ij}= \frac{c_i}{c_j} \frac{(-1)^{i+j}}{x_i-x_j}.

For the diagonal,

D_{ii}=-\sum_{j\ne i}D_{ij}.

Here D maps nodal values U_j=u(x_j) to approximate derivative values:

u'(x_i)\approx \sum_{j=0}^{N}D_{ij}U_j.

Affine derivative scaling

x=\frac{2t-(a+b)}{b-a},

then

\frac{\dd}{\dd t}= \frac{2}{b-a}\frac{\dd}{\dd x}.

Therefore the physical differentiation matrix is

D_t=\frac{2}{b-a}D_x.

Pseudospectral collocation residual

\mathcal L u_N(x_j)=f(x_j), \qquad j=0,\ldots,N.

Here u_N is the spectral approximation, \mathcal L is the differential operator, and x_j are collocation nodes.

Galerkin spectral condition

\langle \mathcal L u_N-f,\phi_j\rangle=0, \qquad j=0,\ldots,N.

Here \phi_j are test functions and the residual is orthogonal to the trial/test space.

Aliasing identity

On an N-point periodic grid,

e^{i(k+mN)x_j}=e^{ikx_j}.

Here m\in\Z, k is a wavenumber, and x_j are the grid points. This causes unresolved high modes to appear as low modes.

Two-thirds dealiasing rule

\widehat u_k=0 \qquad \text{for} \qquad |k|>\frac{N}{3}.

Here N is the number of Fourier grid points and k is the Fourier mode.

Spectral filter

\widehat u_k^{\,F} = \sigma\left(\frac{|k|}{K}\right)\widehat u_k, \qquad \sigma(\eta)=\exp(-\alpha\eta^p).

Here K is the maximum retained wavenumber, p is the filter order, \alpha>0 controls strength, and 0\le\eta\le1.

Chapter 14: Optimization and Nonlinear Systems

Unconstrained optimization problem

\min_{x\in\R^n}f(x).

Here f:\R^n\to\R is the objective function and x is the vector of decision variables.

Gradient descent

x_{k+1}=x_k-\alpha_k\nabla f(x_k).

Here x_k is the current iterate, \alpha_k>0 is the step length, and \nabla f(x_k) is the gradient.

Newton optimization step

\nabla^2 f(x_k)s_k=-\nabla f(x_k), \qquad x_{k+1}=x_k+s_k.

Here \nabla^2 f(x_k) is the Hessian matrix and s_k is the Newton step.

Armijo condition

f(x_k+\alpha p_k) \le f(x_k)+c_1\alpha\nabla f(x_k)^Tp_k.

Here p_k is a descent direction, \alpha>0 is a trial step size, and 0<c_1<1.

BFGS update

B_{k+1} = B_k- \frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k} + \frac{y_ky_k^T}{y_k^Ts_k}.

Here

s_k=x_{k+1}-x_k, \qquad y_k=\nabla f(x_{k+1})-\nabla f(x_k),

and B_k approximates the Hessian.

Inverse BFGS update

H_{k+1} = (I-\rho_ks_ky_k^T)H_k(I-\rho_ky_ks_k^T) + \rho_ks_ks_k^T,

where

\rho_k=\frac{1}{y_k^Ts_k}.

Here H_k approximates the inverse Hessian.

Trust-region subproblem

\min_{\norm{s}\le\Delta_k} m_k(s)=f_k+g_k^Ts+\frac12s^TB_ks.

Here

f_k=f(x_k), \qquad g_k=\nabla f(x_k),

B_k is a Hessian or Hessian approximation, and \Delta_k is the trust-region radius.

Trust-region ratio

\rho_k= \frac{f(x_k)-f(x_k+s_k)} {m_k(0)-m_k(s_k)}.

Here the numerator is actual reduction and the denominator is predicted reduction.

Gauss--Newton step

For

\min_x \frac12\norm{r(x)}_2^2,

the Gauss--Newton step solves

J(x_k)^TJ(x_k)s_k=-J(x_k)^Tr(x_k).

Here r(x) is the residual vector and

J_{ij}(x_k)=\frac{\partial r_i}{\partial x_j}(x_k).

Levenberg--Marquardt step

(J^TJ+\lambda I)s=-J^Tr.

Here J=J(x_k), r=r(x_k), \lambda>0 is the damping parameter, and s is the step.

Nelder--Mead simplex operations

For n-dimensional optimization, use n+1 vertices

x_1,\ldots,x_{n+1}, \qquad f(x_1)\le\cdots\le f(x_{n+1}).

The centroid excluding the worst point is

\bar x=\frac1n\sum_{i=1}^{n}x_i.

Reflection:

x_r=\bar x+\alpha(\bar x-x_{n+1}).

Expansion:

x_e=\bar x+\gamma(x_r-\bar x).

Shrink:

x_i\leftarrow x_1+\sigma(x_i-x_1), \qquad i=2,\ldots,n+1.

Here \alpha,\gamma,\sigma are algorithmic parameters.

KKT conditions

For

\min f(x), \qquad c(x)=0, \qquad g(x)\le0,

define

\mathcal L(x,\lambda,\mu) = f(x)+\lambda^Tc(x)+\mu^Tg(x).

The KKT conditions are

\nabla_x\mathcal L(x,\lambda,\mu)=0,

c(x)=0, \qquad g(x)\le0, \qquad \mu\ge0,

\mu_i g_i(x)=0.

Here \lambda are equality multipliers and \mu are inequality multipliers.

Quadrature-based integral objective

J(\theta)=\int_a^b\Phi(x,\theta)\,\dd x \approx J_Q(\theta)=\sum_{j=1}^{q}w_j^{\rm quad}\Phi(x_j,\theta).

Here \theta is the parameter vector, x_j are quadrature nodes, w_j^{\rm quad} are quadrature weights, and q is the number of nodes. The gradient approximation is

\nabla_\theta J_Q(\theta)= \sum_{j=1}^{q}w_j^{\rm quad}\nabla_\theta\Phi(x_j,\theta).

Chapter 15: Fractional Differential Equations

Riemann--Liouville fractional integral

(I_{a+}^{\alpha}f)(t) = \frac{1}{\Gamma(\alpha)} \int_a^t(t-s)^{\alpha-1}f(s)\,\dd s.

Here \alpha>0 is the fractional order, a is the lower terminal, t>a, and \Gamma is the gamma function.

Caputo derivative

For 0<\alpha<1,

({}^CD_{a+}^{\alpha}f)(t) = \frac{1}{\Gamma(1-\alpha)} \int_a^t(t-s)^{-\alpha}f'(s)\,\dd s.

Here \alpha is the fractional order, f'(s) is the classical derivative, and (t-s)^{-\alpha} is the weakly singular memory kernel.

Caputo derivative of a power

For \beta>0,

{}^CD_{0+}^{\alpha}t^\beta = \frac{\Gamma(\beta+1)} {\Gamma(\beta+1-\alpha)} t^{\beta-\alpha}.

Here 0<\alpha<1, \beta is the power exponent, and t>0. Also,

{}^CD_{0+}^{\alpha}1=0.

Katugampola fractional integral

(I_{a+}^{\alpha,\rho}f)(t) = \frac{\rho^{1-\alpha}}{\Gamma(\alpha)} \int_a^t s^{\rho-1}(t^\rho-s^\rho)^{\alpha-1}f(s)\,\dd s.

Here \alpha>0 is the fractional order and \rho>0 is the Katugampola deformation parameter.

Caputo--Katugampola derivative

For 0<\alpha<1,

({}^CD_{a+}^{\alpha,\rho}f)(t) = \frac{\rho^\alpha}{\Gamma(1-\alpha)} \int_a^t (t^\rho-s^\rho)^{-\alpha}f'(s)\,\dd s.

Here \rho>0, t>a, and the memory kernel is

(t^\rho-s^\rho)^{-\alpha}.

Caputo--Katugampola derivative of a power

For \beta>0,

{}^CD_{0+}^{\alpha,\rho}t^\beta = \rho^\alpha \frac{\Gamma(\beta/\rho+1)} {\Gamma(\beta/\rho+1-\alpha)} t^{\beta-\alpha\rho}.

Here 0<\alpha<1, \rho>0, and \beta is the power exponent. Also,

{}^CD_{0+}^{\alpha,\rho}1=0.

Fractional clocks

Caputo clock:

\tau_\alpha(t)=\frac{t^\alpha}{\Gamma(1+\alpha)}, \qquad {}^CD^\alpha\tau_\alpha(t)=1.

Caputo--Katugampola clock:

\tau_{\alpha,\rho}(t) = \frac{t^{\alpha\rho}}{\rho^\alpha\Gamma(1+\alpha)}, \qquad {}^CD^{\alpha,\rho}\tau_{\alpha,\rho}(t)=1.

Here \tau is a stretched time variable adapted to fractional memory.

Caputo Volterra form

For

{}^CD_{0+}^{\alpha}u(t)=f(t,u(t)), \qquad u(0)=u_0,

the equivalent Volterra equation is

u(t)=u_0+ \frac{1}{\Gamma(\alpha)} \int_0^t(t-s)^{\alpha-1}f(s,u(s))\,\dd s.

Here 0<\alpha<1.

Caputo--Katugampola Volterra form

For

{}^CD_{0+}^{\alpha,\rho}u(t)=f(t,u(t)), \qquad u(0)=u_0,

the Volterra form is

u(t)=u_0+ \frac{\rho^{1-\alpha}}{\Gamma(\alpha)} \int_0^t s^{\rho-1}(t^\rho-s^\rho)^{\alpha-1}f(s,u(s))\,\dd s.

Here 0<\alpha<1 and \rho>0.

L1 weights and L1 formula

Time grid:

t_n=nh, \qquad u_n\approx u(t_n).

Weights:

a_k=(k+1)^{1-\alpha}-k^{1-\alpha}, \qquad k=0,\ldots,n-1.

L1 formula:

{}^CD_{0+}^{\alpha}u(t_n) \approx \frac{1}{h^\alpha\Gamma(2-\alpha)} \sum_{k=0}^{n-1} a_k(u_{n-k}-u_{n-k-1}).

Here h is the uniform time step.

Caputo--Katugampola L1 in stretched time

Use

y=t^\rho, \qquad y_n=nh_y, \qquad t_n=y_n^{1/\rho}.

Then

{}^CD_{0+}^{\alpha,\rho}u(t_n) \approx \frac{\rho^\alpha}{h_y^\alpha\Gamma(2-\alpha)} \sum_{k=0}^{n-1} a_k\bigl(u(t_{n-k})-u(t_{n-k-1})\bigr),

where

a_k=(k+1)^{1-\alpha}-k^{1-\alpha}.

Grunwald--Letnikov weights

g_k^{(\alpha)}=(-1)^k{\alpha\choose k}, \qquad g_0^{(\alpha)}=1.

Recurrence:

g_k^{(\alpha)} = \left(1-\frac{\alpha+1}{k}\right)g_{k-1}^{(\alpha)}.

Approximation:

D_{0+}^{\alpha}u(t_n) \approx h^{-\alpha} \sum_{k=0}^{n}g_k^{(\alpha)}u_{n-k}.

Lubich convolution quadrature weights

Let \delta(\zeta) be a multistep generating polynomial. Define weights \omega_n by

\left(\frac{\delta(\zeta)}{h}\right)^\alpha = \sum_{n=0}^{\infty}\omega_n\zeta^n.

Here h is the time step, \zeta is a generating variable, and \omega_n are convolution weights. Examples:

\delta(\zeta)=1-\zeta \quad\text{for backward Euler},

\delta(\zeta)=\frac32-2\zeta+\frac12\zeta^2 \quad\text{for BDF2}.

Chebyshev--Lobatto panel nodes

On the reference panel [0,1],

\theta_m^{(q)} = \frac{1-\cos(m\pi/q)}{2}, \qquad m=0,\ldots,q.

Physical panel [a_j,b_j]:

s_{j,m}=a_j+(b_j-a_j)\theta_m^{(q)}.

Here q is the polynomial degree and s_{j,m} are physical memory-panel nodes.

Cubic Chebyshev--Lobatto nodes

For q=3,

\theta_0=0, \qquad \theta_1=\frac14, \qquad \theta_2=\frac34, \qquad \theta_3=1.

On [a_j,b_j],

s_{j,m}=a_j+(b_j-a_j)\theta_m.

Quintic Chebyshev--Lobatto nodes

For q=5,

\theta_0=0, \quad \theta_1=\frac{3-\sqrt5}{8}, \quad \theta_2=\frac{5-\sqrt5}{8},

\theta_3=\frac{5+\sqrt5}{8}, \quad \theta_4=\frac{3+\sqrt5}{8}, \quad \theta_5=1.

Physical nodes are s_{j,m}=a_j+(b_j-a_j)\theta_m.

Caputo incomplete-beta moment

\int_a^b(t-s)^{-\alpha}s^\ell\,\dd s = t^{\ell+1-\alpha} \left[ B_{b/t}(\ell+1,1-\alpha) - B_{a/t}(\ell+1,1-\alpha) \right].

Here B_z(p,q) is the incomplete beta function,

B_z(p,q)=\int_0^z r^{p-1}(1-r)^{q-1}\,\dd r.

Caputo--Katugampola incomplete-beta moment

\int_a^b (t^\rho-s^\rho)^{-\alpha}s^\ell\,\dd s = \frac{t^{\ell+1-\alpha\rho}}{\rho} \left[ B_{(b/t)^\rho} \left(\frac{\ell+1}{\rho},1-\alpha\right) - B_{(a/t)^\rho} \left(\frac{\ell+1}{\rho},1-\alpha\right) \right].

Here \ell is a nonnegative integer, 0<a<b<t, 0<\alpha<1, and \rho>0.

Chebyshev--Beta panel weight

(\ell_m^{(q)})'(\theta) = \sum_{r=0}^{q-1}d_{m,r}^{(q)}\theta^r,

then the Caputo panel weight is

W_{j,m}^{C,q}(t) = \sum_{r=0}^{q-1} \frac{d_{m,r}^{(q)}}{h_j^{r+1}} I_r^C(a_j,b_j;t),

where

h_j=b_j-a_j, \qquad I_r^C(a_j,b_j;t)= \int_{a_j}^{b_j}(t-s)^{-\alpha}(s-a_j)^r\,\dd s.

Here \ell_m^{(q)} is the Lagrange basis polynomial on the reference nodes.

Fractional residual minimization

For a fractional model residual

\mathcal R[u_\theta](x,t)=0,

define

J(\theta)= \int_0^T\int_\Omega |\mathcal R[u_\theta](x,t)|^2\,\dd x\,\dd t.

Here u_\theta is a parameterized trial function, \theta is the parameter vector, \Omega is the spatial domain, and T is the final time. A quadrature version is

J_Q(\theta)= \sum_{i=1}^{N_x}\sum_{j=1}^{N_t} w_i^x w_j^t |\mathcal R[u_\theta](x_i,t_j)|^2,

where x_i,t_j are quadrature nodes and w_i^x,w_j^t are quadrature weights.