2024-12-02 | | Total: 28
Traction force microscopy is a method widely used in biophysics and cell biology to determine forces that biological cells apply to their environment. In the experiment, the cells adhere to a soft elastic substrate, which is then deformed in response to cellular traction forces. The inverse problem consists in computing the traction stress applied by the cell from microscopy measurements of the substrate deformations. In this work, we consider a linear model, in which 3D forces are applied at a 2D interface, called 2.5D traction force microscopy, and a nonlinear pure 2D model, from which we directly obtain a linear pure 2D model. All models lead to a linear resp. nonlinear parameter identification problem for a boundary value problem of elasticity. We analyze the respective forward operators and conclude with some numerical experiments for simulated and experimental data.
The randomized Kaczmarz (RK) method is a well-known approach for solving linear least-squares problems with a large number of rows. RK accesses and processes just one row at a time, leading to exponentially fast convergence for consistent linear systems. However, RK fails to converge to the least-squares solution for inconsistent systems. This work presents a simple fix: average the RK iterates produced in the tail part of the algorithm. The proposed tail-averaged randomized Kaczmarz (TARK) converges for both consistent and inconsistent least-squares problems at a polynomial rate, which is known to be optimal for any row-access method. An extension of TARK also leads to efficient solutions for ridge-regularized least-squares problems.
The SE and DE formulas are known as efficient quadrature formulas for integrals with endpoint singularities. Especially for integrals with algebraic singularity, explicit error bounds in a computable form have been given, which are useful for computation with guaranteed accuracy. Such explicit error bounds have also given for integrals with logarithmic singularity. However, the error bounds have two points to be discussed. The first point is on overestimation of divergence speed of logarithmic singularity. The second point is on the case where there exist both logarithmic and algebraic singularity. To remedy these points, this study provides new error bounds for integrals with logarithmic and algebraic singularity. Although existing and new error bounds described above handle integrals over the finite interval, the SE and DE formulas may be applied to integrals over the semi-infinite interval. On the basis of the new results, this study provides new error bounds for integrals over the semi-infinite interval with logarithmic and algebraic singularity at the origin.
In this work we derive higher order error estimates for inverse problems distorted by non-additive noise, in terms of Bregman distances. The results are obtained by means of a novel source condition, inspired by the dual problem. Specifically, we focus on variational regularization having the Kullback-Leibler divergence as data-fidelity, and a convex penalty term. In this framework, we provide an interpretation of the new source condition, and present error estimates also when a variational formulation of the source condition is employed. We show that this approach can be extended to variational regularization that incorporates more general convex data fidelities.
We consider a mixed variational formulation recently proposed for the coupling of the Brinkman--Forchheimer and Darcy equations and develop the first reliable and efficient residual-based a posteriori error estimator for the 2D version of the associated conforming mixed finite element scheme. For the reliability analysis, due to the nonlinear nature of the problem, we make use of the inf-sup condition and the strong monotonicity of the operators involved, along with a stable Helmholtz decomposition in Hilbert spaces and local approximation properties of the Raviart--Thomas and Clément interpolants. On the other hand, inverse inequalities, the localization technique through bubble functions, and known results from previous works are the main tools yielding the efficiency estimate. Finally, several numerical examples confirming the theoretical properties of the estimator and illustrating the performance of the associated adaptive algorithms are reported. In particular, the case of flow through a heterogeneous porous medium is considered.
Physics-Informed Neural Networks (PINNs) are an emerging tool for approximating the solution of Partial Differential Equations (PDEs) in both forward and inverse problems. PINNs minimize a loss function which includes the PDE residual determined for a set of collocation points. Previous work has shown that the number and distribution of these collocation points have a significant influence on the accuracy of the PINN solution. Therefore, the effective placement of these collocation points is an active area of research. Specifically, adaptive collocation point sampling methods have been proposed, which have been reported to scale poorly to higher dimensions. In this work, we address this issue and present the Point Adaptive Collocation Method for Artificial Neural Networks (PACMANN). Inspired by classic optimization problems, this approach incrementally moves collocation points toward regions of higher residuals using gradient-based optimization algorithms guided by the gradient of the squared residual. We apply PACMANN for forward and inverse problems, and demonstrate that this method matches the performance of state-of-the-art methods in terms of the accuracy/efficiency tradeoff for the low-dimensional problems, while outperforming available approaches for high-dimensional problems; the best performance is observed for the Adam optimizer. Key features of the method include its low computational cost and simplicity of integration in existing physics-informed neural network pipelines.
In this paper we propose a novel traffic flow model based on understanding the city as a porous media, this is, streets and building-blocks characterizing the urban landscape are seen now as the fluid-phase and the solid-phase of a porous media, respectively. Moreover, based in the interchange of mass in the porous media models, we can model the interchange of cars between streets and off-street parking-spaces. Therefore, our model is not a standard conservation law, being formulated as the coupling of a non-stationary convection-diffusion-reaction PDE with a Darcy-Brinkman-Forchheimer PDE system. To solve this model, the classical Galerkin P1 finite element method combined with an explicit time marching scheme of strong stability-preserving type was enough to stabilize our numerical solutions. Numerical experiences on an urban-porous domain inspired by the city of Guadalajara (Mexico) allow us to simulate the influence of the porosity terms on the traffic speed, the traffic flow at rush-valley hours, and the streets congestions due to the lack of parking spaces.
This article presents updates to lifex [Africa, SoftwareX (2022)], a C++ library for high-performance finite element simulations of multiphysics, multiscale and multidomain problems. In this release, we introduce an additional intergrid transfer method for non-matching multiphysics coupling on the same domain, significantly optimize nearest-neighbor point searches and interface coupling utilities, extend the support for 2D and mixed-dimensional problems, and provide improved facilities for input/output and simulation serialization and restart. These advancements also propagate to the previously released modules of lifex specifically designed for cardiac modeling and simulation, namely lifex-fiber [Africa et al., BMC Bioinformatics (2023)], lifex-ep [Africa et al., BMC Bioinformatics (2023)] and lifex-cfd [Africa et al., Computer Physics Communications (2024)]. The changes introduced in this release aim at consolidating lifex's position as a valuable and versatile tool for the simulation of multiphysics systems.
In this paper, we consider an elliptic eigenvalue problem with multiscale, randomly perturbed coefficients. For an efficient and accurate approximation of the solutions for many different realizations of the coefficient, we propose a computational multiscale method in the spirit of the Localized Orthogonal Decomposition (LOD) method together with an offline-online strategy similar to [Målqvist, Verfürth, ESIAM Math. Model. Numer. Anal., 56(1):237-260, 2022]. The offline phase computes and stores local contributions to the LOD stiffness matrix for selected defect configurations. Given any perturbed coefficient, the online phase combines the pre-computed quantities in an efficient manner. We further propose a modification in the online phase, for which numerical results indicate enhanced performances for moderate and high defect probabilities. We show rigorous a priori error estimates for eigenfunctions as well as eigenvalues.
We present and analyze a discontinuous Galerkin method for the numerical modeling of a Kelvin-Voigt thermo/poro-viscoelastic problem. We present the derivation of the model, and we develop a stability analysis in the continuous setting that holds both for the full inertial and quasi-static problems and that is robust with respect to most of the physical parameters of the problem. For spatial discretization, we propose an arbitrary-order weighted symmetric interior penalty scheme that supports general polytopal grids and is robust with respect to strong heterogeneities in the model coefficients. For the semi-discrete problem, we prove the extension of the stability result demonstrated in the continuous setting. A wide set of numerical simulations is presented to assess the convergence and robustness properties of the proposed method. Moreover, we test the scheme with literature and physically sound test cases for proof-of-concept applications in the geophysical context.
A centrality measure of the cut-edges of an undirected graph, given in [Altafini et al.~SIMAX 2023] and based on Kemeny's constant, is revisited. A numerically more stable expression is given to compute this measure, and an explicit expression is provided for some classes of graphs, including one-path graphs and trees formed by three or more branches. These results theoretically confirm the good physical behaviour of this centrality measure, experimentally observed in [Altafini et al.~SIMAX 2023]. Numerical tests are reported to check the stability and to confirm the good physical behaviour.
We develop efficient and effective strategies for the update of Katz centralities after node and edge removal in simple graphs. We provide explicit formulas for the ``loss of walks" a network suffers when nodes/edges are removed, and use these to inform our algorithms. The theory builds on the newly introduced concept of $\cF$-avoiding first-passage walks. Further, bounds on the change of total network communicability are also derived. Extensive numerical experiments on synthetic and real-world networks complement our theoretical results.
Stiff systems of ordinary differential equations (ODEs) arise in a wide range of scientific and engineering disciplines and are traditionally solved using implicit integration methods due to their stability and efficiency. However, these methods are computationally expensive, particularly for applications requiring repeated integration, such as parameter estimation, Bayesian inference, neural ODEs, physics-informed neural networks, and MeshGraphNets. Explicit exponential integration methods have been proposed as a potential alternative, leveraging the matrix exponential to address stiffness without requiring nonlinear solvers. This study evaluates several state-of-the-art explicit single-step exponential schemes against classical implicit methods on benchmark stiff ODE problems, analyzing their accuracy, stability, and scalability with step size. Despite their initial appeal, our results reveal that explicit exponential methods significantly lag behind implicit schemes in accuracy and scalability for stiff ODEs. The backward Euler method consistently outperformed higher-order exponential methods in accuracy at small step sizes, with none surpassing the accuracy of the first-order integrating factor Euler method. Exponential methods fail to improve upon first-order accuracy, revealing the integrating factor Euler method as the only reliable choice for repeated, inexpensive integration in applications such as neural ODEs and parameter estimation. This study exposes the limitations of explicit exponential methods and calls for the development of improved algorithms.
In this paper, in order to improve the spatial accuracy, the exponential integrator Fourier Galerkin method (EIFG) is proposed for solving semilinear parabolic equations in rectangular domains. In this proposed method, the spatial discretization is first carried out by the Fourier-based Galerkin approximation, and then the time integration of the resulting semi-discrete system is approximated by the explicit exponential Runge-Kutta approach, which leads to the fully-discrete numerical solution. With certain regularity assumptions on the model problem, error estimate measured in $H^2$-norm is explicitly derived for EIFG method with two RK stages. Several two and three dimensional examples are shown to demonstrate the excellent performance of EIFG method, which are coincident to the theoretical results.
Krylov subspace methods are a powerful tool for efficiently solving high-dimensional linear algebra problems. In this work, we study the approximation quality that a Krylov subspace provides for estimating the numerical range of a matrix. In contrast to prior results, which often depend on the gaps between eigenvalues, our estimates depend only on the dimensions of the matrix and Krylov subspace, and the conditioning of the eigenbasis of the matrix. In addition, we provide nearly matching lower bounds for our estimates, illustrating the tightness of our arguments.
We present a simple universal algorithm for high-dimensional integration which has the optimal error rate (independent of the dimension) in all weighted Korobov classes both in the randomized and the deterministic setting. Our theoretical findings are complemented by numerical tests.
In this paper, we develop monolithic limiting techniques for enforcing nonlinear stability constraints in enriched Galerkin (EG) discretizations of nonlinear scalar hyperbolic equations. To achieve local mass conservation and gain control over the cell averages, the space of continuous (multi-)linear finite element approximations is enriched with piecewise-constant functions. The resulting spatial semi-discretization has the structure of a variational multiscale method. For linear advection equations, it is inherently stable but generally not bound preserving. To satisfy discrete maximum principles and ensure entropy stability in the nonlinear case, we use limiters adapted to the structure of our locally conservative EG method. The cell averages are constrained using a flux limiter, while the nodal values of the continuous component are constrained using a clip-and-scale limiting strategy for antidiffusive element contributions. The design and analysis of our new algorithms build on recent advances in the fields of convex limiting and algebraic entropy fixes for finite element methods. In addition to proving the claimed properties of the proposed approach, we conduct numerical studies for two-dimensional nonlinear hyperbolic problems. The numerical results demonstrate the ability of our limiters to prevent violations of the imposed constraints, while preserving the optimal order of accuracy in experiments with smooth solutions.
We consider a setting in which an evolving surface is implicitly characterized as the zero level of a level set function. Such an implicit surface does not encode any information about the path of a single point on the evolving surface. In the literature different approaches for determining a velocity that induces corresponding paths of points on the surface have been proposed. One of these is based on minimization of the strain energy functional. This then leads to a constrained minimization problem, which has a corresponding equivalent formulation as a saddle point problem. The main topic of this paper is a detailed analysis of this saddle point problem and of a finite element discretization of this problem. We derive well-posedness results for the continuous and discrete problems and optimal error estimates for a finite element discretization that uses standard $H^1$-conforming finite element spaces.
We propose a quadrature-based formula for computing the exponential function of matrices with a non-oscillatory integral on an infinite interval and an oscillatory integral on a finite interval. In the literature, existing quadrature-based formulas are based on the inverse Laplace transform or the Fourier transform. We show these expressions are essentially equivalent in terms of complex integrals and choose the former as a starting point to reduce computational cost. By choosing a simple integral path, we derive an integral expression mentioned above. Then, we can easily apply the double-exponential formula and the Gauss-Legendre formula, which have rigorous error bounds. As numerical experiments show, the proposed formula outperforms the existing formulas when the imaginary parts of the eigenvalues of matrices have large absolute values.
We show that, under certain circumstances, it is possible to automatically compute Jacobian-inverse-vector and Jacobian-inverse-transpose-vector products about as efficiently as Jacobian-vector and Jacobian-transpose-vector products. The key insight is to notice that the Jacobian corresponding to the use of one basis function is of a form whose sparsity is invariant to inversion. The main restriction of the method is a constraint on the number of active variables, which suggests a variety of techniques or generalization to allow the constraint to be enforced or relaxed. This technique has the potential to allow the efficient direct calculation of Newton steps as well as other numeric calculations of interest.
The eigenvalues and eigenvectors of nonnormal matrices can be unstable under perturbations of their entries. This renders an obstacle to the analysis of numerical algorithms for non-Hermitian eigenvalue problems. A recent technique to handle this issue is pseudospectral shattering [BGVKS23], showing that adding a random perturbation to any matrix has a regularizing effect on the stability of the eigenvalues and eigenvectors. Prior work has analyzed the regularizing effect of dense Gaussian perturbations, where independent noise is added to every entry of a given matrix [BVKS20, BGVKS23, BKMS21, JSS21]. We show that the same effect can be achieved by adding a sparse random perturbation. In particular, we show that given any $n\times n$ matrix $M$ of polynomially bounded norm: (a) perturbing $O(n\log^2(n))$ random entries of $M$ by adding i.i.d. complex Gaussians yields $\log\kappa_V(A)=O(\text{poly}\log(n))$ and $\log (1/\eta(A))=O(\text{poly}\log(n))$ with high probability; (b) perturbing $O(n^{1+\alpha})$ random entries of $M$ for any constant $\alpha>0$ yields $\log\kappa_V(A)=O_\alpha(\log(n))$ and $\log(1/\eta(A))=O_\alpha(\log(n))$ with high probability. Here, $\kappa_V(A)$ denotes the condition number of the eigenvectors of the perturbed matrix $A$ and $\eta(A)$ denotes its minimum eigenvalue gap. A key mechanism of the proof is to reduce the study of $\kappa_V(A)$ to control of the pseudospectral area and minimum eigenvalue gap of $A$, which are further reduced to estimates on the least two singular values of shifts of $A$. We obtain the required least singular value estimates via a streamlining of an argument of Tao and Vu [TV07] specialized to the case of sparse complex Gaussian perturbations.
Mutual Information (MI) is a powerful statistical measure that quantifies shared information between random variables, particularly valuable in high-dimensional data analysis across fields like genomics, natural language processing, and network science. However, computing MI becomes computationally prohibitive for large datasets where it is typically required a pairwise computational approach where each column is compared to others. This work introduces a matrix-based algorithm that accelerates MI computation by leveraging vectorized operations and optimized matrix calculations. By transforming traditional pairwise computational approaches into bulk matrix operations, the proposed method enables efficient MI calculation across all variable pairs. Experimental results demonstrate significant performance improvements, with computation times reduced up to 50,000 times in the largest dataset using optimized implementations, particularly when utilizing hardware optimized frameworks. The approach promises to expand MI's applicability in data-driven research by overcoming previous computational limitations.
We introduce a novel method for solving density-based topology optimization problems: \underline{Si}gmoidal \underline{M}irror descent with a \underline{P}rojected \underline{L}atent variable (SiMPL). The SiMPL method (pronounced as "the simple method") optimizes a design using only first-order derivative information of the objective function. The bound constraints on the density field are enforced with the help of the (negative) Fermi--Dirac entropy, which is also used to define a non-symmetric distance function called a Bregman divergence on the set of admissible designs. This Bregman divergence leads to a simple update rule that is further simplified with the help of a so-called latent variable. %Introducing a generalized Barzilai-Borwein step size rule accelerates the convergence of SiMPL. Because the SiMPL method involves discretizing the latent variable, it produces a sequence of pointwise-feasible iterates, even when high-order finite elements are used in the discretization. Numerical experiments demonstrate that the method outperforms other popular first-order optimization algorithms. To outline the general applicability of the technique, we include examples with (self-load) compliance minimization and compliant mechanism optimization problems.
In this work we present a computationally efficient linear optimization approach for estimating the cross--power spectrum of an hidden multivariate stochastic process from that of another observed process. Sparsity in the resulting estimator of the cross--power is induced through $\ell_1$ regularization and the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) is used for computing such an estimator. With respect to a standard implementation, we prove that a proper initialization step is sufficient to guarantee the required symmetric and antisymmetric properties of the involved quantities. Further, we show how structural properties of the forward operator can be exploited within the FISTA update in order to make our approach adequate also for large--scale problems such as those arising in context of brain functional connectivity. The effectiveness of the proposed approach is shown in a practical scenario where we aim at quantifying the statistical relationships between brain regions in the context of non-invasive electromagnetic field recordings. Our results show that our method provide results with an higher specificity that classical approaches based on a two--step procedure where first the hidden process describing the brain activity is estimated through a linear optimization step and then the cortical cross--power spectrum is computed from the estimated time--series.
We present a novel formulation for the dynamics of geometrically exact Timoshenko beams and beam structures made of viscoelastic material featuring complex, arbitrarily curved initial geometries. An $\textrm{SO}(3)$-consistent and second-order accurate time integration scheme for accelerations, velocities and rate-dependent viscoelastic strain measures is adopted. To achieve high efficiency and geometrical flexibility, the spatial discretization is carried out with the isogemetric collocation (IGA-C) method, which permits bypassing elements integration keeping all the advantages of the isogeometric analysis (IGA) in terms of high-order space accuracy and geometry representation. Moreover, a primal formulation guarantees the minimal kinematic unknowns. The generalized Maxwell model is deployed directly to the one-dimensional beam strain and stress measures. This allows to express the internal variables in terms of the same kinematic unknowns, as for the case of linear elastic rate-independent materials bypassing the complexities introduced by the viscoelastic material. As a result, existing $\textrm{SO}(3)$-consistent linearizations of the governing equations in the strong form (and associated updating formulas) can straightforwardly be used. Through a series of numerical tests, the attributes and potentialities of the proposed formulation are demonstrated. In particular, we show the capability to accurately simulate beams and beam systems featuring complex initial geometry and topology, opening interesting perspectives in the inverse design of programmable mechanical meta-materials and objects.