## List of publications (with abstracts)

[> bottom of page (oldest)]

W.Auzinger, K.N.Burdeos, M.Fallahpour, O.Koch, R.G.Mendoza, E.Weinmüller:

**A numerical continuation method for parameter-dependent boundary value problems using**`bvpsuite 2.0`,

Journal of Numerical Analysis, Industrial and Applied Mathematics (JNAIAM) 16(1-2) 1-13, 2022.

__Abstract.__The Matlab package`bvpsuite 2.0`is a numerical collocation code for the approximation of solutions of a broad range of boundary value problems in ordinary differential equations. In this article, its newly implemented pathfollowing module with automated step-length control is presented. Two versions using the pseudo-arclength continuation method, allowing pathfollowing beyond turning points, were developed, both taking advantage of the existing features of`bvpsuite 2.0`such as error estimation and mesh adaptation. The first version is based on the Gauss-Newton method. The second version is now contained in the package`bvpsuite 2.0`and uses its built-in iterative method, the Fast Frozen Newton method. Their operating principles are presented and their performance is compared by means of the computation of some pathfollowing problems. Furthermore, the results of computations with`bvpsuite 2.0`for a problem with path bifurcations are presented.

W.Auzinger, J.Dubois, K.Held, H.Hofstätter, T.Jawecki, A.Kauch, O.Koch, K.Kropielnicka, P.Singh, C.Watzenböck:

**Efficient Magnus-type integrators for solar energy conversion in Hubbard models,**

Journal of Computational Mathematics and Data Science (JCMDS) 2 (2022) 100018. (DOI 10.1016/j.jcmds.2021.100018)

__Abstract.__Strongly interacting electrons in solids are generically described by Hubbard-type models, and the impact of solar light can be modeled by an additional time-dependence. This yields a finite dimensional system of ordinary differential equations (ODE)s of Schrödinger type, which can be solved numerically by exponential time integrators of Magnus type. The efficiency may be enhanced by combining these with operator splittings. We will discuss several different approaches of employing exponential-based methods in conjunction with an adaptive Lanczos method for the evaluation of matrix exponentials and compare their accuracy and efficiency. For each integrator, we use defect-based local error estimators to enable adaptive time-stepping. This serves to reliably control the approximation error and reduce the computational effort.

F.K.Auer, W.Auzinger, J.Burkotová, I.Rachůnková, E.Weinmüller:

**On nonlinear singular BVPs with nonsmooth data. Part 2: Convergence of collocation methods,**

Appl. Numer. Math. 171, 149-175, 2022. (DOI 10.1016/j.apnum.2021.08.016)

__Abstract.__We discuss numerical solution of boundary value problems for systems of nonlinear ordinary differential equations with a time singularity, \[ x'(t) = \frac{M(t)}{t}x(t)+\frac{f(t,x(t))}{t}, \quad t \in (0,1], \quad b(x(0),x(1)) = 0, \] where \( M\colon [0,1] \to {\mathbb R}^{n \times n} \) and \( f\colon [0,1] \times {\mathbb R}^n \to {\mathbb R}^n \) are continuous matrix-valued and vector-valued functions, respectively. Moreover, \( b\colon {\mathbb R}^n \times {\mathbb R}^n \to {\mathbb R}^n \) is a continuous nonlinear mapping which is specified according to a spectrum of the matrix \( M(0)\) to guarantee the BVP to be well-posed. For the case where \( M(0) \) has eigenvalues with nonzero real parts, we prove new convergence results for the collocation method and analytical results about the necessary smoothness of the solution to the problem required in the numerical analysis. We illustrate the theory by means of numerical examples.

W.Auzinger, I.Březinová, A.Grosz, H.Hofstätter, O.Koch, T.Sato:

**Efficient adaptive exponential time integrators for nonlinear Schrödinger equations with nonlocal potential,**

Journal of Computational Mathematics and Data Science (JCMDS) 1 (2021) 100014. (DOI 10.1016/j.jcmds.2021.100014)

__Abstract.__The performance of exponential-based numerical integrators for the time propagation of the equations associated with the multiconfiguration time-dependent Hartree-Fock (MCTDHF) method for the approximation of the multi-particle Schröodinger equation in one space dimension is assessed. Among the most popular integrators such as Runge-Kutta methods, time-splitting, exponential integrators and Lawson methods, exponential Lawson multistep methods with one predictor-corrector step provide the best stability and accuracy at the least effort. This assessment is based on the observation that the evaluation of the nonlocal terms associated with the potential is the computationally most demanding part of such a calculation in our setting. In addition, the predictor step provides an estimator for the local time-stepping error, thus allowing for adaptive time-stepping which reflects the smoothness of the solution and enables to reliably control the accuracy of a computation in a robust way, without the need to guess an optimal stepsize a priori. One-dimensional model examples are studied to compare different time integrators and demonstrate the successful application of our adaptive methods.

W.Auzinger, T.Jawecki, O.Koch, P.Pukach, R.Stolyarchuk, E.Weinmüller:

**Some Aspects on [numerical] Stability of Evolution Equations of Stiff Type; Use of Computer Algebra,**

in: Proceedings of the 2021 IEEE XVIIth International Conference on the Perspective Technologies and Methods in MEMS Design (MEMSTECH), 180-183, 2021. (ISSN 2573-5357)

__Abstract.__We are dealing with stiffness phenomena in initial value problems for systems of ordinary differential equations. Such problems appear in many physical applications, for instance, in the modeling of the evolution in time of mechanical systems or of electronic devices. A general theory of stiffness effects is not given here. Rather, several concrete examples are discussed in order to support understanding such effects. Use of computer algebra is also a relevant issue here.

W.Auzinger, K.Obelovska, I.Dronyuk, K.Pelekh, R.Stolyarchuk:

**A Continuous Model for States in CSMA/CA-based Wireless Local Networks,**

in: Proceedings of the 2nd International Conference on Data Science and Applications (ICDSA 2021). Lecture Notes in Networks and Systems, vol.287, 571-579, Springer, Singapore, 2021. (DOI 10.1007/978-981-16-5348-3_45)

__Abstract.__The effectiveness of wireless local area networks significantly depends on the method of access to the common physical environment. The main access method that is used today in wireless local area networks is Carrier Sense Multiple Access with Collision Avoidance (CSMA/CA), managing the efficient distribution of the common environment between active stations. In this paper as a topic of research we consider the state transition diagram of a CSMA/CA - based network, and we describe it by a system of differential equations. Analytical expressions for the probabilities of being in one of the possible states of a network station are derived.

P.Pukach, A.Slipchuk, W.Auzinger, R.Stolyarchuk, Y.Pukach, A.Kunynets, N.Pabyrivska:

**Asymptotic Method for Studying Mathematical Models of Resonant and Nonresonant Nonlinear Vibrations for Some 1D Moving Bodies,**

in: Proceedings of the 2021 IEEE 16th International Conference on the Experience of Designing and Application of CAD Systems (CADSM), Section 5: 'Models and Methods for Microelectronical Systems', 5/6-5/9, 2021. (DOI 10.1109/CADSM52681.2021.9385221)

__Abstract.__The mathematical model of nonlinear oscillations of an elastic moving one-dimensional body is investigated. The method of research of the specified model in case of resonant and nonresonant modes of fluctuations is developed. The oscillatory systems considered in the work simulate, in particular, dynamic processes in auger electromechanical machines, technological equipment used in drilling works, etc. The equations in standard form, which determine the parameters of the dynamic process for nonlinear oscillations of a moving body, are considered and investigated. Asymptotic approaches to the study of these mathematical models of oscillations make it possible to obtain and analyze the conditions of resonant and non-resonant modes of operation of the above technological systems.

W.Auzinger, H.Hofstätter, O.Koch, M.Quell:

**Adaptive time propagation for time-dependent Schrödinger equations,**

Int. J. Appl. Comput. Math. 7, 6.1-6.14, 2021. (DOI 10.1007/s40819-020-00937-9)

__Abstract.__We compare adaptive time integrators for the numerical solution of linear Schrödinger equations where the Hamiltonian explicitly depends on time. The approximation methods considered are splitting methods, where the time variable is split off and advanced separately, and commutator-free Magnus-type methods. The time-steps are chosen adaptively based on asymptotically correct estimators of the local error in both cases. It is found that splitting methods are more efficient when the Hamiltonian naturally suggests a separation into kinetic and potential part, whereas Magnus-type integrators excel when the structure of the problem only allows to advance the time variable separately.

C.Schattauer, L.Linhart, T.Fabian, T.Jawecki, W.Auzinger, F.Libisch:

**Graphene quantum dot states near defects,**

Phys. Rev. B 102, 155430-1−155430-9, 2020. (DOI 10.1103/PhysRevB.102.155430)

__Abstract.__Smooth confined graphene quantum dots (GQDs) localize Dirac electrons with conserved spin and valley degrees of freedom. Recent experimental realization of such structures using a combination of magnetic fields and a scanning tunneling microscope tip showcased their potential in locally probing and adjusting the valley degree of freedom. The present work models the influence of lattice defects on the level structures of GQDs. We study both the adiabatic level spacing "landscape" - orbital splitting and valley splitting - as well as transition dynamics between GQD-states. The system is modeled using a tight-binding approach with onsite- and hopping parameters in the vicinity of the defect region extracted from density functional theory via Wanner orbitals while time propagation is done using Magnus operators. Different defect types, such as double vacancy, Stone-Walse, flower and Si substitution are considered. We predict tunable valley splittings of the order of 2-20meV. The level structure can thus be tailored at will by engineering appropriate defect geometries.

W.Auzinger, A.A.Kytmanov, S.P.Tsarev:

**A Study of Anomalies in GPS Time Series via Polynomial Filtering,**

in: Proceedings of the IEEE 2020 XV'th International Scientific and Technical Conference on Computer Science and Information Technologies (CSIT 2020), Zbarazh-Lviv, Ukraine, September 23-26, 2020, 175-178.

__Abstract.__For the purpose of analyzing discrete time series, in particular originating from GPS or GLONASS coordinate measurements, we discuss an approach based on linear polynomial filtering. By expanding the given data set in terms of a (discrete) orthogonal polynomial basis it is possible to detect and analyze outliers from an expected smooth progression. Several aspects of this approach are discussed and illustrated via numerical experiments based on available measurements. We also consider a potential alternative approach based on (non-orthogonal) discrete Chebyshev polynomials.

W.Auzinger, K.Obelovska, R.Stolyarchuk:

**A Revised Gomory-Hu Algorithm Taking Account of Physical Unavailability of Network Channels,**

in: 'Computer Networks', P.Gaj, W.Gumiński, A.Kwiecień (Eds.), Springer International Publishing, (2020), 3-13. (DOI 10.1007/978-3-030-50719-0_1)

__Abstract.__The classical Gomory-Hu algorithm aims for finding, for given input flows, a network topology for data transmission and bandwidth of its channels which are optimized subject to minimal bandwidth criteria. In practice, it may occur that some channels between nodes of the network are not active. Ignoring such channels using the topology obtained by the Gomory-Hu algorithm will not lead to an optimal flow-rate. In this paper the focus is on a modified algorithm taking into account deficient channels. While the classical algorithm generates a sequence of ring subnets, in our modified version the use of deficient channels is checked at intermediate stages in each cycle of the algorithm. When forming ring subnets, the availability of new channels to be introduced into the ring subnet is checked and in the case of unavailability another ring closest to the optimal one is selected. The network optimized by this modified algorithm guarantees the transmission of the maximum input stream.

W.Auzinger, Y.Pelekh, M.Ihnatyshyn, R.Stolyarchuk, O.Kozar, S.Mentynskyi:

**The Studying of Hydrogen Diffusion Non-Stationary Processes Near a Crack in the Field of Heterogeneous Mechanical Tensions for the Encapsulated MEMS Devices,**

in: Proceedings of the IEEE 2020 XVI'th International Conference on the Perspective Technologies and Methods in MEMS Design (MEMSTECH 2020), Lviv, Ukraine, April 22-26, 2020, 164-168. (DOI 10.1109/MEMSTECH49584.2020.9109464)

__Abstract.__In this paper we study an elastic-plastic isotropic body weakened by a rectilinear crack, directed along the abscissa axis, under the action of stresses symmetrical about its plane. An approximate (analytical and numerical) solution of this problem is constructed under the condition that the hydrostatic deformation of the split expansion is approximated by a parabola. The concentration of hydrogen near the top of the crack is calculated.

Note that when using the above calculation formulas, at each node we get several approximations to the exact solution. In the proposed methods, only three appeals to the right-hand side of the differential equation are necessary to obtain both the third order accuracy method and the two approximations (upper and lower) of the second order of accuracy.

W.Auzinger, P.Pukach, R.Stolyarchuk, M.Vovk:

**Adaptive Numerics for Linear ODE Systems with Time-Dependent Data; Application in Photovoltaics,**

in: Proceedings of the IEEE 2020 XVI'th International Conference on the Perspective Technologies and Methods in MEMS Design (MEMSTECH 2020), Lviv, Ukraine, April 22-26, 2020, 5-8.

__Abstract.__We give a brief account on the numerical approximation of linear evolution systems with time-dependent data. In particular, exponential integrators of Magnus type are considered. As an application we consider a system modeling the effect of photon impact of electrons in photovoltaic cells, in form of a time-dependent system of ordinary differential equations (ODEs) of spatially discrete Schrödinger type. Adaptive tuning of time steps in course of the numerical integration is an essential point; this is based on recently developed computable local error estimators. Numerical results are presented.

W.Auzinger, K.Obelovska, R.Stolyarchuk:

**A Modified Gomory-Hu Algorithm with DWDM-Oriented Technology,**

in: Lecture Notes in Computer Science, LNCS 11958: 'Large-Scale Scientific Computing', I.Lirkov, S.Margenov (Eds.), 547-554, Springer, 2020. (DOI 10.1007/978-3-030-41032-2_63)

__Abstract.__Optimization of the topology of computer networks based on the classical Gomory-Hu algorithm does not take the specific transfer technology into account. For WDM technology requirements this leads to a redundancy of channel capacities. To reduce the redundancy of allocating network resources, we propose a modification of the Gomory-Hu algorithm which takes account of the specifics of DWDM technology - not at the final stage but already at intermediate stages in the process. The original algorithm proposed by Gomory and Hu involves the decomposition of the graph of the input network into ring subnets of different dimensions. Our modified algorithm takes account of the technical parameters of the DWDM technology for each ring during the decomposition. We illustrate our method by an example. The technique can be extended to large networks, which may lead to a significant economic effect.

W.Auzinger, A.Grosz, H.Hofstätter, O.Koch:

**Adaptive Exponential Integrators for MCTDHF,**

in: Lecture Notes in Computer Science, LNCS 11958: 'Large-Scale Scientific Computing', I.Lirkov, S.Margenov (Eds.), 557-565, Springer, 2020. (DOI 10.1007/978-3-030-41032-2_64)

__Abstract.__We compare exponential-type integrators for the numerical time-propagation of the equations of motion arising in the multiconfiguration time-dependent Hartree-Fock method for the approximation of the high-dimensional multi-particle Schrödinger equation. We find that among the most widely used integrators like Runge-Kutta, exponential splitting, exponential Runge-Kutta, exponential multistep and Lawson methods, exponential Lawson multistep methods with one predictor/corrector step provide optimal stability and accuracy at the least computational cost, taking into account that the evaluation of the nonlocal potential terms is by far the computationally most expensive part of such a calculation. Moreover, the predictor step provides an estimator for the time-stepping error at no additional cost, which enables adaptive time-stepping to reliably control the accuracy of a computation.

T.Jawecki, W.Auzinger, O.Koch:

**Computable upper error bounds for Krylov approximations to matrix exponentials and associated \(\varphi\)-functions,**

BIT Numer. Math. 60(1), 157-197, 2020. [open access] (DOI 10.1007/s10543-019-00771-6)

__Abstract.__An a posteriori bound for the error of a standard Krylov approximation to the exponential of nonexpansive matrices is derived. This applies for instance to skew-Hermitian matrices which appear in the numerical treatment of Schrödinger equations. It is proven that this defect-based bound is strict in contrast to existing approximations of the error, and it can be computed economically in the underlying Krylov space. In view of time-stepping applications, assuming that the given matrix is scaled by a time step, it is shown that the bound is asymptotically correct (with an order related to the dimension of the Krylov space) for the time step tending to zero. Furthermore, this result is extended to Krylov approximations of \( \varphi \)-functions and to improved versions of such approximations. The accuracy of the derived bounds is demonstrated by examples and compared with different variants known from the literature, which are also investigated more closely. Also other error bounds are tested on examples, in particular a version based on the concept of effective order. For the case where the matrix exponential is used in time integration algorithms, a step size selection strategy is proposed and illustrated by experiments.

W.Auzinger, M.Fallahphour, O.Koch, E.Weinmüller:

**Implementation of a pathfollowing strategy with an automatic step-length control: New**`MATLAB`package`bvpsuite 2.0`,

[ASC Report No. 31/2019 (pdf)]

__Abstract.__The main objective of this work was the implementation of a pathfollowing module to complete the`MATLAB`package`bvpsuite 2.0`. In preparation for the implementation, the tangent continuation method was studied, as presented in [7] and [4]. This method was adapted for the already existing code`bvpsuite 2.0`. Also the package was adapted where needed for the new function to fit in. Moreover, some features were implemented in the function, that would allow to tackle a broad range of problems. The implementation is tested on various examples. A secondary objective was the preparation of a manual, which would help first time users and more experienced users of`bvpsuite 2.0`, to use the package as intended. Finally, reports on simulations with`bvpsuite 2.0`that were carried out in collaboration with various researchers are given.

W.Auzinger, O.Koch, E.Weinmüller, S.Wurm:

**Modular version**`bvpsuite 1.2`of the collocation`MATLAB`package`bvpsuite 1.1`,

[ASC Report No. 30/2019 (pdf)]

__Abstract.__This work was devoted to the implementation of an updated version of the open domain`MATLAB`code`bvpsuite 1.1`. The goal was to understand the structure of the old code, identify the drawbacks, and implement missing features. An important task was to implement the new version in a possibly modular form to improve the readability and accessibility of the code. All modules were finally collected in a new software package`bvpsuite 2.0`. These modules were:- Solving BVPs in ODEs on finite and semi-infinite intervals;
- Solving EVPs in ODEs on finite and semi-infinite intervals;
- Solving Index-1 DAEs.

W.Auzinger, R.Stolyarchuk, R.Stoliarchuk:

**Strongly damped Stiff Oscillator under External Force: A Case Study,**

available online in: Proceedings of the 9-th International Youth Science Forum "Litteris et Artibus", Lviv Polytechnic National University, Ukraine; 22-28, 2019.

__Abstract.__This is a numerical study of ODE systems modelling the dynamics of a strongly damped stiff oscillator under periodic force. After reviewing the notion of stiffness, we demonstrate and discuss convergence order reduction for some well-established numerical integrators. We also include some comments on stiffness in nonlinear systems.

P.Pukach, V.Il'kiv, M.Vovk, O.Slyusarchuk, Y.Pukach, Y.Mylyan, W.Auzinger:

**On the mathematical model of nonlinear vibrations of a biologically active rod with consideration of the rheological factor,**

CEUR Workshop Proceedings, vol.2488, 30-42, 2019. (2nd International Workshop on Informatics and Data-Driven Medicine, IDDM 2019, Lviv, Ukraine).

__Abstract.__Qualitative and numerical methods of researching nonlinear vibration systems are used to study the mathematical model of nonlinear vibrations of a biologically active rod. This model is widely used in biomechanics and medical research for designing new materials with biofactor elements that possess certain preset features. Conditions are established for the existence of a unique solution of the boundary value problem for the beam vibration nonlinear differential equation, in which there is an integral summand with the fourth derivative by the spatial variables. This summand models the rheological factor in the system. The existence of classes of nonlinear rheological vibration systems with dissipation that have blow-up regimes is stated theoretically. The relation between nonlinearity indices in such regimes is obtained. The theoretical possibility of using the Runge-Kutta method for numerical solution of the corresponding boundary value problem is shown. The results are illustrated by a model example. The importance of the obtained theoretical assumptions for the practical modeling, analysis, and synthesis of parameters of technological vibration systems is shown.

Z.Nytrebych, V.Il'kiv, O.Malanchuk, W.Auzinger:

**Investigation of mathematical model of acoustic wave propagation through relax environment in ultrasound diagnostics problems,**

CEUR Workshop Proceedings, vol.2488, 280-289, 2019. (2nd International Workshop on Informatics and Data-Driven Medicine, IDDM 2019, Lviv, Ukraine).

__Abstract.__A mathematical model of the process of an acoustic wave propagation in a relax environment has been investigated. This mathematical model is widely used to describe and determine the basic parameters of the wave process in the problems of ultrasound diagnostics. The model is formulated in the form of the Cauchy problem for hyperbolic equation of third order with the initial data, which are analytical functions. The class of entire functions, which is the class of existence and uniqueness of the Cauchy problem solution for the partial differential equation, which describes this wave, is established. In the selected class of functions, the Cauchy problem solution is constructed using the differential-symbol method. Examples of solving problems with specific initial data are given. The obtained results and the indicated methodology allow us to determine the basic parameters of the process of acoustic wave propagation in the problems of ultrasound diagnostics.

I.Dronyiuk, O.Fedevych, R.Stolyarchuk, W.Auzinger:

**OMNET++ and Maple software environments for IT Bachelor studies,**

in: Proceedings of the 6th International Symposium on Emerging Inter-networks, Communication and Mobility, Procedia Computer Science, vol.155, 654-659, 2019.

[open access] (DOI 10.1016/j.procs.2019.08.093)

available online at www.sciencedirect.com

__Abstract.__This article deals with the methodology of computer network modeling for scientific and educational purposes. We choose tools which are useful, easy to understand and enable to model real life communication tasks for IT Bachelor students.

H.Hofstätter, W.Auzinger, O.Koch:

**An algorithm for computing coefficients of words in expressions involving exponentials and its application to the construction of exponential integrators,**

in 'Computer Algebra in Scientific Computing', Lecture Notes in Computer Science 11661, Springer-Verlag, 197-214, 2019. (DOI 978-3-030-26831-2_14)

__Abstract.__This paper discusses an efficient implementation of the generation of order conditions for the construction of exponential integrators like exponential splitting and Magnus-type methods in the computer algebra system Maple. At the core of this implementation is a new algorithm for the computation of coefficients of words in the formal expansion of the local error of the integrator. The underlying theoretical background including an analysis of the structure of the local error is briefly reviewed. As an application the coefficients of all 8th order self-adjoint commutator-free Magnus-type integrators involving the minimum number of 8 exponentials are computed.

W.Auzinger, H.Hofstätter, O.Koch:

**Non-existence of generalized splitting methods with positive coefficients of order higher than four,**

Appl. Math. Lett. 97, 48-52, 2019. (DOI 10.1016/j.aml.2019.05.017)

__Abstract.__We prove that generalized exponential splitting methods making explicit use of commutators of the vector fields are limited to order four when only real coefficients are admitted. This generalizes the restriction to order two for classical splitting methods with only positive coefficients.

W.Auzinger, H.Hofstätter, O.Koch, K.Kropielnicka, P.Singh:

**Time adaptive Zassenhaus splittings for the Schrödinger equation in semiclassical regime,**

Appl. Math. Comput. 362, 124550, 2019. (DOI 10.1016/j.amc.2019.06.064)

[preprint at arXiv.org]

__Abstract.__Time dependent Schrödinger equations with conservative force field commonly constitute a major challenge in the numerical approximation, especially when they are analysed in the semiclassical regime. Extremely high oscillations originate from the semiclassical parameter, and call for appropriate methods. We propose to employ a combination of asymptotic Zassenhaus splitting with time adaptivity. While the former turns the disadvantage of the semiclassical parameter into an advantage, leading to highly efficient methods with low error constants, the latter enables to choose an optimal time step and to speed up the calculations when the oscillations subside. We support the results with numerical examples.

W.Auzinger, I.Dronyuk, O.Fedevych, Y.Klishch, P.Lipinski:

**Information technology of traffic modeling using the Ateb-functions theory for improving service delivery in eHealth systems and networks,**

in: Proceedings of the 1st International Workshop on Informatics and Data-Driven Medicine (IDDM 2018, Lviv, Ukraine), 267-274, 2018.

__Abstract.__This work is devoted to description of developed information technology for e-Health network traffic modeling by using Ateb-functions theory. Prediction model is built based on hypothesis about periodic nature of network traffic pulsations which lead to using a nonlinear oscillating system with single degree-of-freedom under the conditions of small-scale disturbances that in result gives more qualitative solution. Obtained prediction of traffic intensity is used as one of indicators in calculating most optimal routes of packets in a various simulated e-Health network topologies. In addition, this work contains description of developed technology and mechanism of its application while adaptive management decision-making.

W.Auzinger, H.Hofstätter, O.Koch:

**Symmetrized local error estimators for time-reversible one-step methods in nonlinear evolution equations,**

J. Comput. Appl. Math. 356, 339-357, 2019. (DOI 10.1016/j.cam.2019.02.011)

[preprint at arXiv.org]

__Abstract.__Prior work on computable defect-based local error estimators for (linear) time-reversible integrators is extended to nonlinear and nonautonomous evolution equations. We prove that the asymptotic results from the linear case [W.Auzinger and O.Koch, An improved local error estimator for symmetric time-stepping schemes, Appl. Math. Lett.82(2018), 106--110] remain valid, i.e., the modified estimators yield an improved asymptotic order as the step size goes to zero. Typically, the computational effort is only slightly higher than for conventional defect-based estimators, and it may even be lower in some cases. We illustrate this by some examples and present numerical results for evolution equations of Schrödinger type, solved by either time-splitting or Magnus-type integrators.

W.Auzinger, H.Hofstätter, O.Koch, M.Quell, M.Thalhammer:

**A posteriori error estimation for Magnus-type integrators,**

ESAIM: M2AN 53(1), 197-218, 2019. (DOI 10.1051/m2an/2018050)

[preprint: ASC Report No. 01/2018 (pdf)]

__Abstract.__We study high-order Magnus-type exponential integrators for large systems of ordinary differential equations defined by a time-dependent skew-Hermitian matrix. We construct and analyze defect-based local error estimators as the basis for adaptive stepsize selection. The resulting procedures provide a posteriori information on the local error and hence enable the accurate, efficient, and reliable time integration of the model equations. The theoretical results are illustrated on two numerical examples.

W.Auzinger, I.Březinová, H.Hofstätter, O.Koch, M.Quell:

**Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part II: Comparisons of local error estimation and step-selection strategies for nonlinear Schrödinger and wave equations,**

Comput. Phys. Comm. 234, 55-71, 2019. (DOI 10.1016/j.cpc.2018.08.003)

[preprint: ASC Report No. 14/2017 (pdf)]

__Abstract.__We compare the practical performance of adaptive splitting methods for the solution of nonlinear Schrödinger equations. Different methods for local error estimation are assessed with respect to their accuracy and efficiency in conjunction with promising strategies for step-size adaptation. The numerical comparisons comprise the cubic nonlinear Schrödinger equation with a blow-up solution, systems of coupled nonlinear Schrödinger equations, a rotational and a Gross-Pitaevskii equation under an oscillatory potential inducing wave chaos, and a quantum control model with a time-dependent potential. Finally, for nonlinear wave equations we demonstrate the enhanced computational stability ensuing from adaptive step selection strategies close to the border mandated by the CFL condition.

W.Auzinger, O.Koch:

**An improved local error estimator for symmetric time-stepping schemes,**

Appl. Math. Lett. 82, 106-110, 2018. (DOI 10.1016/j.aml.2018.03.001)

[ Archived at u:scholar ]

__Abstract.__We propose a symmetrized version of the defect to be used in the estimation of the local time-stepping error of symmetric one-step methods for the time propagation of linear autonomous evolution equations. Using the anti-commutator of the numerical flow and the right-hand side operator in the definition of the defect of the numerical approximation, a local error estimator is obtained which has higher accuracy asymptotically than an established version using the common defect. This theoretical result is illustrated for a splitting method applied to a linear Schrödinger equation.

W.Auzinger, J.Burkotová, I.Rachůnková, V.Wenin:

**Shooting methods for state-dependent impulsive boundary value problems, with applications,**

Appl. Numer. Math. 128, 217-229, 2018. (DOI 10.1016/j.apnum.2018.02.006)

[preprint: ASC Report No. 02/2018 (pdf)]

__Abstract.__For impulsive boundary value problems whose solutions encounter discontinuities (jumps) at a priori not known positions depending on the solution itself, numerical methods have not been considered so far. We extend the well-known shooting approach to this case, combining Newton iteration with the numerical solution of impulsive initial value problems. We discuss conditions necessary for the procedure to be well-defined, and we present numerical results for several examples obtained with an experimental code realized in`MATLAB`.

W.Auzinger, W.Herfort, O.Koch, M.Thalhammer:

**The BCH-formula and order conditions for splitting methods,**

in "Lie Groups, Differential Equations, and Geometry", 71-83, Springer-Verlag, 1st ed. 2017. (ISBN: 978-3-319-62180-7)

__Abstract.__As an application of the BCH-formula, order conditions for splitting schemes are derived. The same conditions can be obtained by using non-commutative power series techniques and inspecting the coefficients of Lyndon-Shirshov words.

W.Auzinger, T.Kassebacher, O.Koch, M.Thalhammer:

**Convergence of a Strang splitting finite element discretization for the Schrödinger-Poisson equation,**

ESAIM: M2AN 51(4), 1245-1278, 2017. (DOI 10.1051/m2an/2016059)

[preprint: ASC Report No. 40/2015 (pdf)]

__Abstract.__Operator splitting methods combined with finite element spatial discretizations are studied for time-dependent nonlinear Schrödinger equations. In particular, the Schrödinger-Poisson equation under homogeneous Dirichlet boundary conditions on a finite domain is considered. A rigorous stability and error analysis is carried out for the second-order Strang splitting method and conforming polynomial finite element discretizations. For sufficiently regular solutions the classical orders of convergence are retained, that is, second-order convergence in time and polynomial convergence in space is proven. The established convergence result is confirmed and complemented by numerical illustrations.

W.Auzinger, O.Koch, M.Quell:

**Adaptive high-order splitting methods for systems of nonlinear evolution equations with periodic boundary conditions,**

Numer. Algor. 75(1), 261-283, 2017. [open access] (DOI 10.1007/s11075-016-0206-8)

__Abstract.__We assess the applicability and efficiency of time-adaptive high-order splitting methods applied for the numerical solution of (systems of) nonlinear parabolic problems under periodic boundary conditions. We discuss in particular several applications generating intricate patterns and displaying nonsmooth solution dynamics. First we give a general error analysis for splitting methods for parabolic problems under periodic boundary conditions and derive the necessary smoothness requirements on the exact solution in particular for the Gray-Scott equation and the Van der Pol equation. Numerical examples demonstrate the convergence of the methods and serve to compare the efficiency of different time-adaptive splitting schemes and of splitting into either two or three operators, based on appropriately constructed a posteriori local error estimators.

W.Auzinger, H.Hofstätter, D.Ketcheson, O.Koch:

**Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes,**

BIT Numer. Math. 57(1), 55-74, 2017. [open access] (DOI 10.1007/s10543-016-0626-9)

__Abstract.__We present a number of new contributions to the topic of constructing efficient higher-order splitting methods for the numerical integration of evolution equations. Particular schemes are constructed via setup and solution of polynomial systems for the splitting coefficients. To this end we use and modify a recent approach for generating these systems for a large class of splittings. In particular, various types of pairs of schemes intended for use in adaptive integrators are constructed.

W.Auzinger, H. Hofstätter, O.Koch, M.Thalhammer:

**Reduced order of the local error of splitting for parabolic problems,**

in `Numerical Analysis and Applied Mathematics', AIP Conference Proceedings 1863, 160008 (2017). (DOI 10.1063/1.4992342)

__Abstract.__We give a theoretical analysis of the local error of splitting methods applied to parabolic initial-boundary value problems under homogeneous Dirichlet or Neumann boundary conditions. For the Lie-splitting, we provide a theoretical local error analysis that rigorously explains the order reduction observed in numerical experiments.

W.Auzinger, R.Stolyarchuk, M.Tutz:

**Defect correction methods, classic and new,**

Visnyk of the Lviv Univ.Series Mech.Math. 82, 5-19, 2016. [open access]

__Abstract.__Defect correction methods are based on the idea of measuring the quality of an approximate solution to an operator equation by forming the defect, or residual, with respect to the given problem. By an appropriate backsolving procedure, an error estimate is obtained. This process can also be continued in an iterative fashion. One purpose of this overview is the further dissemination of the underlying concepts. Therefore, we first we give a general and consistent review on various types defect correction methods, and its application in the context of discretization schemes for dfferential equations. After describing the general algorithmic templates we discuss some specific techniques used in the solution of ordinary differential equations. Moreover, new results about the application to implicit problems are presented.

W.Auzinger, R.Stolyarchuk:

**A topic in the stability analysis associated with companion matrices,**

Matematychni Studii 46(2), 115-120, 2016. [open access] (DOI 10.15330/ms.46.2.115-120)

__Abstract.__Assume that a family of (non-normal) matrices has stable spectra contained in the complex unit circle. This does not necessarily imply that the 2-norm of such matrices is small, and the question is under what `natural' similarity transformation the transformed matrices will have 2-norm smaller or equal to 1. Already for the 2x2-case this is a nontrivial question, involving the analysis of a function in two complex variables (the eigenvalues) and a positive scaling parameter. We discuss and explain an approach to this problem which has been used before in the analysis of companion matrices. For the 3D case we present a numerical example.

W.Auzinger, W.Herfort, H.Hofstätter, O.Koch:

**Setup of order conditions for splitting methods,**

in 'Computer Algebra in Scientific Computing', Lecture Notes in Computer Science 9890, Springer-Verlag, 2016, 30-42. (DOI 10.1007/978-3-319-45641-6)

[preprint at arXiv.org]

__Abstract.__For operator splitting methods, an approach based on Taylor expansion and the particular structure of its leading term as an element of a free Lie algebra is used for the setup of a system of order conditions. Along with a brief review of the underlying theoretical background, we discuss the implementation of the resulting algorithm in computer algebra, in particular using Maple 18. A parallel version of such a code is described, and its performance on a computational node with 16 threads is documented.

W.Auzinger, H.Hofstätter, O.Koch:

**Symbolic manipulation of flows of nonlinear evolution equations, with application in the analysis of split-step time integrators,**

in 'Computer Algebra in Scientific Computing', Lecture Notes in Computer Science 9890, Springer-Verlag, 2016, 43-57. (DOI 10.1007/978-3-319-45641-6)

[preprint at arXiv.org]

__Abstract.__We describe a package realized in the Julia programming language which performs symbolic manipulations applied to nonlinear evolution equations, their flows, and commutators of such objects. This tool was employed to perform contrived computations arising in the analysis of the local error of operator splitting methods. It enabled the proof of the convergence of the basic method and of the asymptotical correctness of a defect-based error estimator. The performance of our package is illustrated on several examples.

W.Auzinger, T.Kassebacher, O.Koch, M.Thalhammer:

**Adaptive splitting methods for nonlinear Schrödinger equations in the semiclassical regime,**

Numer. Algor. 72(1), 1-35, 2016. (DOI 10.1007/s11075-015-0032-4)

[ Archived at u:scholar ]

__Abstract.__The convergence of exponential splitting methods for nonlinear Schrödinger equations in the semiclassical regime is studied. For the Lie and Strang splitting methods, the exact form of the local error is determined. and the dependence on the semiclassical parameter is identified. This is enabled within a defect-based framework which also suggests asymptotically correct a posteriori error estimators as the basis for adaptive time stepping. Numerical examples, also including higher-order schemes, confirm the theoretical results.

W.Auzinger:

**A note on similarity to contraction for stable 2x2 companion matrices,**

Ukr. Math. J. 68(3), 448-457, 2016.

[Ukr. Mat. Zh. 68(3), 400-407, 2016.] (DOI 10.1007/s11253-016-1233-2)

[preprint: ASC Report No. 01/2013 (pdf)]

__Abstract.__We consider companion matrices of size 2x2 with general complex spectra satisfying a root condition with respect to the closed complex unit circle or the closed left complex half plane. For both cases, smooth and naturally conditioned basis transformations are constructed such that the resulting, transformed matrix is contractive or dissipative, respectively.

W.Auzinger:

**Defect Correction Methods,**

p. 323-332, in Encyclopedia of Applied and Computational Mathematics, Volume 1, Enquist, B. (Ed.), Springer, 2015. (ISBN: 978-3-540-70528-4)

W.Auzinger, O.Koch, M.Schöbinger, E.Weinmüller:

**A new version of the code**`bvpsuite`for singular BVPs in ODEs: Nonlinear solver and its application to \( m \)-Laplacians,

[ASC Report No. 28/2015 (pdf)]

__Abstract.__This report documents the`MATLAB`routine`bvpsuite`, developed at the Technical University of Vienna. It explains the ways`bvpsuite`solves nonlinear boundary value problems. These problems may be singular or posed on a semiinfinite interval.`bvpsuite`can also solve eigenvalue problems and problems containing unknown parameters. Before the`bvpsuite`-specific details an overview of Collocation and the Newton method is given.

This report also shows the changes applied to the routine during a resent update. The use of`bvpsuite`is illustrated by two fully documented examples, a linear and a nonlinear one. An overview of the settings options is also given.

The report ends with the application of`bvpsuite`to \(m\)-Laplacians, both before and after a smoothing transformation. The numerical rates of convergence are evaluated depending on the number of collocation points per collocation interval.

W.Auzinger, O.Koch, M.Quell:

**Splittingverfahren für die Gray-Scott-Gleichung,**

[ASC Report No. 07/2015 (pdf)]

W.Auzinger, W.Herfort:

**An application of Gröbner bases to perturbed polynomial equations,**

[ASC Report No. 01/2015 (pdf)]

__Abstract.__The interplay of symbolic and numerical calculations in polynomial algebra is featured in the seminal monograph by H.J.Stetter (2004). We present a contribution to this field in form of an ideal-theoretic approach for determining solutions of an a perturbed polynomial system which are close to the zeros of the system subject to perturbation. Algorithmically, identifying these zeros is of special interest if the unperturbed system has a nontrivial solution manifold while the perturbed system admits only isolated zeros.

W.Auzinger, O.Koch, M.Thalhammer:

**Representation of the local error for higher-order exponential splitting schemes involving two or three sub-operators,**

in `Numerical Analysis and Applied Mathematics', AIP Conference Proceedings 1648, 150003 (2015). (DOI 10.1063/1.4912433)

[preprint (pdf)]

__Abstract.__We consider the numerical treatment of abstract evolution equations \[ \partial_t u(t) = H u(t) = A u(t) + B u(t)\;[\,+\,C u(t)\,]\,, \quad u(0)~\text{given}\,, \] by higher-order exponential splitting schemes. The main focus is on linear problems. A single step of an exponential splitting scheme with stepsize \( t \) and \( s \) sub-steps comprises a multiplicative composition of sub-flows of the form \( \mathcal{S}_j(t) = [\mathrm{e}^{t c_j C}]\,\,\mathrm{e}^{t b_j B}\,\mathrm{e}^{t a_j A}\,,\;j=1 \ldots s \). We present an algebraic theory of the structure of the local error. The leading term is a linear combination of iterated commutators of the sub-operators \( A,\,B \) [and \( C \)] involved. This fact can be exploited for the automatic setup of order conditions, i.e., systems of polynomial equations for the coefficients \( a_j,\,b_j \) [and \( c_j \)] which have to be satisfied for a desired order \( p \). In view of application to partial differential equations, an explicit, exact representation of the local error is of interest. This can be obtained by performing a multiple variation-of-constants integral expansion involving higher-order defect terms. The latter satisfy a multinomial expansion, and the building blocks in this expansion are determined via yet another recursively defined integral representation. We describe the rich combinatorial structure of this local error expansion which is influenced by iterated commutators of the given sub-operators. A defect-based a posteriori local error estimator is also proposed.

W.Auzinger, O.Koch, M.Thalhammer:

**Defect-based local error estimators for high-order splitting methods involving three linear operators,**

Numer. Algor. 70(1), 61-91, 2015. (DOI 10.1007/s11075-014-9935-8)

[ Archived at u:scholar ]

__Abstract.__Prior work on high-order exponential operator splitting methods is extended to evolution equations defined by three linear operators. A posteriori local error estimators are constructed via a suitable integral representation of the local error involving the defect associated with the splitting solution and quadrature approximation via Hermite interpolation. In order to prove asymptotical correctness, a multiple integral representation involving iterated defects is deduced by repeated application of the variation-of-constant formula. The error analysis within the framework of abstract evolution equations provides the basis for concrete applications. Numerical examples for initial-boundary value problems of Schrödinger and of parabolic type confirm the asymptotical correctness of the proposed a posteriori error estimators.

W.Auzinger, H.Hofstätter, O.Koch, M.Thalhammer:

**Defect-based local error estimators for splitting methods, with application to Schrödinger equations, Part III: The nonlinear case,**

J. Comput. Appl. Math. 273, 182-204, 2015. (DOI 10.1016/j.cam.2014.06.012)

[ Archived at u:scholar ]

__Abstract.__The present work is concerned with the efficient time integration of nonlinear evolution equations by exponential operator splitting methods. Defect-based local error estimators serving as a reliable basis for adaptive stepsize control are constructed and analyzed. In the context of time-dependent nonlinear Schrödinger equations, asymptotical correctness of the local error estimators associated with the first-order Lie-Trotter and second-order Strang splitting methods is proven. Numerical examples confirm the theoretical results and illustrate the performance of adaptive stepsize control.

W.Auzinger, R.Stolyarchuk, M.Tutz:

**Методи корекції дефекту, класичні та нові**(*Defect correction methods, classic and new*, in Ukrainian),

[ASC Report No. 26/2014 (pdf)]

__Abstract.__Методи корекції дефекту базуються на ідеї виміру якості наближеного розв'зку за допомогою формування дефекту, або залишку, по відношенню до даної задачі. За допомогою процедури розв'язку назад отримується оцінка похибки Цей процес може бути продовжений ітеративним чином. Метою цього огляду є подальше поширення концепції що розглядається. Більш того, вперше дано загальний і узгоджений огляд різноманітних типів методів корекції дефекту і їх застосувань в контексті дискритизаційних схем для диференціальних рівнянь. Після опису загального алгоритмічного шаблону ми обговоримо деякі спеціальні технології, що використовуються при в розв'язку звичайних диференціальних рівнянь. Також представлені нові результати стосовно застосування до неявних задач.

W.Auzinger, W.Herfort:

**Local error structures and order conditions in terms of Lie elements for exponential splitting schemes,**

Opuscula Mathematica 34(2), 243-255, 2014. [open access] (DOI 10.7494/OpMath.2014.34.2.243)

__Abstract.__We discuss the structure of the local error of exponential operator splitting methods. In particular, it is shown that the leading error term is a Lie element, i.e., a linear combination of higher-degree commutators of the given operators. This structural assertion can be used to formulate a simple algorithm for the automatic generation of a minimal set of polynomial equations representing the order conditions, for the general case as well as in symmetric settings.

W.Auzinger, O.Koch, A.Saboor Bagherzadeh:

**Error estimation based on locally weighted defect for boundary value problems in second order ordinary differential equations,**

BIT Numer. Math. 54, 873-900, 2014. (DOI 10.1007/s10543-014-0488-y)

[ Archived at reposiTUm ]

__Abstract.__We investigate efficient asymptotically correct a posteriori error estimates for the numerical approximation of two-point boundary value problems for second order ordinary differential equations by polynomial collocation methods. Our error indicators are based on the defect of the collocation solution with respect to an associated exact difference scheme and backsolving using a cheap, low order finite-difference scheme. We prove high asymptotical correctness of this error indicator and illustrate the theoretical analysis by numerical examples.

W.Auzinger, O.Koch, M.Thalhammer:

**Defect-based local error estimators for splitting methods, with application to Schrödinger equations, Part II: Higher-order methods for linear problems,**

J. Comput. Appl. Math. 255, 384-403, 2014. (DOI 10.1016/j.cam.2013.04.043)

[ Archived at u:scholar ]

__Abstract.__In this work, defect-based local error estimators for higher-order exponential operator splitting methods are constructed and analyzed in the context of time-dependent linear Schrödinger equations. The technically involved procedure is carried out in detail for a general three-stage third-order splitting method and then extended to the higher-order case. Asymptotical correctness of the a posteriori local error estimator is proven under natural commutator bounds for the involved operators, and along the way the known (non)stiff order conditions and a priori convergence bounds are recovered. The theoretical error estimates for higher-order splitting methods are confirmed by numerical examples for a test problem of Schrödinger type. Further numerical experiments for a test problem of parabolic type complement the investigations.

I.Reichl, W.Auzinger:

**Identifying Tibio-Femoral Joint Kinematics: Individual Adjustment versus Numerical Robustness,**

IFAC-PapersOnLine, Mathematical Modelling, Volume #7, Part #1, 819-824, 2013. (DOI 10.3182/20120215-3-AT-3016.00145)

[preprint (pdf)]

__Abstract.__The interpretation of joint kinematics data in terms of displacements is sensitive to the type of movement, the measurement technique and the reference axes where rotations and translations are associated with. Misaligned knee joint axes do not only lead to a misinterpretation of knee joint kinematics but, additionally, have a general impact on lower body kinematics. Therefore, several groups independently propose to a posteriori replace empirically palpated axes by functionally calculated axes. Flexion is the most dominant motion of the knee joint. Thus, the flexion axis has the largest effect on the kinematical calculation. The large angular flexion range facilitates the mathematical procedure of axis determination. The second rotation of importance is the internal/external rotation which is sometimes controlled by mathematical optimization techniques. This contribution focuses on the evaluation of the concepts of symmetrical axis of rotation approach (SARA) and finite helical axis (FHA) regarding their applicability in axis reorientation procedures, and to explore which of the two underlying algorithms performs most robust and convenient under typical data perturbations.

W.Auzinger, M.Tutz:

**A review of stability and error theory for collocation methods applied to linear boundary value problems,**

[ASC Report No.35/2012 (pdf)]

__Abstract.__We present a systematic, didactically remastered revision of the convergence proof. The presentation is (only slightly) informal in some details, with emphasis on the structure of the convergence argument.

I.Reichl, W.Auzinger:

**Parameter Identification In Functional Axis Calculation For The Human Knee Joint,**

in Proceedings of the 10th International Symposium on Computer Methods in Biomechanics and Biomedical Engineering (CMBBE-2012), Berlin, Germany, 2012. (ISBN: 978-0-9562121-5-3)

[preprint (pdf)]

__Abstract.__The interpretation of joint kinematics data in terms of displacements is a product of the type of movement, the measurement technique and the underlying model of the joint implemented in parameter optimization procedures. Misaligned knee joint axes do not only lead to a misinterpretation of knee joint kinematics but, additionally, have a general impact on lower body kinematics. Therefore, several groups independently propose different methods for a posteriori replacing the empirical flexion axis (determined from the original gait protocol) by a functionally calculated axis. Among the three main different rotations of the knee joint, flexion is the most dominant motion. Thus, if not determined accurately, it significantly affects any kinematical calculation. On the other hand, the large angular flexion range facilitates the mathematical procedure of axis determination. The second rotation of importance is the internal/external rotation which is sometimes controlled by mathematical optimization techniques. Consecutively to the axis reorientation process, joint displacements are computed via conventional procedures where knee rotation axes and consequently ankle and hip description have been modified. This contribution focuses on the evaluation of methods for the determination of rotation axes via parameter identification. It is the purpose to elaborate which of the available mathematical algorithms performs most robust and convenient in selected clinical settings. Quality parameters for the calculated flexion axis are proposed.

L.Exl, W.Auzinger, S.Bance, M.Gusenbauer, F.Reichel, T.Schrefl:

**Fast stray field computation on tensor grids,**

J. Comput. Phys. 231(7), 2840-2850, 2012. (DOI 10.1016/j.jcp.2011.12.030)

[preprint at arXiv.org)]

__Abstract.__A direct integration algorithm is described to compute the magnetostatic field and energy for given magnetization distributions on not necessarily uniform tensor grids. We use an analytically-based tensor approximation approach for function-related tensors, which reduces calculations to multilinear algebra operations. The algorithm scales with \( N^{4/3} \) for \( N \) computational cells used and with \( N^{2/3} \) (sublinear) when magnetization is given in canonical tensor format. In the final section we confirm our theoretical results concerning computing times and accuracy by means of numerical examples.

W.Auzinger, O.Koch, M.Thalhammer:

**Defect-based local error estimators for splitting methods, with application to Schrödinger equations, Part I: The linear case,**

J. Comput. Appl. Math. 236(10), 2643-2659, 2012. (DOI 10.1016/j.cam.2012.01.001)

[ Archived at u:scholar ]

__Abstract.__We introduce a defect correction principle for exponential operator splitting methods applied to time-dependent linear Schrödinger equations and construct a posteriori local error estimators for the Lie-Trotter and Strang splitting methods. Under natural commutator bounds on the involved operators we prove asymptotical correctness of the local error estimators, and along the way recover the known a priori convergence bounds. Numerical examples illustrate the theoretical local and global error estimates.

W.Auzinger, M.Łapińska:

**Convergence of rational multistep methods of Adams-Padé type,**

BIT Numer. Math. 52(1), 3-20, 2012. (DOI 10.1007/s10543-011-0353-1)

[preprint: ASC Report No. 05/2011 (pdf)]

__Abstract.__Rational generalizations of multistep schemes, where the linear stiff part of a given problem is treated by an A-stable rational approximation, have been proposed by several authors, but a reasonable convergence analysis for stiff problems has not been provided so far. In this paper we directly relate this approach to exponential multistep methods, a subclass of the increasingly popular class of exponential integrators. This natural, but new interpretation of rational multistep methods enables us to prove a convergence result of the same quality as for the exponential version. In particular, we consider schemes of rational Adams type based on A-acceptable Padé approximations to the matrix exponential. A numerical example is also provided.

W.Auzinger, M.Łapińska:

**Rational multistep methods via modified \( \varphi_j \)-functions,**

Proc. Appl. Math. Mech. 11, 757-758, 2011. (DOI 10.1002/pamm.201110368)

[preprint (pdf)]

__Abstract.__A class of rational multistep methods, in particular of Adams-Padé type, is designed via rational modification of the \( \varphi_j \)-functions inherent in exponential integrators. Convergence properties and implementation issues are discussed.

W.Auzinger:

**Error estimation via defect computation and reconstruction: Some particular techniques,**

Journal of Numerical Analysis, Industrial and Applied Mathematics (JNAIAM) 6(1-2), 15-27, 2011.

__Abstract.__The well-known technique of defect correction can be used in various ways for estimating local or global errors of numerical approximations to differential or integral equations. In this paper we describe the general principle in the context of linear and nonlinear problems and indicate the interplay between the auxiliary scheme involved and a correct definition of the defect. Applications discussed include collocation approximations to first and second order boundary value problems for nonlinear ODEs and, in particular, exponential splitting approximations for linear evolution equations. We describe the design of error estimators and their essential properties and give numerical examples. The theoretical tools for the analysis of the asymptotical correctness of such estimators are described, and references to original research papers are given where the complete analysis is provided.

I.Reichl, W.Auzinger, H.-B.Schmiedmayer, E.Weinmüller:

**Model dependence of knee joint kinematics,**

Proceedings of the XXIIIth International Society of Biomechanics Congress (ISB 2011), Brussels, Belgium, 2011. (ISSN: 1476-8259)

[preprint (pdf)]

__Abstract.__Two optimization models for the analysis of kinematics data of the human tibio-femoral joint are presented. It is assumed that the intact joint can be approximated by two rotational axes, only. Model 1 allows for an arbitrary angle between the two axes, model 2 demands an orthogonal constraint. The two alternative optimizations were applied to kinematics data due to a computer simulation. It turned out that model 2 was more robust with respect to noise and is, thus, to be preferred in case of noise typically present in motion capturing with skinmounted markers. Model 1 may possibly lead to more accurate results for medical imaging data.

I.Reichl, W.Auzinger, H.-B.Schmiedmayer, E.Weinmüller:

**Knee joint kinematics: comparison of two optimization models with respect to data noise,**

Proc. Appl. Math. Mech. 11, 111-112, 2011. (DOI 10.1002/pamm.201110047)

__Abstract.__The interpretation of human knee joint kinematics in terms of displacements is an outcome of the underlying model of the joint and the measurement technique. Measurement errors and noise challenge the development of optimization procedures which, based on a reduction in degrees of freedom, aim for the reproducibility of joint displacements by computational techniques. So far, optimization algorithms have been applied which are based on a kinematic model of the healthy human tibio-femoral joint (TFJ) as a compound hinge with two fixed orthogonal axes. On the other hand, empirical studies find non-orthogonal rotational axes. Therefore it was the aim of the present study to investigate the implications of a refined kinematic model on the accuracy of computed joint rotation angles. For the purpose of quantitative comparison, kinematic data of a TFJ with two axes intersecting at an arbitrary angle were simulated. The joint rotations were optimized for the assumption of (a) two orthogonal and intersecting axes (model A), and (b) two axes intersecting at an arbitrary angle (model B). Model B recovers the original input data closer in case of low noise level as encountered in invasive measurement techniques. Skin mounted markers tracking involves non-normally distributed noise which is typically larger by on order of magnitude. In this case, model A exhibits a more favorable performance. These observations motivate the search for alternative kinematic descriptions of the TFJ.

I.Reichl, W.Auzinger, H.-B.Schmiedmayer, E.Weinmüller:

**A study of constrained models for the kinematic analysis of the human knee joint,**

Computer Methods in Biomechanics and Biomedical Engineering, 14(S1), 77-79, 2011. (DOI 10.1080/10255842.2011.592370)

[preprint (pdf)]

__Abstract.__Keywords: tibio-femoral joint; joint rotations; optimisation; error compensation

W.Auzinger, H.Lehner, E.Weinmüller:

**An efficient asymptotically correct error estimator for collocation solutions to singular index-1 DAEs,**

BIT Numer. Math. 51(1), 43-65, 2011. (DOI 10.1007/s10543-011-0321-9)

[preprint: ASC Report No. 18/2010 (pdf)]

__Abstract.__A computationally efficient a posteriori error estimator is introduced and analyzed for collocation solutions to linear index-1 DAEs (differential-algebraic equations) with properly stated leading term exhibiting a singularity of the first kind. The procedure is based on a modified defect correction principle, extending an established technique from the context of ordinary differential equations to the differential-algebraic case. Using recent convergence results for stiffly accurate collocation methods, we prove that the resulting error estimate is asymptotically correct. Numerical examples demonstrate the performance of this approach. To keep the presentation reasonably self-contained, some arguments from the literature on DAEs concerning the decoupling of the problem and its discretization, which is essential for our analysis, are also briefly reviewed. The appendix contains a remark about the interrelation between collocation and implicit Runge-Kutta methods for the DAE case.

I.Reichl, W.Auzinger, H.-B.Schmiedmayer, E.Weinmüller:

**Reconstructing the knee joint mechanism from kinematic data,**

Mathematical and Computer Modelling of Dynamical Systems 16(5), 403-415, 2010. [open access] (DOI 10.1080/13873954.2010.507094)

[preprint: ASC Report No. 48/2009 (pdf)]

__Abstract.__The interpretation of joint kinematics data in terms of displacements is a product of the type of movement, the measurement technique, and the underlying model of the joint implemented in optimization procedures. Kinematic constraints reducing the number of degrees of freedom are expected to compensate for measurement errors and noise, thus, increasing the reproducibility of joint angles. One approach already successfully applied by several groups approximates the healthy human knee joint as a compound hinge joint with minimal varus/valgus rotation. Most of these optimizations involve an orthogonality constraint. This contribution compares the effect of a model with and without orthogonality constraint on the obtained joint rotation angles. For this purpose kinematic data is simulated without noise and with normally distributed noise of varying size. For small noise the unconstrained model provides more accurate results while for larger noise this is the case for the constrained model. This can be attributed to the shape of the objective function of the unconstrained model near its minimum.

W.Auzinger, F.Kramer:

**On the stability and error structure of BDF schemes applied to linear parabolic evolution equations,**

BIT Numer. Math. 50(3), 455-480, 2010. (DOI 10.1007/s10543-010-0269-1)

[preprint: ASC Report No.29/2009 (pdf)]

__Abstract.__We continue the work of various authors on the stability and error structure of multistep schemes applied to linear evolution equations. BDF schemes are considered, and, as far as reasonable, explicit expressions for all occurring bounds are specified, exploiting prior work on the location of characteristic roots. The 2-step BDF scheme is considered in particular detail, and for problems of sectorial type, an asymptotic error expansion is derived based on damping properties of the scheme.

W.Auzinger:

**Normal forms for companion matrices and contractivity in inner product norms,**

[ASC Report No.24/2009 (pdf)]

__Abstract.__We study the problem of finding a inner product norm in which a given companion matrix with a [weakly] stable spectrum becomes contractive (or dissipative), via a preferably well-conditioned change of basis. To this end we use a basis transformation related to a rescaled LQ decomposition of the associated Vandermonde matrix which is robust to w.r.t. confluent or non-confluent spectra. For n=2 we give an explicit construction. The transformed, contractive matrix is non-normal in general, and it depends on the distribution of the spectrum in a nonlinear way. This analysis cannot be directly generalized to higher dimension, but it suggests an algebraic/numerical algorithm for a numerically given spectrum. This has been tested for small values of \( n \) and appears to be successful.

W.Auzinger, H.Lehner, E.Weinmüller:

**Defect-based a posteriori error estimation for index-1 DAEs with a singularity of the first kind,**

in `Numerical Analysis and Applied Mathematics', AIP Conference Proceedings 1168, vol.2, 1009-1012. (DOI 10.1063/1.3241219)

[preprint (pdf)]

__Abstract.__We show how the QDeC estimator, an efficient and asymptotically correct a posteriori estimator for collocation solutions to ODEs can be extended to differential algebraic equations (DAEs), where a singularity of the first kind is admitted.

W.Auzinger, H.Lehner, E.Weinmüller:

**Defect-based a-posteriori error estimation for differential-algebraic equations,**

PAMM 7, 1023101-1023102, 2007.

[preprint (pdf)]

__Abstract.__We show how the QDeC estimator, an efficient and asymptotically correct a-posteriori error estimator for collocation solutions to ODE systems, can be extended to differential-algebraic equations of index 1.

W.Auzinger, H.Lehner, E.Weinmüller:

**Defect-based a posteriori error estimation for index-1 DAEs,**

[ASC Report No.20/2007 (pdf)]

__Abstract.__A computationally efficient a posteriori error estimator is introduced and analyzed for collocation solutions to linear index-1 DAEs with properly stated leading term. The procedure is based on a modified defect correction principle, extending an established technique from the ODE context to the DAE case. We prove that the resulting error estimate is asymptotically correct, and illustrate the method by means of a numerical example.

To keep the presentation reasonably self-contained, we also briefly review some arguments from the literature on DAEs concerning the decoupling of the problem and its discretization, which is essential for our analysis.

W.Auzinger, W.Herfort:

**A uniform quantitative stiff stability estimate for BDF schemes,**

Opuscula Mathematica 26(2), 203-227, 2006. [open access]

__Abstract.__The concepts of stability regions, A- and A(alpha)-stability - albeit based on scalar models - turned out to be essential for the identification of implicit methods suitable for the integration of stiff ODEs. However, for multistep methods, knowledge of the stability region provides no information about the quantitative stability behavior of the scheme.

In this paper we fill this gap for the important class of Backward Differentiation Formulas (BDF). Quantitative stability bounds are derived which are uniformly valid in the stability region of the method. Our analysis is based on a study of the separation of the characteristic roots and a special similarity decomposition of the associated companion matrix.

W.Auzinger, E.Karner, O.Koch, E.Weinmüller:

**Collocation methods for the solution of eigenvalue problems for singular ordinary differential equations,**

Opuscula Mathematica 26(2), 229-241, 2006. [open access]

__Abstract.__We demonstrate that eigenvalue problems in ordinary differential equations can be recast in a formulation suitable for the solution by polynomial collocation. It is shown that the well-posedness of the two formulations is equivalent in the regular as well as in the singular case. Thus, a collocation code equipped with asymptotically correct error estimation and adaptive mesh selection can be successfully applied to compute the eigenvalues and eigenfunctions efficiently and with reliable control of the accuracy. Numerical examples illustrate this claim.

W.Auzinger, O.Koch, E.Weinmüller:

**Efficient mesh selection for collocation methods applied to singular BVPs,**

J. Comput. Appl. Math., 180(1), 213-227, 2005. (DOI 10.1016/j.cam.2004.10.013)

[ANUM Preprint No.7/04 (pdf)]

__Abstract.__We describe a mesh selection strategy for the numerical solution of boundary value problems for singular ordinary differential equations. This mesh adaptation procedure is implemented in ourcode`MATLAB``sbvp`which is based on polynomial collocation. We prove that under realistic assumptions our mesh selection strategy serves to approximately equidistribute the global error of the collocation solution, thus enabling to reach prescribed tolerances efficiently. Moreover, we demonstrate that this strategy yields a favorable performance of the code and compare its computational effort with other implementations of polynomial collocation.

W.Auzinger, O.Koch, D.Praetorius, E.Weinmüller:

**New a posteriori error estimates for singular boundary value problems,**

Numer. Algor. 40(1), 79-100, 2005. (DOI 10.1007/s11075-005-3791-5)

[preprint (pdf)]

__Abstract.__In this paper, we discuss the asymptotic properties and efficiency of several a posteriori estimates for the global error of collocation methods. Proofs of the asymptotic correctness are given for regular problems and for problems with a singularity of the first kind. We were also strongly interested in finding out which of our error estimates can be applied for the efficient solution of boundary value problems in ordinary differential equations with an essential singularity. Particularly, we compare estimates based on the defect correction principle with a strategy based on mesh halving.

W.Auzinger, O.Koch, E.Weinmüller:

**Analysis of a new error estimate for collocation methods applied to singular boundary value problems,**

SIAM J. Numer. Anal. 42(6), 2366-2386, 2005. (DOI 10.1137/S0036142902418928)

[ANUM Preprint No.13/02 (pdf)]

__Abstract.__We discuss an a-posteriori error estimate for the numerical solution of boundary value problems for nonlinear systems of ordinary differential equations with a singularity of the first kind. The estimate for the global error of an approximation obtained by collocation with piecewise polynomial functions is based on the defect correction principle. We prove that for collocation methods based on an even number of equidistant collocation nodes, the error estimate is asymptotically correct. As an essential prerequisite we derive convergence results for collocation methods applied to nonlinear singular problems.

W.Auzinger, H.Hofstätter, W.Kreuzer, E.Weinmüller:

**Modified defect correction algorithms for ODEs. Part II: Stiff initial value problems,**

Numer. Algor. 40(3), 285-303, 2005. (DOI 10.1007/s11075-005-5327-4)

[ANUM Preprint No.2/03 (pdf)]

__Abstract.__As shown in Part I of this paper and references therein, the classical method of Iterated Defect Correction (IDeC) can be modified in several nontrivial ways, extending the flexibility and range and applications of this type of iterative schemes. The essential point is an adequate definition of the defect, leading to a significantly more robust convergence behavior of the IDeC iteration, in particular for nonequidistant grids.

The present Part II is devoted to the efficient high-order integration of stiff initial value problems. By means of model problem investigation and systematic numerical experiment with a set of stiff test problems featuring different degrees of difficulty, our new versions of defect correction are systematically evaluated, and further algorithmic measures are proposed for the stiff case. The performance of the different variants under consideration is compared, and it is shown how strong coupling between non-stiff and stiff components can be successfully handled.

W.Auzinger, H.Hofstätter, W.Kreuzer, E.Weinmüller:

**Modified defect correction algorithms for ODEs. Part I: General theory,**

Numer. Algor. 36(2), 135-156, 2004. (DOI 10.1023/B:NUMA.0000033129.73715.7f)

Extended version: [ANUM Preprint No.1/03 (pdf)]

__Abstract.__The well-known method of Iterated Defect Correction (IDeC) is based on the following idea: Compute a simple, basic approximation and form its defect w.r.t. the given ODE via a piecewise interpolant. This defect is used to define an auxiliary, neighboring problem whose exact solution is known. Solving the neighboring problem with the basic discretization scheme yields a global error estimate. This can be used to construct an improved approximation, and the procedure can be iterated. The fixed point of such an iterative process corresponds to a certain collocating solution.

We present a variety of modifications to this algorithm. Some of these have been proposed only recently, and together they form a family of iterative techniques, each with its particular advantages. These modifications are based on techniques like defect quadrature (IQDeC), defect interpolation (IPDeC), and combinations thereof. We investigate the convergence on locally equidistant and nonequidistant grids and show how superconvergent approximations can be obtained. Numerical examples illustrate our considerations.

The application to stiff initial value problems will be discussed in Part II of this paper.

W.Auzinger, H.Hofstätter, O.Koch, W.Kreuzer, E.Weinmüller:

**Superconvergent defect correction algorithms,**

WSEAS Transactions on Systems 4, 1378-1383, 2004.

[preprint (pdf)]

__Abstract.__In this paper we discuss several variants of the acceleration technique known as Iterated Defect Correction (IDeC) for the numerical solution of initial value problems for ODEs. A first approximation, computed by a low order basic method, is iteratively improved to obtain higher order solutions. We propose new versions of the IDeC algorithm with maximal achievable (super-)convergence order twice as high as in the classical setting. Moreover, if the basic numerical method is designed for a special type of ODE only, as is the case for many geometric integrators, the idea of classical IDeC is not applicable in a straightforward way. Our approach enables the application of the defect correction principle in such cases as well.

W.Auzinger, O.Koch, E.Weinmüller:

**Collocation methods for boundary value problems with an essential singularity,**

in*Large-Scale Scientific Computing*(I.Lirkov, S.Margenov, J.Wasniewski and P.Yalamov, eds.), Springer Lecture Notes in Computer Science 2907, 347-354, 2004. (DOI 10.1007/978-3-540-24588-9_39)

[preprint (pdf)]

__Abstract.__We investigate collocation methods for the efficient solution of singular boundary value problems with an essential singularity. We give numerical evidence that this approach indeed yields high order solutions. Moreover, we discuss the issue of a posteriori error estimation for the collocation solution. An estimate based on the defect correction principle, which has been successfully applied to problems with a singularity of the first kind, is less robust with respect to an essential singularity than a classical strategy based on mesh halving.

W.Auzinger, C.Fabianek:

**Iterative solution of large linear systems arising in the 3-dimensional modelling of an electric field in human thigh,**

[Technical Report / ANUM Preprint No.12/04 (pdf)]

__Abstract.__Functional electrical stimulation (FES) of denervated skeletal muscles found in recent years an entry into rehabilitation of paraplegics. The aim of current research is the optimization of stimulation parameters as well as studying the ramifications of stimulating the muscles in the thigh.

To simulate the distribution of the electric field a 3-dimensional model of the thigh is created which is based on the theory of activation functions from Rattay 1990 and the muscle model from Reichel 1999. With this model it is possible to simulate the electrical activity of muscle fibers. In the current implementation conductivity values based on grey values obtained from CT data are used in the Poisson equation. Via discretization by the method of finite differences and solution of the arising systems this leads to the voltage respectively current distribution. This work aims to implement and integrate an efficient solver for the large linear systems arising and shorten the time of the solution process.

After a thorough analysis of the numerical properties of these systems, various solution strategies have been investigated. Beginning with direct solvers which take the symmetric and sparse nature into account, iterative solvers and finally Krylov subspace methods have been implemented and tested. In the sequel we focused on the preconditioned Conjugate Gradient method.

Because of the bad condition of the given systems the solvers mentioned above led to unsatisfactory results and special preconditioning methods became mandatory. Again various methods have been implemented and tested until multigrid methods finally led to outstanding results in terms of convergence speed. However, the quite high memory requirements could not be satisfied on the target computer and therefore tools for remote computing were used to leverage external machines. Now these multigrid methods are performed remotely on high-end computers and the solution time from previously 9 hours shortens to about an hour.

In course of the evaluation of the methods described above it became necessary to work with various model problems, and a new simple and especially small thigh model was developed. Using this model existing results could be verified and new insight was gained. In particular, this allows to evaluate new methods easier and quicker than before.

W.Auzinger, O.Koch, D.Praetorius, G.Pulverer, E.Weinmüller:

**Performance of collocation software for singular BVPs,**

[Technical Report ANUM Preprint No.4/04 (pdf)]

__Abstract.__In this study we investigate the performance of collocation codes applied to approximate the solution of systems of ordinary differential equations of the first order,

\( z'(t) = f(t,z(t)), t \in (a,b), B(z(a),z(b)) = 0 \).

The right-hand side of the differential equation may contain a singularity of the first kind, which means that it can be of the form

\( f(t,z(t)) = (1/t)M(t)z(t) + g(t,z(t)), t \in (0,1) \)

where \( M(t) \) is a matrix that depends continuously on \( t \) and \( g(t,z(t)) \) is a smooth vector-valued function.

The main purpose of this work is to develop the second version`sbvp 2.0`of our firstCode`MATLAB``sbvp 1.0`published in 2002. Therefore various tests are executed to find out about the performance of the code with different settings of some of the parameters.

In Chapter 2 the improved version of the Newton method is discussed. Here the old version of the code`sbvp 1.0`is used. In Chapter 3 the test runs for`sbvp 2.0`are recorded. The parameters are first set as in the beta-version of`sbvp`and then updated in accordance to the results of the tests. We first investigate the performance of the code for different parameter settings in context of the Jacobian update, further tests cover the questions dealing with the number of meshpoints on the starting grid, the setting for the monitor function or some possibilities of mesh distribution by varying the ratio between the largest and the smallest stepsize.

After the completion of this main phase of this work, in Chapter 4 we compare the new code to other collocation solvers. These other codes are`sbvp 1.0`, the old version of our code as well as thestandard code`MATLAB``BVP4C`and the`Fortran 90`routine`COLNEW`.

Main criterion for most of the observations is the runtime required to solve the boundary value problem. As`COLNEW`is not asolver, this criterion is to be exchanged with the`MATLAB``fcounts`and`dfcounts`and the number of meshpoints on the final grid.

The boundary value problems that are used as test models throughout this study can be found in Chapter 5.

W.Auzinger, O.Koch, E.Weinmüller:

**New variants of defect correction for boundary value problems in ordinary differential equations,**

in*Current Trends in Scientific Computing*, Z.Chen, R.Glowinski, K.Li (eds), AMS Series in Contemporary Mathematics 329, 43-50, 2003.

[preprint (pdf)]

__Abstract.__In this paper we discuss new variants of the acceleration technique known as Iterated Defect Correction (IDeC) for the numerical solution of boundary value problems in ordinary differential equations. A first approximation computed by the backward Euler scheme, is iteratively improved to obtain a high order solution. Typically, the maximal attainable accuracy is limited by the smoothness of the exact solution and by technical details of the procedure. However, for our new version of the IDeC algorithm the maximal achievable order is higher than in the classical setting. Moreover, our procedure can be shown to be convergent on arbitrary grids, while the classical IDeC iteration is restricted to the equidistant case. Finally, the performance of this new algorithm for singular boundary value problems is discussed.

W.Auzinger:

**Sectorial operators and normalized numerical range,**

Appl. Numer. Math. 45(4), 367-388, 2003. (DOI 10.1016/S0168-9274(02)00254-4)

[article (pdf)]

__Abstract.__We study the notion of sectorial operator in a Hilbert space. According to the classical definition, the numerical range \( R(A) \) of a sectorial operator A is contained in a sector \( S_\sigma = { z \in {\mathbb C}: |Arg z| \leq \sigma } \), and this is equivalent to a certain inverse estimate valid outside \( S_\sigma \). In this paper we show that the validity of the same estimate, but with a factor >1 , is equivalent to the validity of a certain strengthened Cauchy-Schwarz inequality for all pairs \( w, Aw \). This extends the original characterization in terms of \( R(A) \) by a more general characterization based on a normalized numerical range \( R_N(A) \). We also show how \( R_N(A) \) can be computed numerically.

W.Auzinger, G.Kneisl, O.Koch, E.Weinmüller:

**A collocation code for boundary value problems in ordinary differential equations,**

Numer. Algor. 33, 27-39, 2003. (DOI 10.1023/A:1025531130904)

[ANUM Preprint No.18/01 (pdf)]

__Abstract.__We present apackage for boundary value problems in ordinary differential equations. Our aim is the efficient numerical solution of systems of ODEs with a singularity of the first kind, but the solver is also well suited for regular problems. The basic solution is computed using collocation methods and a new, efficient estimate of the global error is used for adaptive mesh selection. Here, we analyze some of the numerical aspects relevant for the implementation, describe measures to increase the efficiency of the code and compare its performance with the performance of established standard codes for boundary value problems.`MATLAB`

W.Auzinger, E.Karner, O.Koch, D.Praetorius, E.Weinmüller:

**Globale Fehlerschätzer für Randwertprobleme mit einer Singularität zweiter Art,**

[Technical Report / ANUM Preprint No.6/03 (pdf)]

__Abstract.__Singuläre Randwertprobleme (singular boundary value problems, BVPs) mit einer Singularität der ersten Art treten oft in Anwendungen auf, z.B. wenn eine partielle Differentialgleichung (partial differential equation, PDE) bei Vorhandensein von Symmetrien zu einer gewöhnlichen Differentialgleichung (ordinary differential equation, ODE) reduziert wird. Schwierigere Probleme mit einer Singularität der zweiten Art (essentielle Singularität) treten zum Beispiel dann auf, wenn man die Randwertaufgabe für eine gewöhnliche Differentialgleichung auf einem unendlichen Intervall zu einer auf einem endlichen Intervall umformt. Weiters finden wir sie in Gebieten der Quantenphysik, Mechanik, usw.

Der-Code`MATLAB``sbvp`löst singuläre Randwertprobleme mittels eines Kollokationsverfahrens unter Zuhilfenahme eines a posteriori Fehlerschätzers. Dieser Code war ursprünglich zur Lösung von regulären Randwertproblemen und Randwertproblemen mit einer Singularität der ersten Art gedacht. Diese Arbeit untersucht verschiedene a posteriori Fehlerschätzer bei Beispielen mit einer Singularität zweiter Art auf deren Effizienz.

W.Auzinger, O.Koch, J.Petrickovic, E.Weinmüller:

**Numerical solution of boundary value problems with an essential singularity,**

[Technical Report / ANUM Preprint No.3/03 (pdf)]

__Abstract.__Singular boundary value problems with a singularity of the first kind frequently occur in applications, e.g. when a partial differential equation PDE is reduced to an ordinary differential equation in the presence of symmetries. More difficult problems featuring a singularity of the second kind (essential singularity) arise, for instance, when a boundary value ODE is transformed from an infinite interval to a finite one. They are also very common in quantum physics, mechanics, investigations of ferromagnetic systems (e.g. Ginzburg-Landau equations) etc.

In this work we present a study of certain numerical techniques - which are well known to perform efficiently and accurately in the case of a singularity of the first kind - applied to problems with an essential singularity. In particular, collocation methods are considered together with a posteriori error estimates intended to provide the basis for a grid selection strategy.

W.Auzinger, O.Koch, E.Weinmüller:

**Efficient collocation schemes for singular boundary value problems,**

Numer. Algor. 31, 5-25, 2002. (DOI 10.1023/A:1021151821275)

[article (pdf)]

__Abstract.__We discuss an error estimation procedure for the global error of collocation schemes applied to solve singular boundary value problems with a singularity of the first kind. This estimate of the global error was proposed by Stetter in 1978 and is based on the idea of Defect Correction, originally due to Zadunaisky. A carefully designed modification of this estimate presented here, not only results in less computational work but also proves robust with respect to the singularity. This estimate also provides a basis for a grid selection routine in which the grid is modified with the aim to equidistribute the global error. This procedure yields meshes suitable for an efficient numerical solution - most importantly, we observe that the grid is refined in a way reflecting only the behavior of the solution and remains unaffected by the unsmooth direction field close to the singular point.

W.Poth, M.Matzl, W.Auzinger, A.Steindl, H.Troger:

**Comparison of displacement versus natural variables for the numerical simulation of a string pendulum,**

Nonlinear Dynamics, Chaos, Control and their Applications to Engineering Sciences, Vol.5: Chaos, Control and Time Series (J.M.Balthazar, P.B.Goncalves, R.F.Brasil, I.L.Caldas, F.B. Rizatto, Eds.) 127-136, 2002.

[preprint (pdf)]

__Abstract.__Simulation of the dynamic behavior of tethered satellite systems requires the numerical solution of a system of coupled nonlinear partial and ordinary differential equations. Due to the design of tethered satellite systems this set of equations of motion is a stiff system. To integrate this system of equations of motion numerically creates a number of problems. In order to address and solve some of these problems we consider the simpler mechanical model of the string pendulum which, however, displays all relevant features of a tethered satellite system.

In this paper a comparison of two different choices of variables and of various different choices of time integrators is presented in their influence on the efficient and accurate numerical integration of large amplitude motions of the string pendulum. our results show that for stiff strings a description of the deformation of the string alternative to the usual description by displacements may result in a substantial increase both of accuracy and reduction of computation time. Further the proper choice of the time integrator has great influence as well.

W.Auzinger, O.Koch, E.Weinmüller:

**Theory and solution techniques for singular boundary value problems in ordinary differential equations,**

in `Parallel Processing and Applied Mathematics', Proceedings, 4th International Conference, PPAM 2001, 851-861. (Springer-Verlag, Berlin, Heidelberg, 2002.)

[article (pdf)]

__Abstract.__In this paper we present an overview of analytical results and numerical methods for singular boundary value problems for ordinary differential equations with a singularity of the first kind. Special attention is paid to the analysis of shooting methods, where the associated initial value problems are solved by the acceleration technique known as Iterated Defect Correction (IDeC) based on the backward Euler method, and on direct discretization using collocation schemes. Convergence, error estimation and mesh selection are discussed for both approaches. Moreover, we study the fixed point convergence of the IDeC iteration, where the fixed point corresponds to a collocation solution.

We show how to implement the collocation scheme and the error estimation routine, provide a roundoff-error analysis of the solution process, devise some numerical tests that illustrate the performance of the solver and error estimate and finally describe usage and development of the solver-package.

W.Auzinger, O.Koch, S.Lammer, E.Weinmüller:

**Variationen der Defektkorrektur zur effizienten numerischen Lösung gewöhnlicher Differentialgleichungen,**

[Technical Report / ANUM Preprint No.8/02 (pdf)]

__Abstract.__

W.Auzinger, G.Kneisl, O.Koch, E.Weinmüller:

**sbvp 1.0 - A**`MATLAB`solver for singular boundary value problems,

[Technical Report / ANUM Preprint No.2/02 (pdf)]

__Abstract.__

W.Auzinger, G.Kneisl, O.Koch, E.Weinmüller:

**A solution routine for singular boundary value problems,**

[Technical Report / ANUM Preprint No.1/02 (pdf)]

__Abstract.__In this report, we discuss the implementation and numerical aspects of the`MATLAB`solver`sbvp`designed for the solution of two-point boundary value problems, which may include a singularity of the first kind, \( z'(t)=f(t,z(t)) := \frac{1}{(t-a)} M(t)z(t) + g(t,z(t)), t \in (a,b), R(z(a),z(b)) = 0 \). The code is based on collocation at either equidistant or Gaussian collocation points. For singular problems, the choice of equidistant nodes does not imply a loss of efficiency since no convergence order higher than the stage order can be expected in general.

An estimate of the global error of the approximate solution is also provided. This estimate is obtained by a modification of the Defect Correction idea originally proposed by Zadunaisky in 1976. This estimate has been proven to be asymptotically correct and provides the basis for an adaptive mesh selection strategy. Here the grid is modified with the aim to equidistribute the global error. Most importantly, we observe that the grid is refined in a way reflecting merely the smoothness of the solution, unaffected by the singularity in \( f \).

We discuss details of the efficient implementation in`MATLAB`6 and illustrate the performance of the code by comparing it with the standard`MATLAB`solver`bvp4c`and the Fortran 90 code COLNEW.

The`sbvp`package is available at MATLAB Central.

W.Auzinger, A.Eder:

**A note on Lyapunov transformation and exponential decay in linear ODE systems,**

Mathematical Models and Methods in Applied Sciences (M3AS), Vol.11, No.1, 23-31, 2001. (DOI 10.1142/S0218202501000714)

[preprint (pdf)]

__Abstract.__In this paper we consider a class of matrices \( A \) where all eigenvalues have negative real part and are of a common magnitude \( O(1) \). Concerning the behavior of \( e^{tA} \) we provide a necessary and sufficient condition, via Lyapunov transformation, for an estimate of the form \( \|e^{tA}\|_2 \leq K_0 \exp(-t/K_1) \) to be valid uniformly for \( t \gt 0 \) with moderate-sized constants \( K_0 \) and \( K_1 \). All relevant relations are quantitatively specified.

W.Auzinger, G.Kneisl, O.Koch, W.Polster, E.Weinmüller:

**Implementation of a solution routine for singular boundary value problems,**

ARGESIM - Working Group Simulation News S70 (2001), Vienna University of Technology, Austria.

[ANUM Preprint No.24/01 (pdf)]

__Abstract.__This paper deals with the implementation of a`MATLAB`-solver for singular boundary value problems with a singularity of the first kind. This solver includes a mesh adaption routine that selects meshes according to the local smoothness of the solution rather than to the smoothness of the direction field. The mesh selection is based on a new global error estimator for collocation schemes, which has recently been developed at the Institute for Applied Mathematics and Numerical Analysis and proves robust with respect to the singularity.

W.Auzinger, A.Eder, R.Frank:

**Extending convergence theory for nonlinear stiff problems, part II,**

[Technical Report / ANUM Preprint No.22/01 (pdf)]

__Abstract.__This paper deals with the convergence properties of high-order implicit Runge-Kutta methods applied to nonlinear stiff initial value problems. Earlier convergence concepts like the theory of B-convergence or singular perturbation analysis are not successfully applicable to stiff problems with a general 'geometry'. To overcome this drawback to a certain extent, in Part I of this paper a problem class was introduced, where the stiffness is axiomatically characterized in natural geometric terms, with a nonlinearly varying 'stiff eigendirection' corresponding to a 'stiff eigenvalue'. Furthermore, a convergence analysis for the implicit Euler scheme was presented.

In the meantime, some work has been done to extend this theory to general Runge-Kutta schemes. In particular, for the class of stiff problems considered in Part I, a complete analysis of Radau and Gauss schemes has been developed. In this paper these results are presented, commented, and the essential facts are proved.

Due to their excellent stability properties, the Radau IIa schemes are of particular interest. Therefore we will work out the theory for these methods in some details, and content ourselves with reporting the corresponding results for Radau Ia and Gauss schemes, which are obtained in an analogous way.

Lack of space prevents us from presenting all the technical details, which can be found in an underlying report. In the proofs given, straightforward but lengthy algebraic manipulations are not always worked out in full detail.

Our error estimates show that for certain stepsize ranges, the classical order of superconvergence can even be obtained in the stiff case. This is a further improvement compared to the B-convergence theory, where only the stage order is maintained in the error bounds.

W.Auzinger, O.Koch, W.Polster, E.Weinmüller:

**Ein Algorithmus zur Gittersteuerung bei Kollokationsverfahren für singuläre Randwertprobleme,**

[Technical Report / ANUM Preprint No.21/01 (pdf)]

__Abstract.__

W.Auzinger, R.Frank, W.Kreuzer, E.Weinmüller:

**Computergestützte Analyse verschiedener Varianten der Iterierten Defektkorrektur unter besonderer Berücksichtigung steifer Differentialgleichungen,**

[Technical Report / ANUM Preprint No.20/01 (pdf)]

__Abstract.__Diese Arbeit stellt eine Fortsetzung der Untersuchungen von Report Nr.134/00 dar. Im Schwerpunkt stehen die experimentelle Erprobung verschiedener Basisverfahren für die Defektkorrektur und der Vergleich der interpolierenden Defektkorrektur mit ähnlichen, von anderer Seite vorgeschlagenen Varianten.

W.Auzinger, R.Frank, H.Hofstätter, E.Weinmüller:

**Defektkorrektur zur numerischen Lösung steifer Anfangswertprobleme,**

[Report No.131/00, Institut für Angewandte und Numerische Mathematik, TU Wien, 2000 (pdf)]

__Abstract.__Zur effizienten Lösung hoch-impliziter Runge-Kutta Verfahren, wie sie für die numerische Integration steifer Anfangswertprobleme verwendet werden, bietet sich der iterative Zugang mittels Defektkorrektur an. In dieser Arbeit werden verschiedenste Varianten der Defektkorrektur theoretisch und experimentell untersucht, insbesondere auch ihre Kombination mit verschiedenen Methoden zur Konvergenzbeschleunigung.

W.Auzinger, A.Eder, R.Frank:

**Eine erweiterte Konvergenztheorie impliziter Runge-Kutta Verfahren für nichtlineare steife Differentialgleichungen,**

[Report No.130/00, Institut für Angewandte und Numerische Mathematik, TU Wien, 2000 (pdf)]

__Abstract.__Dieser Report enthält eine detaillierte Konvergenzanalyse für implizite Runge-Kutta Verfahren vom Typ Gauss- und Radau-Typ bei Anwendung auf eine Klasse nichtlinearer steifer Differentialgleichungen, bei der die Steifheit mittels einer natürlichen nichtlinearen Verallgemeinereung der Begriffe `steifer Eigenwert' und `steifer Eigenvektor' charakterisiert ist. Es werden präzise Ordnungsaussagen für den globalen Fehler bzw. für seine glatten und steifen Anteile hergeleitet.

W.Auzinger, O.Koch, P.Kofler, E.Weinmüller:

**Acceleration techniques for singular initial value problems,**

[Report No.129/00, Institut für Angewandte und Numerische Mathematik, TU Wien, 2000 (pdf)]

__Abstract.__We consider the numerical solution of singular initial and terminal value problems using various low-order Runge-Kutta methods. With these basic solutions, we investigate the acceleration technique known as Iterated Defect Correction (IDeC). We show that the performance depends crucially on the asymptotic expansions of the global error. The results are compared with the asymptotic properties of extrapolation.

W.Auzinger, A.Eder, R.Frank:

**Convergence theory for implicit Runge-Kutta methods applied to a one-parameter family of stiff autonomous differential equations,**

[Report No.123/00, Institut für Angewandte und Numerische Mathematik, TU Wien, 2000 (pdf)]

__Abstract.__The intention of this paper is to extend the convergence concepts for discretization methods applied to nonlinear stiff problems. First a one-paremeter family of stiff autonomous differential equations is introduced, where stiffness is axiomatically characterized in geometric terms. Then the discretization methods are analyzed, where we restrict our considerations to the implicit Runge-Kutta methods of the type Radau Ia, IIa and Gauss. For these methods we prove the solvability of the algebraic equations and derive global error bounds.

W.Auzinger, O.Koch, P.Kofler, E.Weinmüller:

**The application of shooting to singular boundary value problems,**

[Report No.126/99, Institut für Angewandte und Numerische Mathematik, TU Wien, 2000 (pdf)]

__Abstract.__A certain class of singular boundary value problems is solved numerically using the shooting method. The associated initial value problems are solved using Iterated Defect Correction, where the implicit Euler rule serves as the basic method.

W.Auzinger, P.Kofler, E.Weinmüller:

**Steuerungsmaßnahmen bei numerischer Lösung singulärer Anfangswertaufgaben,**

[Report No.124/98, Institut für Angewandte und Numerische Mathematik, TU Wien, 1998 (pdf)]

__Abstract.__Wir betrachten nichtlineare Anfangswertaufgaben für Systeme gewöhnlicher Differentialgleichungen erster Ordnung mit einer Singularität erster Art. Basierend auf einer Arbeit von Koch, Kofler und Weinmüller, welche das analytische Hintergrundwissen über die betrachtete Problemklasse liefert, werden verschiedene Methoden der Schätzung des lokalen und globalen Diskretisierungsfehlers von Runge-Kutta Verfahren sowie verschiedene Gitterwahlstrategien experimentell untersucht.

Es zeigt sich, dass in der Nähe der Singularität die Ordnung des lokalen Fehlers fast immer reduziert ist, wodurch die Schrittweitenwahl basierend auf der Schätzung des lokalen Fehlers nicht sehr zuverlässig ist. Im Gegensatz dazu funktionieren Methoden zur Schätzung des globalen Diskretisierungsfehlers ausgesprochen gut. Der globale Fehler wird gut wiedergegeben, und man erhält bereits auf groben Gittern eine relativ genaue Schätzung, woraus ein zur Integration singulärer Anfangswertaufgaben geeigneter Algorithmus hergeleitet werden kann.

W.Auzinger, R.Frank:

**B-Convergence,**

in*Encyclopedia of Mathematics Supplement Volume I*, p.71, Kluwer Academic Publishers, 1997.

[preprint (pdf)]

__Abstract.__This is a short encyclopedical overview on the theory of B-convergence, which predicts the convergence properties of discretization methods, implicit Runge-Kutta methods in particular, applied to initial value problems for systems of non-linear ordinary differential equations \( y'=f(t,y) \) where \( f \) satisfies a one-sided Lipschitz condition.

W.Auzinger, R.Frank, H.J.Stetter:

**Vienna contributions to the development of RK-methods,**

Appl. Numer. Math. 22, 35-49, 1996. (DOI 10.1016/S0168-9274(96)00024-4)

[preprint (pdf)]

__Abstract.__In 1965, H.J.Stetter moved from his home in Munich to Vienna, to take a newly created chair for Numerical Mathematics at the Technical University of Vienna. At the Technical University of Munich, since 1953, he had grown up in one of the first major research groups in the new ara of Numerical Analysis, with one of the first electronic computers available since 1955. His own interests had been strongly tied to the computational solution of initial value problems in both ordinary and partial differential equations. He took these interests with him to Vienna where they spread to his students and collaborators. Thus, a research group in ODEs originated at the Technical University in Vienna and RK-methods were an important objective. In the following, we sketch the related activities and achievements of this group from 1965 to these days.

The main emphasis has been given to the work related with stiff systems. In particular, section 5 is devoted to a rather detailed discussion of the nature of stiffness and its implications on the convergence theory of RK-methods for nonlinear stiff problems from a point of view which has originated from the Vienna group since the late Seventies.

W.Auzinger, R.Frank, G.Kirlinger:

**Extending convergence theory for nonlinear stiff problems, part I,**

BIT Numer. Math. 36(4), 635-652, 1996. (DOI 10.1007/BF01733784)

[preprint (pdf)]

__Abstract.__Existing convergence concepts for the analysis of discretizations of nonlinear stiff problems suffer from considerable drawbacks. Our intention is to extend the convergence theory to a relevant class of nonlinear problems, where stiffness is axiomatically characterized in natural geometric terms.

Our results will be presented in a series of papers. In the present paper (Part I) we motivate the need for such an extension of the existing theory, and our approach is illustrated by means of a convergence argument for the Implicit Euler scheme.

W.Auzinger, G.Kirlinger:

**Kreiss resolvent conditions and strengthened Cauchy-Schwarz inequalities,**

Appl. Numer. Math. 18, 57-67, 1995. (DOI 10.1016/0168-9274(95)00043-T)

[article (pdf)]

__Abstract.__In the present paper it is shown how the resolvent conditions in the Kreiss Matrix Theorem for \( e^{tA} \) and \( A^\nu \), respectively, can be reformulated as certain Strengthened Cauchy-Schwarz Inequalities (SCSIs) to be satisfied for all pairs \( w,Aw \) (\( w \in {\mathbb C}^n \)). This yields certain generalizations of the notions `logarithmic norm' and `numerical radius', respectively.

W.Auzinger, P.Trieb, E.Weinmüller:

**Beschleunigte Algorithmen zur Lösung von singulären Randwertaufgaben,**

Report No.113/94, Institut für Angewandte und Numerische Mathematik, TU Wien, 1994. [outline (pdf)]

__Abstract.__Betrachtet werden lineare Randwertaufgaben zweiter Ordnung mit einer Singularität erster Art. Im Falle hoher Genauigkeitsforderungen sind Verfahren höherer Ordnung wie z.B. Kollokationsverfahren meist mit sehr hohem Aufwand verbunden, eine Tatsache, die speziell in den nichtlinearen Anwendungsbeispielen aus der Mechanik und der Chemie deutlich zutage tritt. Demgegenüber besitzen Beschleunigte Algorithmen den erheblichen Vorteil, eine sehr genaue Lösung bei vernünftigem Aufwand zu liefern. In dieser Arbeit betrachten wir eine einfache Variante der sogenannten Iterierten Defektkorrektur (IDeC), die auf einer Idee von Zadunaisky basiert. In dieser Arbeit wird das Verhalten dieser Methode bei Anwendung auf singuläre Randwertprobleme experimentell untersucht. Die erzielten Ergebnisse werden mit theoretischen Resultaten für nichtsinguläre Randwertprobleme verglichen; signifikante Abweichungen werden diskutiert und dienen als Anlass für Modifikationen des Verfahrens.

W.Auzinger, R.Frank, G.Kirlinger:

**Modern convergence theory for stiff initial-value problems,**

J. Comput. Appl. Math. 45, 5-16, 1993. (DOI 10.1016/0377-0427(93)90260-I)

[preprint (pdf)]

__Abstract.__In this paper we give a brief review of available theoretical results about convergence and error structures for discretizations of stiff initial value problems. We point out limitations of the various approaches and discuss some recent developments.

W.Auzinger, R.Frank, G.Kirlinger:

**An extension of B-convergence for Runge-Kutta methods,**

Appl. Numer. Math. 9(2), 91-109, 1992. (DOI 10.1016/0168-9274(92)90008-2)

[preprint (pdf)]

__Abstract.__The well-known concepts of B-stability and B-convergence for the analysis of one-step methods applied to stiff initial value problems are based on the notion of one-sided Lipschitz continuity. In a recent paper the authors have pointed out that the one-sided Lipschitz constant \( m \) must often be expected to be very large (positive and of the order of magnitude of the stiff eigenvalues) despite a (globally) well-conditioned behavior of the underlying problem. As a consequence, the existing B-theory suffers from considerable restrictions; e.g., not even linear systems with time dependent coefficients are satisfactorily covered. The purpose of the present paper is to fill this gap; for implicit Runge-Kutta methods we extend the B-convergence theory such as to be valid for a class of non-autonomous weakly nonlinear stiff systems; reference to the (potentially large) one-sided Lipschitz constant is avoided. Unique solvability of the system of algebraic equations is shown, and global error bounds are derived.

W.Auzinger, R.Frank, E.Weinmüller:

**On the efficient numerical integration of conservation laws with stiff source terms,**

[Report No.94/92, Institut für Angewandte und Numerische Mathematik, TU Wien, 1992 (pdf)]

__Abstract.__We consider scalar 1D-conservation laws with stiff source terms, \( u_t + (f(u))_x = \psi(u) \). Such problems arise e.g. as simple models in the study of chemically reacting gas dynamics in combustion processes. While for pure conservation laws a variety of reliable numerical methods are known, the presence of a stiff source term causes severe additional difficulties. Several standard approaches (e.g. splitting into pure conservation laws and solution of stiff ODEs) fail unless extremely fine meshes are used. We propose an alternative method which successfully and efficiently solves certain model problems. The generalization of these ideas to more general problem classes remains to be studied.

W.Auzinger, R.Frank, G.Kirlinger:

**The state of the art in the convergence theory for stiff ODEs,**

Proceedings, IMACS 91, Volume 1, 280-284 (Eds. R.Vichnevetsky, J.J.H.Miller, Dublin, 1991).

[preprint (pdf)]

__Abstract.__In this paper we give a brief review on available theoretical results about convergence and error structures for discretizations of stiff initial value problems. We point out limitations of the various approaches and discuss some recent developments.

W.Auzinger:

**Error structures and extrapolation in time for parabolic problems,**

Bonner Mathematische Schriften, Nr.228, 41-56, Bonn, 1991.

[preprint (pdf)]

__Abstract.__We study asymptotic expansions of the time discretization error for some time stepping schemes applied to (semi-discretizations of) linear parabolic initial/boundary value problems. For this type of problems, the `classical' theory of asymptotic error expansions cannot be applied directly; a more careful analysis is necessary to understand the structure of the time discretization error. We recall and discuss our respective results for the smooth case and present numerical experiments with extrapolation confirming these results. In the final section we discuss the case of nonsmooth solutions caused by inconsistent initial/boundary data.

W.Auzinger, R.Frank, G.Kirlinger:

**Error structures for stiff initial value problems,**

in `Numerical Treatment of Differential Equations', Proceedings, NUMDIFF-5 (Ed. K.Strehmel), Halle, 1989. (Teubner-Texte zur Mathematik 121, 1991, 16-25.)

[preprint (pdf)]

__Abstract.__Within the convergence theory of discretization methods, the existence of asymptotic error expansions in powers of the stepsize \( h \) is the classical prerequisite for the theoretical justification of stepsize control mechanisms and acceleration techniques like extrapolation or defect correction. In the present paper we discuss the existence of asymptotic error expansions and their algorithmic utilization in the context of stiff nonlinear initial value problems. Most of the material presented is discussed in more detail in a recent series of papers on this subject.

W.Auzinger:

**On error structures and extrapolation for stiff systems, with application in the method of lines,**

Computing 44, 331-356, 1990. (DOI 10.1007/BF02241272)

[preprint (pdf)]

__Abstract.__In this paper, the structure of the global error is studied for some time discretization schemes, applied to a class of stiff initial value problems as they typically arise from the semi-discretization of parabolic initial/boundary value problems (method of lines). The implicit Euler and trapezoidal schemes and a locally one-dimensional splitting method are considered, and `perturbed' asymptotic error expansions are derived which are valid independent of the stiffness (independent of the meshwidth in space). The key point within the analysis is a careful, quantitative description of the remainder term in such an expansion. The results are applicable in the method of lines setting and enable the prediction of the behavior of extrapolation algorithms for the class of problems under consideration. These theoretical considerations are illustrated by numerical examples.

W.Auzinger, R.Frank, G.Kirlinger:

**A note on convergence concepts for stiff problems,**

Computing 44, 197-208, 1990. (DOI 10.1007/BF02262216)

[preprint (pdf)]

__Abstract.__Most convergence concepts for discretizations of nonlinear stiff initial value problems are based on one-sided Lipschitz continuity. Therefore only those stiff problems that admit moderately sized one-sided Lipschitz constants are covered in a satisfactory way by the respective theory. In the present note we show that the assumption of moderately sized one-sided Lipschitz constants is violated for many stiff problems. We recall some convergence results that are not based on one-sided Lipschitz constants; the concept of singular perturbations is one of the key issues. Numerical experience with stiff problems that are not covered by available convergence results is reported.

W.Auzinger, R.Frank, G.Kirlinger:

**Asymptotic error expansions for stiff equations: applications,**

Computing 43, 223-253, 1990. (DOI 10.1007/BF02242919)

[preprint (pdf)]

__Abstract.__In a series of foregoing papers we have studied the structure of the global discretization error for the implicit Euler scheme and the implicit midpoint and trapezoidal rules applied to a general class of nonlinear stiff initial value problems. Full asymptotic error expansions (in the conventional sense) exist only in special situations; for the general case, asymptotic expansions in a weaker sense have been derived. In the present paper we demonstrate how these results can be used for an analysis of acceleration techniques applied to stiff problems. In particular, extrapolation and defect correction algorithms are considered. Various numerical results are presented and discussed.

W.Auzinger, R.Frank, F.Macsek:

**Asymptotic error expansions for stiff equations: The implicit Euler scheme,**

SIAM J. Numer. Anal. 27(1), 67-104, 1990. (DOI 10.1137/0727005)

[preprint (pdf)]

__Abstract.__This paper discusses the existence of an asymptotic expansion for the global error of the implicit Euler scheme applied to stiff systems of ordinary differential equations. It is shown that in strongly stiff situations, a full asymptotic expansion exists at all gridpoints. For the mildly stiff case it is shown that the full order of the remainder term, which inevitably breaks down at the first gridpoints after a stepsize change, reappears at the subsequent gridpoints. Our analysis is based on singular perturbation techniques.

W.Auzinger, H.J.Stetter:

**A study of numerical elimination for the solution of multivariate polynomial systems,**

Report, Institut für Angewandte und Numerische Mathematik, TU Wien, 1989.

[report (pdf)]

__Abstract.__In an earlier paper we had motivated and described am algorithm for the computation of the zeros of multivariate polynomial systems by means of multiplication tables and solution of an associated eigenvalue problem. In this paper we give a more transparent description of this approach, including degenerate situations. Both the theoretical results and the algorithmic procedures are applied to examples, and the results are displayed and discussed.

W.Auzinger:

**On the error structure of the implicit Euler scheme applied to stiff systems of differential equations,**

Computing 43, 115-131, 1989. (DOI 10.1007/BF02241856)

[preprint (pdf)]

__Abstract.__In this paper we investigate the structure of the global discretization error of the implicit Euler scheme applied to systems of stiff differential equations, extending earlier work on this subject. We restrain our considerations to the linear, self-adjoint, constant coefficient case but we make no assumptions about the nature of the stiff spectrum; the stiff eigenvalues may be arbitrarily distributed on the negative real axis.

Our main result says that the global error of the implicit Euler scheme admits an asymptotic expansion in powers of the stepsize \( \tau \) which is not asymptotically correct in the conventional sense: Near the initial point \( t=0 \) the expansion is spoiled at the \( O(\tau^2) \)-level by `irregular' error components which are, however, (algebraically) damped, such that away from \( t=0 \) the `pure' asymptotic expansion reappears. We present numerical experiments confirming this result. Our considerations should be particularly helpful for a rigorous, quantitative analysis of the structure of the full (space & time) discretization error in the PDE (method of lines) context, and thus for a sound theoretical justification of extrapolation techniques for this important class of stiff problems.

W.Auzinger, R.Frank:

**Asymptotic expansions of the global discretization error for stiff problems,**

SIAM J. Sci. Stat. Comput. 10(5), 950-963, 1989. (DOI 10.1137/0910055)

[preprint (pdf)]

__Abstract.__We discuss the existence of asymptotic expansions of the global discretization error for a general class of nonlinear stiff differential equations \( y'(t) = A(t)y(t) + \phi(t,y(t)) \) where \( A(t) \) has a `stiff spectrum' characterized by a small parameter \( \epsilon \) and where \( \phi(t,y) \) is smooth. The following methods are considered: Implicit Euler, implicit midpoint and trapezoidal rules. In strongly stiff situations (\( \epsilon \) significantly smaller than the stepsize \( h \)) the implicit Euler scheme admits a full asymptotic expansion; the same is true for the midpoint rule and for the trapezoidal rule under certain coupling conditions. In those strongly stiff cases where a full expansion does not exist for the midpoint or trapezoidal rule, the remainder term is of a reduced order but shows a regular, oscillating behavior which we describe in detail. In mildly stiff situations, order reductions of the remainder term inevitably occur in any case after the start or after the change of stepsize but - as can be shown by discrete singular perturbation techniques - these order reductions are rapidly damped out as the integration proceeds. Our results are illustrated by various numerical examples; in particular, numerical experience with extrapolation and Defect Correction is reported.

W.Auzinger, R.Frank:

**Asymptotic error expansions for stiff equations: An analysis for the implicit midpoint and trapezoidal rules in the strongly stiff case,**

Numer. Math. 56(5), 469-499, 1989. (DOI 10.1007/BF01396649)

[preprint (pdf)]

__Abstract.__The structure of the global discretization error is studied for the implicit midpoint and trapezoidal rules applied to nonlinear*stiff*initial value problems. The point is that, in general, the global error contains nonsmooth (oscillating) terms at the dominant \( O(h^2) \)-level. However, it is shown in the present paper that for special classes of stiff problems these nonsmooth terms contain an additional factor \( \epsilon \) (where \( -1/\epsilon \) is the magnitude of the stiff eigenvalues). In these cases a `full' asymptotic error expansion exists in the*strongly stiff case*(\( \epsilon \) sufficiently small compared to the stepsize \( h \)).

W.Auzinger, R.Frank:

**Asymptotische Fehlerentwicklungen bei steifen Differentialgleichungen: Ergebnisse und numerische Verifikation,**

ZAMM 69(4), T273-T274, 1989.

[preprint (pdf)]

__Abstract.__Diese Arbeit gibt einen kurzen überblick über neuere Resultate bezüglich der Fehlerstruktur der impliziten Mittelpunktsregel und der impliziten Trapezregel bei Anwendung auf steife Anfangswertprobleme. Es ergeben sich hier `gestörte' asymptotische Fehlerentwicklungen, deren Restglied im allgemeinen nicht die volle Ordnung aufweist und typischerweise einen oszillierenden Charakter besitzt. Diese Resultate liefern Information darüber, wie konventionelle Beschleunigte Algorithmen zu modifizieren sind, um auch im steifen Fall weitgehend effizient zu funktionieren.

W.Auzinger, R.Frank, F.Macsek:

**Fehlerstrukturen bei steifen Differentialgleichungen,**

ZAMM 69(4), T133-T134, 1989.

[preprint (pdf)]

__Abstract.__Diese Arbeit gibt einen kurzen überblick über neuere Resultate bezüglich der Fehlerstruktur des impliziten Euler-Verfahrens bei Anwendung auf steife Anfangswertprobleme. Es ergeben sich hier `gestörte' asymptotische Fehlerentwicklungen, deren Restglied speziell in mittelsteifen Situationen erst nach einer übergangsphase die volle Ordnung erreicht.

W.Auzinger, R.Frank:

**Asymptotic error expansions for stiff equations: The implicit midpoint rule,**

[Report No.77/88, Institut für Angewandte und Numerische Mathematik, TU Wien, 1988 (pdf)]

__Abstract.__This is a study of the structure of the global discretization error of the implicit midpoint rule applied to nonlinear stiff initial value problems. Those special cases where an asymptotically correct expansion in even powers of the stepsize \( h \) exists have been described in an earlier paper. In the general stiff case the global error does not admit a "pure" asymptotic expansion but there occur dominant oscillating terms at the \( h^2 \)-level. The amplitude of these oscillations shows a very smooth behavior and, consequently, the error structure is sufficiently regular to guarantee a satisfactory efficiency of acceleration algorithms like extrapolation or defect correction. These algorithmic applications are discussed in a separate paper.

W.Auzinger, H.J.Stetter:

**An elimination algorithm for the computation of all zeros of a system of multivariate polynomial equations,**

Proceedings, `International Conference on Numerical Mathematics', Singapore, 1988. (Birkhäuser ISNM 86, 1988, 11-30.)

[preprint (pdf)]

__Abstract.__A direct numerical method is proposed for the determination of all isolated zeros of a system of multivariate polynomial equations. By "polynomial combination", the system is reduced to a special form which may be interpreted as a multiplication table for power products modulo the system. The zeros are then formed from an ordinary eigenvalue problem for the matrix of the multiplication table. Degenerate situations may be handled by perturbing them into general form and reaching the zeros of the unperturbed system via a homotopy method.

W.Auzinger:

**A quantitative discrete \( H^2 \)-regularity estimate for the Shortley-Weller scheme in convex domains,**

Numer. Math. 52(5) 523-537, 1988. (DOI 10.1007/BF01400890)

[preprint (pdf)]

__Abstract.__We investigate the discrete \( H^2 \)-regularity properties of the Shortley-Weller discretization of Poisson's equation (with homogeneous Dirichlet boundary conditions) in bounded convex domains \( \Omega \in {\mathbb R}^2 \). It is shown that the regularity constant is \(1\) independent of the mesh size \( h \) if the \( H^2 \)-seminorm is defined in a way assigning less weight to the (unsymmetric) differences near the boundary.

W.Auzinger, R.Frank, F.Macsek:

**On asymptotic expansions of the discretization error for stiff equations,**

in `Numerical Treatment of Differential Equations', Proceedings, NUMDIFF-4 (Ed. K.Strehmel), Halle, 1987. (Teubner-Texte zur Mathematik 104, 1988, 31-39.)

[preprint (pdf)]

__Abstract.__We discuss asymptotic expansions of the global discretization error for the implicit Euler scheme applied to stiff nonlinear initial value problems, giving a survey on recent results.

W.Auzinger, J.P.Monnet:

**IDeC-convergence independent of error asymptotics,**

BIT Numer. Math. 27(3), 350-367, 1987. (DOI 10.1007/BF01933730)

[preprint (pdf)]

__Abstract.__The present paper is concerned with the method of Iterated Defect Correction (IDeC) for two-point boundary value problems. We investigate the contractive behavior of the IDeC iteration in a completely discrete setting. Our results (which are a generalization of "classical" results based on asymptotic expansions of the discretization error) imply the stability of the collocation method which defines the fixed point of the IDeC iteration.

W.Auzinger:

**Defect correction for nonlinear elliptic difference equations,**

Numer. Math. 51(2), 199-208, 1987. (DOI 10.1007/BF01396749)

[preprint (pdf)]

__Abstract.__The present paper is concerned with the study of a high-order defect correction technique for discretizations of nonlinear elliptic boundary value problems. The convergence of the method is analyzed in general and, in more detail, for a specific example. The algorithmic combination of defect correction and multigrid techniques is also discussed.

W.Auzinger:

**Defect corrections for multigrid solutions of the Dirichlet problem in general domains,**

Math. Comp. 48, 471-484, 1987. (DOI 10.1090/S0025-5718-1987-0878685-7)

[article (pdf)]

__Abstract.__Recently, the technique of defect correction for the refinement of discrete solutions to elliptic boundary value problems has gained new acceptance in connection with the multigrid approach. In the present paper we give an analysis for a specific application, namely to finite-difference analogues of the Dirichlet problem for Helmholtz' equation, emphasizing the case of nonrectangular domains. A quantitative convergence proof is presented for a class of convex polygonal domains.

W.Auzinger, R.Frank, F.Macsek:

**Asymptotic error expansions for stiff equations. Part 2: The mildly stiff case,**

[Report No.68/86, Institut für Angewandte und Numerische Mathematik, TU Wien, 1986 (pdf)]

__Abstract.__We derive an asymptotic expansion for the global error of the implicit Euler scheme applied to stiff nonlinear systems of ordinary differential equations. Part 2, in particular, is devoted to the mildly stiff case. It is shown by a discrete singular perturbation analysis of the remainder term equation that the desired order \( O(h^{q+1}) \) of the remainder term, which inevitably breaks down at the first grid points, reappears at the subsequent grid points.

W.Auzinger, R.Frank, F.Macsek:

**Asymptotic error expansions for stiff equations. Part 1: The strongly stiff case,**

[Report No.67/86, Institut für Angewandte und Numerische Mathematik, TU Wien, 1986 (pdf)]

__Abstract.__In this paper we derive an asymptotic expansion for the global error of the implicit Euler scheme applied to stiff nonlinear systems of ordinary differential equations. Part 1 is devoted to the strongly stiff case. It is shown that in strongly stiff situations a full asymptotic error expansion exists at all grid points. Our analysis is based on singular perturbation techniques.

W.Auzinger, H.J.Stetter:

**Accurate arithmetic results for decimal data on non-decimal computers,**

Computing 35, 141-151, 1985. (DOI 10.1007/BF02260501)

[preprint (pdf)]

__Abstract.__Recently, techniques have been devised and implemented which permit the computation of smallest enclosing machine number intervals for the exact results of a good number of highly composite operations. These exact results refer, however, to the data as they are represented in the computer. This note shows how the conversion of decimal data into non-decimal representation may be joined with the mathematical operation on the data into one high-accuracy algorithm. Such an algorithm is explicitly presented for the solution of systems of linear equations.

W.Auzinger:

**DCMG01: A multigrid code with defect correction to solve \( \Delta U - c(x,y)U = f(x,y) \) on \( \Omega \), \( U = g(x,y) \) on \( \partial\Omega \), on nonrectangular bounded domains \( \Omega \) with high accuracy,**

Arbeitspapiere der GMD, Nr. 127, 1985.

[report (pdf)]

__Abstract.__DCMG01 is an extension of the Multigrid program MG01 to solve the Helmholtz equation with Dirichlet boundary conditions, \( \Delta U - c(x,y)U = f(x,y) \) (on \( \Omega \)), \( U = g(x,y) \) (on \( \partial\Omega \)), by difference methods in general bounded domains \( \Omega \in {\mathbb R}^2 \) with high accuracy. The extension is based on a*defect correction*principle. This report describes the properties of DCMG01 and how to use it. Numerical examples are also presented.

W.Auzinger:

**Defektkorrektur für Diskretisierungen des Dirichlet-Problems in allgemeinen Gebieten,**

Dissertation, TU Wien, 1984.

[text (pdf)]

__Abstract.__Defektkorrekturmethoden für die numerische Lösung von Randwertaufgaben bei partiellen Differentialgleichungen sind des öfteren vorgeschlagen und auch verwendet worden. Bei der Analyse dieser Verfahren ist man aus verschiedenen Gründen gezwungen, das Konvergenzverhalten in einem sehr weiten Rahmen (insbesondere ohne asymptotische Fehlerentwicklungen) zu untersuchen. Damit treten eine Reihe neuer Probleme auf, die bisher kaum zufriedenstellend behandelt worden sind.

In der vorliegenden Arbeit wird ein Defektkorrektur-Algorithmus für die Lösung der Helmholtz-Gleichung in allgemeinen ebenen Gebieten (mit Dirichlet-Randbedingung) vorgestellt. Der Diskretisierungsfehler ist \( O(h^4) \) trotz geringerer lokaler Konsistenz in der Nähe des Randes. Die Frage nach der Kontraktivität wird für eine recht grosse Klasse von Gebieten positiv beantwortet: Die Defektkorrektur kontrahiert gleichmäßig in \( h \) bezüglich einer diskreten Sobolev-Norm.

Diskutiert werden weiters die Bedeutung der Fehlerglättung und die Kombination der Defektkorrektur mit einem Multigrid-Solver.

W.Auzinger, H.J.Stetter:

**Defect corrections and multigrid iterations,**

in `Multigrid Methods', Proceedings, Köln-Porz, Nov. 1981, Lecture Notes in Mathematics 960, Springer, Berlin 1982.

[preprint (pdf)]

__Abstract.__At first we demonstrate how any defect correction iteration may be combined with an iterative procedure for its effective solution operator and derive a bound on the convergence factor for the combined iterative process. Then we consider the situation when the iterative solution procedure consists of multigrid cycles. Various implementation versions of a defect correction multigrid cycle are discussed, particularly for the case where defect correction is to achieve a higher order approximation to the differential equation. We also regard some algorithmic and quantitative details in the context of model problems. Finally we present some numerical examples.

[> top of page (recent)]