Software

Tsunamis are long-wavelength ocean waves generated by large-scale disturbances such as underwater earthquakes, volcanic eruptions, or landslides. Unlike ordinary wind-driven waves, which affect only the ocean surface, tsunamis involve the motion of the entire water column from the seabed to the surface. Because of this, they can travel across entire ocean basins with very little loss of energy. In the open ocean, a tsunami may pass almost unnoticed, with wave heights of less than a meter but wavelengths that can exceed hundreds of kilometers.

As a tsunami approaches shallower coastal waters, its behavior changes dramatically. The decreasing depth causes the wave to slow down while its height increases, often leading to powerful surges that can inundate coastal regions. This process, known as wave shoaling, concentrates the wave's energy and can produce devastating effects upon landfall. In some cases, the sea may first recede significantly before the arrival of the main wave, exposing the seabed and serving as a natural warning sign.

The prediction and analysis of tsunami propagation rely on mathematical models derived from the equations of fluid dynamics. In many practical applications, tsunamis are approximated as shallow-water waves, allowing their speed to be described in terms of the local ocean depth. This approximation enables the efficient computation of travel times and wavefront propagation across complex ocean basins. More advanced models incorporate dispersion and nonlinearity to capture finer details of wave evolution, particularly near coastlines.

Despite advances in modeling and monitoring, tsunami hazards remain difficult to predict with complete accuracy. Early warning systems combine seismic data, ocean buoys, and numerical simulations to estimate arrival times and potential impact. These systems play a crucial role in mitigating risk, but public awareness and preparedness are equally important in reducing the loss of life and property in vulnerable coastal regions.

Tsunami travel time describes how long it takes for a wave to propagate from its source to different locations across the ocean. This quantity is central to early warning systems, as it allows for rapid estimation of arrival times without resolving the full wave dynamics. Instead, one focuses on the motion of the leading wavefront, treating the ocean as a medium where wave speed depends on the local depth.

This propagation is governed by the eikonal equation, which relates the spatial variation of travel time to the wave speed, typically approximated by the shallow-water formula. In practice, this equation is solved numerically using efficient methods such as the Fast Marching Method, implemented for example in the Python library scikit-fmm. This approach makes it possible to compute tsunami arrival times quickly and accurately over large, complex ocean basins, providing a practical and widely used tool for large-scale tsunami modeling. An implementation for estimating tsunami travel time can be downloaded from my GitHub page. The pictures below present the estimated tsunami travel times for the 2004 Indian Ocean tsunami (Boxing Day tsunami) and the 2025 Kamchatka tsunami, and have been generated with the previously mentioned Python code.

The Complex-Step Newton Method is a numerical root-finding method that combines Newton iterations with complex-step derivative approximations. By avoiding the numerical cancellation errors that appear in finite-difference methods, the approach provides highly accurate derivative estimates and improved computational reliability. The complex-step iteration can be expressed as the equation or system of equations of the form:

where the product between the Jacobian matrix and the vector is approximated by the complex-step Jacobian approximation

for a small complex-step parameter \(h > 0\).

Therefore, the corresponding linear system can be solved using Krylov methods without forming the Jacobian but using instead this matrix-vector product approximation.

An important feature of the method is its convergence behavior. For systems of nonlinear equations, the method retains quadratic convergence even for relatively large values of the complex-step parameter. In the scalar case, the convergence is generally first order for finite values of the complex step, but approaches quadratic convergence as the complex-step parameter \(h\) tends to zero. In practice, this does not present a significant limitation, since the complex-step approximation remains numerically stable even for extremely small values of \(h\), unlike classical finite-difference methods which suffer from subtractive cancellation errors.

My work on this method includes theoretical analysis, numerical implementation, and the development of open-source software tools for scientific computing applications. The project is motivated by the broader goal of designing accurate and efficient numerical methods for nonlinear problems arising in applied mathematics and engineering.

An open-source implementation of the method has also been developed in C++ with a Python interface and can be accessed on PyPI and on GitHub.

References

• D. Mitsotakis, The complex-step Newton method and its convergence, Numerische Mathematik, 157 (2025), 993–1021, DOI:10.1007/s00211-025-01471-w (PDF) (Appendix)
• D. Mitsotakis, Computational Mathematics: An Introduction to Numerical Analysis and Scientific Computing with Python, Chapman and Hall/CRC, New York, 2023, DOI:10.1201/9781003287292

FEMd is a finite element library for partial differential equations in one and two space dimensions, especially designed for running efficiently simulations on modern personal computers. Its core is a header-only C++17 library that is used from Python, and can be integrated with OpenMP for parallel computing. The user provides a domain, a polynomial degree and boundary conditions, writes the equations in weak form, and FEMd assembles and solves the discrete problem. The form language follows the FEniCS project, so a weak formulation is written in code as it reads on paper.

In one dimension the basis functions can be B-splines of any degree and continuity, classical Lagrange elements or discontinuous Galerkin elements. Dirichlet, Neumann, clamped, periodic and Robin boundary conditions are built into the basis, and banded, cyclic and FFT solvers exploit the structure of the resulting matrices.

In two dimensions FEMd generates its own meshes. Its Delaunay mesh generator guarantees a minimum angle and handles holes and boundary markers, and the meshes can be refined locally, smoothed and coarsened. Structured meshes and meshes imported from Gmsh are also supported. On these meshes the library provides continuous \(P_k\) elements on triangles and \(Q_k\) elements on quadrilaterals of any degree, vector and mixed spaces for problems such as the Stokes and Navier-Stokes equations, Raviart-Thomas and Nédélec elements, and discontinuous Galerkin methods for conservation laws such as the shallow water equations. Solutions can be written to VTK files and visualized in ParaView.

Nonlinear problems are solved with Newton's method using symbolic Jacobians, or with Newton-Krylov methods based on the complex-step approximation of the Jacobian described above. The library also includes sparse direct and iterative linear solvers, implicit and explicit Runge-Kutta methods for time-dependent problems, and sparse eigenvalue solvers.

As an example, the boundary value problem \(-u'' + u = f\) in \((0, 1)\) with \(u(0) = 0\) and \(u(1) = 1\) is solved with cubic splines on 20 elements by

import numpy as np
import femd as fd
from femd import D, dx

V = fd.SplineSpace(np.linspace(0, 1, 21), 3, bc="dirichlet", left=0.0, right=1.0)
u, v = fd.TrialFunction(V), fd.TestFunction(V)
f = (np.pi**2 + 1) * fd.sin(np.pi * fd.x) + fd.x
uh = fd.solve(D(u) * D(v) * dx + u * v * dx, f * v * dx)

The library, its manual, the API reference and example notebooks can be downloaded from my GitHub page.