Use the ordered displacement vector $q=(x_1,y_1,x_2,y_2,x_3,y_3,x_4,y_4{)}^T$, and define $K$ by $V=\frac{1}{2}{q}^{T}Kq$. For an unstretched spring, the extension to first order is the dot product of the relative displacement with its equilibrium unit vector. Summing its squared extensions gives the central spring stiffness matrix. The four side springs and two diagonal springs giveThereforeThe kinetic energy is $\frac{m}{2}{\dot{q}}^T\dot{q}$, so a normal mode satisfies $Kq=m{\omega }^{2}q$. Keeping the factor $1/2$ in the definition of $V$ fixes the normalization of $K$ and the resulting frequencies.

Write the symmetry group of a square as $D_4=⟨r,s\mid r^4=s^2=e,srs=r^{-1}⟩$, with $r$ a counterclockwise quarter-turn and $s$ reflection across the horizontal axis. The conjugacy classes of the symmetry group of a square areHere $r^2$ is the half-turn and $r^3$ the clockwise quarter-turn. The two members of the fourth class reflect across the horizontal and vertical axes. Those of the fifth reflect across $y=x$ and $y=-x$. The identity and half-turn are central elements, giving singleton classes. The relations $srs=r^3$, $rs{r}^{-1}=r^{2}s$, and $r(rs)r^{-1}=r^{3}s$ show that each displayed pair consists of conjugate elements. More generally, conjugating $r^{k}s$ by $r$ changes $k$ to $k+2$, and conjugating by $s$ changes $k$ to $-k$; parity is preserved, so the two reflection pairs cannot merge.

In a one-dimensional representation, the relation $srs=r^{-1}$ forces the image of $r$ to equal its inverse. Thus $r$ and $s$ can be assigned independent signs $\alpha ,\beta \in \{1,-1\}$. These four choices give four irreducible characters. The ordinary two-dimensional action on the plane gives one more character: a rotation through $\theta$ has trace $2\cos \theta$, and a reflection has trace zero. Therefore the character table of the symmetry group of a square isFor instance the last row has squared character norm $(4+4)/8=1$, so the planar representation is irreducible; its inner product with each one-dimensional row is zero. All five rows obey character orthogonality. Their squared dimensions sum to $1+1+1+1+4=8$, so they exhaust the irreducible representations of this order-eight group.

If $Hv=\lambda v$ for a Hermitian matrix and a nonzero eigenvector, then $v^{†}Hv=\lambda v^{†}v$. The scalar on the left equals its complex conjugate because $H^{†}=H$, while $v^{†}v>0$. Thus $\lambda$ is real. If also $Hw=\mu w$, thenConsequently distinct eigenvalues force $v^{†}w=0$, proving orthogonality of their eigenvectors.

For a unitary matrix, $U^{†}U=I$ preserves the norm induced by the standard Hermitian inner product. If $Uv=\zeta v$ with $v\ne 0$, then $\Vert v{\Vert }^2=\Vert Uv{\Vert }^2=|\zeta {|}^2\Vert v{\Vert }^2$, so $\boxed{|\zeta |=1}$.

A normal mode is a fixed displacement pattern whose components share a common harmonic time dependence in the equations for small oscillations. Its normal frequency is the associated oscillation frequency; below, $\omega$ denotes the angular frequency, so the frequency in cycles per unit time is $\omega /(2\pi )$. A normal coordinate is the amplitude of one such pattern, chosen so that the linear equations decouple. For mass matrix $M$ and stiffness matrix $K$, the patterns solve $Kv={\omega }^{2}Mv$. Zero eigenvalues describe free rigid motions rather than periodic restoring oscillations.

For a piecewise continuous function of exponential order on $t\ge 0$, its unilateral Laplace transform isThe Bromwich inversion formula, at a point of continuity, isThe Bromwich contour is a straight vertical line, traversed upwards, with $\gamma$ to the right of the abscissa of convergence and hence of all singularities of the transform. Under the usual piecewise smoothness assumptions, at a jump the symmetric inversion limit gives the mean of the two one-sided values.

In the transform in the question, put $\lambda =e^{-t}$. Then $d\lambda =-e^{-t}dt$ andThus the restricted Mellin transform on the unit interval is the Laplace transform of $f(t)=h(e^{-t})$. Apply the Bromwich inversion formula with $t=-\log \lambda$, where $0<\lambda <1$, to obtainHere $\gamma$ lies in the convergence half-plane of $H$. At a discontinuity this formula has the same mean-value interpretation. The real natural logarithm defines ${\lambda }^{-s}$ unambiguously.

Commuting inverse. Assume $I-M$ is invertible, as required by the displayed inverse. Multiplication givesMultiply on both sides by $(I-M{)}^{-1}$ to obtainMore generally, a matrix commuting with an invertible matrix also commutes with its inverse.

Orthogonality. For a real skew-symmetric matrix $A$, $x^{T}Ax=0$. If $(I-A)x=0$ for a real vector $x$, multiply by $x^T$ to obtain $\Vert x{\Vert }^2=0$. Thus $I-A$ is invertible; the same argument applies to $I+A$. For the Cayley transform of a real skew-symmetric matrix,The factors commute by the first result, soTherefore $\boxed{N\text{ is an orthogonal matrix}.}$

Eigenvector mapping. If $Av=av$, then $(I-A{)}^{-1}v=(1-a{)}^{-1}v$. ConsequentlyThe denominator cannot vanish. In fact the eigenvalues of a real skew-symmetric matrix are purely imaginary: $v^{†}Av$ is the negative of its complex conjugate, so $a=-\overline{a}$. For $a=i\beta$ the corresponding eigenvalue of $N$ is $(1+i\beta )/(1-i\beta )=e^{2i\arctan \beta }$, consistent with orthogonality.

Spectrum in three dimensions. Let $Q$ be any real three-by-three orthogonal matrix. For a complex eigenvector $v$, $\Vert Qv\Vert =\Vert v\Vert$ implies $|\lambda |=1$. A real cubic characteristic polynomial has a real root, necessarily $+1$ or $-1$. Any nonreal roots form a complex-conjugate pair of unit modulus. If all three roots are real, two of the three signs agree and can be grouped as a pair with angle zero or $\pi$. Thus in every case the eigenvalues can be writtenTheir product is $\epsilon =\det Q$.

Determinant and fixed axis for the Cayley transform. Since $A^T=-A$,so $\boxed{\det N=1}$. This determinant argument works in every dimension. In dimension three the preceding spectrum result then guarantees an eigenvector with $Nx=x$. Set $z=(I-A{)}^{-1}x$. The fixed-vector equation says $(I+A)z=(I-A)z$, so $Az=0$ and $x=(I-A)z=z$. HenceThis is the axis fixed by the three-dimensional rotation. Equivalently, a real odd-dimensional skew-symmetric matrix has zero determinant, because $\det A=\det A^T=\det (-A)=-\det A$.

Local errors of the Euler schemes. Let $h=\Delta t>0$, and evaluate $f$ and its partial derivatives at an exact starting point $(y(t),t)$. DefineThe exact Taylor expansion is $y(t+h)=y+hf+\frac{1}{2}{h}^{2}A+O(h^3)$. The forward Euler method retains only $y+hf$, so its one-step defect, exact minus computed value, isFor the Backward Euler method, substitution of the exact solution into the scheme gives the residualbecause $f(y(t+h),t+h)=f+hA+O(h^2)$. Under local smoothness, solving the implicit step changes this residual by a factor $1+O(h)$, so the one-step defect has the same leading term:These are unnormalized local truncation errors, as in the question. A convention dividing by $h$ instead gives $O(h)$. Both schemes have first-order global error under the usual smoothness and stability hypotheses.

Euler stability. For the Dahlquist test equation $y^{\prime }=-\Lambda y$, put $z=h\Lambda$, with $\mathrm{Re}z>0$. The numerical amplification factors areStrict decay requires $|R|<1$; the boundary $|R|=1$ gives bounded but generally nondecaying values. For Forward Euler method,For positive real $\Lambda$, this reduces to $0<h\Lambda <2$. The absolute stability region in the $z$ plane is the disk $|z-1|\le 1$. For Backward Euler method, $|1+z|>1$ whenever $\mathrm{Re}z>0$, so it decays for every step size. It is A-stable and also L-stable, since $R_{BE}(z)\to 0$ for large $z$ in this half-plane. Here $z$ has the opposite sign from the alternative convention $z=h\lambda$ for $y^{\prime }=\lambda y$.

Predictor-corrector accuracy. The printed two-stage construction is Heun's method, equivalentlyAlthough the second stage resembles a Backward Euler method evaluation, it uses the predicted value, so the combined method is explicit. Expand that stage around the exact starting point, puttingThenand the numerical step is $y+hf+\frac{1}{2}{h}^{2}A+\frac{1}{4}{h}^{3}B+O(h^4)$. Along the exact solution, $y^{\prime \prime \prime }=B+f_{y}A$, so comparison with the exact Taylor expansion givesThus the local defect improves by one power of $h$, giving second-order global error under the usual hypotheses. For the linear test equation the defect reduces to $-h^3{\Lambda }^{3}y/6+O(h^4)$.

Predictor-corrector stability. Applied to the test equation, the predictor is $y^{(1)}=(1-z)y_n$ and the second value is $y^{(2)}=(1-z+z^2)y_n$. Their average givesFor real positive $z$, $R_H(z)=\frac{1}{2}[(z-1{)}^2+1]$, so strict decay holds exactly for $0<z<2$, the same real-axis step bound as Forward Euler method. For complex $z=x+iy$ the precise condition isThe region is bounded and does not contain the full right half-plane: for example $z=3$ gives $R_H=5/2$. It is therefore not A-stable. Its imaginary-axis values satisfy $|R_H(iy){|}^2=1+y^4/4>1$ for $y\ne 0$. Better accuracy has not inherited the unconditional stability of an actual implicit Backward Euler method step.