Multiply the complex Fourier series by $e^{-i\omega rt}$ and integrate over one period. Fourier orthogonality removes every term except the one with index $r$, givingFor general square-integrable functions, this is an orthogonal projection onto the corresponding Fourier mode. It remains valid when the expansion converges only in L2 space: integrate finite Fourier partial sums first, then pass to the limit using the Cauchy-Schwarz inequality. No unjustified pointwise exchange of an infinite series and an integral is needed.

Set $\omega =2\pi /T$ and $q=m+n$. If $q=0$, the integrand is one and its integral is $T$. If $q$ is a nonzero integer,ThusThis is Fourier orthogonality in its unconjugated form. With a complex conjugate on the second exponential, the nonzero case would instead be $m=n$.

Choose a fixed time increment $h=\Delta t>0$, let $t_n=nh$, and approximate $\mathbf{c}(t_n)$ by ${\mathbf{c}}_n$. The Backward Euler method evaluates the transport term at the new time:Rearrange to the linear systemAlthough $A$ is a singular matrix, $B$ is an invertible matrix for every $h>0$. In each row, its diagonal entry is $1+2h$ and the sum of the absolute values of its off-diagonal entries is $2h$; hence it is a strictly diagonally dominant matrix. Its Gershgorin discs have centre $1+2h$ and radius $2h$, excluding zero. In particular, the zero eigenvalue of $A$ becomes an eigenvalue $1$ of $B$.

The algorithm uses Reusing LU factors for multiple right-hand sides:

Both triangular solves use the factors of $I-hA$. The LU factorization of the singular transport matrix $A$ from part (a) is not a solver for this new linear system. No matrix inverse needs to be formed. Since $A$ and $h$ are fixed, the same LU factorization serves every step; only the right-hand side changes. A varying $\mathbf{b}(t)$ would give $r_n={\mathbf{c}}_n-h\mathbf{b}(t_{n+1})$ without changing these factors.

As a check, multiplying the update by $1^T$ givesexactly reproducing the inventory balance at the sampled times, apart from roundoff error. The Backward Euler method has temporal global error $O(h)$ for smooth solutions, so its A-stability does not remove the need to select $h$ for accuracy.

A steady state ${\mathbf{c}}_{\ast }$ satisfies $0=A{\mathbf{c}}_{\ast }-\mathbf{b}$, henceThe sign is positive on the right after rearranging the supplied equation. A solution exists exactly when $\mathbf{b}$ lies in the column space of $A$, which above was found to be the zero-sum vector subspace. Equivalently, multiply by the left null space vector $1^T$:This condition is necessary and sufficient for an algebraic steady state, since both the column space and that zero-sum vector subspace have dimension $3$. It expresses global conservation of mass: in the closed box model, total tracer production must balance total tracer removal. With the paper's convention $-\mathbf{b}$, a positive $b_i$ is a sink and a negative $b_i$ is a source. If the sum is nonzero, $dM/dt=-{\sum }_i{b}_i$ and the inventory changes linearly in time, preventing a steady state.

When the condition holds, all algebraic steady states form the affine solution space of a linear equation ${\mathbf{c}}_p+\alpha 1$, for any particular solution $A{\mathbf{c}}_p=\mathbf{b}$. The initial inventory selects one member through $1^T{\mathbf{c}}_{\ast }=1^T\mathbf{c}(0)$. For a physically admissible solution, its four concentrations must additionally be nonnegative; the compatibility condition alone imposes no such inequalities.

The null space consists of concentration patterns on which transport produces no change. Here it is the familyEvery box has the same concentration. Ocean circulation may continue, but inflow and outflow of the chemical tracer balance in each box. With $\mathbf{b}=0$ these are exactly the steady states; physically admissible members have $\alpha \ge 0$. With nonzero $\mathbf{b}$, adding $\alpha 1$ to any particular steady state leaves its transport balance unchanged, but a null space vector alone need not balance the sources and sinks.

To identify the tracer inventory, apply conserved inventory from a left null vector. Directly summing each column givesConsequently, for the supplied ordinary differential equations,For $\mathbf{b}=0$, $M$ is a conserved quantity. The connected box model equalizes concentrations, approaching the unique member of the null space with the initial inventory:The arithmetic mean is appropriate here because the supplied matrix has zero column sums. In a general box model the conserved quantity is $w^T\mathbf{c}$, where $w$ is a positive left null space vector; for physical concentrations expressed per unit volume, these weights represent the box volumes up to a common factor. One must establish the weights from the equations, rather than infer equal volumes from a schematic drawing.

Nonnegative initial concentrations remain nonnegative under transport alone. Indeed $A+2I$ is a nonnegative matrix, sohas nonnegative entries for $t\ge 0$: it is a positive matrix semigroup. Its column sums are one by $1^{T}A=0$. To see why equalization occurs, set $d=\mathbf{c}-(M/4)1$ for $\mathbf{b}=0$; then $1^{T}d=0$ andThe right side vanishes only for constant $d$, which in the zero-sum vector subspace means $d=0$. Its strict negative bound on the unit sphere in that vector subspace gives decay to zero. Thus the matrix exponential redistributes inventory without creating negative concentrations; arbitrary prescribed sinks can still exhaust the tracer if they are allowed to act indefinitely.

Regard $A$ as a linear map ${\mathbb{R}}^4\to {\mathbb{R}}^4$. By the rank-nullity theorem and matrix rank $3$,All three are vector subspaces of ${\mathbb{R}}^4$, with the row space identified with four-component row vectors in the usual way. The null space lies in the input vector space of concentrations; the column space lies in the output vector space of concentration-change rates. Both have four coordinates, although they play different physical roles.

In this particular matrix, the row sums and column sums both vanish. Put $1=(1,1,1,1{)}^T$. Then $A1=0$ and $1^{T}A=0$, and the two three-dimensional vector subspaces coincide:This equality follows here from both sum identities and the dimensions of a vector space; it does not make $A$ a symmetric matrix and is not true for every square matrix.

Since $L$ is an invertible matrix, $Av=0$ is equivalent to $Uv=0$. The last nontrivial row gives $v_3=v_4$; the second gives $v_2=v_4$; and the first gives $2v_1=v_2+v_4$, so $v_1=v_4$. A basis of the null space is thereforeIts single basis vector has Euclidean norm $2$, so the unit null space basis vector is $(1,1,1,1{)}^T/2$.

The pivot elements of $U$ lie in columns $1,2,3$. Take these columns from the original matrix $A$ to obtain a basis of its column space:The fourth column is minus the sum of these three columns, and the leading three-by-three matrix minor has nonzero determinant $-4$, confirming linear independence.

Elementary row operations preserve the row space. Thus the nonzero rows of $U$ give a basis of the row space of $A$:These basis vectors need not have unit Euclidean norm or be orthogonal. Their Euclidean norms are $\sqrt{6},\sqrt{2},2\sqrt{2}$, respectively; the three column basis vectors have Euclidean norms $2\sqrt{2},\sqrt{6},2\sqrt{2}$.

Apply Gaussian elimination to the supplied nonsymmetric matrix,The first pivot element is $-2$, with elimination multiplier $l_{21}=2/(-2)=-1$. Perform $R_2\leftarrow R_2+R_1$, givingThe second pivot element is $-1$, with elimination multiplier $l_{32}=-1$. Perform $R_3\leftarrow R_3+R_2$, then use the third pivot element $-2$ and elimination multiplier $l_{43}=-1$ to perform $R_4\leftarrow R_4+R_3$. The LU factorization isHere $L$ is a unit lower triangular matrix, hence an invertible matrix, and $U$ is an upper triangular matrix. No row interchange is needed: all three pivot elements used in elimination are nonzero. The final diagonal entry is zero, with an entirely zero last row; there is no further division by it. The three nonzero rows of $U$ are linearly independent, so the matrix rank is $3$ and $A$ is a singular matrix. Direct matrix multiplication gives $LU=A$.

Each supplied eigenvector gives relative displacements at $x=1,2,3$ in a standing wave. Insert zero at $x=0,4$ to describe the whole sampled string. A real normal mode is $u_j(t)=C{v}_j\cos (\omega t+\varphi )$: its shape stays fixed while its amplitude oscillates.

These semidiscrete fixed-string normal modes are samples of $\sin (k\pi x/4)$, with $k=3,2,1$ in the order listed. More precisely, their interior vectors are $v_j^{(k)}=2^{-1/2}\sin (k\pi j/4)$, which have unit Euclidean norm and form an orthonormal basis of ${\mathbb{R}}^3$. Their eigenvalues and angular frequencies obeyThus larger absolute value of a negative eigenvalue means faster oscillation, rather than faster decay. The numerical angular frequencies approximate the continuous fixed-string values $k\pi c/L$; the finite difference spatial approximation changes them. Only three independent interior normal modes exist on this grid, and any interior displacement or initial velocity vector can be expanded in this orthonormal basis.

For each row, a Gershgorin disc has centre $a_{jj}=-2$ and radius $r_j={\sum }_{k\ne j}|a_{jk}|$. Here $r_1=r_3=1$ and $r_2=2$, so the Gershgorin circle theorem givesSince $A$ is a real symmetric matrix, its eigenvalues are real. Intersecting the union of these Gershgorin discs with the real axis yieldsThis bound permits zero but does not assert a zero eigenvalue. The fixed-end Dirichlet boundary conditions exclude it: for any nonzero real vector $v=(v_1,v_2,v_3{)}^T$, extend by $v_0=v_4=0$ and calculateIndeed, equality would force every adjacent value to agree with the zero endpoint. Thus $-A$ is a symmetric positive-definite matrix, $A$ is a negative-definite matrix, and the eigenvalues satisfy $-4\le \lambda <0$. The same argument works for any number of interior grid points, with prefactor $c^2/h^2>0$ in the sum. The fixed string has no zero-frequency rigid translation.

The interior grid is $x=1,2,3$, with fixed-end displacements $u_0=u_4=0$. HenceUse the monic characteristic polynomial convention $p(\lambda )=\det (\lambda I-A)$, where $I$ is the identity matrix. Expansion of the determinant along the first row givesThe eigenvalues are the roots of a polynomial $p(\lambda )=0$. With the alternative determinant convention $\det (A-\lambda I)$, the characteristic polynomial is $-p(\lambda )$ because the matrix has odd size; its roots of a polynomial are unchanged.

For a nonzero, time-independent vector $\hat{\mathbf{u}}$, substitute $\mathbf{u}(t)=\hat{\mathbf{u}}{e}^{i\omega t}$. Taking two time derivatives givesCancel the nonzero exponential to obtain the eigenvalue problemThus an eigenvector specifies a normal mode shape and its negative eigenvalue specifies the square of its angular frequency, $\omega =\sqrt{-\lambda }$. Real displacements are obtained by taking the real part, or by combining the two temporal solutions as $\hat{\mathbf{u}}(a\cos \omega t+b\sin \omega t)$.

Put $h=\Delta x$ and suppose $L=(n+1)h$, so there are $n$ interior grid points $x_j=jh$, $j=1,\dots ,n$. Write $u_j(t)\simeq u(x_j,t)$. The fixed ends impose homogeneous Dirichlet boundary conditions,For a sufficiently smooth displacement field, the two Taylor series at fixed time givewhere all displayed derivatives on the right are evaluated at $(x_j,t)$. Adding cancels the odd terms, leaving the second-order central differenceSubstitute this finite difference into the wave equation and retain continuous time. The resulting method of lines gives coupled second-order ordinary differential equations,In particular, the first and last equations are ${\ddot{u}}_1=(c^2/h^2)(-2u_1+u_2)$ and ${\ddot{u}}_n=(c^2/h^2)(u_{n-1}-2u_n)$ for $n\ge 2$. With the interior displacement vector only,Equivalently, $A_{jj}=-2c^2/h^2$, $A_{j,j+1}=A_{j+1,j}=c^2/h^2$, and all other entries vanish. This entrywise definition also covers small $n$, including $A=(-2c^2/h^2)$ when $n=1$. The matrix is a symmetric matrix and a tridiagonal matrix. The spatial approximation has error $O(h^2)$ in Big O notation; time has not yet been discretized.

Excavate across the mapped fault scarp, approximately perpendicular to its strike, at a site where young sediments record and preserve deformation. On cleaned trench walls log bedding, soil horizons, shear zones, fissure fills, buried surfaces and colluvial wedges. Correlate the same marker beds across the fault before reconstructing individual displacements.

Look for event horizons, not merely a stack of sediments:

In the drawing, A was faulted in event 1; W1 accumulated; B buried and stabilized that surface. Event 2 displaced A, W1 and B, after which W2 accumulated and C sealed the fault tip. Thus A records two events, B one and C none. This combination of offset and sealing relationships is stronger evidence than identifying two wedge-shaped units alone.

Bracket event ages with multiple samples just below and above the event horizons. Suitable short-lived organic material can be dated by radiocarbon dating; optically stimulated luminescence dating can constrain deposition of appropriate mineral grains. Reworked old charcoal, younger roots, incomplete bleaching, erosional gaps and time elapsed between rupture and wedge deposition must enter the age interpretation. For calendar dates, $t_{\mathrm{youngest}\mathrm{cut}}<t_{\mathrm{event}}<t_{\mathrm{oldest}\mathrm{seal}}$ only when those deposits genuinely bracket the rupture; a radiocarbon age in years before present reverses the numerical ordering.^[ref]

One earthquake can generate several debris pulses and wedges; repeated earthquakes may reuse one strand without leaving separable deposits. Erosion, bioturbation, nondeposition, lateral facies changes and upward-dying splays can hide events or mimic terminations. Confirm relationships on both trench walls and, when possible, in adjacent exposures. The result is a minimum supported event sequence with age intervals, not an assumption of perfectly periodic recurrence or one earthquake per wedge.

A surface-reaching normal fault earthquake commonly produces a down-facing step towards the hanging wall, accompanied by extension across the break. A fresh free face can expose bedrock, polish or striations, roots or offset soil. The actual rupture is usually irregular: fault slip changes along strike, separate strands overlap or step, and tension fissures, small grabens and secondary cracks distribute part of the extension.

Exposed competent rock or very thin cover: a relatively narrow, sharp break can closely follow the principal fault plane. A steep free face can be recognizable immediately, then weather and shed debris. The face need not be vertical or maintain the same height along strike.

Loose granular or poorly consolidated cover: the main shear can splay upwards through alluvium or regolith. Open fissures, step-like scarps, rotated blocks, synthetic and antithetic normal faults, and a small intervening graben can form across a wider zone. Collapse and raveling quickly build colluvium at the toe. A single rupture event can therefore create several parallel breaks; counting them does not count earthquakes.

Thick weak or cohesive cover: deeper fault slip can be transferred into distributed shear strain, tilted beds, broad warping or monoclines, with tension cracks above a propagating tip. A discrete break may be attenuated or absent if the fault remains blind. Cover thickness alone is not enough to predict the response: cohesion, grain size, layering, saturation, rupture displacement and fault depth also matter. Saturated loose sediment can additionally develop soil liquefaction, sand blows or lateral spreading; these are secondary shaking effects and do not independently locate the tectonic rupture.

The observed vertical step, crack aperture and total width of disturbed ground should be measured separately. Reconstruct the total displacement across all strands and account for the original ground slope. Landslides or gravitational spreading can imitate normal-fault scarps, so map connections to the main fault and displaced markers rather than assuming every open crack is primary tectonic slip.