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)
Project
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.