Numerical Accuracy

Each solution method in gFlex has a distinct accuracy profile.

Analytical solutions (SAS and SAS_NG)

The superposition-of-analytical-solutions methods are exact for a plate of constant elastic thickness under the NoOutsideLoads boundary condition (zero deflection at infinity). In one dimension the Green’s function is the exponential sinusoid of Eq. (3) in Wickert (2016); in two dimensions it is the zeroth-order Kelvin function \(\mathrm{kei}\), evaluated via scipy.special.kei(). Floating-point errors in these special-function evaluations are at machine precision and negligible in practice.

The only meaningful source of error is a domain that is too small: if the load produces non-trivial deflection at the domain boundary, the NoOutsideLoads assumption is violated and the solution is physically wrong — not because of numerics, but because the boundary condition does not match the situation. Use gflex.flexural_wavelengths() to estimate the flexural parameter \(\alpha\) and ensure the domain extends several \(\alpha\) beyond the loaded region. Unlike the finite difference solver, there is no grid-spacing error; output sampling density affects only the display, not the solution.

FFT spectral solver

The FFT solver is spectrally accurate for uniform elastic thickness: for a smooth load field the solution error decays faster than any finite power of the grid spacing, limited in practice only by floating-point arithmetic. For periodic boundary conditions (Periodic on all sides) the solution is exact to machine precision. For all other boundary conditions the load is zero-padded by \(4\alpha\) on each side before the transform, which approximates the NoOutsideLoads condition; the padding introduces a small error near the domain edges that decreases as the pad width increases relative to the load’s flexural footprint. For typical geoscience applications the \(4\alpha\) default padding is more than sufficient.

Because the FFT assumes a single, spatially uniform flexural rigidity \(D\), it cannot represent variable-\(T_e\) problems; those require the finite difference solver.

Finite difference solver (FD)

Comparison of numerical (FD) and analytical (SAS) solutions in one dimension (a) and two dimensions (c), and their differences (b, d), for a 100 km central line load / circular load. The 1–2 m offset in (b) is due primarily to the NoOutsideLoads BC of the analytical solution versus the 0Displacement0Slope BC of the FD solution; the cross-shaped residual in (d) reflects boundary effects along the longer diagonal boundaries. Reproduced from Wickert (2016), Fig. 3; CC BY 3.0.

---

2-D finite-difference solver: second-order convergence

The 2-D finite-difference solver (Method = 'FD', PlateSolutionType = 'vWC1994') is second-order accurate in space: halving the grid spacing \(\Delta x\) reduces the numerical error by a factor of approximately four (\(\mathcal{O}(\Delta x^2)\)).

This is verified two ways in gFlex:

1. Method of Manufactured Solutions (MMS) — a formal unit test (tests.test_2D_FD.test_2d_fd_convergence_order()) that runs with every CI build. A cosine deflection field is chosen as the exact solution; the corresponding load is derived analytically and fed to the solver. Three grid refinements (N = 30, 60, 120 cells per side) on a periodic domain confirm a convergence rate > 1.8 at all refinement levels.

2. Kelvin-function benchmark — a grid-convergence study comparing the FD solver to the analytical Kelvin-function solution for an infinite plate under a point load. Results are summarised in the table below.

Kelvin-function benchmark setup

Young’s modulus \(E\)	65 GPa
Elastic thickness \(T_e\)	10 km
Poisson’s ratio \(\nu\)	0.25
Mantle density \(\rho_m\)	3300 kg m⁻³
Infill density \(\rho_\text{fill}\)	0 kg m⁻³ (air)
\(g\)	9.81 m s⁻²
Flexural parameter \(\alpha\)	≈ 21 km
Domain	600 km × 600 km, point load at centre
Boundary conditions	0Moment0Shear on all four sides
Grid spacings tested	\(\Delta x\) = 10 000, 5 000, 2 500, 1 250 m

Convergence orders (coarse → fine)

The table gives the local convergence order between successive grid refinements at several distances from the load, expressed as multiples of the flexural parameter α.

\(r/\alpha\)	\(\Delta x\) 10→5 km	\(\Delta x\) 5→2.5 km	\(\Delta x\) 2.5→1.25 km
0.97	2.20	2.07	2.01
1.46	2.68	2.18	2.04
1.95	1.07	1.88	1.97
2.92	−0.07	1.77	1.95

At the finest spacing (\(\Delta x = 1250\) m, \(\Delta x/\alpha \approx 0.06\)) relative errors are 0.02–0.1 % everywhere except very close to the singularity of the point load.

The sub-second-order rates at \(r/\alpha\) = 1.95 and 2.92 on the coarsest grids reflect pre-asymptotic effects near the forebulge zero-crossing, where the solution changes sign and the absolute error temporarily dominates the relative error. All rates converge to ≈ 2 at finer resolution, confirming asymptotic second-order behaviour.

Practical guidance

For most geoscience applications, \(\Delta x / \alpha \leq 0.1\) (ten cells per flexural parameter) is sufficient to keep errors below 1 %. At \(\Delta x / \alpha \approx 0.5\) (two cells per flexural parameter) errors can be several percent, particularly near the load centre and the forebulge. Use gflex.flexural_wavelengths() to estimate \(\alpha\) before choosing a grid spacing.