Adaptive Step Size Control Research Articles

On the adaptive solution of phase‐field problems with A‐stable explicit last‐stage diagonally implicit Runge–Kutta (ELDIRK) methods

AbstractRunge–Kutta algorithms with adaptive step size control provide reliable tools for the solution of initial value problems with diagonally implicit Runge‐Kutta (DIRK) methods as the most common approach. In this contribution, the new low‐order explicit last‐stage diagonally implicit Runge–Kutta (ELDIRK) methods are investigated, combining implicit schemes with an additional explicit evaluation as an explicit last stage. This results in Butcher tableaus with two solutions of different convergence orders suitable for embedded methods, where the higher‐order solution is achieved by additional explicit evaluations. Thus, the iterative solution of non‐linear systems is omitted for the additional stage, presenting a major reduction in computational cost for the determination of a local error estimate. The key contribution is the application of the novel Butcher tableaux to phase‐field problems, solved with the finite‐element method, leading to substantial numerical investigations with an efficient approach for DIRK schemes. The most important aspects are the investigations of stability properties which lead to the novel class of A‐stable ELDIRK methods. In accordance, the study of the convergence orders is presented. A local error estimator is presented, such that adaptive step size control for the new low‐order embedded schemes based on an empirical approach for error estimation is achieved. A suitable parallel algorithm is presented with conclusive two‐dimensional phase‐field simulations based on a Kobayashi‐Warren‐Carter model. The higher‐order convergence suggested by the novel schemes is confirmed, and their effective results are demonstrated, resulting in a valuable semi‐explicit addition to the family of Runge–Kutta time integration schemes.

FLARE: field line analysis and reconstruction for 3D boundary plasma modeling

Abstract The FLARE code is a magnetic mesh generator that is integrated within a suite of tools for the analysis of the magnetic geometry in toroidal fusion devices. A magnetic mesh is constructed from field line segments and permits fast reconstruction of field lines in 3D boundary plasma codes such as EMC3-EIRENE. Both intrinsically non-axisymmetric configurations (stellarators) and those with symmetry breaking perturbations of an axisymmetric equilibrium (tokamaks) are supported. The code itself is written in Modern Fortran with MPI support for parallel computing, and it incorporates object-oriented programming for the definition of the magnetic field and the material surface geometry. Extended derived types for a number of different magnetohydrodynamic equilibrium and plasma response models are implemented. The core element of FLARE is a field line tracer with adaptive step-size control, and this is integrated into tools for the construction of Poincaré maps and invariant manifolds of X-points. A collection of high-level procedures that generate output files for visualization is build on top of that. The analysis modules are build with Python frontends that facilitate customization of tasks and/or scripting of parameter scans.

A novel chaos control method based on conjugate direction and adaptive step size

The Hasofer-Lind-Rackwitz-Fiessler method (HL-RF method) with the negative gradient direction as the iteration direction often oscillates or even does not converge during the iteration process due to the unreasonable control of the negative gradient direction. This greatly reduces the convergence and robustness of the HL-RF method. Chaos control method (CC method) is also a method with negative gradient direction as the iterative direction. However, in order to alleviate the oscillation problem caused by unreasonable control of negative gradient direction, it introduces a smaller iterative step size in the iterative process. Because its iterative step size is a fixed small value, the iterative efficiency of the CC method is not high. In view of the above problems, in order to reduce the iterative oscillation phenomenon and increase the iterative efficiency at the same time, it proposes a novel chaos control method based on conjugate direction and adaptive step size (CDCC method). Firstly, it uses the conjugate direction instead of the negative gradient direction, so as to fundamentally avoid the problem that the negative gradient direction is easy to cause oscillation. It effectively increases the convergence and robustness of the method. Secondly, it uses iteration information to adaptively control the iteration step size, increasing the computational efficiency of the method effectively. In addition, it also discusses the optimal value of the modification coefficient in the adaptive step size control. Finally, eight examples are given to show that the proposed method has obvious advantages in computational efficiency, convergence and robustness.

Numerical Analysis of Free Play-Induced Aeroelastic Phenomena: A Numerical Approach With Adaptive Step Size Control

This study presents a detailed numerical analysis of nonlinear aeroelastic behavior in a two degree of freedom (DOF) model, focusing on plunge and pitch motions and employing the continuation method (CM) with an adaptive step size control algorithm. The research incorporates free-play nonlinearity at the plunge hinge, a common structural nonlinearity in aeronautics that can induce detrimental limit cycle oscillations (LCOs) during flight. By examining three scenarios—linear response, unhindered plunge motion, and nonlinear stiffness behavior—the study assesses the effects of free play on flutter and LCO phenomena, including discontinuity-induced bifurcations like grazing bifurcation. Additionally, the study explores parameter variation for nonlinear flutter analysis, revealing the dynamics of grazing bifurcation and its impact on LCO behavior. The research also demonstrates the method’s superior accuracy in flutter speed estimation and mode-switching identification, despite higher computational demands. The findings underscore the diminishing influence of nonlinear free-play behavior on LCO amplitude, providing insights with significant implications for aeroelastic design and aircraft safety.

Adaptive step-size control for global approximation of SDEs driven by countably dimensional Wiener process

In this paper, we deal with the global approximation of solutions of stochastic differential equations (SDEs) driven by countably dimensional Wiener process. Under certain regularity conditions imposed on the coefficients, we show lower bounds for exact asymptotic error behaviour. For that reason, we analyse separately two classes of admissible algorithms: based on equidistant, and possibly not equidistant meshes. Our results indicate that in both cases, decrease of any method error requires significant increase of the cost term, which is illustrated by the product of cost and error diverging to infinity. This is, however, not visible in the finite-dimensional case. In addition, we propose an implementable, path-independent Euler algorithm with adaptive step-size control, which is asymptotically optimal among algorithms using specified truncation levels of the underlying Wiener process. Our theoretical findings are supported by numerical simulation in Python language.

On the accuracy, stability and computational efficiency of explicit last-stage diagonally implicit Runge–Kutta methods (ELDIRK) for the adaptive solution of phase-field problems

Runge–Kutta algorithms with adaptive step size control provide reliable tools for the solution of initial value problems. Most widespread within the vast family of Runge–Kutta methods are the diagonally implicit Runge–Kutta (DIRK) methods. In particular, the new low-order explicit last-stage diagonally implicit Runge–Kutta (ELDIRK) methods are investigated, combining implicit schemes with an additional explicit evaluation as an explicit last stage. This results in Butcher tableaus with two solutions of different convergence orders suitable for embedded methods, where the higher-order solution is achieved by additional explicit evaluations. Thus, the iterative solution of non-linear systems is omitted for the additional stage, presenting a major reduction in computational cost for the determination of a local error estimate. The key contribution is the application of the novel Butcher tableaus to phase-field problems, solved with the finite-element method, leading to substantial numerical investigations with an efficient approach for diagonally implicit Runge–Kutta schemes. The most important aspects are the extension towards fourth-order methods, the study of the convergence orders for the ELDIRK schemes, their respective stability regions and computational efficiency. A local error estimator is presented capturing the evolution of the phase-field problems, such that adaptive step size control for the new low-order embedded schemes based on an empirical approach for error estimation is achieved. A suitable parallel algorithm is presented with conclusive two-dimensional phase-field simulations based on a Kobayashi–Warren–Carter model including benchmarks for computational efficiency. The higher-order convergence suggested by the novel schemes is confirmed, and their effective results are demonstrated, resulting in a valuable semi-explicit addition to the family of Runge–Kutta time integration schemes.

Self-Adaptive PDLC Control Strategy With Smart Light Intensity Adjustment Using Photovoltaic-Thermoelectric Hybrid Energy Supply Technology

As an electronically controlled glass whose transmittance varies with voltage, polymer dispersed liquid crystal (PDLC) has the advantages of good privacy, high transmittance, fast response speed, and good heat shielding performance. It can be applied to electronically controlled building windows. However, PDLC needs to consume extra power, and electronically controlled PDLC has not yet been integrated with intelligent home environment control. In this paper, we proposed an adaptive variable voltage step-size fast control algorithm in Field Programmable Gate Array (FPGA)-based home environment control system, which is powered by photovoltaic-thermoelectric (PV-TEG) hybrid sources. The power supply voltage of the PDLC can be dynamically changed through the adaptive Pulse width modulation (PWM) duty cycle adjustment technology to realize intelligent transmittance control, thereby improving the indoor lighting environment. Experimental results show that by applying the proposed algorithm, the indoor temperature is reduced by 2°C, and the adaptive and stable regulation of the optimal indoor illumination of 500 Lux can be completed in less than 2.5 seconds. This paper provides a feasible solution for the low-energy intelligent integrated control and monitoring of the modern home environment.

Numerical investigations of new low‐order explicit last stage diagonal implicit Runge–Kutta schemes with the finite‐element method

AbstractInitial value problems can be solved efficiently by means of Runge–Kutta algorithms with adaptive step size control. Diagonally implicit Runge–Kutta (DIRK) methods are the most popular class among the diverse family of Runge–Kutta algorithms. In this paper, the novel class of low‐order explicit last‐stage diagonally implicit Runge–Kutta (ELDIRK) methods are explored, which combine implicit schemes with an additional explicit evaluation as an explicit last stage. ELDIRK Butcher tableaus are used to control embedded RK methods to obtain solutions of different orders. The lower‐order solution is obtained by classical implicit RK stages and the higher‐order solution is obtained by additional explicit evaluation. As a result, a significant reduction in computational cost is achieved by skipping the iterative solution of nonlinear systems for the additional step. The examination of the heat problem and the use of the innovative Butcher tableau in the finite‐element method are the main contributions of this work. Thus, it is possible to establish adaptive step size control for the new low‐order embedded methods based on an empirical method for error estimation. Two‐dimensional simulations are used to show an appropriate algorithm for the ELDIRK schemes. The new Runge–Kutta schemes' predictions of higher‐order convergence are confirmed, and their successful outcomes are illustrated.

Predicting wing-pylon-nacelle configuration flutter characteristics using adaptive continuation method

In this study, the aeroelastic response of a wing-pylon-nacelle system in subsonic and low supersonic flow regimes is analyzed using the continuation method in conjunction with an adaptive step size control algorithm. Idealizing the pylon and nacelle as a point mass, the computed effects of a standard structural analysis of the wing together with the pylon and nacelle are compared with those of a clean wing to build a reduced-order model for analysis. The aerodynamic forces relating to different reduced frequencies are assessed using the Doublet Lattice Method (DLM) in the subsonic flow regime and supersonic lifting surface theory relying on the unsteady linearized small-disturbance potential flow model in the low supersonic flow regime. The Rational Function Approximation (RFA) method is then utilized for the state-space formulation of the system equations, appended with the continuation method for flutter prediction. Thereafter, the linearized aeroelastic equations are resolved using the continuation method with adaptive step size, the results of which are matched with those obtained from the traditional p-k method to emphasize that the continuation method exhibits a distinct advantage in achieving better accuracy in estimating the flutter speed and identifying the “mode switching” phenomenon.

High fidelity transient solver in STREAM based on multigroup coarse-mesh finite difference method

This study incorporates a high-fidelity transient analysis solver based on multigroup CMFD in the MOC code STREAM. Transport modeling with heterogeneous geometries of the reactor core increases computational cost in terms of memory and time, whereas the multigroup CMFD reduces the computational cost. The reactor condition does not change at every time step, which is a vital point for the utilization of CMFD. CMFD correction factors are updated from the transport solution whenever the reactor core condition changes, and the simulation continues until the end. The transport solution is adjusted once CMFD achieves the solution. The flux-weighted method is used for rod decusping to update the partially inserted control rod cell material, which maintains the solution's stability. A smaller time-step size is needed to obtain an accurate solution, which increases the computational cost. The adaptive step-size control algorithm is robust for controlling the time step size. This algorithm is based on local errors and has the potential capability to accept or reject the solution. Several numerical problems are selected to analyze the performance and numerical accuracy of parallel computing, rod decusping, and adaptive time step control. Lastly, a typical pressurized LWR was chosen to study the rod-ejection accident.

Dynamic Modeling and Simulation of Integrated Electricity and Gas Systems

The dynamic simulation of the integrated electricity and gas system (IEGS) is effective for accurate operational guidance. However, the IEGS model governed by partial and ordinary differential-algebraic equations (PDAE) is challenging for analysis, because it poses heterogeneous properties, high dimensionality, and strong nonlinearity. To eliminate the nonlinear terms, this paper firstly derives the expression of average flow velocity (AFV) and reformulates the gas flow dynamics into the variable-coefficient partial differential equations (VC-PDE), whose consistency, stability, and convergence can be theoretically guaranteed. Then, a novel adaptive step size control (ASC) strategy is developed for VC-PDE solvers to improve computational efficiency and accuracy. Finally, this paper builds a framework for the simulation of IEGS considering the dynamics in the natural gas network (NGN) and coupling units, where the bidirectional interaction between the electric power system (EPS) and NGN is investigated. Case studies verify the superiority of the proposed VC-PDE solver with ASC regarding accuracy, stability, convergence, and efficiency in different systems.

An extended and improved particle‐spectral method for analysis of unsteady inviscid incompressible flows through a channel of finite length

AbstractThis paper addresses the development of a method to simulate and analyze unsteady flows of incompressible and inviscid fluids in rectangular channels. Mathematically, the flow problem formulates as Euler equations in terms of the stream function and vorticity. The normal velocity and the vorticity are given at the inlet of the channel, and only the normal velocity is specified at the outlet. The particle‐spectral method describes the vorticity field change in a Lagrangian sense by dynamics of marker particles and an Eulerian description of the velocity field using the global Galerkin approximation. The method involves the subdivision of flow domain on rectangular cells; the approximation of vorticity distribution in cells by least‐squares method; Galerkin coefficients calculation using analytical integration on space variables. For time integration, the five‐stage Runge‐Kutta method of the third accuracy order and the sixth one of symplecticity with adaptive step size control uses. The extension of the particle‐spectral approach comprises algorithms for fluid flows qualitative analysis jointly with the time simulation of unsteady flows. This includes the trajectories of marker particles calculation, computing the particle's residence time in the channel, kinetic energy dynamics, and algorithms for analyzing the dynamics and mixing of vortex patches. Test calculations based on reference solutions showed the effectiveness and good accuracy of the particle‐spectral approach. Using the developed method the dynamics of one and two initial vortex patches in the channel were investigated. We found the final state of their long‐time dynamics is time‐periodic or quasi‐periodic movements. The application of qualitative analysis algorithms demonstrates these final states produce a chaotic scattering of fluid particles in the flow‐through zone of the channel.

A C++ code for simulations of RHEED intensity oscillations within the kinematical approximation

This paper proposes a new software implementation for the rapid development and deployment of computer simulations of reflection high-energy electron diffraction (RHEED) intensity oscillations during the growth of thin epitaxial films. The calculations are based on the kinematical approximation. In this model the RHEED intensity oscillations from growing layers is considered as an effect of changes of surface roughness. New version program summaryProgram Title: growth22CPC Library link to program files:https://doi.org/10.17632/zkpwhrrfng.1Licensing provisions: GNU General Public License 3Programming language: C++Journal reference of previous version: Computer Physics Communications 170 (2005) 265–286Does the new version supersede the previous version?: YesReasons for the new version:The previous version of the program was written in Object Pascal [1]. This programming language is no longer popular in the scientific community. Responding to user's feedback, the current distribution of the software has been implemented in C++. This version implements improvements for ergonomics, computational performances and code readability. In addition, the output generated by the growth22 (saved in the coverageProfiles.dat file) can be used directly by program [2].Summary of revisions:The growth22 is written in standard C++, therefore one can use the GNU g++ compiler or any ANSI C++ compliant compiler to build the executable program. The code is distributed in the form of one growth22.cpp source file. A Table 1 lists the retired Object Pascal procedures and shows replacements by C++ functions. The present version of the program enables the user to easily manage the input data. Input data are loaded from a manually prepared inputData.dat text file where the user provides information on the parameters of the simulation, i.e. the growth model, the number of growing layers, the upper limit of growth time interval, the values of the phenomenological parameters that measure the net rate of transfer from one layer to the next, and the values of the growth rate for each layer. The user can edit this file as per the requirement of the problem domain. Details of the parameters in the sample inputData.dat file are described in the readme.pdf file provided with the software.The nature of problem:Generally, in order to explain the experimental data two theories—kinematical and dynamical—associated with high-energy electron diffraction in epitaxial layers are used. The dynamical theories take into account multiple scattering within the crystal and are used whenever diffraction from a perfect epitaxial layers is considered [3]. On the other hand, the kinematical theory treats the scattering from each plane in the crystal sample as being independent of others [1,4].In the kinematical theory of RHEED it is assumed that electron wave scatters from atom in the surface layer only once and then interferes with waves scattered on other atoms. This model can be used to describe diffraction effects for small angles of incidence, when the electron beam does not penetrate the deeper layers of the crystal. In the kinematical model the intensity oscillations of specularly reflected electron beam from growing layer is considered as an effect of changes of surface roughness.In this paper we present a time-efficient computer model for kinematical calculations of RHEED intensity oscillations during the growth of thin epitaxial films.The method of solution:According to kinematical approximation the diffracted-beam intensity is given by [4]:(1)I(S‖→=0,Sz)=|∑n=0∞(θn−θn+1)exp(−iSznd)|2 where S→ is the scattering vector, d is the interlayer spacing and n is an integer. θn−θn+1 is the exposed layer coverage which is treated as the atomic distribution function gives the position of surface atoms. For Sz=π/d Eq. (1) simplifies to:(2)I=|∑n=0∞(θn−θn+1)cos⁡(nπ)|2.Variation of the coverage ratio θn of the nth surface layer is modelled by a set of non-linear differential equations [1,4]. The Cash-Karp fifth-order Runge-Kutta method [5] with error estimation and adaptive stepsize control was used for solving initial value problem for non-linear differential equations [6]. The absolute and relative error tolerances and number of integration intervals are fully user-configurable.The proposed code approaches the modelling of RHEED from two different points of view, i.e. the data and functional perspectives. The data perspective (i.e. calculating the layer coverages as a function of their growth time) is used to characterize the data sources of the program whereas the functional perspective includes services provided by the program. It should be noted that the perspective of the data can be expressed by a completely different model introduced by the user.The coverage evolution curves for the diffusive growth model, corresponding kinematical diffracted intensity and total roughness of the surface are shown in Fig. 1. During the numerical calculations showed in Fig. 1 it was assumed that the number of growing layers (numLayers) equals 15, the upper limit of growth time interval (tMax) equals 12, the growth rate equals 1 MLs−1 (for all layers), and the values of the parameter kn (that determines the amount of interlayer diffusion) changing according to the linear relation:(3)kn=k1−(k1−kN)n−1N−1, where N is the total number of growing layers.

New low order Runge–Kutta schemes for asymptotically exact global error estimation of embedded methods without order reduction

Embedded Runge–Kutta (RK) methods provide an attractive methodology for solution of initial value problems (IVP) by means of adaptive step size control. According to Fehlbergs suggestion, only one extra function calculation is required to estimate the local error of the embedded method. In the present paper, this methodology is applied to three prominent low order implicit RK-schemes. To this end, new embedded Butcher-tableaus are presented for the implicit Euler-method, the trapezoidal-rule and the midpoint-rule by means of corresponding higher order RK-schemes. In order to account for the global error, an asymptotically exact global error estimator is derived for the so-called reversed embedded RK-method. Moreover, a cause of order reduction is illustrated by aid of a weak discontinuous formulation. It is prevented in a new hybrid integration algorithm by means of a Hermite interpolation for the source term of the IVP, and its implementation renders a convenient “node-to-node” stepping scheme for a laser beam simulation. Two numerical examples are presented with the following main results: The higher order convergence behaviors of the proposed new RK-schemes are verified, its successful performances for asymptotically exact global error estimation of the reversed embedded RK-method are shown, the proposed methodology is illustrated for two different global error functions, and the desired effect of the hybrid algorithm to prevent order reduction is illustrated.

New Step Size Control Algorithm for Semi-Implicit Composition ODE Solvers

Composition is a powerful and simple approach for obtaining numerical integration methods of high accuracy order while preserving the geometric properties of a basic integrator. Adaptive step size control allows one to significantly increase the performance of numerical integration methods. However, there is a lack of efficient step size control algorithms for composition solvers due to some known difficulties in constructing a low-cost embedded local error estimator. In this paper, we propose a novel local error estimator based on a difference between the semi-implicit CD method and semi-explicit midpoint methods within a common composition scheme. We evaluate the performance of adaptive composition schemes with the proposed local error estimator, comparing it with the other state-of-the-art approaches. We show that composition ODE solvers with the proposed step size control algorithm possess higher numerical efficiency than known methods, by using a comprehensive set of nonlinear test problems.

An adaptive homotopy tracking algorithm for solving nonlinear parametric systems with applications in nonlinear ODEs

Homotopy continuation has played an important role in solving nonlinear parametric systems. The homotopy tracking algorithm is the key part of this method. In this paper, we develop an adaptive stepsize control strategy for the homotopy tracking algorithm. We also extend this algorithm to the predictor/corrector methods for solving systems of nonlinear ODEs. The adaptive strategy for adjusting the stepsize is to guarantee the convergence of Newton’s method at each iteration along the solution path. Several numerical examples are used to show the efficiency of this new adaptive strategy by comparing it with other existing methods.

WKB-based scheme with adaptive step size control for the Schrödinger equation in the highly oscillatory regime

This paper is concerned with an efficient numerical method for solving the 1D stationary Schrödinger equation in the highly oscillatory regime. Being a hybrid, analytical–numerical approach it does not have to resolve each oscillation, in contrast to standard schemes for ODEs. We build upon the WKB-based (named after the physicists Wentzel, Kramers, Brillouin) marching method from Arnold et al. (2011) and extend it in two ways: By comparing the O(h) and O(h2) methods from Arnold et al. (2011) we design an adaptive step size controller for the WKB method. While this WKB method is very efficient in the highly oscillatory regime, it cannot be used close to turning points. Hence, we introduce for such regions an automated methods switching, choosing between the WKB method for the oscillatory region and a standard Runge–Kutta–Fehlberg 4(5) method in smooth regions.A similar approach was proposed recently in [Handley et al. (2016), Agocs et al. (2020)], however, only for an O(h)-method. Hence, we compare our new strategy to their method on two examples (Airy function on the spatial interval [0,108] with one turning point at x=0 and on a parabolic cylinder function having two turning points), and illustrate the advantages of the new approach w.r.t. accuracy and efficiency.

Mixed precision path tracking for polynomial homotopy continuation

This article develops a new predictor-corrector algorithm for numerical path tracking in the context of polynomial homotopy continuation. In the corrector step, it uses a newly developed Newton corrector algorithm which rejects an initial guess if it is not an approximate zero. The algorithm also uses an adaptive step size control that builds on a local understanding of the region of convergence of Newton’s method and the distance to the closest singularity following Telen, Van Barel, and Verschelde. To handle numerically challenging situations, the algorithm uses mixed precision arithmetic. The efficiency and robustness are demonstrated in several numerical examples.

Availability analysis of Tube Light Manufacturing industry using numerical technique

Continuous improvement in the performance of any industrial system is highly desirable. But failure is an inevitable phenomenon. In order to run industrial systems without failure reliability and availability of each subsystem needs to be maintained. So the maintenance of the working units of the system plays a very important role to make system available for a pre define period of time. The system performance can be increased to the maximum extent by adopting the suitable maintenance strategy. In this paper a computational technique is used for the analysis of Time Dependant System Availability (TDSA) of Tube Light manufacturing system has been addressed.. This system comprises of seven subsystems namely End Sealing cum Coating Machine, Baker Machine, Wiper Machine, Sealing Machine and Exhaust Machine, Capping Machine and Activator machine. Mathematical formulation of the system is done using Markov Approach (MA). The model has been solved TDSA with the help of Adaptive step size control Runge-Kutta Method. After carrying out the TDSA analysis it is found that subsystem- II i.e Baker machine is critical subsystem having more effect on TDSA of the system on basis of repair rate on as compared to other subsystems. Further, present research is also helpful in deciding maintenance repair priorities of different subsystems.

A shooting approach to the scaled boundary finite element equations of elastodynamics in the frequency domain

Dealing with numerical analysis of problems, especially ones with semi-infinite boundaries, scaled boundary finite element method has emerged as one of the efficient tools for the task. Combining the exactness of strong forms with the flexibility of weak formulations makes the method an improvement to its predecessors. Problem with the method arises when the analytical solution of the semi-discretized system is not available, which is the case for numerous problems. In the most recent attempt to solve the issue, a shooting method was proposed for elastostatic problems. Generality of the method removes any concerns regarding the type of governing equations since it no longer needs any analytical solutions. In this study, the method is further developed to deal with elastodynamic problems of viscoelastic media in the frequency domain, for which the closed-form solution is not available. Adopting the shooting method, the two-point boundary value problem of each subdomain is treated as a second-order initial value problem. Boundary conditions of the bounded subdomains are formulated by enforcing equilibrium at both ends. On the other hand, in contrast to linear elasticity, neither forces nor displacements vanish at infinity. As a result, a special treatment is required to deal with the boundary conditions of the unbounded subdomains. A Runge–Kutta–Nyström method, which consists of an embedded pair of adjacent order, is utilized to integrate the strong form of oscillatory equations enabling local error estimation and adaptive step-size control. This guarantees efficiency since the integration proceeds not by equal step-sizes, but by equal user-defined local errors. In the first step, analytical solution of a viscoelastic half-space subjected to a harmonic strip load is used for the verification process. First, the half-space is discretized by Voronoi cells and each truncated edge is defined using an unbounded subdomain. Next, one bounded subdomain is used to define the whole bounded area, and one unbounded subdomain resembles the radiation condition of the far-field boundaries. For the cases involving only bounded subdomains, three other numerical tests, including a perforated plate, an L-shaped panel, and a bi-material specimen are employed, and the effect of discretization method on the results is analyzed. The method is a huge step forward since the size of subdomains are not restricted, far-field boundaries can be arbitrary-shaped, hysteretic damping can be naturally included, no additional unknowns are introduced, no compound matrices are defined, no preconditioning matrices are needed, and a simple unified approach is used for both bounded and unbounded subdomains.