Discretization Of Partial Differential Equations Research Articles

Discrete and continuous mathematical models of sharp-fronted collective cell migration and invasion

Mathematical models describing the spatial spreading and invasion of populations of biological cells are often developed in a continuum modelling framework using reaction–diffusion equations. While continuum models based on linear diffusion are routinely employed and known to capture key experimental observations, linear diffusion fails to predict well-defined sharp fronts that are often observed experimentally. This observation has motivated the use of nonlinear degenerate diffusion; however, these nonlinear models and the associated parameters lack a clear biological motivation and interpretation. Here, we take a different approach by developing a stochastic discrete lattice-based model incorporating biologically inspired mechanisms and then deriving the reaction–diffusion continuum limit. Inspired by experimental observations, agents in the simulation deposit extracellular material, which we call a substrate , locally onto the lattice, and the motility of agents is taken to be proportional to the substrate density. Discrete simulations that mimic a two-dimensional circular barrier assay illustrate how the discrete model supports both smooth and sharp-fronted density profiles depending on the rate of substrate deposition. Coarse-graining the discrete model leads to a novel partial differential equation (PDE) model whose solution accurately approximates averaged data from the discrete model. The new discrete model and PDE approximation provide a simple, biologically motivated framework for modelling the spreading, growth and invasion of cell populations with well-defined sharp fronts. Open-source Julia code to replicate all results in this work is available on GitHub .

Solving inverse problems in physics by optimizing a discrete loss: Fast and accurate learning without neural networks.

In recent years, advances in computing hardware and computational methods have prompted a wealth of activities for solving inverse problems in physics. These problems are often described by systems of partial differential equations (PDEs). The advent of machine learning has reinvigorated the interest in solving inverse problems using neural networks (NNs). In these efforts, the solution of the PDEs is expressed as NNs trained through the minimization of a loss function involving the PDE. Here, we show how to accelerate this approach by five orders of magnitude by deploying, instead of NNs, conventional PDE approximations. The framework of optimizing a discrete loss (ODIL) minimizes a cost function for discrete approximations of the PDEs using gradient-based and Newton's methods. The framework relies on grid-based discretizations of PDEs and inherits their accuracy, convergence, and conservation properties. The implementation of the method is facilitated by adopting machine-learning tools for automatic differentiation. We also propose a multigrid technique to accelerate the convergence of gradient-based optimizers. We present applications to PDE-constrained optimization, optical flow, system identification, and data assimilation. We compare ODIL with the popular method of physics-informed neural networks and show that it outperforms it by several orders of magnitude in computational speed while having better accuracy and convergence rates. We evaluate ODIL on inverse problems involving linear and nonlinear PDEs including the Navier-Stokes equations for flow reconstruction problems. ODIL bridges numerical methods and machine learning and presents a powerful tool for solving challenging, inverse problems across scientific domains.

Symmetrized Two-Scale Finite Element Discretizations for Partial Differential Equations with Symmetric Solutions

Abstract In this paper, some symmetrized two-scale finite element methods are proposed for a class of partial differential equations with symmetric solutions. With these methods, the finite element approximation on a fine tensor-product grid is reduced to the finite element approximations on a much coarser grid and a univariant fine grid. It is shown by both theory and numerics including electronic structure calculations that the resulting approximations still maintain an asymptotically optimal accuracy. By symmetrized two-scale finite element methods, the computational cost can be reduced further by a factor of 𝑑 approximately compared with two-scale finite element methods when Ω = ( 0 , 1 ) d \Omega=(0,1)^{d} . Consequently, symmetrized two-scale finite element methods reduce computational cost significantly.

Reduced order model predictive control for parametrized parabolic partial differential equations

Model Predictive Control (MPC) is a well-established approach to solve infinite horizon optimal control problems. Since optimization over an infinite time horizon is generally infeasible, MPC determines a suboptimal feedback control by repeatedly solving finite time optimal control problems. Although MPC has been successfully used in many applications, applying MPC to large-scale systems – arising, e.g., through discretization of partial differential equations – requires the solution of high-dimensional optimization problems and thus poses immense computational effort. We consider systems governed by parametrized parabolic partial differential equations and employ the reduced basis method as a low-dimensional surrogate model for the finite time optimal control problem. The reduced order optimal control serves as feedback control for the original large-scale system. We analyze the proposed RB-MPC approach by first developing a posteriori error bounds for the errors in the optimal control and associated cost functional. These bounds can be evaluated efficiently in an offline-online computational procedure and allow us to guarantee asymptotic stability of the closed-loop system using the RB-MPC approach in several practical scenarios. We also propose an adaptive strategy to choose the prediction horizon of the finite time optimal control problem. Numerical results are presented to illustrate the theoretical properties of our approach.

Parsimonious physics-informed random projection neural networks for initial value problems of ODEs and index-1 DAEs.

We present a numerical method based on random projections with Gaussian kernels and physics-informed neural networks for the numerical solution of initial value problems (IVPs) of nonlinear stiff ordinary differential equations (ODEs) and index-1 differential algebraic equations (DAEs), which may also arise from spatial discretization of partial differential equations (PDEs). The internal weights are fixed to ones while the unknown weights between the hidden and output layer are computed with Newton's iterations using the Moore-Penrose pseudo-inverse for low to medium scale and sparse QR decomposition with L 2 regularization for medium- to large-scale systems. Building on previous works on random projections, we also prove its approximation accuracy. To deal with stiffness and sharp gradients, we propose an adaptive step-size scheme and address a continuation method for providing good initial guesses for Newton iterations. The "optimal" bounds of the uniform distribution from which the values of the shape parameters of the Gaussian kernels are sampled and the number of basis functions are "parsimoniously" chosen based on bias-variance trade-off decomposition. To assess the performance of the scheme in terms of both numerical approximation accuracy and computational cost, we used eight benchmark problems (three index-1 DAEs problems, and five stiff ODEs problems including the Hindmarsh-Rose neuronal model of chaotic dynamics and the Allen-Cahn phase-field PDE). The efficiency of the scheme was compared against two stiff ODEs/DAEs solvers, namely, ode15s and ode23t solvers of the MATLAB ODE suite as well as against deep learning as implemented in the DeepXDE library for scientific machine learning and physics-informed learning for the solution of the Lotka-Volterra ODEs included in the demos of the library. A software/toolbox in Matlab (that we call RanDiffNet) with demos is also provided.

Space and chaos‐expansion Galerkin proper orthogonal decomposition low‐order discretization of partial differential equations for uncertainty quantification

AbstractThe quantification of multivariate uncertainties in partial differential equations can easily exceed any computing capacity unless proper measures are taken to reduce the complexity of the model. In this work, we propose a multidimensional Galerkin proper orthogonal decomposition that optimally reduces each dimension of a tensorized product space. We provide the analytical framework and results that define and quantify the low‐dimensional approximation. We illustrate its application for uncertainty modeling with polynomial chaos expansions and show its efficiency in a numerical example.

Parallel block preconditioners for virtual element discretizations of the time-dependent Maxwell equations

The focus of this study is the construction and numerical validation of parallel block preconditioners for low order virtual element discretizations of the three-dimensional Maxwell equations. The virtual element method (VEM) is a recent technology for the numerical approximation of partial differential equations (PDEs), that generalizes finite elements to polytopal computational grids. So far, VEM has been extended to several problems described by PDEs, and recently also to the time-dependent Maxwell equations. When the time discretization of PDEs is performed implicitly, at each time-step a large-scale and ill-conditioned linear system must be solved, that, in case of Maxwell equations, is particularly challenging, because of the presence of both H(div) and H(curl) discretization spaces. We propose here a parallel preconditioner, that exploits the Schur complement block factorization of the linear system matrix and consists of a Jacobi preconditioner for the H(div) block and an auxiliary space preconditioner for the H(curl) block. Several parallel numerical tests have been performed to study the robustness of the solver with respect to mesh refinement, shape of polyhedral elements, time step size and the VEM stabilization parameter.

Isogeometric Analysis Approximation of Linear Elliptic Equations with L1 Data

Isogeometric Analysis (IgA) is a recent technique for the discretization of Partial Differential Equations (PDEs). The main feature of the method is the ability to maintain the same exact description of the computational domain geometry throughout the analysis process, including refinement. In the present paper, we consider, in dimension d >= 2 the Isogeometric Analysis approximation of second order elliptic equations in divergence form with right-hand side in L1 . We assume that the family of meshes is shape regular and satisfies the discrete maximum principle. When the right-hand side belongs to L1(\Omega), we prove that the unique solution of the discrete problem converges to the unique renormalized solution in W01,q(\Omega), 1 <= q < d/(d-1) . We also prove some error estimates and include numerical tests for data with low smoothness.

A mathematical model that integrates cardiac electrophysiology, mechanics, and fluid dynamics: Application to the human left heart.

We propose a mathematical and numerical model for the simulation of the heart function that couples cardiac electrophysiology, active and passive mechanics and hemodynamics, and includes reduced models for cardiac valves and the circulatory system. Our model accounts for the major feedback effects among the different processes that characterize the heart function, including electro-mechanical and mechano-electrical feedback as well as force-strain and force-velocity relationships. Moreover, it provides a three-dimensional representation of both the cardiac muscle and the hemodynamics, coupled in a fluid-structure interaction (FSI) model. By leveraging the multiphysics nature of the problem, we discretize it in time with a segregated electrophysiology-force generation-FSI approach, allowing for efficiency and flexibility in the numerical solution. We employ a monolithic approach for the numerical discretization of the FSI problem. We use finite elements for the spatial discretization of partial differential equations. We carry out a numerical simulation on a realistic human left heart model, obtaining results that are qualitatively and quantitatively in agreement with physiological ranges and medical images.

Differential Evolution Based Numerical Variable Speed Limit Control Method with a Non-Equilibrium Traffic Model

This paper introduces a numerical variable speed limit (VSL) control method on a motorway, modeled by the system of partial differential equations (PDEs) of a non- equilibrium continuum traffic model. The method consists of a macroscopic simulation (i.e., numerical solution of the system of PDEs of the continuum model), introduction of the solution-based cost function and numerical optimization with a differential evolution algorithm (DE). Due to the numerical solution scheme, the method enables application of a wide range of continuum traffic models without prior discretization of PDEs. In this way, the method overcomes the limitations of the basic continuum models and represents a step towards more accurate traffic modelling in control strategies. In this paper, we determine optimal variable speed limits with the DE algorithm on a motorway section modeled by the modified switching curve model, which is a non-equilibrium continuum model consistent with the three-phase traffic flow theory. The effectiveness of the determined variable speed limits is validated using microsimulations of the test section, which show promising reductions of queue lengths and number of stops.

Accelerating Algebraic Multigrid Methods via Artificial Neural Networks

We present a novel deep learning-based algorithm to accelerate—through the use of Artificial Neural Networks (ANNs)—the convergence of Algebraic Multigrid (AMG) methods for the iterative solution of the linear systems of equations stemming from finite element discretizations of Partial Differential Equations (PDE). We show that ANNs can be successfully used to predict the strong connection parameter that enters in the construction of the sequence of increasingly smaller matrix problems standing at the basis of the AMG algorithm, so as to maximize the corresponding convergence factor of the AMG scheme. To demonstrate the practical capabilities of the proposed algorithm, which we call AMG-ANN, we consider the iterative solution of the algebraic system of equations stemming from finite element discretizations of two-dimensional model problems. First, we consider an elliptic equation with a highly heterogeneous diffusion coefficient and then a stationary Stokes problem. We train (off-line) our ANN with a rich dataset and present an in-depth analysis of the effects of tuning the strong threshold parameter on the convergence factor of the resulting AMG iterative scheme.

Secant Acceleration of Sequential Residual Methods for Solving Large-Scale Nonlinear Systems of Equations

Sequential Residual Methods try to solve nonlinear systems of equations $F(x)=0$ by iteratively updating the current approximate solution along a residual-related direction. Therefore, memory requirements are minimal and, consequently, these methods are attractive for solving large-scale nonlinear systems. However, the convergence of these algorithms may be slow in critical cases; therefore, acceleration procedures are welcome. In this paper, we suggest to employ a variation of the Sequential Secant Method in order to accelerate Sequential Residual Methods. The performance of the resulting algorithm is illustrated by applying it to the solution of very large problems coming from the discretization of partial differential equations.

BAMPHI: Matrix-free and transpose-free action of linear combinations of [formula omitted]-functions from exponential integrators

The time integration of stiff systems of Ordinary Differential Equations (ODEs), usually arising from the spatial discretization of Partial Differential Equations (PDEs), constitutes a hot topic in numerical analysis. In particular, exponential-type integrators have attracted much attention for their capacity to effectively handle the stiffness, allowing integration with large time steps. The efficiency of exponential-type integrators strongly relies on the fast computation of the action of the exponential of matrices which often change only slightly during the integration. The authors exploited this characteristic of exponential integrators to develop a backward stable algorithm, bamphi, which is designed to reuse the information gathered through the exponential integration steps, reaching unmatched levels of speed and accuracy on a variety of numerical experiments.

A Bootstrap Multigrid Eigensolver

This paper introduces bootstrap multigrid methods for solving eigenvalue problems arising from the discretization of partial differential equations. Inspired by the full bootstrap algebraic multigrid (BAMG) setup algorithm that includes an AMG eigensolver, it is illustrated how the algorithm can be simplified for the case of a discretized partial differential equation (PDE), thereby developing a bootstrap geometric multigrid (BMG) approach. We illustrate numerically the efficacy of the BMG method for: (1) recovering eigenvalues having large multiplicity, (2) computing interior eigenvalues, and (3) approximating shifted indefinite eigenvalue problems. Numerical experiments are presented to illustrate the basic components and ideas behind the success of the overall bootstrap multigrid approach. For completeness, we present a simplified error analysis of a two-grid bootstrap algorithm for the Laplace-Beltrami eigenvalue problem.

Acceleration of RBF-FD meshless phase-field modelling of dendritic solidification by space-time adaptive approach

• A novel meshless approach for phase-field modelling of dendritic growth. • Reduction of discretisation-induced anisotropy by the scattered nodes. • Increase in computational efficiency by a space-time adaptive approach. A novel adaptive numerical approach is developed for an accurate and computationally efficient phase-field modelling of dendritic solidification. The adaptivity is based on the dynamic quadtree domain decomposition. A quadtree decomposes the computational domain into rectangular sub-domains of different sizes. Each sub-domain is extended to ensure overlap communication between neighbouring sub-domains. In each sub-domain, uniform distribution of computational nodes is generated. The product between the node density and the sub-domain area is fixed to ensure the h-adaptivity. The adaptive approach provides the highest density of computational nodes at the solid-liquid interface and the lowest density in the bulk of the phases. The meshless radial basis function generated finite difference (RBF-FD) method is applied for the spatial discretisation of the partial differential equations which arise from the phase-field model. The RBF-FD method is especially appealing since it allows straightforward spatial discretisation of partial differential equations on scattered node distributions. The use of scattered node distribution reduces the discretisation-induced anisotropy in the phase-field modelling of dendritic growth. The forward Euler scheme is used for temporal discretisation. The adaptive time-stepping is employed to speed up the calculations further. The performance of the novel numerical approach is tested for dendritic solidification of supercooled pure melts and supersaturated dilute binary alloys at arbitrary preferential growth directions. The impact of the numerical parameters on the accuracy and computational efficiency is thoroughly analysed. It is shown that the RBF-FD method, defined on scattered node distribution, together with the space-time adaptive approach, represents an accurate and efficient technique for solving the phase-field models of dendritic solidification.

Adaptive Solution of Initial Value Problems by a Dynamical Galerkin Scheme

We study dynamical Galerkin schemes for evolutionary partial differential equations (PDEs), where the projection operator changes over time. When selecting a subset of basis functions, the projection operator is non-differentiable in time and an integral formulation has to be used. We analyze the projected equations with respect to existence and uniqueness of the solution and prove that non-smooth projection operators introduce dissipation, a result which is crucial for adaptive discretizations of PDEs, e.g., adaptive wavelet methods. For the Burgers equation we illustrate numerically that thresholding the wavelet coefficients, and thus changing the projection space, will indeed introduce dissipation of energy. We discuss consequences for the so-called `pseudo-adaptive' simulations, where time evolution and dealiasing are done in Fourier space, whilst thresholding is carried out in wavelet space. Numerical examples are given for the inviscid Burgers equation in 1D and the incompressible Euler equations in 2D and 3D.

Reducing errors caused by geometrical inaccuracy to solve partial differential equations with moving frames on curvilinear domain

In the numerical discretization of partial differential equations (PDEs) with moving frames on curved surfaces, the discretization error does not converge for a high p≥5. Moreover, the conservation error remains significant even in a refined mesh and does not converge as the polynomial order p increases. We postulate that the inaccurate location of the internal grid points of curved elements causes this problem; this is called the internal point error. This bottleneck of convergence persists even when the vertices of the elements are located exactly on the curved domain. We first theoretically explain the effects of the internal point error on numerical accuracy. We then computationally validate the postulation by using a high-order mesh with a negligible internal point error, even for a high-order polynomial of order p. For computational validations, four differential operators and four PDEs are numerically solved using moving frames on the unit sphere to demonstrate the effect of the internal point error.

F-actin bending facilitates net actomyosin contraction By inhibiting expansion with plus-end-located myosin motors

Contraction of actomyosin networks underpins important cellular processes including motility and division. The mechanical origin of actomyosin contraction is not fully-understood. We investigate whether contraction arises on the scale of individual filaments, without needing to invoke network-scale interactions. We derive discrete force-balance and continuum partial differential equations for two symmetric, semi-flexible actin filaments with an attached myosin motor. Assuming the system exists within a homogeneous background material, our method enables computation of the stress tensor, providing a measure of contractility. After deriving the model, we use a combination of asymptotic analysis and numerical solutions to show how F-actin bending facilitates contraction on the scale of two filaments. Rigid filaments exhibit polarity-reversal symmetry as the motor travels from the minus to plus-ends, such that contractile and expansive components cancel. Filament bending induces a geometric asymmetry that brings the filaments closer to parallel as a myosin motor approaches their plus-ends, decreasing the effective spring force opposing motor motion. The reduced spring force enables the motor to move faster close to filament plus-ends, which reduces expansive stress and gives rise to net contraction. Bending-induced geometric asymmetry provides both new understanding of actomyosin contraction mechanics, and a hypothesis that can be tested in experiments.

On the spectrum of the finite element approximation of a three field formulation for linear elasticity

We continue the investigation on the spectrum of operators arising from the discretization of partial differential equations. In this paper we consider a three field formulation recently introduced for the finite element least-squares approximation of linear elasticity. We discuss in particular the distribution of the discrete eigenvalues in the complex plane and how they approximate the positive real eigenvalues of the continuous problem. The dependence of the spectrum on the Lamé parameters is considered as well and its behavior when approaching the incompressible limit.

A Decision-Making Machine Learning Approach in Hermite Spectral Approximations of Partial Differential Equations

The accuracy and effectiveness of Hermite spectral methods for the numerical discretization of partial differential equations on unbounded domains are strongly affected by the amplitude of the Gaussian weight function employed to describe the approximation space. This is particularly true if the problem is under-resolved, i.e., there are no enough degrees of freedom. The issue becomes even more crucial when the equation under study is time-dependent, forcing in this way the choice of Hermite functions where the corresponding weight depends on time. In order to adapt dynamically the approximation space, it is here proposed an automatic decision-making process that relies on machine learning techniques, such as deep neural networks and support vector machines. The algorithm is numerically tested with success on a simple 1D problem, but the main goal is its exportability in the context of more serious applications. As a matter of fact we also show at the end an application in the framework of plasma physics.