The talk aims at a generalization of the known Gauss-Weierstrass integral operators introduced two decades ago by Eugeniusz Wachnicki. It is intimately connected by a generalization of the heat equation.

Our new results concern two aspects: obtaining an asymptotic expansion for this generalization transform when applied to a certain class of signals and getting the asymptotic expansion of derivatives of any order of Wachnicki’s operators. All coefficients are explicitly calculated and distinct expressions are provided for analytical functions.

Many engineering control and optimization problems occupy an extensive area of space and time. These include production cost models in energy or the control of central plants. Model predictive control, one of the favorite approaches for physics-centered control, generally relies on direct numerical solvers, which tend to run out of memory for such problems. Approaches relying on domain decomposition are then key to fit in memory, but theory for such second-order methods tends to be lacking in multicomponent systems.

To address this issue, we prove that certain classes of graph-indexed optimization (GIO) problems exhibit exponential decay of sensitivity with respect to perturbation in the data. GIOs include dynamic optimization (for the linear graph), or distributed, including network control (for the mesh/network-time product graph). This feature allows for very efficient approximation, solutions, or policies based on domain decomposition with overlap relative to centralized or monolithic approaches. In particular, we prove that the proper efficiency metric increases exponentially fast with the overlap size. Immediate consequences of such behavior are that distributed control policies with overlap approach the performance of centralized policies exponentially fast and that Schwarz-type algorithms exhibit, in addition to exceptional parallelism and reduced memory footprint per subproblem, a linear rate of convergence that tends exponentially fast to zero.

We consider  a model convection-diffusion problem  and  present our recent numerical and analysis results regarding mixed finite element formulation and discretization in the singular perturbed case when the convection term dominates the problem. Using  the concepts of optimal norm and saddle point reformulation, we found new error  estimates  for the case of uniform meshes.  We compare  the standard linear Galerkin discretization to a saddle point least square discretization that uses quadratic test functions, and explain the non-physical oscillations of the discrete solutions. We also relate  a known upwinding Petrov-Galerkin method and  the stream-line diffusion discretization method, by emphasizing the resulting  linear  systems and by comparing appropriate error norms.  The results can be extended to the multidimensional case in order to find efficient approximations for more general singular perturbed problems including convection dominated models.

Stiff ordinary differential equations (ODEs) are challenging to solve numerically, but artificial intelligence (AI) offers a promising new approach. Physics-informed neural networks (PINNs) can be trained to learn the dynamics of a system from both data and physical laws. PINNs have been shown to be effective at solving stiff ODEs with high accuracy and efficiency. This presentation will discuss the potential of AI to aid in solving stiff ODEs, introduce PINNs, and exemplify their application to a stiff ODE.

We analyze the strict superlinear (SSL) convergence order, which is clearly much faster than linear, but not necessarily slower than, say, quadratic.

We show that the set of the SSL sequences can be partitioned in four distinct classes: weak, medium, strong (in increasing speed) and mixed (which cannot be assessed).

The sequences from the medium and weak classes are slower than the sequences with high classical C-orders $p>1$, while an example shows that certain sequences from the mixed class may be term-by-term faster than some sequences with infinite C-order.

One can evaluate numerically to which class a SSL sequence belongs.

Some recent results of Rodomanov and Nesterov (2022), resp. Ye et al. (2023) show that certain classical quasi-Newton methods (DFP, BFGS and SR1) belong to the weak class.

In this talk we will present the fitting of the theoretical predictions of the rotation curves in conformal gravity with data for several galaxies included in the SPARC database by using the Multi Start and Global Search methods. In the total expression of the tangential velocity we also include the effects of the baryonic matter, and  the mass to luminosity ratio. The preliminary findings suggest that the solution offered by Weyl geometric gravity effectively explains a significant portion of the rotation curves within the SPARC sample. (joint work with Tiberiu Harko)

The aim of this research is to develop efficient multigrid preconditioners for a classic linear-quadratic optimization problem constrained by an elliptic equation with stochastic coefficients. Using a discretize-then-optimize approach, in previous work we have shown that strategies inherited from the associated deterministic optimal control problem extend to the stochastic version when a stochastic Galerkin discretization is employed. In this talk we show that similar strategies succeed when discretizing the elliptic equation using a sparse grid stochastic collocation approach.
Authors: Sumaya Alzuhairy, Andrei Draganescu and Bedrich Sousedik.

The most important example of the pseudo parabolic problem we deal with is a Cauchy one attached to the Camassa-Holm equation formulated on an unbounded domain. To remove the inconveniences related to the domain not being bounded, we will use the collocation based on sinc functions. These make it possible to introduce the behavior of the solution to infinity. To simulate the evolution in time, various schemes with finite differences are used.
The treated problem has the property of being bi-Hamiltonian and the methods used failed to satisfactorily conserve these invariants.  We will focus our future efforts in this direction.
The solutions already obtained are presented and commented on the accuracy.

Over the past decade, several stabilization strategies have been developed to tackle the reduced order model (ROM) inaccuracy in the convection-dominated, under-resolved regime, i.e., when the number of degrees of freedom is too small to capture the complex underlying dynamics.  In this talk, I will survey regularized reduced order models (Reg-ROMs), which are simple, modular stabilizations that employ ROM spatial filtering of various terms in the Navier-Stokes equations (NSE) to alleviate the spurious numerical oscillations generally produced by standard ROMs in the convection-dominated, under-resolved regime.  I will focus on two different types of Reg-ROM strategies: (i) the evolve-filter-relax ROM (EFR-ROM), which first filters an intermediate velocity approximation, and then relaxes it; and (ii) the Leray-ROM (L-ROM), which filters the convective term in the NSE.  Throughout my talk, I will highlight the impact made by ROM spatial filtering on the Reg-ROM development.  Specifically, I will talk about the two main types of ROM spatial filters: (i) the ROM differential filter; and (ii) the ROM projection.  I will also propose two novel higher-order ROM differential filters.  An important role played in ROM spatial filters and Reg-ROMs is the ROM lengthscale.

In my talk, I will put forth a novel ROM lengthscale, which is constructed by leveraging energy balancing arguments.  I emphasize that this novel energy-based lengthscale is fundamentally different from the standard ROM lengthscale introduced decades ago, which is based on simple dimensional arguments.  Finally, I will illustrate the success achieved by ROM spatial filters and Reg-ROMs in under-resolved numerical simulations of the turbulent channel flow.

Throughout my talk, I will discuss numerical analysis results proved for the Reg-ROMs that we propose, including fundamental properties, e.g., stability, convergence, and parameter scalings. I will also present some of the challenges and open questions in the development of rigorous numerical analysis foundations for ROM closures and stabilizations.

Surface reconstruction from scattered data is a 2D interpolation problem. There are many theoretical results and more many numerical algorithms related to the construction of a surface from sampled points. If the measurements are affected by random errors, the constructed surface exhibits oscillations and the numerical algorithms must cure the errors to produce a smooth surface. If one is only interested in the shape of the surface, then the penalized least squares method generally furnishes good results. The problem becomes more difficult if one wishes to generate surfaces with smooth curvatures. In this talk we will analyze a class of interpolation methods able to solve such problem.

Joint work with Doerin Marinescu and Stefan Cruceanu.

In computer-aided mathematical proofs, a basic, yet critical, building block is the problem of actually obtaining numerical values. In practice, one strives to achieve precise and/or guaranteed results without compromising efficiency. We briefly introduce basic tools like interval arithmetic and then focus on effectively computing polynomial approximations (expressed as truncated Taylor, Chebyshev or trigonometric series) together with validated error bounds. In particular, we discuss about a posteriori validation with Newton-like operators. Finally, all these techniques are applied to the efficient finite precision evaluation of numerical functions. We provide examples ranging from robust space mission design to computer-aided proofs in dynamical systems.

Stochastic partial differential equations are widely used to model the evolution of uncertain dynamical systems in geophysical fluid dynamics. For a judicious modelling of the evolution of a fluid flow, the noise term needs to be properly calibrated. Lagrangian methods have been developed in this sense, where particle trajectories are simulated starting from each point on both the physical grid and its refined version, then the differences between the particle positions are used to calibrate the noise. This is computationally expensive and not fully justified from a theoretical perspective. We propose an Eulerian alternative which can be applied to a general class of stochastic models. In Part I we will introduce the theoretical justification for this new calibration method, while in Part II we will present the calibration results for a stochastic rotating shallow water model.

This is joint work with Alexander Lobbe (Imperial College London) and the results are being published in “Noise calibration for SPDEs: a case study for the rotating shallow water model” (Foundations of Data Science), arXiv:2305.03548.

The main aim of this talk is to present a method for reconstructing damaged black and white or color images, considering both a global and a local approach. The method is based on the bivariate Shepard operator combined with the inverse quadratic and the inverse multiquadric radial basis functions. We evaluate the results by studying the approximation errors.

We investigate, from both the theoretical and the numerical viewpoints, the acceleration of the iterative algorithms of Kozlov et al. (Comput. Math. Math. Phys. 31 45-52, 1991) for the accurate, convergent and stable reconstruction of the missing data (both the temperature and the normal heat flux) on an inaccessible boundary of the domain occupied by a solid from the knowledge of Cauchy data on the remaining and accessible boundary in the framework of stationary anisotropic heat conduction without heat sources. For each of the two algorithms with relaxation considered, this inverse Cauchy problem in anisotropic heat conduction with exact data is transformed into an equivalent fixed-point problem for an associated operator that is defined on and takes values in a suitable function space, and accounts for the relaxation parameter. Consequently, the convergence of each relaxation algorithm reduces to investigating the properties of the corresponding operator. This enables one to determine the admissible range for the relaxation parameter along with a criterion for selecting its optimal value at each iteration for each iterative algorithm and exact Cauchy data. The numerical implementation is realised for 2D and 3D homogeneous anisotropic solids via the meshless method of fundamental solutions and confirms a significant reduction in the number of iterations and hence the CPU time required for the two relaxation algorithms proposed to achieve convergence, provided that the dynamical selection of the optimal value for the relaxation parameter is employed.

Biological tissue exhibits significant variability in mechanical properties that is usually attributed to random microscopic realizations. A computational approach to nonlinear mechanics of tissue with an underlying stochastic component is introduced. At the bulk continuum level, a hyperelastic formulation is combined with linear viscous models to provide the mean stress and strain fields for microscale constituent elements of the tissue. The hyperelastic behavior is expressed as a system of conservation laws, $\boldsymbol{q}_{, t} + \nabla_{\boldsymbol{a}} \boldsymbol{f}(\boldsymbol{q}) =\boldsymbol{\psi} (\boldsymbol{q}, \boldsymbol{\xi})$, for the velocity and deformation gradient $\boldsymbol{q}= (\boldsymbol{v}, \boldsymbol{F})$ in a Lagrangian framework. The flux $\boldsymbol{f}= (-\boldsymbol{P}(\boldsymbol{F}),-\boldsymbol{Z} (\boldsymbol{v}))$ requires definition of the dependence of the first Piola-Kirhhoff tensor $\boldsymbol{P}$ in terms of the deformation gradient $\boldsymbol{F}$. Viscous effects are captured through the sink term $\boldsymbol{\psi}$, that incorporates history dependence through the auxilliary variables $\boldsymbol{\xi}$.

Constitutive laws of known analytical form utilize model parameters $\boldsymbol{y} (\boldsymbol{\omega})$ that are assumed here to be random functions of the microstate $\omega$, and can be determined by carrying out integration over all possible microstates $\boldsymbol{P} (\boldsymbol{F}) = \int_{\Omega} \boldsymbol{G} (\boldsymbol{F}, \omega) p (\omega) d \omega$, where $\boldsymbol{G}$ is the generating function for $\boldsymbol{P}$. The probability density function (PDF) $p (\omega)$ is extracted from microscopic simulations carried out under known mean fields; Brownian Dynamics and Markov Chain Monte Carlo methods are considered. Probability distribution functions characterizing the microscale configuration are organized as a statistical manifold $\mathcal{M}= \{ p (\omega ; \lambda)\}$, where $\lambda$ is a parametrization of the manifold. Singular value decomposition of data from successive microscopic states defines the localcurvature of $\mathcal{M}$, and computational geometry provides a local reconstruction of the manifold. Low-parameter PDFs are obtained by geodesic transport on $\mathcal{M}$ applying concepts of information geometry. The construction is updated during time evolution, forming a macro-microscale interaction loop. The procedure is applied to propagation of shear waves in the brain, using realistic data obtained from histology and compared to ultrasonic shear wave imaging.

Pole placement is a well-known and classical method for controlling finite-dimensional linear time-invariant (LTI) systems. Roughly speaking, under appropriate controllability conditions, the method consists in assigning poles of the closed-loop system to specified locations by an appropriate choice of the controller gains guaranteeing the exponential stability in closed-loop with a prescribed decay rate of the corresponding system’s solution. The construction makes use of the characteristic polynomial degree and the structural properties (controllability) of the system. Its extension for infinite-dimensional linear systems is far from being trivial. In the case of dynamical systems governed by retarded and/or neutral delay-differential equations (DDEs), there exist only two effective extensions – continuous pole placement and partial pole placement.

This talk provides some insights into the so-called partial pole placement method. By exploiting the notion of degree of a quasipolynomial, and by assigning a real root with a maximal admissible multiplicity, we try to answer to a simple question: Does such an assignment (of real roots) guarantee the closed-loop stability of the closed-loop system? In other words, we are interested to explore the situation when the assigned (multiple real) root «govern » the location of the remaining (infinitely many?!) roots of the characteristic function. More precisely, in some cases, we can prove that such a multiple root is « dominant » and explicitly defines the spectral abscissa of the dynamical system. Such a property, called « multiplicity-induced-dominancy » (MID), has no a « natural » (non-trivial) equivalent in finite-dimension and it is specific to DDEs. To perform such an analysis we exploit the analytical properties of Kummer/Whittaker confluent hypergeometric functions. Finally, such ideas open interesting perspectives in interpreting Padé approximations.

(joint work with Islam BOUSSAADA & Guilherme MAZANTI)

Incremental learning refers to a class of machine learning models in which the learning process goes through a continuous model adaptation based on a constantly arriving data stream. This type of learning is needed when access to past data is limited or impossible, but is affected by catastrophic forgetting. This phenomenon consists in a drastic performance drop for previously learned information when ingesting new data. We tackle class-incremental learning without memory by adapting prediction bias correction, a method which makes predictions of past and new classes more comparable. In the same context, we also propose a method to control the dark knowledge values also known as soft targets, with the purpose of improving the training by transfer of knowledge.

In this paper, we prove that, if $\inf \limits_{x\in A}|f(x)|=m>0$, then the
partial differential operator $D$ defined by $D(u)= {\textstyle \sum \limits_{k=1}^{n}} f_{k}\frac{\partial u}{\partial x_{k}}-fu$, where $f,f_{i}\in C(A,\mathbb{R}),u\in C^{1}(A,X),\ i=1,\ldots,n,I\subset \mathbb{R}$ is an interval, $A=I\times \mathbb{R}^{n-1}$ and $X$ is a Banach space, is Ulam stable with the Ulam constant $K=\frac{1}{m}$. Moreover, if $\inf \limits_{x\in A}|f(x)|=0$, we prove that $D$ is not generally Ulam stable.

Joint work with Adela Novac and Dorian Popa.

References

[1] S. M. Jung, Hyers-Ulam stability of linear partial differential equations of first order. Appl. Math. Lett. 2009, 22, 70-74.

[2] N. Lungu, D. Popa, Hyers-Ulam stability of a first order partial differential equation. J. Math. Anal. Appl.  2012, 385, 86–91.

[3] D. Popa, I. Rasa, Hyers-Ulam stability of the Laplace operator. Fixed Point Theory, 2018, 19(1), 379-382.

Model-based projections of complex systems will play a central role in prediction and decision-making, e.g., anticipate ice sheet contribution to sea level rise. However, models are typically subject to uncertainties stemming from uncertain inputs in the model such as coefficient fields, source terms, etc. Refining our understanding of such models requires solving a (Bayesian) inverse problem that integrates measurement data and the model to estimate and quantify the uncertainties in the parameters. Such problems face a number of challenges, e.g., high-dimensionality of the inversion parameters, expensive-to-evaluate models, and model uncertainty additional to the uncertainty in inversion parameters. In this talk, I will focus on identifying and exploiting low-dimensional structure in such inverse problems stemming for example from local sensitivity of the data with respect to parameters, diffusive models, sparse data, etc. The key idea is to look at the Hessian operator arising in these inverse problems, which often exhibits a low rank and even a lower off-diagonal structure, and build approximations of this operator. This approximation can then be used as a preconditioner for the Newton-Krylov systems arising when solving the inverse problem, or to build more-informative proposals in an MCMC context. Combining an efficient matrix-free point spread function (PSF) method for approximating operators with fast hierarchical (H) matrix methods, we show that Hessian information can be obtained using only a small number of operator applications.

Co-authors: Noemi Petra (UC Merced), Nick Alger (UT Austin), Tucker Hartland (LLNL) and Omar Ghattas (UT Austin)

In this talk I will present a geometric multilevel optimization approach choosing as case study a regularised inverse problem. In particular, the approach is motivated by variational models that arise as the discretization of some underlying infinite dimensional problem. Such problems naturally lead to a hierarchy of discretized models. We employ multilevel optimization to take advantage of this hierarchy: while working at the fine level we compute the search direction based on a coarse model. By utilising concepts of information geometry in our formulation, we propose a smoothing operator that only uses first-order information and incorporates constraints smoothly. We show that the proposed algorithm is well suited for ill-posed reconstruction problems and demonstrate its efficiency on several large-scale examples.

We consider a nonlinear, possibly degenerated parabolic equation, modelling e.g. unsaturated flow in a porous medium, or biofilm growth. Implicit time stepping methods are popular due to their stability, allowing to avoid severe restrictions on the time step. This leads to nonlinear, time-discrete elliptic equations, for which linear iterative schemes are needed for approximating the solution. Due to the degeneracy, the Newton scheme may fail to converge, unless very small time steps are used, possibly after applying a regularisation step.

In this talk, we present an iterative scheme that does not require any regularization, for which the linear convergence can be proved rigorously under a mild restriction (if any) on the time step. This convergence is obtained at the level of the elliptic problem, so it is not restricted to any spatial discretisation or mesh. Moreover, it features an improved the convergence behavior, being linearly convergent in cases when the Newton diverges, and requiring at most a comparable number of iterations whenever the Newton scheme does converge.

Modern supercomputers vary in compute power and network capabilities. For example, Summit and Frontier make use of GPUs for accelerating computation, while relying on either a fat-tree or dragonfly topology for transferring data between the nodes. On the other hand, Fugaku has thousands of CPUs and uses a six-dimensional torus for communication. In this work, we want to study these differences in the context of scaling eigenvalue solvers used in planewave DFT calculations. More specifically, we will focus on the components that appear in nonlinear eigenvalue solvers like Conjugate Gradient, Jacobi Davidson and Unconstrained. We will look at the distributed Fourier computation and some of the linear algebra operations. We will provide an in-depth analysis and offer solutions for scaling the computation to a large number of GPUs. We will emphasize that systematic approaches can be derived to guide the parallelization such that the computation can effectively use the compute and network resources.

We discuss the localization of velocity for a problem of the type
$$\begin{equation} \left\{ \begin{array}{l} -\text{div }\left( A\left( x\right) \nabla u\right) +\eta _{0}\left( x\right) u+\kappa _{0}\left( x\right) \left( u\cdot \nabla \right) u+\nabla p=\Phi \left( x,u\right) \ \ \ \text{in }\Omega \\ \text{div }u=0\ \ \ \text{in }\Omega \\ u=0\ \ \text{on }\partial \Omega , \end{array}% \right. \end{equation}$$
where $\Phi$ is a reaction term dependent on velocity. First we obtain the
localization of the enstrophy, namely $\ r\leq \left\vert u\right\vert _{H_{0}^{1}\left( \Omega \right) }\leq R,\ $and then, the localization of the kinetic energy, that is $\ r\leq \left\vert u\right\vert _{L^{2}\left( \Omega \right) }\leq R.$ The bounds $r$ and $R$ are estimated in terms of the reaction force $\Phi$ and of system coefficients. The proofs are based on the fixed point formulation of the problem and on the fixed point index.
The results come from a joint work with Mirela Kohr, in continuation of the paper: M. Kohr and R. Precup. Analysis of Navier-Stokes models for flows in bidisperse porous media. J. Math. Fluid Mech. (2023)25:38.

Simulation of complex systems, composed of multiple subsystems, is often done by resolving individual susystems, and approximating their interaction with the other. This approach is called co-simulation. In this talk we will review popular co-simulation techniques, and latest developments in the field.

We deal with a class of elliptic quasivariational  inequalities with constraints in a reflexive Banach space. We use arguments of monotonicity, convexity and compactness in order to prove a convergence criterion for such inequalities. This criterion allows us to consider a new well-posedness concept in the study of the corresponding inequalities, which extends the classical Tykhonov and Levitin-Polyak well-posedness concepts used in the literature. Then, we introduce a new penalty method, for which we provide a convergence result. Finally, we consider a variational inequality which describes the equilibrium of a spring-rods system, under the action of external forces.
We apply our abstract results in the study of this inequality and provide the corresponding mechanical interpretations. We also present numerical simulations which  validate our convergence results.

Using fixed points arguments together with critical point techniques we obtain a hybrid existence result for a system with multiple operator equations, where only a part of them admit a variational structure.
The part of the solution of the equations which admit a variational structure are a Nash-type equilibrium for the corresponding energy functionals, and moreover they are localized in a bounded convex conical set defined by a norm and a concave upper semicontinuous functional.
This is achieved by an hybrid iterative scheme, partially based on Minty-Browder theorem and partially based on Ekeland variational principle. An abstract example of a system of differential equations with Dirichlet boundary conditions is given.

The Richards’ equation is modelling saturated/unsaturated flow porous media. It is a highly nonlinear, degenerate elliptic-parabolic equation, which is known to be very challenging to solve. In this talk we present a new solver for Richards´equation based on an adaptive switching between Newton and L-scheme. The theoretical properties of the scheme will be discussed. Illustrative numerical examples will be shown.

We prove that discrete systems of particles with piece-wise analytical trajectories can be described macroscopically by almost everywhere continuous fields. The transition from the microscopic to the macroscopic description is achieved through space-time averages on d-dimensional cubes and symmetric time intervals. The fields obtained in this way verify non-closed balance equations that are further used to verify the averaging procedure by estimating the intrinsic diffusion coefficients. The approach is applied to space-time upscaling of the reactive transport in variably saturated porous media. Processes at the microscopic level are modeled by as many random walkers as molecules involved in chemical reactions. The space-time scales of the averaging procedure may be chosen as representative for the macroscopic observations and the experimental measurements.

I will present a recent reconstruction of the full planar vector fields from its Doppler and its first moment Doppler transform data in fan beam coordinates. The method of proof is constructive  and is based on the theory of A-analytic maps.

We propose and analyze a second-order partitioned time-stepping method for a two phase flow problem in porous media. The algorithm is based on a refactorization of Cauchy’s one-legged θ-method.

The main part consists of the implicit backward Euler method, while part two uses a linear extrapolation. In the backward Euler step, the decoupled equations are solved iteratively. We prove that the iterations converge linearly to the solution of the coupled problem, under some conditions on the data. When θ=1/2, the algorithm is equivalent to the symplectic midpoint method. Similar to the continuous case, we also prove a discrete Gibbs free energy balance, without numerical dissipation. We compare this midpoint method with the classic backward Euler method, and two implicit-explicit time-lagging schemes. The midpoint method outperforms the other schemes in terms of rates of convergence, long-time behaviour and energy approximation, for small and large values of the time step.

We present high-order Convolution Quadratures (CQ) methods for the solution of the wave equation in unbounded domains in two and three dimensions that rely on Nyström discretizations for the Boundary Integral Equation formulations of the ensemble of associated Laplace domain modified Helmholtz problems. Both Dirichlet and Neumann boundary conditions, imposed on open-arc/open surfaces as well  as Lipschitz closed scatterers, are considered. A variety of accuracy tests are presented that showcase the high-order in time convergence (up to and including fifth order) that the Nyström CQ discretizations are capable of delivering and we compare to numerical results in the literature pertaining to time-domain multiple scattering problems solved with other methods.

We consider infinite dimensional time invariant linear systems which are exactly controllable in every positive time. We focus on the time optimal control problem. After showing the existence of solutions we show that they satisfy a maximum principle. We then conclude that under a particular approximate type controllability assumption, these controls are bang-bang. We finally explain how these abstract results apply to systems described by the Schrodinger and plate equations.