Refine
Year of publication
Language
- English (97)
Keywords
- optimal control (11)
- differential-algebraic equation (4)
- finite elements (4)
- metastability (4)
- trajectory storage (4)
- Lyapunov exponent (3)
- Sacker-Sell spectrum (3)
- exponential dichotomy (3)
- geodesic finite elements (3)
- rare events (3)
Automatic, or algorithmic, differentiation addresses the need for the accurate
and efficient calculation of derivative values in scientific computing. To this
end procedural programs for the evaluation of problem-specific functions are
transformed into programs that also compute the required derivative values
at the same numerical arguments in floating point arithmetic. Disregarding
many important implementation issues, we examine in this article complexity
bounds and other more mathematical aspects of the program transformation
task sketched above.
Quasi-Newton methods based on least change secant updating
formulas that solve linear equations $Ax=b$ in $n=\dim(x)=\dim(b)$ steps
can be expected to solve corresponding smooth nonlinear
systems $n$-step quadratically, i.e. with an $r$-order
of $\rho = 2^{1/n} = 1 + 1/n +O(1/n^2)$. The best rate one can
possibly expect on general problems is given by the positive root
$\rho_n$ of $\rho^n(\rho -1)=1$, for which
$\rho_n-1 = \ln(n)/n + O(1/n^2)$. To show that this upper bound is
actually achieved one usually has to impose a priori some kind of
linear independence condition on the sequence of steps taken by the
quasi-Newton iteration in question. Without any such assumptions we
establish in this paper the convergence order $\rho_n$ for the
two-sided rank one formula proposed by Schlenkrich et al in \cite{SGW06}.
It requires the evaluation of adjoint vectors, is invariant with respect
to linear transformations on the variable domain and combines the
properties of bounded deterioration and heredity.
We present a new extension of the well-known
Perron-Frobenius theorem to regular matrix pairs $(E,A)$.
The new extension is based on projector chains and is motivated from
the solution of positive differential-algebraic systems or descriptor
systems. We present several examples where the new condition holds, whereas conditions
in previous literature are not satisfied.
We present globally convergent multigrid methods for the nonsymmetric
obstacle problems as arising from the discretization of Black–Scholes models of
American options with local volatilities and discrete data. No tuning or regularization
parameters occur. Our approach relies on symmetrization by transformation
and data recovery by superconvergence.
We present a domain decomposition approach for the computation of the
electromagnetic field within periodic structures. We use a
Schwarz method with transparent boundary conditions at the interfaces of
the domains. Transparent boundary conditions are approximated by the
perfectly matched layer method (PML). To cope with Wood anomalies
appearing in periodic structures an adaptive strategy to determine
optimal PML parameters is developed. \\ We focus on the application to
typical EUV lithography line masks. Light propagation within the
multi-layer stack of the EUV mask is treated analytically. This results
in a drastic reduction of the computational costs and allows for the
simulation of next generation lithography masks
on a standard personal computer.
We present a new solver for large-scale two-body contact problems in nonlinear elasticity. It is based on an SQP-trust-region approach.
This guarantees global convergence to a first-order critical point of
the energy functional. The linearized contact conditions are
discretized using mortar elements. A
special basis transformation known from linear contact problems
allows to use a monotone multigrid solver for the inner quadratic programs.
They can thus be solved with multigrid complexity. Our algorithm
does not contain any regularization or penalization parameters,
and can be used for all hyperelastic material models.
The purpose of the paper is to apply monotone multigrid methods
to static and dynamic biomechanical contact problems.
In space, a finite element method involving a mortar
discretization of the contact conditions is used.
In time, a new contact--stabilized Newmark scheme is presented.
Numerical experiments for a two body Hertzian contact problem
and a biomechanical knee problem are reported.
For the solution of nonlinear equation systems
quasi-Newton methods based on low-rank updates are of particular interest. We analyze a class
of TR1 update formulas to approximate the system Jacobian. The local q-superlinear convergence for nonlinear problems is proved for a particular subclass of updates. Moreover, we give an estimate of the r-order of convergence. Numerical results comparing the TR1 method to Newton's and other quasi-Newton methods atr presented.
A new concept is introduced for the adaptive finite element discretization of partial differential equations that have a sparsely
representable solution. Motivated by recent work on compressed sensing, a recursive mesh refinement procedure is presented that uses linear programming to find a good approximation to the sparse solution on a given refinement level. Then only those parts of the mesh are refined that belong to nonzero expansion coefficients. Error estimates for this procedure are refined and the behavior of the procedure is demonstrated via some simple elliptic model problems.
Interior Point Methods in Function Space for State Constraints - Inexact Newton and Adaptivity
(2009)
We consider an interior point method in function space for PDE constrained optimal control problems with state constraints. Our emphasis is on the construction and analysis of an algorithm that integrates a Newton path-following method with adaptive grid refinement. This is done in the framework of inexact Newton methods in function space, where the discretization error of each Newton step is controlled by adaptive grid refinement in the innermost loop. This allows to perform most of the required Newton steps on coarse grids, such that the overall computational time is dominated by the last few steps. For this purpose we propose an a-posteriori error estimator for a problem suited norm.
A continuity result for Nemyckii Operators and some applications in PDE constrained optimal control
(2008)
This work explores two applications of a classical result on the continuity of Nemyckii operators to optimal control with PDEs. First, we present an alternative approach to the analysis of Newton's method for function space problems involving semi-smooth Nemyckii operators. A concise proof for superlinear convergence is presented, and sharpened bounds on the rate of convergence are derived. Second, we derive second order sufficient conditions for problems, where the underlying PDE has poor regularity properties. We point out that the analytical structure in both topics is essentially the same.
We propose and analyse an interior point path-following method in function space
for state constrained optimal control. Our emphasis is on proving convergence in
function space and on constructing a practical path-following algorithm. In particular, the introduction of a pointwise damping step leads to a very efficient method, as verified by numerical experiments.
We consider hybrid systems of differential-algebraic equations and present
a general framework for general nonlinear over- and underdetermined hybrid
systems that allows the
analysis of existence and uniqueness and the application of index reduction
methods for hybrid differential-algebraic systems.
A particular difficulty in the numerical simulation of hybrid systems is
(numerical) chattering, i.e., fast oscillations between modes of operations.
A regularization technique using sliding modes allows to regularize the
system behavior in the case of chattering.
Further, we show how chattering behavior during the numerical solution can
be prevented using sliding mode simulation. The advantage of the sliding mode
simulation is illustrated by numerical examples.
A generalization of the method of Chu, Liu and Mehrmann
for the computation of the Hamiltonian real Schur form is presented.
The new method avoids some of the difficulties that may arise when
a Hamiltonian matrix has tightly clustered groups of eigenvalues.
A detailed analysis of the method is presented and several numerical examples demonstrate the superior behavior of the method.
We present and analyze novel hierarchical a posteriori error estimates
for self-adjoint elliptic obstacle problems.
Our approach differs from straightforward, but non-reliable estimators~\cite{RHWHoppe_RKornhuber_1994a}
by an additional extra term accounting for the deviation
of the discrete free boundary in the localization step.
We prove efficiency and reliability
on a saturation assumption and a regularity condition on the underlying grid.
Heuristic arguments suggest
that the extra term is of higher order and preserves full locality.
Numerical computations confirm our theoretical findings.
The main focus of this paper is on an a-posteriori analysis for the method of proper orthogonal decomposition (POD) applied to optimal control problems governed by
parabolic and elliptic PDEs. Based on a perturbation method it is deduced how far the suboptimal
control, computed on the basis of the POD model, is from the (unknown)
exact one. Numerical examples illustrate the realization of the proposed approach for linear-quadratic problems governed by parabolic and elliptic partial differential equations.
The complexity of molecular kinetics can be reduced significantly by a restriction to metastable conformations which are almost invariant sets of molecular dynamical systems. With the Robust Perron Cl uster Analysis PCCA+, developed by Weber and Deuflhard, we have a tool available which can be used to identify these conformations from a transition probability matrix. This method can also be applied to the corresponding transition rate matrix which provides important information concerning transition pathways of single molecules. In the present paper, we explain the relationship between these tw o concepts and the extraction of conformation kinetics from transition rates. Moreover, we show how transition rates can be approximated and conclude with numerical examples.
Wigner functions are functions on classical phase space, which are in one-to-one correspondence to square integrable functions on configuration space. For molecular quantum systems, classical transport of Wigner functions provides the basis of asymptotic approximation methods in the high energy regime. The article addresses the sampling of Wigner functions by Monte Carlo techniques. The approximation step is realized by an adaption of the Metropolis algorithm for real-valued functions with disconnected support. The quadrature, which computes values of the Wigner function, uses importance sampling with a Gaussian weight function. The numerical experiments combine the sampling with a surface hopping algorithm for non-adiabatic quantum dynamics. In agreement with theoretical considerations, the obtained results show an accuracy of two to four percent.
In this review, we intend to clarify the underlying ideas and the relations
between various multigrid methods ranging from subset decomposition,
to projected subspace decomposition and truncated multigrid.
In addition, we present a novel globally convergent inexact active set method
which is closely related to truncated multigrid. The numerical properties
of algorithms are carefully assessed by means of a degenerate problem and
a problem with a complicated coincidence set.
We present a new inexact nonsmooth Newton method for the solution
of convex minimization problems with piecewise smooth, pointwise
nonlinearities. The algorithm consists of a nonlinear smoothing
step on the fine level and a linear coarse correction.
Suitable postprocessing guarantees global convergence
even in the case of a single multigrid step for each linear subproblem.
Numerical examples show that the overall efficiency
is comparable to multigrid for similar linear problems.
In this paper, we propose and investigate numerical methods based on QR factorization for computing all or some Lyapunov or Sacker-Sell spectral intervals for
linear differential-algebraic equations.
Furthermore, a perturbation and error analysis for these methods is presented. We
investigate how errors in the data and in the numerical integration affect the
accuracy of the approximate spectral intervals. Although we need to integrate
numerically some differential-algebraic systems on usually very long
time-intervals, under certain assumptions, it is shown that the error of the
computed spectral intervals can be controlled by the local error of numerical
integration and the error in solving the algebraic constraint.
Some numerical examples are presented to illustrate the theoretical results.
Motivated by the analysis of passive control systems, we undertake a detailed perturbation analysis of Hamiltonian matrices that have eigenvalues on the imaginary axis. We construct minimal Hamiltonian perturbations that move and coalesce eigenvalues of opposite sign characteristic to form multiple eigenvalues with mixed sign characteristics, which are then moved from the imaginary axis to specific locations in the complex plane by small Hamiltonian perturbations. We also present a numerical method to compute upper bounds for the minimal perturbations that move all eigenvalues of a given Hamiltonian matrix outside a vertical strip along the imaginary axis.
The PSurface Library
(2010)
We describe psurface, a C++ library that allows to store and access piecewise linear mappings between simplicial surfaces in $\R^2$ and $\R^3$. These mappings are stored in a graph data structure and can be constructed explicitly, by projection, or by surface simplification. Piecewise linear maps can be used, e.g., to construct boundary
approximations for finite element grids, and grid intersections for domain decomposition methods. In computer graphics the mappings allow to build level-of-detail representations as well as texture- and bump maps. We document the data structures and algorithms used and show how \psurface is used in the numerical analysis framework Dune
and the visualization software Amira.
This paper is devoted to the numerical approximation of Lyapunov and Sacker-Sell spectral intervals for linear differential-algebraic equations (DAEs). The spectral analysis for DAEs is improved and the concepts of leading directions and solution subspaces associated with spectral intervals are extended to DAEs. Numerical methods
based on smooth singular value decompositions are introduced for computing all or only some spectral intervals and their associated leading directions. The numerical algorithms as well as implementation issues are discussed in detail and numerical examples are presented to illustrate the theoretical results.
This paper presents three different adaptive algorithms for eigenvalue problems
associated with non-selfadjoint partial differential
operators. The basis for the developed algorithms is a homotopy method.
The homotopy method starts from a well-understood selfadjoint problem,
for which well-established adaptive methods are available.
Apart from the adaptive grid refinement, the progress of the homotopy as
well as the solution of the iterative method are adapted to balance the contributions
of the different error sources.
The first algorithm balances the homotopy, discretization and approximation errors with respect
to a fixed step-size $\tau$ in the homotopy.
The second algorithm combines the adaptive step-size control for the homotopy
with an adaptation in space that ensures an error below a fixed tolerance $\varepsilon$.
The third algorithm allows the complete adaptivity in space,
homotopy step-size as well as the iterative algebraic eigenvalue solver.
All three algorithms are compared in numerical examples.
In this work we propose a general framework for the structured perturbation
analysis of several classes of structured matrix polynomials in homogeneous
form, including complex symmetric, skew-symmetric, even and odd matrix polynomials. We introduce structured backward errors for approximate eigenvalues and eigenvectors and we construct minimal structured perturbations such that an approximate eigenpair is an exact eigenpair of an appropriately perturbed matrix polynomial. This work extends previous work for the non-homogeneous case (we include infinite eigenvalues) and we show that the structured backward errors improve the known unstructured backward errors.
The paper proposes goal-oriented error estimation and mesh refinement
for optimal control problems with elliptic PDE constraints using the value
of the reduced cost functional as quantity of interest. Error representation,
hierarchical error estimators, and greedy-style error indicators are derived and
compared to their counterparts when using the all-at-once cost functional as
quantity of interest. Finally, the efficiency of the error estimator and generated
meshes are demonstrated on numerical examples.
We investigate geodesic finite elements for functions with values in a space of zero curvature, like a torus or the M\"obius strip. Unlike in the general case, a closed-form expression for geodesic finite element functions is then available. This simplifies computations, and allows us to prove optimal estimates for the interpolation error in 1d and 2d. We also show the somewhat surprising result that the discretization by Kirchhoff transformation of the Richards equation proposed by Berninger et al. is a discretization by geodesic finite elements in the manifold $\mathbb{R}$ with a special metric.
We propose a generalization of the Structured Doubling Algorithm (SDA) to compute invariant subspaces
of structured matrix pencils
that arise in the context of solving linear quadratic optimal control problems.
The new algorithm is
designed to attain better accuracy when the classical Riccati equation approach for the solution of the optimal control problem is not well suited because
the stable and unstable invariant subspaces are not well separated (due to eigenvalues near or on the imaginary
axis) or in the case when the Riccati solution does not exist at all. We analyze the convergence
of the method and compare the new method with the classical SDA algorithm as well as some recent structured QR-methods.
This paper presents concepts and implementation of the finite element toolbox Kaskade 7, a flexible C++ code for solving elliptic and parabolic PDE systems. Issues such as problem formulation, assembly and adaptivity are discussed at the example of optimal control problems. Trajectory compression for parabolic optimization problems is considered as a case study.
Deflated and augmented Krylov subspace methods: Basic Facts and a Breakdown-free deflated MINRES
(2011)
In this paper we consider deflation and augmentation techniques for accelerating
the convergence of Krylov subspace methods for the solution of nonsingular linear
algebraic systems. The two techniques are conceptually different from
preconditioning. Deflation "removes" certain parts from the operator, while
augmentation adds a subspace to the Krylov subspace. Both approaches have been
used in a variety of methods and settings. For Krylov subspace methods that
satisfy a (Petrov-) Galerkin condition we show that augmentation can in general
be achieved implicitly by projecting the residuals appropriately and correcting
the approximate solutions in a final step. In this context, we analyze known
methods to deflate CG, GMRes and MinRes. Our analysis reveals that the recently
proposed RMinRes method can break down. We show how such breakdowns can be
avoided by choosing a special initial guess, and we derive a breakdown-free
deflated MinRes method. In numerical experiments we study the properties of
different variants of MinRes analyzed in this paper.
One of the most challenging problems in dynamic concurrent multiscale
simulations is the reflectionless transfer of physical quantities between the
different scales. In particular, when coupling molecular dynamics and finite
element discretizations in solid body mechanics, often spurious wave reflections
are introduced by the applied coupling technique. The reflected waves are
typically of high frequency and are arguably of little importance in the domain
where the finite element discretization drives the simulation.
In this work, we provide an analysis of this phenomenon.
Based on the gained
insight, we derive a new coupling approach, which neatly separates high and low
frequency waves. Whereas low frequency waves are permitted to
bridge the scales, high frequency waves can be removed by applying damping techniques without affecting the coupled share of the solution. As a consequence, our new method almost completely eliminates unphysical wave reflections and deals in a consistent way with waves of arbitrary frequencies. The separation of
wavelengths is achieved by employing a discrete $L^2$-projection, which acts as a
low pass filter. Our coupling constraints enforce matching in the range of this projection. With respect to the numerical realization this approach
has the advantage of a small number of constraints, which is computationally
efficient. Numerical results in one and two dimensions confirm our theoretical
findings and illustrate the performance of our new weak coupling approach.
We investigate optimal elliptic
regularity (within the scale of Sobolev spaces) of anisotropic
div--grad operators in three dimensions at a multi-material vertex on
the Neumann boundary part of a polyhedral spatial domain. The
gradient of a solution to the corresponding elliptic PDE (in a
neighbourhood of the vertex) is integrable to an index greater than
three.
We consider Large Deformation Diffeomorphic Metric Mapping of general $m$-currents. After stating an optimization algorithm in the function space of admissable morph generating velocity fields, two innovative aspects in this framework are presented and numerically investigated: First, we spatially discretize the velocity field with conforming adaptive finite elements and discuss advantages of this new approach. Second, we directly compute the temporal evolution of discrete $m$-current attributes.
We consider the mechanical coupling of a geometrically exact Cosserat rod to a linear elastic continuum. The coupling conditions are formulated in the nonlinear rod configuration space. We describe a Dirichlet--Neumann algorithm for the coupled system, and use it to simulate the static stresses in a human knee joint, where the Cosserat rods are models for the ligaments.
In optimal control problems with nonlinear time-dependent 3D PDEs, full 4D discretizations are usually prohibitive due to the storage requirement. For this reason gradient and Newton type methods working on the reduced functional are often employed. The computation of the reduced gradient requires one solve of the state equation forward in time, and one backward solve of the adjoint equation. The state enters into the adjoint equation, again requiring the storage of a full 4D data set. We propose a lossy compression algorithm using an inexact but cheap predictor for the state data, with additional entropy coding of prediction errors. As the data is used inside a discretized, iterative algorithm, lossy compression
maintaining a certain error bound turns out to be sufficient.
We introduce geodesic finite elements as a conforming way to discretize partial differential equations for functions $v : \Omega \to M$, where $\Omega$ is an open subset of $\R^d$ and $M$ is a Riemannian manifold. These geodesic finite elements naturally generalize standard first-order
finite elements for Euclidean spaces. They also generalize the geodesic finite elements proposed for $d=1$ by the author. Our formulation is equivariant under isometries of $M$, and hence preserves objectivity of continuous problem formulations. We concentrate on partial differential equations that can be formulated as minimization problems. Discretization leads to algebraic minimization problems on product manifolds $M^n$. These can be solved efficiently
using a Riemannian trust-region method. We propose a monotone multigrid method to solve the constrained inner problems with linear multigrid speed. As an example we numerically compute harmonic maps from a domain in $\R^3$
to $S^2$.
Logical modeling of biological regulatory networks gives rise to a representation of the system's dynamics as a so-called state transition graph. Analysis of such a graph in its entirety allows for a comprehensive understanding of the functionalities and behavior of the modeled system. However, the size of the vertex set of the graph is exponential in the number of the network components making analysis costly, motivating development of reduction methods. In this paper, we present results allowing for a complete description of an asynchronous state transition graph of a Thomas network solely based on the analysis of the subgraph induced by certain extremal states. Utilizing this notion, we compare the behavior of a simple multi-valued network and a corresponding Boolean network and analyze the conservation of dynamical properties between them. Understanding the relation between such coarser and finer models is a necessary step towards meaningful network reduction as well as model refinement methods.
Mathematical modeling often helps to provide a systems perspective on gene regulatory networks. In particular, qualitative approaches are useful when detailed kinetic information is lacking. Multiple methods have been developed that implement qualitative information in different ways, e.g., in purely discrete or hybrid discrete/continuous models. In this paper, we compare the discrete asynchronous logical modeling formalism for gene regulatory networks due to R. Thomas with piecewise affine differential equation models.
We provide a local characterization of the qualitative dynamics of a piecewise affine differential equation model using the discrete dynamics of a corresponding Thomas model. Based on this result, we investigate the consistency of higher-level dynamical properties such as attractor characteristics and reachability. We show that although the two approaches are based on equivalent information, the resulting qualitative dynamics are different. In particular, the dynamics of the piecewise affine differential equation model is not a simple refinement of the dynamics of the Thomas model.
The authors propose a recycling MINRES scheme for a solution of subsequent self-adjoint linear systems as appearing, for example, in the Newton process for solving nonlinear equations. Ritz vectors are automatically extracted from one MINRES run and then used for self-adjoint deflation in the next. The method is designed to work with a preconditioner and arbitrary inner products. Numerical experiments with nonlinear Schrödinger equations indicate a substantial decrease in computation time when recycling is used.
We consider a shape implant design problem that arises in the context of facial surgery. We introduce a reformulation as an optimal control problem, where the control acts as a boundary force. The state is modelled as a minimizer of a polyconvex hyperelastic energy functional. We show existence of optimal solutions and derive - on a formal level - first order optimality conditions. Finally, preliminary numerical results are presented.
The pole condition approach for deriving transparent boundary conditions is extended to the time-dependent, two-dimensional case. Non-physical modes of the solution are identified by the position of poles of the solution's spatial Laplace transform in the complex plane. By requiring the Laplace transform to be analytic on some
problem dependent complex half-plane, these modes can be
suppressed. The resulting algorithm computes a finite number of coefficients of a series expansion of the Laplace transform, thereby providing an approximation to the exact boundary condition. The resulting error decays super-algebraically with the number of coefficients, so relatively few additional degrees of freedom are
sufficient to reduce the error to the level of the discretization error in the interior of the computational domain. The approach shows good results for the Schroedinger and the drift-diffusion equation
but, in contrast to the one-dimensional case, exhibits instabilities for the wave and Klein-Gordon equation. Numerical examples are shown that demonstrate the good performance in the former and the instabilities in the latter case.
Flux variability analysis (FVA) is an important tool to further analyze the results obtained by flux balance analysis (FBA) on genome-scale metabolic networks. Standard FVA may predict unbounded fluxes through some reactions in the network even if the nutrient uptake rate is bounded. These fluxes violate the second law of thermodynamics. They may be eliminated by extending flux variability analysis with thermodynamic constraints.
We present a new algorithm for efficient flux variability (and flux balance) analysis with thermodynamic constraints, suitable for analyzing genome-scale metabolic networks. We first show that flux balance analysis with thermodynamic constraints is NP-hard. Then we derive a theoretical tractability result, which can be applied to metabolic networks in practice. We use this result to develop a new constraint programming algorithm Fast-tFVA for fast flux variability analysis with thermodynamic constraints (tFVA). Computational comparisons with previous methods demonstrate the efficiency of the new method. For tFVA, a speed-up of factor 30-300 is achieved.
In an analysis of genome-scale metabolic networks in the BioModels database, we found that in 485 out of 716 networks additional irreversible or fixed reactions could be detected.
This article describes Fortran 77 subroutines for computing eigenvalues and invariant subspaces
of Hamiltonian and skew-Hamiltonian matrices. The implemented algorithms are based on orthogonal
symplectic decompositions, implying numerical backward stability as well as symmetry
preservation for the computed eigenvalues. These algorithms are supplemented with balancing and
block algorithms, which can lead to considerable accuracy and performance improvements. As a
by-product, an efficient implementation for computing symplectic QR decompositions is provided.
We demonstrate the usefulness of the subroutines for several, practically relevant examples.
Stewart's recently introduced Krylov-Schur algorithm
is a modification of the implicitly restarted Arnoldi algorithm which
employs reordered Schur decompositions to perform restarts and de-
ations in a numerically reliable manner. This paper describes a variant
of the Krylov-Schur algorithm suitable for addressing eigenvalue
problems associated with products of large and sparse matrices. It
performs restarts and de
ations via reordered periodic Schur decompositions
and, by taking the product structure into account, it is
capable to achieve qualitatively better approximations to the eigenvalues
of small magnitude.
The role of larger bulges in the QR algorithm is controversial. Large bulges are infamous
for having a strong, negative influence on the convergence of the implicit shifted QR algorithm.
This paper provides a new explanation of this shift blurring effect by connecting the computation of
the first column of the shift polynomial to the notoriously ill-conditioned pole placement problem.
To avoid shift blurring, modern variants of the QR algorithm employ chains of tightly coupled tiny
bulges instead of one large bulge. It turns out that larger bulges still play a positive role in these
variants; a slight increase of the bulge sizes often results in considerable performance improvements.
We investigate the condition number for a complex eigenvalue of a real matrix under
real perturbations. Based on an explicit formula, it is shown that this number is never
smaller than 1/
p
2 times the corresponding condition number with respect to complex
perturbations. This result can be generalized to the condition number of an arbitrary
complex-valued function under real perturbations. This extends to related condition
numbers.
We give an algorithm to compute N steps of a convolution quadrature approximation
to a continuous temporal convolution using only O(N logN) multiplications and O(logN) active
memory. The method does not require evaluations of the convolution kernel, but instead O(logN)
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 integrodifferential
equations of convolution type. In a numerical example we apply it to solve a subdiffusion
equation with transparent boundary conditions.
In this paper we propose a new finite element realization of the Perfectly Matched
Layer method (PML-method). Our approach allows to deal with a wide class of
polygonal domains and with certain types of inhomogeneous exterior domains.
Among the covered inhomogeneities are open waveguide structures playing an essential
role in integrated optics. We give a detailed insight into implementation
aspects. Numerical examples show exponential convergence behavior to the exact
solution with the thickness of the PML sponge layer.
Lyapunov and exponential dichotomy spectral theory is extended
from ordinary differential equations (ODEs) to nonautonomous
differential-algebraic equations (DAEs). By using orthogonal
changes of variables, the original DAE system is transformed into
appropriate condensed forms, for which concepts such as Lyapunov
exponents, Bohl exponents, exponential dichotomy and spectral
intervals of various kinds can be analyzed via the resulting
underlying ODE. Some essential differences between the spectral
theory for ODEs and that for DAEs are pointed out. Numerical
methods for computing the spectral intervals associated with
Lyapunov and Sacker-Sell (exponential dichotomy) spectra are
derived by modifying and extending those methods proposed for ODEs. Perturbation theory and error analysis are discussed, as
well. Finally, some numerical examples are presented to illustrate
the theoretical results and the properties of the numerical
methods.
A Generic Grid Interface for Parallel and Adaptive Scientific Computing. Part I: Abstract Framework
(2007)
We give a mathematically rigorous definition of a grid for algorithms solving
partial differential equations. Unlike previous approaches, our grids have a
hierarchical structure. This makes them suitable for geometric multigrid
algorithms and hierarchical local grid refinement. The description is also
general enough to include geometrically nonconforming grids. The definitions
in this article serve as the basis for an implementation of an abstract grid
interface as C++ classes in the DUNE.
In a companion paper [Matheon-Preprint 403] we introduced an abstract definition of a parallel and adaptive hierarchical grid for scientific computing. Based on this
definition we derive an efficient interface specification as a set of C++ classes.
This interface separates the applications from the grid data structures.
Thus, user implementations become independent of the underlying grid
implementation. Modern C++ template techniques are used to provide an
interface implementation without big performance losses.
The implementation is realized as part of the
software environment DUNE.
Numerical tests demonstrate the flexibility and the efficiency of our approach.
In many applications such as data compression, imaging or
genomic data analysis,
it is important to approximate a given $m\times n$ matrix $A$
by a matrix $B$ of rank at most $k$ which is much smaller than $m$ and $n$.
The best rank $k$ approximation can be determined via
the singular value decomposition
which, however, has prohibitively
high computational complexity and storage requirements
for very large $m$ and $n$.
We present an optimal least squares algorithm for computing a rank $k$
approximation to an $m\times n$ matrix $A$ by reading
only a limited number of rows and columns of $A$.
The algorithm has complexity $\mathcal O(k^2\max(m,n))$ and
allows to iteratively improve given rank $k$
approximations by reading additional rows and
columns of $A$. We also show how this approach can be extended
to tensors and present numerical results.
In this paper, we discuss stability properties of positive descriptor systems in the continuous-time as well as in the discrete-time case. We present different characterisations of positivity and establish generalised stability criteria for the case of positive descriptor systems. We show that if the spectral projector onto the right finite deflating subspace of the matrix pair $(E,A)$ is non-negative, then all stability criteria for standard positive systems take a comparably simple form in the positive descriptor case. Furthermore, we provide sufficient conditions that guarantee entry-wise non-negativity along with positive semi-definiteness of solutions of generalised projected Lyapunov equations. As an application of the framework established throughout this paper, we exemplarily generalise two criteria for the stability of two switched standard positive systems under arbitrary switching to the descriptor case.
Adjoint Broyden a la GMRES
(2007)
It is shown here that a compact storage implementation of a quasi-Newton
method based on the adjoint Broyden update reduces in the affine
case exactly to the well established GMRES procedure. Generally,
storage and linear algebra effort per step are small multiples of $n\cdot k$,
where $n$ is the number of variables and $k$ the number of steps taken
in the current cycle. In the affine case the storage is exactly $(n+k)\cdot k$
and in the nonlinear case the same bound can be achieved if adjoints,
i.e. transposed Jacobian-vector products are available. A
transposed-free variant that relies exclusively on Jacobian-vector
products (or possibly their approximation by divided differences)
requires roughly twice the storage and turns out to be somewhat slower
in our numerical experiments reported at the end.
Structured eigenvalue conditioning and backward error of a class of polynomial eigenvalue problems
(2007)
Characterisations of simple eigenvalues of complex matrix polynomials with *-even/odd and *-palindromic/antipalindromic structures that have the same normwise condition number with respect to structure preserving and arbitrary perturbations are obtained. Here * denotes either the transpose T or the conjugate tranpose *. In the process we obtain formulae for the normwise structured condition number of simple eigenvalues of T-palindromic/antipalindromic and *-even/odd polynomials. Moreover, conditions under which the normwise structured backward error of approximate eigenvalues of such polynomials is equal to the unstructured error are also derived. These lead to complete characterisations of approximate eigenvalues that have the same structured and unstructured backward errors for the *-even/odd and T-even/odd polynomials.
The paper considers the time integration of frictionless dynamical contact problems between viscoelastic bodies in the frame of the Signorini condition. Among the numerical integrators, interest focuses on the contact-stabilized Newmark method recently suggested by Deuflhard et al., which is compared to the classical Newmark method and an improved energy dissipative version due to Kane et al. In the absence of contact, any such variant is equivalent to the Störmer-Verlet scheme, which is well-known to have consistency order 2. In the presence of contact, however, the classical approach to discretization errors would not show consistency at all because of the discontinuity at the contact. Surprisingly, the question of consistency in the constrained situation has not been solved yet. The present paper fills this gap by means of a novel proof technique using specific norms based on earlier perturbation results due to the authors. The corresponding estimation of the local discretization error requires the bounded total variation of the solution. The results have consequences for the construction of an adaptive timestep control, which will be worked out subsequently in a forthcoming paper.
We present an extension module for the Dune system. This module, called dune-subgrid, allows to mark elements of another Dune hierarchical grid. The set of marked elements can then be accessed as a Dune grid in its own right. dune-subgrid is free software and is available for download. We describe the functionality and use of dune-subgrid, comment on its implementation, and give two example applications.
First, we show how dune-subgrid can be used for micro-FE simulations of trabecular bone. Then we present an algorithm that allows to use exact residuals for the adaptive solution of the spatial problems of time-discretized evolution equations.
We present a hierarchical a~posteriori error analysis for the
minimum value of the energy functional in symmetric obstacle
problems. The main result is that the energy of the exact solution
is, up to data oscillation, equivalent to an appropriate
hierarchical estimator. The proof of the main result does not
invoke any saturation assumption. Moreover, we prove an a
posteriori error estimate indicating that the estimator from
\cite{RHWHoppe_RKornhuber_1994a} is asymptotically reliable and we
give sufficient conditions for the validity of a saturation
assumption. Finally, we corroborate and complement our theoretical
results with numerical experiments.
We introduce geodesic finite elements as a new way to discretize
the nonlinear configuration space of a geometrically exact Cosserat rod.
These geodesic finite elements naturally generalize standard one-dimensional
finite elements to spaces of functions with values in a Riemannian manifold.
For the special orthogonal group, our approach reproduces the
interpolation formulas of [Crisfield/Jelenic:1999].
Geodesic finite elements are
conforming and lead to objective and path-independent problem formulations.
We introduce geodesic finite elements for general Riemannian manifolds,
discuss the relationship between geodesic finite elements and
coefficient vectors, and estimate the interpolation error.
Then we use them to find static equilibria of hyperelastic Cosserat rods.
Using the Riemannian trust-region algorithm of [Absil/Mahony/Sepulchre:2008]
we show numerically that the discretization error depends optimally on
the mesh size.
The identification of metastable conformations of molecules plays an
important role in computational drug design. One main difficulty is the
fact that the underlying dynamic processes take place in high dimensional
spaces. Although the restriction of degrees of freedom to a few dihedral
angles significantly reduces the complexity of the problem, the existing
algorithms are time-consuming. They are mainly based on the approximation
of a transfer operator by an extensive sampling of states according
to the Boltzmann distribution and short-time Hamiltonian dynamics simulations.
We present a method which can identify metastable conformations
without sampling the complete distribution. Our algorithm is based
on local transition rates and uses only pointwise information about the
potential energy surface. In order to apply the cluster algorithm PCCA+,
we compute a few eigenvectors of the rate matrix by the Jacobi-Davidson
method. Interpolation techniques are applied to approximate the thermodynamical
weights of the clusters. The concluding example illustrates
our approach for epigallocatechine, a molecule which can be described by
seven dihedral angles.
In order to compute the thermodynamic weights of the different metastable conformations
of a molecule, we want to approximate the molecule’s Boltzmann distribution in a reasonable
time. This is an essential issue in computational drug design. The energy landscape of active
biomolecules is generally very rough with a lot of high barriers and low regions. Many of the
algorithms that perform such samplings (e.g. the hybrid Monte Carlo method) have difficulties
with such landscapes. They are trapped in low-energy regions for a very long time and cannot
overcome high barriers. Moving from one low-energy region to another is a very rare event. For
these reasons, the distribution of the generated sampling points converges very slowly against
the thermodynamically correct distribution of the molecule.
The idea of ConfJump is to use a priori knowledge of the localization of low-energy regions
to enhance the sampling with artificial jumps between these low-energy regions. The artificial
jumps are combined with the hybrid Monte Carlo method. This allows the computation of
some dynamical properties of the molecule. In ConfJump, the detailed balance condition is
satisfied and the mathematically correct molecular distribution is sampled.
Biochemical interactions are determined by the 3D-structure of the involved components –
thus the identification of conformations is a key for many applications in rational drug design.
ConFlow is a new multilevel approach to conformational analysis with main focus on
completeness in investigation of conformational space.
In contrast to known conformational analysis, the starting point for design is a space-based
description of conformational areas. A tight integration of sampling and analysis leads to an
identification of conformational areas simultaneously during sampling. An incremental
decomposition of high-dimensional conformational space is used to guide the analysis. A new
concept for the description of conformations and their path connected components based on
convex hulls and Hypercubes is developed. The first results of the ConFlow application
constitute a ‘proof of concept’ and are further more highly encouraging. In comparison to
conventional industrial applications, ConFlow achieves higher accuracy and a specified
degree of completeness with comparable effort.
In this paper, we investigate the interconversion processes of the major flame retardant - 1,2,5,6,9,10-hexabromocyclododecane (HBCD) - by the means of statistical thermodynamics based on classical force-fields. Three ideas will be presented. First, the application of classical hybrid Monte-Carlo simulations for quantum mechanical processes will be justified. Second, the problem of insufficient convergence properties of hybrid Monte-Carlo methods for the generation of low temperature canonical ensembles will be solved by an interpolation approach. Furthermore, it will be shown how free energy differences can be used for a rate matrix computation. The results of our numerical simulations will be compared to experimental results.
Keywords: Markov process, molecular dynamics, rate matrix
This article deals with an efficient sampling of the stationary distribution
of dynamical systems in the presence of metastabilities. For such
systems, standard sampling schemes suffer from trapping problems and
critical slowing down. Starting multiple trajectories in different regions of
the sampling space is a promising way out. The different samplings represent
the stationary distribution locally very well, but are still far away
from ergodicity or from the global stationary distribution. We will show
how these samplings can be joined together in order to get one global
sampling of the stationary distribution.
For the treatment of equilibrated molecular systems in a heat bath we
propose a transition state theory that is based on conformation dynamics.
In general, a set-based discretization of a Markov operator P does not
preserve the Markov property. In this article, we propose a discretization
method which is based on a Galerkin approach. This discretization
method preserves the Markov property of the operator and can be interpreted
as a decomposition of the state space into (fuzzy) sets. The
conformation-based transition state theory presented here can be seen as
a first step in conformation dynamics towards the computation of essential
dynamical properties of molecular systems without time-consuming
molecular dynamics simulations.
For an analysis of a molecular system from a computational statistical thermodynamics point of view, extensive
molecular dynamics simulations are very inefficient. During this procedure, at lot of redundant data is generated. Whereas
the algorithms spend most of the computing time for a sampling of configurations within the basins of the potential energy
landscape of the molecular system, the important information about the long-time behaviour of the molecules is given by
transition regions and barriers between the basins, which are sampled rarely only. Thinking of molecular dynamics trajectories,
researchers try to figure out which kind of dynamical model is suitable for an efficient simulation. This article suggests to
change the point of view from extensive simulation of molecular dynamics trajectories to more efficient sampling strategies
of the conformation dynamics approach.
Deformable surface models are often represented as triangular meshes in image segmentation applications. For a fast and easily regularized deformation onto the target object boundary, the vertices of the mesh are commonly moved along line segments (typically surface normals). However, in case of high mesh curvature, these lines may intersect with the target boundary at “non-corresponding” positions, or may not intersect at all. Consequently, certain deformations cannot be achieved. We propose omnidirectional displacements for deformable surfaces (ODDS) to overcome this limitation. ODDS allow each vertex to move not only along a line segment but within a surrounding sphere, and achieve globally optimal deformations subject to local regularization con-
straints. However, allowing a ball-shaped instead of a linear range of motion per vertex significantly increases runtime and memory. To alleviate this drawback, we propose a hybrid approach, fastODDS, with improved runtime and reduced memory requirements. Furthermore, fastODDS can also cope with simultaneous segmentation of multiple objects. We show the theoretical benefits of ODDS with experiments on synthetic data, and evaluate ODDS and fastODDS quantitatively on clinical image data of the mandible and the hip bones. There, we assess both the global segmentation accuracy as well as local accuracy in high curvature regions, such as the tip-shaped mandibular coronoid processes and the ridge-shaped acetabular rims of
the hip bones.
The mixed regularity of electronic wave functions in fractional order and weighted Sobolev spaces
(2012)
The paper continues the study of the regularity of electronic wave functions in Hilbert spaces of mixed derivatives. It is shown that the eigenfunctions of electronic Schr\"odinger operators and their
exponentially weighted counterparts possess, roughly speaking, square integrable mixed weak derivatives of fractional order $\vartheta$ for $\vartheta<3/4$. The bound $3/4$ is best possible and can neither be reached nor surpassed. Such results are important for the study
of sparse grid-like expansions of the wave functions and show that their asymptotic convergence rate measured in terms of the number of ansatz functions involved does not deteriorate with the number of electrons.
We study the expansion of the eigenfunctions of Schrödinger operators with smooth confinement potentials in Hermite functions; confinement potentials are potentials that become unbounded at infinity. The key result is that such eigenfunctions and all their derivatives decay more rapidly than any exponential function under some mild growth
conditions to the potential and its derivatives. Their expansion in Hermite functions converges therefore very fast, super-algebraically.
Particle methods have become indispensible in conformation dynamics to compute transition rates in protein folding, binding processes and molecular design, to mention a few. Conformation dynamics requires at a decomposition of a molecule's position space into metastable conformations. In this paper, we show how this decomposition can be obtained via the design of either ``soft'' or ``hard'' molecular conformations. We show, that the soft approach results in a larger metastabilitiy of the decomposition and is thus more advantegous. This is illustrated by a simulation of Alanine Dipeptide.
Convergence Analysis of Smoothing Methods for Optimal Control of Stationary Variational Inequalities
(2012)
In the article an optimal control problem subject to a stationary variational inequality is investigated. The optimal control problem is complemented with pointwise
control constraints. The convergence of a smoothing scheme is analyzed. There, the variational inequality is replaced by a semilinear elliptic equation. It is shown that solutions of the regularized optimal control problem converge to solutions of the original one. Passing to
the limit in the optimality system of the regularized problem allows to prove C-stationarity of local solutions of the original problem. Moreover, convergence rates with respect to the regularization parameter for the error in the control are obtained. These rates coincide with
rates obtained by numerical experiments, which are included in the paper.
We derive a-priori estimates on the length of the primal-dual path that results from a Moreau-Yosida approximation of the feasible set for state constrained optimal control problems. These bounds depend on the regularity of the state and the dimension of the problem. Comparison with numerical results indicates that these bounds are sharp and are attained for the case of a single active point.
We present a unified computational framework for matching 3d geometric objects (points, lines, surfaces, volumes) of highly varying shape. Our approach is based on the Large Deformation Diffeomorphic Metric Mapping (LDDMM) method acting on $m$-currents. After stating an optimization algorithm in the function space of admissible morph generating velocity fields, two innovative aspects in this framework are presented: First, we spatially discretize the velocity field with conforming adaptive finite elements and discuss advantages of this new approach. Secondly, we directly compute the temporal evolution of discrete $m$-current attributes. Several numerical experiments demonstrate the effectiveness of this approach.
Genome-scale metabolic networks are useful tools for achieving a system-level understanding of metabolism. However, due to their large size, analysis of such networks may be difficult and algorithms can be very slow. Therefore, some authors have suggested to analyze subsystems instead of the original genome-scale models.
Flux coupling analysis (FCA) is a well-known method for detecting functionally related reactions in metabolic networks.
In this paper, we study how flux coupling relations may change if we analyze a subsystem instead of the original network. We show mathematically that a pair of fully, partially or directionally coupled reactions may be detected as uncoupled in certain subsystems. Interestingly, this behavior is the opposite of the flux coupling changes that may occur due to missing reactions, or equivalently, deletion of reactions.
Computational experiments suggest that the analysis of plastid (but not mitochondrial) subsystems may significantly influence the results of FCA. Therefore, the results of FCA for subsystems, especially plastid subsystems, should be interpreted with care.
We formulate the static mechanical coupling of a geometrically exact Cosserat rod to an elastic continuum. The coupling conditions accommodate for the difference in dimension between the two models. Also, the Cosserat rod model incorporates director variables, which are not present in the elastic continuum model. Two alternative coupling conditions are proposed, which correspond to two different configuration trace spaces. For both we show existence of solutions of the coupled problems. We also derive the corresponding conditions for the dual variables and interpret them in mechanical terms.
We introduce nonsmooth Schur-Newton methods for the solution of the nonlinear discrete saddle-point problems arising from discretized vector-valued Cahn-Hilliard equations with logarithmic and obstacle potentials. The discrete problems are obtained by semi-implicit discretization in time and a first order finite element discretization in space. We incorporate the linear constraints that enforce solutions to stay on the Gibbs simplex using Lagrangian multipliers and prove existence of these multipliers under the assumption of a non-trivial initial condition for the order parameters.
Supercomputers can simulate complex molecular systems. However, there is a very large gap between the fastest oscillations of covalent bonds of a molecule and the time-scale of the dominant processes. In order to extract the dominant time-scales and to identify the dominant processes, a clustering of information is needed. This thesis shows that only the subspace-based Robust Perron Cluster Analysis (PCCA+) can solve this problem correctly by the construction of a Markov State Model. PCCA+ allows for time-extrapolation in molecular kinetics. This thesis shows the difference between molecular dynamics and molecular kinetics. Only in the molecular kinetics framework a definition of transition rates is possible. In this context, the existence of an infinitesimal generator of the dynamical processes is discussed. If the existence is assumed, the Theorem of Gauß can be applied in order to compute transition rates efficiently. Molecular dynamics, however, is not able to provide a suitable statistical basis for the determination of the transition pattern.
Markov state models have become very popular for the description of conformation dynamics of molecules over
long timescales. The construction of such models requires a partitioning of the configuration space such that the discretization
can serve as an approximation of metastable conformations. Since the computational complexity for the construction of a
Markov state model increases quadratically with the number of sets, it is desirable to obtain as few sets as necessary. In this
paper we propose an algorithm for the adaptive refinement of an initial coarse partitioning. A spectral clustering method is
applied to the final partitioning to detect the metastable conformations. We apply this method to the conformation analysis of
a model tri-peptide molecule, where metastable beta- and gamma-turn conformations can be identified.
In high accuracy numerical simulations and optimal control of time-dependent processes, often both many time steps and fine spatial discretizations are needed. Adjoint gradient computation, or post-processing of simulation results, requires the storage of the solution trajectories over the whole time, if necessary together with the adaptively refined spatial grids. In this paper we discuss various techniques to reduce the memory requirements, focusing first on the storage of the solution data, which typically are double precision floating point values. We highlight advantages and disadvantages of the different approaches. Moreover, we present an algorithm for the efficient storage of adaptively refined, hierarchic grids, and the integration with the compressed storage of solution data.
Coarsening in solder alloys is a widely accepted indicator for
possible failure of joints in electronic devices.
Based on the well-established Cahn--Larch\'e model
with logarithmic chemical energy density (see, e.g., Dreyer and Mueller 2001),
we present a numerical environment for the efficient
and reliable simulation of coarsening in binary alloys.
Main features are adaptive mesh refinement based on hierarchical error estimates,
fast and reliable algebraic solution by multigrid and Schur--Newton multigrid methods,
and the quantification of the coarsening speed by the temporal growth of mean phase radii.
We provide a detailed description and a numerical assessment of the algorithm and its different components,
together with a practical application to a eutectic AgCu brazing alloy.
Elementary flux mode (EM) analysis is an important tool in the constraint-based analysis of genome-scale metabolic networks. Due to the combinatorial complexity of these networks, as well as the advances in the level of detail to which they can be reconstructed, an exhaustive enumeration of all EMs is often not practical. Therefore, in recent years interest has shifted towards searching EMs with specific properties. We present a novel method that allows computing EMs containing a given set of target reactions. This generalizes previous algorithms where the set of target reactions consists of a single reaction. In the one-reaction case, our method compares favorably to the previous approaches. In addition, we present several applications of our algorithm for computing EMs containing two target reactions in genome-scale metabolic networks.
For the solution of optimal control problems governed by nonlinear parabolic PDEs, methods working on the reduced objective functional are often employed to avoid a full
spatio-temporal discretization of the problem. The evaluation of the reduced gradient requires one solve of the state equation forward in time, and one backward solve of
the adjoint equation. The state enters into the adjoint equation, requiring the storage of a full 4D data set. If Newton-CG methods are used, two additional trajectories
have to be stored. To get numerical results which are accurate enough, in many case very fine discretizations in time and space are necessary, which leads to a significant
amount of data to be stored and transmitted to mass storage. Lossy compression methods were developed to overcome the storage problem by reducing the accuracy of the stored
trajectories. The inexact data induces errors in the reduced gradient and reduced Hessian. In this paper, we analyze the influence of such a lossy trajectory compression
method on Newton-CG methods for optimal control of parabolic PDEs and design an adaptive strategy for choosing appropriate quantization tolerances.
Spectral deferred correction methods for solving stiff ODEs are known to
converge rapidly towards the collocation limit solution on equidistant grids,
but show a much less favourable contraction on non-equidistant grids such as
Radau-IIa points. We interprete SDC methods as fixed point iterations for the
collocation system and propose new DIRK-type sweeps for stiff problems based
on purely linear algebraic considerations. Good convergence is recovered also
on non-equidistant grids. The properties of different variants are explored on
a couple of numerical examples.
Constraint-based analysis of metabolic networks has become a widely used approach in computational systems biology.
In the simplest form, a metabolic network is represented by a stoichiometric matrix and thermodynamic information on the irreversibility of certain reactions.
Then one studies the set of all steady-state flux vectors satisfying these stoichiometric and thermodynamic constraints.
We introduce a new lattice-theoretic framework for the computational analysis of metabolic networks, which focuses on the support of the flux vectors,
i.e., we consider only the qualitative information whether or not a certain reaction is active, but not its specific flux rate.
Our lattice-theoretic view includes classical metabolic pathway analysis as a special case, but turns out to be much more flexible and general,
with a wide range of possible applications.
We show how important concepts from metabolic pathway analysis, such as blocked reactions, flux coupling, or elementary modes, can be generalized to arbitrary lattice-based models. We develop corresponding general algorithms and present a number of computational results.
The huge number of elementary flux modes (EFM) in genome-scale metabolic networks makes analysis based on elementary flux modes intrinsically difficult. However, it has been shown that the elementary flux modes with optimal yield often contain highly redundant information. The set of optimal-yield elementary flux modes can be compressed using modules. Up to now, this compression was only possible by first enumerating the whole set of all optimal-yield elementary flux modes.
We present a direct method for computing modules of the thermodynamically constrained optimal flux space of a metabolic network. This method can be used to decompose the set of optimal-yield elementary flux modes in a modular way and to speed up their computation. In addition, it provides a new form of coupling information that is not obtained by classical flux coupling analysis. We illustrate our approach on a set of model organisms.
The chemical master equation is a fundamental equation in chemical kinetics. It is an adequate substitute for the classical reaction-rate equations whenever stochastic effects become relevant. In the present paper we give a simple argument showing that the solutions of a large class of chemical master equations, including all those in which elementary reactions between two and more molecules do not generate a larger number of molecules than existed before,
are bounded in weighted $\ell_1$-spaces. As an illustration for the implications of this kind of regularity we analyze the effect of truncating the state space. This leads to an error analysis of the finite state projection of the chemical master equation, an approximation that underlies many numerical methods.
The long standing problem is discussed of how to deflate the part associated with the eigenvalue infinity in a structured matrix pencil using structure preserving unitary transformations. We derive such a deflation procedure and apply this new technique to symmetric, Hermitian or alternating pencils and in a modified form to (anti)-palindromic pencils. We present a detailed error and perturbation analysis of this and other deflation procedures and demonstrate the properties of the new algorithm with several numerical examples.
Flux coupling analysis (FCA) has become a useful tool for aiding metabolic reconstructions and guiding genetic manipulations. Originally, it was introduced for constraint-based models of metabolic networks that are based on the steady-state assumption. Recently, we have shown that the steady-state assumption can be replaced by a much weaker lattice-theoretic property related to the supports of metabolic fluxes. In this paper, we further extend our approach and delevelop an efficient algorithm for general qualitative flux coupling analysis (QFCA). We illustrate our method by thermodynamic flux coupling analysis (tFCA), which allows studying steady-state metabolic models with loop-law thermodynamic constraints. These models do not satisfy the lattice-theoretic properties required in our previous work. For a selection of genome-scale metabolic network reconstructions, we discuss both theoretically and practically, how thermodynamic constraints strengthen the coupling results that can be obtained with classical FCA.
Constraint-based modeling of genome-scale metabolic network reconstructions has become a widely used approach in computational biology. Flux coupling analysis is a constraint-based method that analyses the impact of single reaction knockouts on other reactions in the network.We present an extension of flux coupling analysis for double and multiple gene or reaction knockouts, and develop corresponding algorithms for an in silico simulation. To evaluate our method, we perform a full single and double knockout analysis on a selection of genome-scale metabolic network reconstructions and compare the results.
Asymptotic behavior is often of particular interest when analyzing asynchronous Boolean networks representing biological systems such as signal transduction or gene regulatory networks. Methods based on a generalization of the steady state notion, the so-called symbolic steady states, can be exploited to investigate attractor properties as well as for model reduction techniques conserving attractors. In this paper, we propose a novel optimization-based method for computing all maximal symbolic steady states and motivate their use. In particular, we add a new result yielding a lower bound for the number of cyclic attractors and illustrate the methods with a short study of a MAPK pathway model.
This paper presents efficient computational techniques for solving an optimization problem in cardiac defibrillation governed by the monodomain equations. Time-dependent electrical currents injected at different spatial positions act as the control. Inexact Newton-CG methods are used, with reduced gradient computation by adjoint solves. In order to reduce the computational complexity, adaptive mesh refinement for state and adjoint equations is performed. To reduce the high storage and bandwidth demand imposed by adjoint gradient and Hessian-vector evaluations, a lossy compression technique for storing trajectory data is applied. An adaptive choice of quantization tolerance based on error estimates is developed in order to ensure convergence. The efficiency of the proposed approach is demonstrated on numerical examples.