Search

It is well known that various dynamical systems preserve volume in phase space such as all Hamiltonian systems. This qualitative geometrical property of the analytical solution should be respected in the sense of Geometric Integration. This paper studies the volume-preserving property of exponential integrators for different vector fields. For exponential integrators, we first derive a necessary and sufficient condition of volume preservation. Then based on this condition, volume-preserving exponential integrators are discussed in detail for four kinds of vector fields. It is shown that symplectic exponential integrators can be volume preserving for a much larger class of vector fields than Hamiltonian systems. On the basis of the analysis, some applications of volume-preserving exponential integrators are discussed. For solving highly oscillatory second-order systems, novel volume-preserving exponential integrators are derived, and for separable partitioned systems, extended Runge–Kutta–Nystrom (ERKN) integrators of volume preservation are presented. Moreover, the volume preservation of Runge–Kutta–Nystrom (RKN) methods is also discussed. Five illustrative numerical experiments are carried out to demonstrate the notable superiority of volume-preserving exponential integrators in comparison with volume-preserving Runge-Kutta methods.

Motivated by the recent paper (Wang and Wu, 2018 [12] ), where the authors proposed a class of functionally-fitted energy-preserving integrators for Poisson systems, in this paper we aim at developing exponential integrators with quadratic energy preservation for the important class of linear Poisson systems. By rewriting the linear Poisson system in a new formulation, we apply the symplectic exponential Runge-Kutta methods and thus obtain a novel class of exponential integrators for the linear Poisson system. It is proved that the obtained exponential integrator can exactly preserve the quadratic energy. In comparison with the functionally-fitted energy-preserving methods whose numerical performance depends on the fitted frequency, the new exponential integrators are independent of the fitted frequency and convenient to use. Furthermore, we investigate the convergence of the fixed-point iteration used in the implementation of the proposed exponential integrator, and show that the convergence region of the stepsize is independent of the norm of the skew-symmetric matrix. Finally, numerical results with the Euler equations of a free rigid body soundly support the remarkable performance of the exponential integrator in both the global errors and the energy preservation.

In this paper, we consider exponential integrators for semilinear Poisson systems. Two types of exponential integrators are constructed, one preserves the Poisson structure, and the other preserves energy. Numerical experiments for semilinear Possion systems obtained by semi-discretizing Hamiltonian PDEs are presented. These geometric exponential integrators exhibit better long time stability properties as compared to non-geometric integrators, and are computationally more efficient than traditional symplectic integrators and energy-preserving methods based on the discrete gradient method., Comment: 18 pages, 11 figures

This paper presents a new algorithm KIOPS for computing linear combinations of φ -functions that appear in exponential integrators. This algorithm is suitable for large-scale problems in computational physics where little or no information about the spectrum or norm of the Jacobian matrix is known a priori . We first show that such problems can be solved efficiently by computing a single exponential of a modified matrix. Then our approach is to compute an appropriate basis for the Krylov subspace using the incomplete orthogonalization procedure and project the matrix exponential on this subspace. We also present a novel adaptive procedure that significantly reduces the computational complexity of exponential integrators. Our numerical experiments demonstrate that KIOPS outperforms the current state-of-the-art adaptive Krylov algorithm phipm .

Rational exponential integrators (REXI) are a class of numerical methods that are well suited for the time integration of linear partial differential equations with imaginary eigenvalues. Since these methods can be parallelized in time (in addition to the spatial parallelization that is commonly performed) they are well suited to exploit modern high performance computing systems. In this paper, we propose a novel REXI scheme that drastically improves accuracy and efficiency. The chosen approach will also allow us to easily determine how many terms are required in the approximation in order to obtain accurate results. We provide comparative numerical simulations for a shallow water equation that highlight the efficiency of our approach and demonstrate that REXI schemes can be efficiently implemented on graphic processing units., 25 pages, 5 figures, submitted to Journal of Computational Physics

The nonlinear Schrodinger equation (NLSE) is widely used in telecommunication applications, since it allows one to describe the propagation of pulses in an optical fiber. Recently some new approaches based on the nonlinear Fourier transform (NFT) have been actively explored to compensate for fiber nonlinearity and to exceed the limitations of nonlinearity-imposed limits of linear transmission methods. Despite the fact that the numerical solution of NLSE is a general problem, nevertheless, the optical community has been focusing on this issue. Improving the accuracy of the NFT algorithms remains an urgent problem in optics. In particular, it is important to increase the approximation order of the methods, especially in problems where it is necessary to analyze the structure of complex waveforms. To correctly describe them and their spectral parameters, more accurate and fast numerical methods are needed. We propose a novel general approach for constructing sixth-order (with respect to an integration step) finite-difference schemes for first-order linear differential systems. These schemes are based on the generalized Cayley transform and include exponential integrators as a special case. If the system has a time-dependent skew-hermitian matrix then the schemes conserve the quadratic first integral automatically. Then we apply our method to solve the direct spectral problem for the Zakharov-Shabat system. New schemes with fractional rational transition matrix allow the use of fast algorithms to solve the initial problem for a large number of values of the spectral parameter.

We present a new symmetric linearly implicit exponential integrator that preserves the polynomial first integrals or the Lyapunov functions for the conservative and dissipative stiff equations, respectively. The method is tested by both oscillated ordinary differential equations and partial differential equations, e.g., an averaged system in wind-induced oscillation, the Fermi–Pasta–Ulam systems, and the polynomial pendulum oscillators. The numerical simulations confirm the conservative properties of the proposed method and demonstrate its good behavior in superior running speed when compared with fully implicit schemes for long-time simulations.

The efficient numerical solution of many kinetic models in plasma physics is impeded by the stiffness of these systems. Exponential integrators are attractive in this context as they remove the CFL condition induced by the linear part of the system, which in practice is often the most stringent stability constraint. In the literature, these schemes have been found to perform well, e.g., for drift-kinetic problems. Despite their overall efficiency and their many favorable properties, most of the commonly used exponential integrators behave rather erratically in terms of the allowed time step size in some situations. This severely limits their utility and robustness. Our goal in this paper is to explain the observed behavior and suggest exponential methods that do not suffer from the stated deficiencies. To accomplish this we study the stability of exponential integrators for a linearized problem. This analysis shows that classic exponential integrators exhibit severe deficiencies in that regard. Based on the analysis conducted we propose to use Lawson methods, which can be shown not to suffer from the same stability issues. We confirm these results and demonstrate the efficiency of Lawson methods by performing numerical simulations for both the Vlasov–Poisson system and a drift-kinetic model of a ion temperature gradient instability.

Taking into account the limited accuracy of the energy-preserving exponential integrator of order two (Li and Wu, 2016 [29] ) for conservative or dissipative systems with highly oscillatory solutions, this paper is devoted to presenting a uniform framework to design energy-preserving exponential integrators of arbitrarily high order based on the modifying integrator theory. To this end, we first show that the second-order energy-preserving exponential integrator is a B-series method. Using the adapted substitution law, we then prove that there exist arbitrary order energy-preserving exponential integrators and show how to design arbitrarily high-order integrators by finding the truncated modified differential equations. As an example, the fourth-order energy-preserving exponential integrator is constructed in detail. The stability and convergence of the proposed integrators are analyzed as well. Finally, numerical experiments are accompanied, including both ODEs and PDEs, and the numerical results demonstrate the remarkable superiority over the existing energy-preserving integrators for highly oscillatory systems in the literature.

Exponential and IMEX integrators are considered for wave turbulence problems, which are stiff systems that often require long-time simulations. Applying a multi-scale asymptotic analysis to the difference equations for a simple wave turbulence problem reveals that Exponential Runge-Kutta methods incorrectly model interactions between resonant waves unless they are specifically designed otherwise, and IMplicit-EXplicit (IMEX) integrators introduce dispersive and dissipative errors that severely degrade the accuracy of slow energy exchange between resonant waves. The same analysis shows Integrating Factor (IF) methods to be the most suitable. Numerical experiments performed on the resonant triad model and the Majda-McLaughlin-Tabak (MMT) model verify the discrete-asymptotic analysis. In the MMT problem, Exponential Time Differencing methods require a time step more than two orders of magnitude smaller than IF methods to achieve comparable accuracy, and IMEX methods require more than six orders of magnitude smaller time steps. A novel near-minimax rational approximation of the matrix exponential is proposed to be used for integrating factor methods, and is demonstrated to be as accurate as the IF method but without the need to compute a matrix exponential.

We propose two new classes of time integrators for stiff DEs: the implicit exponential (IMEXP) and the hybrid exponential methods. In contrast to the existing exponential schemes, the new methods offer significant computational advantages when used with preconditioners. Any preconditioner can be used with any of these new schemes. This leads to a broader applicability of exponential methods. The proof of convergence of these integrators and numerical demonstration of their efficiency are presented.

The exponential propagation methods were applied in the past for accurate integration of the shallow water equations on the sphere. Despite obvious advantages related to the exact solution of the linear part of the system, their use for the solution of practical problems in geophysics has been limited because efficiency of the traditional algorithm for evaluating the exponential of Jacobian matrix is inadequate. In order to circumvent this limitation, we modify the existing scheme by using the Incomplete Orthogonalization Method instead of the Arnoldi iteration. We also propose a simple strategy to determine the initial size of the Krylov space using information from previous time instants. This strategy is ideally suited for the integration of fluid equations where the structure of the system Jacobian does not change rapidly between the subsequent time steps. A series of standard numerical tests performed with the shallow water model on a geodesic icosahedral grid shows that the new scheme achieves efficiency comparable to the semi-implicit methods. This fact, combined with the accuracy and the mass conservation of the exponential propagation scheme, makes the presented method a good candidate for solving many practical problems, including numerical weather prediction.

The structural flexibility of the exponential propagation iterative methods of Runge-Kutta type (EPIRK) enables construction of particularly efficient exponential time integrators. While the EPIRK methods have been shown to perform well on stiff problems, all of the schemes proposed up to now have been derived using classical order conditions. In this paper we extend the stiff order conditions and the convergence theory developed for the exponential Rosenbrock methods to the EPIRK integrators. We derive stiff order conditions for the EPIRK methods and develop algorithms to solve them to obtain specific schemes. Moreover, we propose a new approach to constructing particularly efficient EPIRK integrators that are optimized to work with an adaptive Krylov algorithm. We use a set of numerical examples to illustrate the computational advantages that the newly constructed EPIRK methods offer compared to previously proposed exponential integrators.

In this paper, new Locally Linearized Runge-Kutta formulas of Dormand and Prince are introduced with the purpose of integrating large systems of initial value problems. The new formulas are obtained by replacing the Pade approximation for matrix exponential by a Krylov-Pade approximation in the embedded formulas proposed in a previous work. Unlike other high order exponential integrators, these new formulas involve the approximation of just a single phi-function times vector, which makes them particularly attractive from a computational viewpoint. With these new formulas, adaptive schemes are constructed with novelties in the approximation of phi-functions times vectors by Krylov-Pade approximations, and in the strategies for controlling the approximation errors, for estimating the dimension of the Krylov subspaces, and for the reuse of the Jacobians. In addition, numerical experiments are performed to show the potential of the new embedded formulas and adaptive codes in the integration of known physical, biophysical, and physical-chemistry models.

In this paper, we propose a semi-Lagrangian discontinuous Galerkin method coupled with Runge-Kutta exponential integrators (SLDG-RKEI) for nonlinear Vlasov dynamics. The commutator-free Runge-Kutta (RK) exponential integrators (EI) were proposed by Celledoni, et al. (FGCS, 2003). In the nonlinear transport setting, the RKEI can be used to decompose the evolution of the nonlinear transport into a composition of a sequence of linearized dynamics. The resulting linearized transport equations can be solved by the semi-Lagrangian (SL) discontinuous Galerkin (DG) method proposed in Cai, et al. (JSC, 2017). The proposed method can achieve high order spatial accuracy via the SLDG framework, and high order temporal accuracy via the RK EI. Due to the SL nature, the proposed SLDG-RKEI method is not subject to the CFL condition, thus they have the potential in using larger time-stepping sizes than those in the Eulerian approach. Inheriting advantages from the SLDG method, the proposed SLDG-RKEI schemes are mass conservative, positivity-preserving, have no dimensional splitting error, perform well in resolving complex solution structures, and can be evolved with adaptive time stepping sizes. We show the performance of the SLDG-RKEI algorithm by classical test problems for the nonlinear Vlasov-Poisson system, as well as the Guiding center Vlasov model. Though that it is not our focus of this paper to explore the SLDG-RKEI scheme for nonlinear hyperbolic conservation laws that develop shocks, we show some preliminary results on schemes' performance on the Burgers' equation.

形態: カラー図版あり, Physical characteristics: Original contains color illustrations, Accepted: 2016-06-14, 資料番号: PA1710033000

We have introduced a new explicit numerical method, based on a discrete stochastic process, for solving a class of fractional partial differential equations that model reaction subdiffusion. The scheme is derived from the master equations for the evolution of the probability density of a sum of discrete time random walks. We show that the diffusion limit of the master equations recovers the fractional partial differential equation of interest. This limiting procedure guarantees the consistency of the numerical scheme. The positivity of the solution and stability results are simply obtained, provided that the underlying process is well posed. We also show that the method can be applied to standard reaction-diffusion equations. This work highlights the broader applicability of using discrete stochastic processes to provide numerical schemes for partial differential equations, including fractional partial differential equations.

In this paper we investigate the solution of boundary value problems on polygonal domains for elliptic partial differential equations. We observe that when the problems are formulated as the boundary integral equations of classical potential theory, the solutions are representable by series of elementary functions. In addition to being analytically perspicuous, the resulting expressions lend themselves to the construction of accurate and efficient numerical algorithms. The results are illustrated by a number of numerical examples.

Despite numerous computational advances over the last few decades, molecular dynamics still favors explicit (and thus easily-parallelizable) time integrators for large scale numerical simulation. As a consequence, computational efficiency in solving its typically stiff oscillatory equations of motion is hampered by stringent stability requirements on the time step size. In this paper, we present a semi-analytical integration scheme that offers a total speedup of a factor 30 compared to the Verlet method on typical MD simulation by allowing over three orders of magnitude larger step sizes. By efficiently approximating the exact integration of the strong (harmonic) forces of covalent bonds through matrix functions, far improved stability with respect to time step size is achieved without sacrificing the explicit, symplectic, time-reversible, or fine-grained parallelizable nature of the integration scheme. We demonstrate the efficiency and scalability of our integrator on simulations ranging from DNA strand unbinding and protein folding to nanotube resonators.

This article presents explicit exponential integrators for stochastic Maxwell's equations driven by both multiplicative and additive noises. By utilizing the regularity estimate of the mild solution, we first prove that the strong order of the numerical approximation is $\frac 12$ for general multiplicative noise. Combing a proper decomposition with the stochastic Fubini's theorem, the strong order of the proposed scheme is shown to be $1$ for additive noise. Moreover, for linear stochastic Maxwell's equation with additive noise, the proposed time integrator is shown to preserve exactly the symplectic structure, the evolution of the energy as well as the evolution of the divergence in the sense of expectation. Several numerical experiments are presented in order to verify our theoretical findings., Comment: 21 Pages

In this paper we propose a space-time adaptive finite element method for the phase field model of pitting corrosion, which is a parabolic partial differential equation system consisting of a phase variable and a concentration variable. A major challenge in solving this phase field model is that the problem is very stiff, which makes the time step size extremely small for standard temporal discretizations. Another difficulty is that a high spatial resolution is required to capture the steep gradients within the diffused interface, which results in very large number of degrees of freedom for uniform meshes. To overcome the stiffness of this model, we combine the Rosenbrock–Euler exponential integrator with Crank–Nicolson scheme for the temporal discretization. Moreover, by exploiting the fact that the speed of the corroding interface decreases with time, we derive an adaptive time stepping formula. For the spatial approximation, we propose a simple and efficient strategy to generate adaptive meshes that reduce the computational cost significantly. Thus, the proposed method utilizes local adaptivity and mesh refinement for efficient simulation of the corrosive dissolution over long times in heterogeneous media with complex microstructures. We also present an extensive set of numerical experiments in both two and three dimensional spaces to demonstrate efficiency and robustness of the proposed method.

We develop a new matrix-based approach to the numerical solution of partial differential equations (PDE) and apply it to the numerical solution of partial fractional differential equations (PFDE). The proposed method is to discretize a given PFDE as a Sylvester Equation, and parameterize the integral surface using matrix algebra. The combination of these two notions results in an algorithm which can solve a general class of PFDE efficiently and accurately by means of an O ( n 3 ) algorithm for solving the Sylvester Matrix Equation (over an m i? n grid with m ~ n ). The proposed parametrization of the integral surface allows for the solution with the more general Robin boundary conditions, and allows for high-order approximations to derivative boundary conditions. To achieve our ends, we also develop a new matrix-based approximation to fractional order derivatives. The proposed method is demonstrated by the numerical solution of the fractional diffusion equation with fractional derivatives in both the temporal and spatial directions.

This paper proposes an approach for high-order time integration within a multi-domain setting for time- fractional differential equations. Since the kernel is singular or nearly singular, two main difficulties arise after the domain decomposition: how to properly account for the history/memory part and how to perform the integration accurately. To address these issues, we propose a novel hybrid approach for the numerical integration based on the combination of three-term-recurrence relations of Jacobi polynomials and high-order Gauss quadrature. The different approximations used in the hybrid approach are justified theoretically and through numerical examples. Based on this, we propose a new multi-domain spectral method for high-order accurate time integrations and study its stability properties by identifying the method as a generalized linear method. Numerical experiments confirm hp-convergence for both time-fractional differential equations and time-fractional partial differential equations.

Multiple time-stepping (MTS) algorithms allow to efficiently integrate large systems of ordinary differential equations, where a few stiff terms restrict the timestep of an otherwise non-stiff system. In this work, we discuss a flexible class of MTS techniques, based on multistep methods. Our approach contains several popular methods as special cases and it allows for the easy construction of novel and efficient higher-order MTS schemes. In addition, we demonstrate how to adapt the stability contour of the non-stiff time-integration to the physical system at hand. This allows significantly larger timesteps when compared to previously known multistep MTS approaches. As an example, we derive novel predictor-corrector (PCMTS) schemes specifically optimized for the time-integration of damped wave equations on locally refined meshes. In a set of numerical experiments, we demonstrate the performance of our scheme on discontinuous Galerkin time-domain (DGTD) simulations of Maxwell's equations.

We propose a new, energy conserving, spectral element, discontinuous Galerkin method for the approximation of the Vlasov–Poisson system in arbitrary dimension, using Cartesian grids. The method is derived from the one proposed in [4] , with two modifications: energy conservation is obtained by a suitable projection operator acting on the solution of the Poisson problem, rather than by solving multiple Poisson problems, and all the integrals appearing in the finite element formulation are approximated with Gauss–Lobatto quadrature, thereby yielding a spectral element formulation. The resulting method has the following properties: exact energy conservation (up to errors introduced by the time discretization), stability (thanks to the use of upwind numerical fluxes), high order accuracy and high locality. For the time discretization, we consider both Runge–Kutta methods and exponential integrators, and show results for 1D and 2D cases (2D and 4D in phase space, respectively).

In this paper the framework of the exponential propagation iterative methods of Runge–Kutta type (EPIRK) is extended to construct split EPIRK (sEPIRK) integrators for semilinear systems of ODEs. The structure of the sEPIRK methods possesses the flexibility and generality that allows construction of very efficient schemes. We demonstrate how bicolored trees-based B-series can be used to derive the order conditions for the new integrators. An algorithm is developed to solve the order conditions up to order three and several new schemes with desirable properties are proposed. The numerical results illustrate the advantages offered by the new class of integrators. The experiments also address the comparative performance of split vs. non-split EPIRK methods and the question of improving efficiency by optimizing coefficients of the sEPIRK schemes. It is shown that specific schemes can be custom-built to improve computational efficiency for a particular problem.

In this paper a method for the numerical solution of distributed order FDEs (fractional differential equations) of a general form is presented. The method applies to both linear and nonlinear equations. The Caputo type fractional derivative is employed. The distributed order FDE is approximated with a multi-term FDE, which is then solved by adjusting appropriately the numerical method developed for multi-term FDEs by Katsikadelis. Several example equations are solved and the response of mechanical systems described by such equations is studied. The convergence and the accuracy of the method for linear and nonlinear equations are demonstrated through well corroborated numerical results.

Partial differential equations describing the patterning of thin crystalline films are typically of fourth or sixth order, they are quasi- or semilinear and they are mostly defined on simple geometries such as rectangular domains. For the numerical simulation of these kinds of problems spectral methods are an efficient approach. We apply several implicit-explicit schemes to one recently derived PDE that we express in terms of coefficients of trigonometric interpolants. While the simplest IMEX scheme turns out to have the mildest step-size restriction, higher order SBDF schemes tend to be more unstable and exponential time integrators are fastest for the calculation of very accurate solutions. We implemented a reduced model in the EXPINT package syntax [3] and compared various exponential schemes. A convexity splitting approach was employed to stabilize the SBDF1 scheme. We show that accuracy control is crucial when using this idea, therefore we present a time-adaptive SBDF1/SBDF1-2-step method that yields convincing results reflecting the change in timescales during topological changes of the nanostructures. The implementation of all presented methods is carried out in Matlab. We used the open source GPUmat package to gain up to 5-fold runtime benefits by carrying out calculations on a low-cost GPU without having to prescribe any knowledge in low-level programming or CUDA implementations and found comparable speedups as with Matlab's PCT or with GPUmat run on Octave.

Overlapping domain decomposition methods, otherwise known as overset grid or chimera methods, are useful for simplifying the discretization of partial differential equations in or around complex geometries. Though in wide use, such methods are prone to numerical instability unless numerical diffusion or some other form of regularization is used, especially for higher-order methods. To address this shortcoming, high-order, provably energy stable, overlapping domain decomposition methods are derived for hyperbolic initial boundary value problems. The overlap is treated by splitting the domain into pieces and using generalized summation-by-parts derivative operators and polynomial interpolation. New implicit and explicit operators are derived that do not require regularization for stability in the linear limit. Applications to linear and nonlinear problems in one and two dimensions are presented, where it is found the explicit operators are preferred to the implicit ones.

We propose a new class of the exponential propagation iterative methods of Runge-Kutta-type (EPIRK). The EPIRK schemes are exponential integrators that can be competitive with explicit and implicit methods for integration of large stiff systems of ODEs. Introducing the new, more general than previously proposed, ansatz for EPIRK schemes allows for more flexibility in deriving computationally efficient high-order integrators. Recent extension of the theory of B-series to exponential integrators [1] is used to derive classical order conditions for schemes up to order five. An algorithm to systematically solve these conditions is presented and several new fifth order schemes are constructed. Several numerical examples are used to verify the order of the methods and to illustrate the performance of the new schemes.

Implicit-explicit (IMEX) multistep methods are very useful for the time discretization of convection diffusion PDE problems such as the Burgers equations and the incompressible Navier-Stokes equations. In the latter as well as in PDE models of plasma physics and of electromechanical systems, semi-discretization in space gives rise to differential-algebraic (DAE) system of equations often of index higher than 1. In this paper we propose a new class of exponential integrators for index 2 DAEs arising from the semi-discretization of PDEs with a dominating and typically nonlinear convection term. This class of problems includes the incompressible Navier-Stokes equations. The integration methods are based on the backward differentiation formulae (BDF) and they can be applied without modifications in the semi-Lagrangian integration of convection diffusion problems. The approach gives improved performance at low viscosity regimes.

We present an exponential time integrator in conjunction with a finite volume discretisation in space for simulating transport by advection and diffusion including chemical reactions in highly heterogeneous porous media representative of geological reservoirs. These numerical integrators are based on the variation of constants solution and solving the linear system exactly. This is at the expense of computing the exponential of the stiff matrix comprising the finite volume discretisation. Using real Leja points or a Krylov subspace technique compared to standard finite difference-based time integrators. We observe for a variety of example applications that numerical solutions with exponential methods are generally more accurate and require less computational cost. They hence comprise an efficient and accurate method for simulating non-linear advection-dominated transport in geological formations.

A two-dimensional fluid laser-plasma integrator for stratified plasma-vacuum systems is presented. Inside a plasma, a laser pulse can be longitudinally compressed from ten or more wave-lengths to one or two cycles. However, for physically realistic simulations, transversal effects have to be included, because transversal instabilities can destroy the pulse and transversal compression in the plasma as well as focusing in vacuum allows much higher intensities to be reached. In contrast to the one-dimensional case, where a two-step implementation of the Gautschi-type exponential integrator with constant step-size turned out to be sufficient, it is essential to enable changes of the time step-size for the two-dimensional case. The use of a one-step version of the Gautschi-type integrator, being accurate of second order independent of the highest frequencies arising in the system, is proposed. In vacuum this allows to take arbitrarily large time-steps. To optimize runtime and memory requirements within the plasma, a splitting of the Laplacian is suggested. This splitting allows to evaluate the matrix functions arising in the Gautschi-type method by one-dimensional Fourier transforms. It is also demonstrated how the different variants of the scheme can be parallelized. Numerical experiments illustrate the superior performance and accuracy of the integrator compared to the standard leap-frog method. Finally, we discuss the simulation of a layered plasma-vacuum structure using the new method.

Fractional order partial differential equations, as generalizations of classical integer order partial differential equations, are increasingly used to model problems in fluid flow, finance and other areas of application. In this paper we discuss a practical alternating directions implicit method to solve a class of two-dimensional initial-boundary value fractional partial differential equations with variable coefficients on a finite domain. First-order consistency, unconditional stability, and (therefore) first-order convergence of the method are proven using a novel shifted version of the classical Grunwald finite difference approximation for the fractional derivatives. A numerical example with known exact solution is also presented, and the behavior of the error is examined to verify the order of convergence.

A multidomain Legendre pseudospectral method is developed for the solution of linear hyperbolic initial boundary value problems, with mixed boundary conditions, in general two-dimensional and axisymmetric geometries. A weak collocation spectral method is utilized for the spatial approximation of a generic wave evolution equation over multiple nonoverlapping subdomains. The system of ordinary differential equations that stems from the above procedure is integrated in time by implicit as well as explicit high order temporal approximation algorithms. The weak formalism of the influence matrix method is combined with the implicit approximation, so as to efficiently solve the coupled system of linear equations after the full discretization, while a novel technique for avoiding the amplification of roundoff error at high temporal resolution simulations with the implicit temporal integration methods, is also studied. An innovative method for the treatment of Dirichlet boundary conditions is proposed, in order to avoid the order reduction which usually arises with the utilization of the explicit time integrator. Furthermore, appropriate modifications are reported, for dealing with the pole singularity problem faced by the weak formulation of axisymmetric problems. Finally, numerical simulations of a variety of wave problems on curvilinear geometries and unstructured subdomain configurations are presented in order to assess the capabilities of the proposed methodology in handling efficiently general hyperbolic differential operators.

Anomalous diffusion is a possible mechanism underlying plasma transport in magnetically confined plasmas. To model this transport mechanism, fractional order space derivative operators can be used. Here, the numerical properties of partial differential equations of fractional order α, 1 ≤ α ≤ 2, are studied. Two numerical schemes, an explicit and a semi-implicit one, are used in solving these equations. Two different discretization methods of the fractional derivative operator have also been used. The accuracy and stability of these methods are investigated for several standard types of problems involving partial differential equations of fractional order.

We develop a class of numerical methods for stiff systems, based on the method of exponential time differencing. We describe schemes with second- and higher-order accuracy, introduce new Runge?Kutta versions of these schemes, and extend the method to show how it may be applied to systems whose linear part is nondiagonal. We test the method against other common schemes, including integrating factor and linearly implicit methods, and show how it is more accurate in a number of applications. We apply the method to both dissipative and dispersive partial differential equations, after illustrating its behavior using forced ordinary differential equations with stiff linear parts.

A quasi-steady-state method is presented that integrates stiff differential equations arising from reaction kinetics. This predictor–corrector method is A-stable for linear equations and second-order accurate. The method is used for all species regardless of the time scales of the individual equations, and it works well for problems typical of hydrocarbon combustion. Start-up costs are low, making the method ideal for use in process-split reacting-flow simulations which require the solution of an initial-value problem in every computational cell for every global time step. The algorithm is described, and error analysis and linear stability analysis are included. The algorithm is also applied to several test problems, and the results are compared to those of the stiff integrator CHEMEQ. The method, which we call α-QSS, is more stable, more accurate, and less costly than CHEMEQ.

Based on a genuine multidimensional numerical scheme, called the Method of Transport, we derive a form of the compressible Euler equations, capable of a linearization for any space dimension. This form enables a rigorous error analysis of the linearization error without the knowledge of the numerical method used to solve the linear equations. The generated error can be eliminated by special correction terms in the linear equations. Hence, existing scalar high order methods can be used to solve the linear equations and obtain high order accuracy in space and time for the non-linear conservation law. In this approach, the scalar version of the method of transport is used to solve the linear equations. This method is multidimensional and reduces the solution of the partial differential equation to an integration process. Convergence histories presented at the end of the paper show that the numerical results agree with the theoretical predictions.

Linear mode conversion, an important issue in the physics of plasma waves, involves an ordinary differential equation of at least fourth order. Attempts to integrate this differential equation using conventional numerical methods typically fail to provide a physically sensible result. This failure occurs because small, unavoidable truncation and differencing errors excite nonphysical, mathematically allowed exponentially growing solutions which quickly overwhelm the desired physical solution. We present here a two-point boundary numerical method which avoids exciting these unwanted nonphysical solutions and so provides solutions that are physically significant. The numerical algorithm is a generalization of the standard tridiagonal method and provides single-pass solutions which satisfy the original difference equations to within the numerical accuracy of the computer.

An unconventional numerical method for solving a restrictive and yet often-encountered class of ordinary differential equations is proposed. The method has a crucial, what we callreflexive, property and requires solving one linear system per time-step, but is second-order accurate. A systematical and easily implementable scheme is proposed to enhance the computational efficiency of such methods whenever needed. Applications are reported on how the idea can be applied to solve the Korteweg?de Vries Equation discretized in space.

An eighth-order P-stable two-step method for the numerical integration of second-order periodic initial-value problems is developed in this paper. This method has seven stages per iteration and an interval of periodicity equal to (0, ∞); i.e., it is P-stable. Numerical and theoretical results obtained for several well-known problems show the efficiency of the new method.

A multilevel wavelet collection method for the solution of partial differential equations is developed. Two different approaches of treating general boundary conditions are suggested. Both are based on the wavelet interpolation technique developed in the present research. The first approach uses wavelets as a basis and results in a differential-algebraic system of equations, where the algebraic part arises from boundary conditions. The second approach utilizes extended wavelets, which satisfy boundary conditions exactly. This approach results in a system of coupled ordinary differential equations. The method is tested on the one-dimensional Burgers equation with small viscosity. The solutions are compared with those resulting from the use of other numerical algorithms. The present results indicate that the method is competitive with well-established numerical algorithms.