Source author record

Christian Lubich

Christian Lubich appears in the imported research catalog. Authorship, coauthor and topic links are available while profile ownership is still unclaimed.

ResearcherUnclaimed source record

Catalog footprint

What is connected

28works
12topics
4close collaborators

Actions

Connect this record

Log in to claim

Research graph

See the researcher in context

Open full explorer

Inspect adjacent papers, topics, institutions and collaborators without losing the researcher page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

Published work

28 published item(s)

preprint2022arXiv

On a large-stepsize integrator for charged-particle dynamics

Xiao and Qin [Computer Physics Comm., 265:107981, 2021] recently proposed a remarkably simple modification of the Boris algorithm to compute the guiding centre of the highly oscillatory motion of a charged particle with step sizes that are much larger than the period of gyrorotations. They gave strong numerical evidence but no error analysis. This paper provides an analysis of the large-stepsize modified Boris method in a setting that has a strong non-uniform magnetic field and moderately bounded velocities, considered over a fixed finite time interval. The error analysis is based on comparing the modulated Fourier expansions of the exact and numerical solutions, for which the differential equations of the dominant terms are derived explicitly. Numerical experiments illustrate and complement the theoretical results.

preprint2022arXiv

Polarized high-frequency wave propagation beyond the nonlinear Schrödinger approximation

This paper studies highly oscillatory solutions to a class of systems of semilinear hyperbolic equations with a small parameter, in a setting that includes Klein--Gordon equations and the Maxwell--Lorentz system. The interest here is in solutions that are polarized in the sense that up to a small error, the oscillations in the solution depend on only one of the frequencies that satisfy the dispersion relation with a given wave vector appearing in the initial wave packet. The construction and analysis of such polarized solutions is done using modulated Fourier expansions. This approach includes higher harmonics and yields approximations to polarized solutions that are of arbitrary order in the small parameter, going well beyond the known first-order approximation via a nonlinear Schrödinger equation. The given construction of polarized solutions is explicit, uses in addition a linear Schrödinger equation for each further order of approximation, and is accessible to direct numerical approximation.

preprint2022arXiv

Rank-$1$ matrix differential equations for structured eigenvalue optimization

A new approach to solving eigenvalue optimization problems for large structured matrices is proposed and studied. The class of optimization problems considered is related to computing structured pseudospectra and their extremal points, and to structured matrix nearness problems such as computing the structured distance to instability or to singularity. The structure can be a general linear structure and includes, for example, large matrices with a given sparsity pattern, matrices with given range and co-range, and Hamiltonian matrices. Remarkably, the eigenvalue optimization can be performed on the manifold of complex (or real) rank-1 matrices, which yields a significant reduction of storage and in some cases of the computational cost. The method relies on a constrained gradient system and the projection of the gradient onto the tangent space of the manifold of complex rank-$1$ matrices. It is shown that near a local minimizer this projection is very close to the identity map, and so the computationally favorable rank-1 projected system behaves locally like the %computationally expensive gradient system.

preprint2022arXiv

Rank-adaptive time integration of tree tensor networks

A rank-adaptive integrator for the approximate solution of high-order tensor differential equations by tree tensor networks is proposed and analyzed. In a recursion from the leaves to the root, the integrator updates bases and then evolves connection tensors by a Galerkin method in the augmented subspace spanned by the new and old bases. This is followed by rank truncation within a specified error tolerance. The memory requirements are linear in the order of the tensor and linear in the maximal mode dimension. The integrator is robust to small singular values of matricizations of the connection tensors. Up to the rank truncation error, which is controlled by the given error tolerance, the integrator preserves norm and energy for Schrodinger equations, and it dissipates the energy in gradient systems. Numerical experiments with a basic quantum spin system illustrate the behavior of the proposed algorithm.

preprint2022arXiv

Time-dependent electromagnetic scattering from thin layers

The scattering of electromagnetic waves from obstacles with wave-material interaction in thin layers on the surface is described by generalized impedance boundary conditions, which provide effective approximate models. In particular, this includes a thin coating around a perfect conductor and the skin effect of a highly conducting material. The approach taken in this work is to derive, analyse and discretize a system of time-dependent boundary integral equations that determines the tangential traces of the scattered electric and magnetic fields. In a familiar second step, the fields are evaluated in the exterior domain by a representation formula, which uses the time-dependent potential operators of Maxwell's equations. The time-dependent boundary integral equationis discretized with Runge--Kutta based convolution quadrature in time and Raviart--Thomas boundary elements in space. Using the frequency-explicit bounds from the well-posedness analysis given here together with known approximation properties of the numerical methods, the full discretization is proved to be stable and convergent, with explicitly given rates in the case of sufficient regularity. Taking the same Runge--Kutta based convolution quadrature for discretizing the time-dependent representation formulas, the optimal order of convergence is obtained away from the scattering boundary, whereas an order reduction occurs close to the boundary. The theoretical results are illustrated by numerical experiments.

preprint2021arXiv

Finding the nearest passive or non-passive system via Hamiltonian eigenvalue optimization

We propose and study an algorithm for computing a nearest passive system to a given non-passive linear time-invariant system (with much freedom in the choice of the metric defining `nearest', which may be restricted to structured perturbations), and also a closely related algorithm for computing the structured distance of a given passive system to non-passivity. Both problems are addressed by solving eigenvalue optimization problems for Hamiltonian matrices that are constructed from perturbed system matrices. The proposed algorithms are two-level methods that optimize the Hamiltonian eigenvalue of smallest positive real part over perturbations of a fixed size in the inner iteration, using a constrained gradient flow. They optimize over the perturbation size in the outer iteration, which is shown to converge quadratically in the typical case of a defective coalescence of simple eigenvalues approaching the imaginary axis. For large systems, we propose a variant of the algorithm that takes advantage of the inherent low-rank structure of the problem. Numerical experiments illustrate the behavior of the proposed algorithms.

preprint2021arXiv

Large-stepsize integrators for charged-particle dynamics over multiple time scales

The Boris algorithm, a closely related variational integrator and a newly proposed filtered variational integrator are studied when they are used to numerically integrate the equations of motion of a charged particle in a non-uniform strong magnetic field, taking step sizes that are much larger than the period of the Larmor rotations. For the Boris algorithm and the standard (unfiltered) variational integrator, satisfactory behaviour is only obtained when the component of the initial velocity orthogonal to the magnetic field is filtered out. The particle motion shows varying behaviour over multiple time scales: fast Larmor rotation, guiding centre motion, slow perpendicular drift, near-conservation of the magnetic moment over very long times and conservation of energy for all times. Using modulated Fourier expansions of the exact and numerical solutions, it is analysed to which extent this behaviour is reproduced by the three numerical integrators used with large step sizes.

preprint2020arXiv

A convergent algorithm for forced mean curvature flow driven by diffusion on the surface

The evolution of a closed two-dimensional surface driven by both mean curvature flow and a reaction--diffusion process on the surface is formulated into a system, which couples the velocity law not only to the surface partial differential equation but also to the evolution equations for the geometric quantities, namely the normal vector and the mean curvature on the surface. Two algorithms are considered for the obtained system. Both methods combine surface finite elements as a space discretisation and linearly implicit backward difference formulae for time integration. Based on our recent results for mean curvature flow, one of the algorithms directly admits a convergence proof for its full discretisation in the case of finite elements of polynomial degree at least two and backward difference formulae of orders two to five. Numerical examples are provided to support and complement the theoretical convergence results (demonstrating the convergence properties of the method without error estimate), and demonstrate the effectiveness of the methods in simulating a three-dimensional tumour growth model.

preprint2020arXiv

A convergent evolving finite element algorithm for Willmore flow of closed surfaces

A proof of convergence is given for a novel evolving surface finite element semi-discretization of Willmore flow of closed two-dimensional surfaces, and also of surface diffusion flow. The numerical method proposed and studied here discretizes fourth-order evolution equations for the normal vector and mean curvature, reformulated as a system of second-order equations, and uses these evolving geometric quantities in the velocity law interpolated to the finite element space. This numerical method admits a convergence analysis in the case of continuous finite elements of polynomial degree at least two. The error analysis combines stability estimates and consistency estimates to yield optimal-order $H^1$-norm error bounds for the computed surface position, velocity, normal vector and mean curvature. The stability analysis is based on the matrix--vector formulation of the finite element method and does not use geometric arguments. The geometry enters only into the consistency estimates. Numerical experiments illustrate and complement the theoretical results.

preprint2020arXiv

Higher-order linearly implicit full discretization of the Landau--Lifshitz--Gilbert equation

For the Landau--Lifshitz--Gilbert (LLG) equation of micromagnetics we study linearly implicit backward difference formula (BDF) time discretizations up to order $5$ combined with higher-order non-conforming finite element space discretizations, which are based on the weak formulation due to Alouges but use approximate tangent spaces that are defined by $L^2$-averaged instead of nodal orthogonality constraints. We prove stability and optimal-order error bounds in the situation of a sufficiently regular solution. For the BDF methods of orders $3$ to~$5$, this requires %a mild time step restriction $τ\leqslant ch$ and that the damping parameter in the LLG equations be above a positive threshold; this condition is not needed for the A-stable methods of orders $1$ and $2$, for which furthermore a discrete energy inequality irrespective of solution regularity is proved.

preprint2020arXiv

Measuring the stability of spectral clustering

As an indicator of the stability of spectral clustering of an undirected weighted graph into $k$ clusters, the $k$th spectral gap of the graph Laplacian is often considered. The $k$th spectral gap is characterized in this paper as an unstructured distance to ambiguity, namely as the minimal distance of the Laplacian to arbitrary symmetric matrices with vanishing $k$th spectral gap. As a conceptually more appropriate measure of stability, the structured distance to ambiguity of the $k$-clustering is introduced as the minimal distance of the Laplacian to Laplacians of graphs with the same vertices and edges but with weights that are perturbed such that the $k$th spectral gap vanishes. To compute a solution to this matrix nearness problem, a two-level iterative algorithm is proposed that uses a constrained gradient system of matrix differential equations in the inner iteration and a one-dimensional optimization of the perturbation size in the outer iteration. The structured and unstructured distances to ambiguity are compared on some example graphs. The numerical experiments show, in particular, that selecting the number $k$ of clusters according to the criterion of maximal stability can lead to different results for the structured and unstructured stability indicators.

preprint2020arXiv

Time integration of tree tensor networks

Dynamical low-rank approximation by tree tensor networks is studied for the data-sparse approximation to large time-dependent data tensors and unknown solutions of tensor differential equations. A time integration method for tree tensor networks of prescribed tree rank is presented and analyzed. It extends the known projector-splitting integrators for dynamical low-rank approximation by matrices and Tucker tensors and is shown to inherit their favorable properties. The integrator is based on recursively applying the Tucker tensor integrator. In every time step, the integrator climbs up and down the tree: it uses a recursion that passes from the root to the leaves of the tree for the construction of initial value problems on subtree tensor networks using appropriate restrictions and prolongations, and another recursion that passes from the leaves to the root for the update of the factors in the tree tensor network. The integrator reproduces given time-dependent tree tensor networks of the specified tree rank exactly and is robust to the typical presence of small singular values in matricizations of the connection tensors, in contrast to standard integrators applied to the differential equations for the factors in the dynamical low-rank approximation by tree tensor networks.

preprint2020arXiv

Time-dependent acoustic scattering from generalized impedance boundary conditions via boundary elements and convolution quadrature

Generalized impedance boundary conditions are effective, approximate boundary conditions that describe scattering of waves in situations where the wave interaction with the material involves multiple scales. In particular, this includes materials with a thin coating (with the thickness of the coating as the small scale) and strongly absorbing materials. For the acoustic scattering from generalized impedance boundary conditions, the approach taken here first determines the Dirichlet and Neumann boundary data from a system of time-dependent boundary integral equations with the usual boundary integral operators, and then the scattered wave is obtained from the Kirchhoff representation. The system of time-dependent boundary integral equations is discretized by boundary elements in space and convolution quadrature in time. The well-posedness of the problem and the stability of the numerical discretization rely on the coercivity of the Calderón operator for the Helmholtz equation with frequencies in a complex half-plane. Convergence of optimal order in the natural norms is proved for the full discretization. Numerical experiments illustrate the behaviour of the proposed numerical method.

preprint2016arXiv

Long-term analysis of semilinear wave equations with slowly varying wave speed

A semilinear wave equation with slowly varying wave speed is considered in one to three space dimensions on a bounded interval, a rectangle or a box, respectively. It is shown that the action, which is the harmonic energy divided by the wave speed and multiplied with the diameter of the spatial domain, is an adiabatic invariant: it remains nearly conserved over long times, longer than any fixed power of the time scale of changes in the wave speed in the case of one space dimension, and longer than can be attained by standard perturbation arguments in the two- and three-dimensional cases. The long-time near-conservation of the action yields long-time existence of the solution. The proofs use modulated Fourier expansions in time.

preprint2015arXiv

A-stable time discretizations preserve maximal parabolic regularity

It is shown that for a parabolic problem with maximal $L^p$-regularity (for $1<p<\infty$), the time discretization by a linear multistep method or Runge--Kutta method has maximal $\ell^p$-regularity uniformly in the stepsize if the method is A-stable (and satisfies minor additional conditions). In particular, the implicit Euler method, the Crank-Nicolson method, the second-order backward difference formula (BDF), and the Radau IIA and Gauss Runge--Kutta methods of all orders preserve maximal regularity. The proof uses Weis' characterization of maximal $L^p$-regularity in terms of $R$-boundedness of the resolvent, a discrete operator-valued Fourier multiplier theorem by Blunck, and generating function techniques that have been familiar in the stability analysis of time discretization methods since the work of Dahlquist. The A($α$)-stable higher-order BDF methods have maximal $\ell^p$-regularity under an $R$-boundedness condition in a larger sector. As an illustration of the use of maximal regularity in the error analysis of discretized nonlinear parabolic equations, it is shown how error bounds are obtained without using any growth condition on the nonlinearity or for nonlinearities having singularities.

preprint2015arXiv

Numerical analysis of parabolic problems with dynamic boundary conditions

Space and time discretizations of parabolic differential equations with dynamic boundary conditions are studied in a weak formulation that fits into the standard abstract formulation of parabolic problems, just that the usual L^2(Ω) inner product is replaced by an L^2(Ω) \oplus L^2(Γ) inner product. The class of parabolic equations considered includes linear problems with time- and space-dependent coefficients and semilinear problems such as reaction-diffusion on a surface coupled to diffusion in the bulk. The spatial discretization by finite elements is studied in the proposed framework, with particular attention to the error analysis of the Ritz map for the elliptic bilinear form in relation to the inner product, both of which contain boundary integrals. The error analysis is done for both polygonal and smooth domains. We further consider mass lumping, which enables us to use exponential integrators and bulk-surface splitting for time integration.

preprint2015arXiv

Time integration of tensor trains

A robust and efficient time integrator for dynamical tensor approximation in the tensor train or matrix product state format is presented. The method is based on splitting the projector onto the tangent space of the tensor manifold. The algorithm can be used for updating time-dependent tensors in the given data-sparse tensor train / matrix product state format and for computing an approximate solution to high-dimensional tensor differential equations within this data-sparse format. The formulation, implementation and theoretical properties of the proposed integrator are studied, and numerical experiments with problems from quantum molecular dynamics and with iterative processes in the tensor train format are included.

preprint2015arXiv

Unifying time evolution and optimization with matrix product states

We show that the time-dependent variational principle provides a unifying framework for time-evolution methods and optimisation methods in the context of matrix product states. In particular, we introduce a new integration scheme for studying time-evolution, which can cope with arbitrary Hamiltonians, including those with long-range interactions. Rather than a Suzuki-Trotter splitting of the Hamiltonian, which is the idea behind the adaptive time-dependent density matrix renormalization group method or time-evolving block decimation, our method is based on splitting the projector onto the matrix product state tangent space as it appears in the Dirac-Frenkel time-dependent variational principle. We discuss how the resulting algorithm resembles the density matrix renormalization group (DMRG) algorithm for finding ground states so closely that it can be implemented by changing just a few lines of code and it inherits the same stability and efficiency. In particular, our method is compatible with any Hamiltonian for which DMRG can be implemented efficiently and DMRG is obtained as a special case of imaginary time evolution with infinite time step.

preprint2014arXiv

Numerical integrators for motion under a strong constraining force

This paper deals with the numerical integration of Hamiltonian systems in which a stiff anharmonic potential causes highly oscillatory solution behavior with solution-dependent frequencies. The impulse method, which uses micro- and macro-steps for the integration of fast and slow parts, respectively, does not work satisfactorily on such problems. Here it is shown that variants of the impulse method with suitable projection preserve the actions as adiabatic invariants and yield accurate approximations, with macro-stepsizes that are not restricted by the stiffness parameter.

preprint2013arXiv

A projector-splitting integrator for dynamical low-rank approximation

The dynamical low-rank approximation of time-dependent matrices is a low-rank factorization updating technique. It leads to differential equations for factors of the matrices, which need to be solved numerically. We propose and analyze a fully ex- plicit, computationally inexpensive integrator that is based on splitting the orthogonal projector onto the tangent space of the low-rank manifold. As is shown by theory and illustrated by numerical experiments, the integrator enjoys robustness properties that are not shared by any standard numerical integrator. This robustness can be exploited to change the rank adaptively. Another application is in optimization algorithms for low-rank matrices where truncation back to the given low rank can be done efficiently by applying a step of the integrator proposed here.

preprint2013arXiv

Plane wave stability of the split-step Fourier method for the nonlinear Schrödinger equation

Plane wave solutions to the cubic nonlinear Schrödinger equation on a torus have recently been shown to behave orbitally stable. Under generic perturbations of the initial data that are small in a high-order Sobolev norm, plane waves are stable over long times that extend to arbitrary negative powers of the smallness parameter. The present paper studies the question as to whether numerical discretizations by the split-step Fourier method inherit such a generic long-time stability property. This can indeed be shown under a condition of linear stability and a non-resonance condition. They can both be verified if the time step-size is restricted by a CFL condition in the case of a constant plane wave. The proof first uses a Hamiltonian reduction and transformation and then modulated Fourier expansions in time. It provides detailed insight into the structure of the numerical solution.

preprint2012arXiv

Energy separation in oscillatory Hamiltonian systems without any non-resonance condition

We consider multiscale Hamiltonian systems in which harmonic oscillators with several high frequencies are coupled to a slow system. It is shown that the oscillatory energy is nearly preserved over long times eps^{-N} for arbitrary N>1, where eps^{-1} is the size of the smallest high frequency. The result is uniform in the frequencies and does not require non-resonance conditions.

preprint2012arXiv

Sobolev stability of plane wave solutions to the cubic nonlinear Schrödinger equation on a torus

It is shown that plane wave solutions to the cubic nonlinear Schrödinger equation on a torus behave orbitally stable under generic perturbations of the initial data that are small in a high-order Sobolev norm, over long times that extend to arbitrary negative powers of the smallness parameter. The perturbation stays small in the same Sobolev norm over such long times. The proof uses a Hamiltonian reduction and transformation and, alternatively, Birkhoff normal forms or modulated Fourier expansions in time.

preprint2011arXiv

Operator splitting for partial differential equations with Burgers nonlinearity

We provide a new analytical approach to operator splitting for equations of the type $u_t=Au+u u_x$ where $A$ is a linear differential operator such that the equation is well-posed. Particular examples include the viscous Burgers' equation, the Korteweg-de Vries (KdV) equation, the Benney-Lin equation, and the Kawahara equation. We show that the Strang splitting method converges with the expected rate if the initial data are sufficiently regular. In particular, for the KdV equation we obtain second-order convergence in $H^r$ for initial data in $H^{r+5}$ with arbitrary $r\ge 1$.

preprint2010arXiv

Symplectic Integration of Post-Newtonian Equations of Motion with Spin

We present a non-canonically symplectic integration scheme tailored to numerically computing the post-Newtonian motion of a spinning black-hole binary. Using a splitting approach we combine the flows of orbital and spin contributions. In the context of the splitting, it is possible to integrate the individual terms of the spin-orbit and spin-spin Hamiltonians analytically, exploiting the special structure of the underlying equations of motion. The outcome is a symplectic, time-reversible integrator, which can be raised to arbitrary order by composition. A fourth-order version is shown to give excellent behavior concerning error growth and conservation of energy and angular momentum in long-term simulations. Favorable properties of the integrator are retained in the presence of weak dissipative forces due to radiation damping in the full post-Newtonian equations.

preprint2005arXiv

Fast and oblivious convolution quadrature

We give an algorithm to compute $N$ steps of a convolution quadrature approximation to a continuous temporal convolution using only $O(N \log N)$ multiplications and $O(\log N)$ active memory. The method does not require evaluations of the convolution kernel, but instead $O(\log N)$ evaluations of its Laplace transform, which is assumed sectorial. The algorithm can be used for the stable numerical solution with quasi-optimal complexity of linear and nonlinear integral and integro-differential equations of convolution type. In a numerical example we apply it to solve a subdiffusion equation with transparent boundary conditions.