Mattheij, R.M.M.

Savcenco, V.

English

Pure OAI Repository

 Multirate numerical integration for stiff ODEsCitation for published version (APA):Savcenco, V., & Mattheij, R. M. M. (2008). Multirate numerical integration for stiff ODEs. (CASA-report; Vol.0839). Technische Universiteit Eindhoven.Document status and date:Published: 01/01/2008Document Version:Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)Please check the document version of this publication:• A submitted manuscript is the version of the article upon submission and before peer-review. There can beimportant differences between the submitted version and the official published version of record. Peopleinterested in the research are advised to contact the author for the final version of the publication, or visit theDOI to the publisher's website.• The final author version and the galley proof are versions of the publication after peer review.• The final published version features the final layout of the paper including the volume, issue and pagenumbers.Link to publicationGeneral rightsCopyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright ownersand it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.            • Users may download and print one copy of any publication from the public portal for the purpose of private study or research.            • You may not further distribute the material or use it for any profit-making activity or commercial gain            • You may freely distribute the URL identifying the publication in the public portal.If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, pleasefollow below link for the End User Agreement:www.tue.nl/taverneTake down policyIf you believe that this document breaches copyright please contact us at:openaccess@tue.nlproviding details and we will investigate your claim.Download date: 04. Oct. 2023Multirate Numerical Integration for Stiff ODEsV. Savcenco and R.M.M. MattheijTechnische Universiteit Eindhoven, P.O. Box 513, 5600 MB Eindhoven, TheNetherlandsSummary. This paper contains an overview of a self-adjusting multirate method.A simple extension which allows the improvement of the efficiency of the method isintroduced. The performance of the extended and the original method is comparedfor a test problem.1 IntroductionFor the numerical solution of systems of ODEs there are many methods avail-able. These methods use time steps that are varying in time, but are constantover the components. However, there are many problems of practical inter-est, where the temporal variations have different time scales for different setsof the components. For example, cellular phones consist of coupled digitaland analogue sub-circuits, which operate in nano- and micro-seconds, respec-tively. The motion of the particles around a star, which attracts mass froma secondary star, in astrophysics is described by a large system of ordinarydifferential equations. In this system the components, that correspond to theparticles near the center, are much faster than those corresponding to thedistant ones. To exploit these local time scale variations, one needs multi-rate methods that use different, local time steps over the components. In thesemethods big time steps are used for the slow components and small time stepsare used for the fast ones.Various multirate methods were developed for solving systems with differ-ent time scales. The first descriptions of multirate schemes were given by Gearand Wells [4] for multistep methods. Sand and Skelboe [9] studied the sta-bility of backward Euler multirate methods. Multirate methods for non-stiffproblems have been examined by Engstler and Lubich [3]. A multirate schemebased on the partitioned Runge-Kutta methods was introduced by Günther,Kværnø and Rentrop [5]. In [1, 13, 14] multirate methods have been appliedto the modeling of electrical networks. Multirate methods for hyperbolic con-servation laws were studied by Constantinescu and Sandu [2].2 V. Savcenco and R.M.M. MattheijA multirate method based on the Rosenbrock methods, together with aself-adjusting partitioning strategy was introduced and analyzed in [7, 10, 11,12]. In this paper we present an overview of this method and suggest a way toimprove it. The comparison of the numerical results obtained with the originaland extended strategy is presented.The paper is organized as follows. In Section 2 we will introduce the Rosen-brock methods which will be used as our basic numerical integration methodsand describe the multirate time stepping strategy. The performance of theextended version of the multirate strategy for a test problem is discussed inSection 3. Finally, Section 4 contains the conclusions.2 A Multirate time stepping strategyWe will consider multirate methods for solving systems of ODEsw′(t) = F (t, w(t)), w(0) = w0, (1)with given initial solution w0 ∈ Rm. The approximations at the global timelevels tn will be denoted by wn.Our multirate time stepping strategy is based on local temporal errorestimation. For a given time step ∆tn = tn−tn−1, we compute a first, tentativeapproximation at the new time level for all components. For those componentsfor which the error estimator indicates that the local temporal error is largerthan a given tolerance Tol, the computation is redone with smaller steps. Therefinement is recursively continued until the error estimator is below Tol for allcomponents. Schematically, with components horizontally and time vertically,the multirate time stepping is displayed in Figure 1.tn−1tnFig. 1. Multirate time stepping for a time slab [tn−1, tn].In the original strategy [12], the refinement is performed by recalculatingthe required components with halved steps. For many problems, the timesteps needed for the active components are much smaller than those neededfor the slow ones. In such cases it is more efficient to immediately recomputethe active components with more than two smaller steps instead of doingseveral halving recursive refinements. Therefore in this paper we will extendthe original strategy and will assume that the number of smaller time stepsat the refinement stage can be also larger than two. Using ideas from [12], itis possible to design an adaptive procedure of choosing the size of the timeslabs for the extended strategy.Multirate Numerical Integration for Stiff ODEs 32.1 Main time integration methodsIn this paper we will use the Rosenbrock methods [6] as our basic numericalintegration methods. To proceed from tn−1 to tn = tn−1 + τ , an s-stageRosenbrock method calculateswn = wn−1 +s∑i=1biki , (2)ki = τFtn−1 + αiτ, wn−1 +i−1∑j=1αijkj + τ∂F∂w(tn−1, wn−1)i∑j=1γijkj+γiτ2 ∂F∂t(tn−1, wn−1), i = 1, . . . , s , (3)where αij , γij , αi, γi, bi are real parameters defining the method and τ denotesthe step size. For the local error estimation withing the variable step sizecontrol we use the embedded formulawn = wn−1 +s∑i=1biki , (4)which uses the same ki-values as (2), but has different weights.2.2 Interface treatmentDuring the refinement stage, values at the intermediate time levels of compo-nents which are not refined might be needed. These values can be obtainedby use of dense output built in the time integration methodwI(tn−1 + θτ) = wn−1 +s∑i=1bi(θ)ki, 0 ≤ θ ≤ 1 . (5)Proper interface treatment is very important for multirate schemes. Use ofdense output of order lower than the order of the main time integrationmethod can lead to order reduction.It is well known that use of Rosenbrock methods for problems with stiffsource terms can lead to order reduction. During the refinement step, subproblems with stiff source terms have to be solved. An easy applicable tech-nique, to avoid the order reduction, was proposed in [11]. Assuming that g(t)is a component of F (t, w(t)), this technique suggests that the source termsg(tn−1 + αiτ) + γiτg′(tn−1) in a Rosenbrock method (3) of order p shall bereplaced by gn−1,i with gn−1 = [gn−1,i] chosen asgn−1 =p∑k=0Bkeτkg(k)(tn−1) , (6)where B = [αij + γij ] ∈ Rs×s and e = [1] ∈ Rs.4 V. Savcenco and R.M.M. Mattheij3 Numerical testIn this section we will present numerical results for a test problem. We con-sider the behavior of both strategies: original [12] and extended (where refine-ment with more than two steps is possible). The results are compared to thesingle-rate approach, also using the same basic time integration method. As ameasure for the amount of work we consider the total number of componentsat which solutions are computed over the complete time integration interval,multiplied by the number of stages of the method. The fact that with themultirate approach some solution components are computed several times atcertain time levels is taken into account. As the basic time integration method,for solving this problem, we use the two-stage second-order Rosenbrock ROS2method [8].3.1 An Inverter Chain ProblemAn inverter is an electrical sub-circuit which transforms a logical input signalto its negation. The inverter chain is a concatenation of several inverters,where the output of an inverter serves as input for the succeeding one. Aninverter chain with an even number of inverters will delay a given input signaland will also provide some smoothing of the signal.The model for m inverters consists of the equations{w′1(t) = Uop − w1(t) − Υg(uin(t), w1(t)),w′j(t) = Uop − wj(t) − Υg(wj−1(t), wj(t)), j = 2, . . . ,m ,(7)whereg(u, v) =(max(u − Uthres, 0))2−(max(u − v − Uthres, 0))2. (8)The coefficient Υ serves as stiffness parameter. We solve the problem for achain of m = 500 inverters with Υ = 100, Uthres = 1 and Uop = 5. The initialcondition iswj(0) = 6.247 · 10−3 for j even, wj(0) = 5 for j odd. (9)The input signal is given byuin(t) =t − 5 for 5 ≤ t ≤ 10 ,5 for 10 ≤ t ≤ 15 ,52 (17 − t) for 15 ≤ t ≤ 17 ,0 otherwise.(10)An illustration of the solution for some of the even components is given inFigure 2.Multirate Numerical Integration for Stiff ODEs 50 30 60 90 12002.55w2 w126 w250 w374 w498Fig. 2. Solution components wj(t), j = 2, 126, 250, 374, 498, for problem (7)–(10).In Table 1 the errors at output time T = 130 (measured in the maximumnorm with respect to an accurate reference solution) together with the amountof work are presented for several tolerances for the single-rate method andthe extended multirate strategy. The speedup for both original and extendedstrategies is calculated.It is seen from the table that a substantial improvement in amount of workis obtained for this problem. For the single-rate scheme, the amount of work isalmost 18 times larger than for the extended multirate scheme. Moreover, theerror behavior of the multirate scheme is very good. We can also see that forthis problem we get a considerably larger speedup for the extended strategycompared to the original strategy.Table 1. Errors and work amount for problem (7)–(9).Single-rate Multirate (extended) Speeduptol error work error work original extended5 · 10−4 1.44 · 10−1 45649872 1.47 · 10−2 2846068 8.7 16.010−4 3.94 · 10−2 94524592 7.16 · 10−3 5512400 13.0 17.15 · 10−5 1.37 · 10−2 131413560 3.24 · 10−3 6980676 13.5 18.810−5 2.04 · 10−3 287207252 9.22 · 10−4 14332486 11.1 20.04 ConclusionsIn this paper we made on overview and extended the multirate time steppingstrategy introduced in [7, 10, 11, 12].As seen from the numerical tests, the efficiency of time integration methodscan be significantly improved by using large time steps for inactive compo-nents, without sacrificing accuracy. Comparing the results obtained for theoriginal and the extended strategies, we do have preference for the extended6 V. Savcenco and R.M.M. Mattheijapproach, where the values of the active components can be recalculated bythe use of more than two smaller time steps.References1. A. Bartel, M. Günther, A multirate W-method for electrical networks in statespace formulation. Journal Comp. Appl. Math., 147, 411–425 (2002)2. E.M. Constantinescu, A. Sandu, Multirate time stepping methods for hyper-bolic conservation laws. Journal on Scientific Computing, 33, 239–278 (2007)3. C. Engstler, C. Lubich, Multirate extrapolation methods for differential equa-tions with different time scales. Computing, 58, 173–185 (1997)4. C. Gear, D. Wells, Multirate linear multistep methods. BIT, 24, 484–502 (1984)5. M. Günther, A. Kværnø, P. Rentrop, Multirate partitioned Runge-Kutta meth-ods. BIT, 41, 504–514 (2001)6. E. Hairer, G. Wanner, Solving Ordinary Differential Equations II – Stiff andDifferential-Algebraic Problems. Springer, Berlin (1996)7. W. Hundsdorfer, V. Savcenco, Analysis of a multirate theta-method for stiffODEs. To appear in Journal Appl. Num. Math.8. W. Hundsdorfer, J.G. Verwer, Numerical Solution of Time-DependentAdvection-Diffusion-Reaction Equations. Springer, Berlin (2003)9. J. Sand, S. Skelboe, Stability of backward Euler multirate methods and con-vergence of waveform relaxation. BIT, 32, 350–366 (1992)10. V. Savcenco, Comparison of the asymptotic stability properties for two multi-rate strategies. Journal Comp. Appl. Math., 220, 508–524 (2008)11. V. Savcenco, Construction of a multirate RODAS method for stiff ODEs. Toappear in Journal Comp. Appl. Math.12. V. Savcenco, W. Hundsdorfer, J.G. Verwer, A multirate time stepping strategyfor stiff ordinary differential equations. BIT, 47, 137–155 (2007)13. M. Striebel, M. Günther, A charge oriented mixed multirate method for aspecial class of index-1 network equations in chip design. Journal Appl. Numer.Math, 53, 489–507 (2005)14. A. Verhoeven, A.El. Guennouni, E.J.W. ter Maten, R.M.M. Mattheij, Ageneral compound multirate method for circuit simulation problems. InA.M. Anile, G. Ali, G. Mascali, editors, Scientific Computing in ElectricalEngineering. Springer, Berlin Heidelberg (2006)

Multirate numerical integration for stiff ODEs

Abstract

Similar works

Full text

Available Versions

Pure OAI Repository

NARCIS

Pure OAI Repository