A self adjusting multirate algorithm based on the TR-BDF2 method
Luca Bonaventura Francesco Casella Ludovica Delpopolo Akshay Ranade
math.NA





A self adjusting multirate algorithm based on the TR-BDF2  method
    Luca Bonaventura^(1), Francesco Casella^(2) 

Ludovica Delpopolo^(1), Akshay Ranade^(3) 
    December 30, 2023
============================================================================================





^(1) MOX – Modelling and Scientific Computing 

Dipartimento di Matematica, Politecnico di Milano 

Via Bonardi 9, 20133 Milano, Italy

 luca.bonaventura@polimi.it, ludovica.delpopolo@polimi.it






^(2)  Dipartimento Elettronica, Informazione e Bioingegneria

 Politecnico di Milano 

Via Ponzio 34/35, 20133 Milano, Italy

francesco.casella@polimi.it






^(3)  Department of Computer Science

University College Cork

Western Road

Cork, Ireland

akshay.ranade@cs.ucc.ie






Keywords:  Multirate methods, stiff ODE problems, Implicit Runge Kutta methods, SDIRK methods




AMS Subject Classification:  65L04, 65L05, 65L07, 65M12, 65M20







We propose a self adjusting multirate method based on the TR-BDF2 solver. The potential advantages
of using  TR-BDF2  as the key component of a 
 multirate framework are highlighted. A linear stability analysis of the resulting approach is 
presented and the stability features of the resulting  algorithm are analysed. 
The analysis framework is completely general and allows to study along the same lines the stability of self adjusting multirate methods based on a generic one step solver. A number of numerical experiments demonstrate the efficiency and accuracy of the resulting approach also the time discretization of hyperbolic partial differential equations.









§ INTRODUCTION

  

 
 The concept of  multirate methods was first proposed in <cit.> and a great attention
 has been devoted to these methods, see e.g. <cit.>, <cit.> and the more comprehensive review  reported in <cit.>.  The basic idea of multirate methods is to integrate each component using a time step that is the most appropriate to the timescale of that component. Slowly varying components are integrated with larger time steps, while smaller time steps are used for fast components only. Thus, multirate methods can avoid a significant amount of the computation that is necessary in the standard, single rate approach. In other words, in the multirate approach,  the most appropriate time resolution is employed for each  state variable of the system. In the context of multirate algorithms, the faster components  are often referred to as  active and the slow ones are called instead  latent. 
  It is also to be remarked that conceptually similar approaches
 have been introduced   in the  literature on numerical weather prediction, see e.g. <cit.>,
  for the purpose of avoiding the time step stability restrictions related  to fast propagating acoustic waves.
 
 
 The specific way in which the partitioning is performed and the degree of coupling
 allowed between latent and active components  determines the efficiency and the stability of the resulting methods. The stability of  multirate methods has been investigated, among others,  in <cit.> and <cit.>, while the stability  of split-explicit methods has been analyzed in   <cit.>, <cit.>.  
  In general, it is seen that the stability region of the multirate methods decreases as the strength of the coupling between the partitions  or the stiffness of the system increases. 
   In spite of this, there is ample evidence in the literature that multirate algorithms can offer a
    promising path for the numerical integration of large systems with widely varying time scales 
    or with highly localised periods of activity.
 In order to address the potential problems of multirate methods,
 the so-called compound step multirate methods have been proposed in  <cit.>. 
A self-adjusting, recursive time stepping strategy has been   proposed in  <cit.>.
In this kind of approach, a tentative  global  step is first taken for all components, using a robust, unconditionally stable method. The time step is then reduced
only  for those components for which some local error estimator is greater than the specified tolerance.
In this way, automatic detection of the latent components is achieved.
 In <cit.>, an algorithm based on this recursive strategy was shown to be more stable than the original compound step algorithm, while a stability analysis of a self-adjusting  algorithm based on the θ-method  has been proposed in  <cit.>.
 
 In this paper, we propose a self-adjusting multirate method based on the same concept as
 <cit.>, that relies on the TR-BDF2 method as fundamental single rate solver.  The TR-BDF2 method has been originally introduced in <cit.> and more thoroughly analyzed in  <cit.>. 
 It is a second order, one step,
 L-stable implicit method endowed with a number  of interesting properties, as discussed in  <cit.>.
  For example, it allows for cheap error estimation
 via an embedded third order method, it can provide  continuous output via a cubic Hermite interpolant  and 
 it has an explicit second order
 companion with which it can be combined to form a second order, implicit-explicit additive Runge Kutta method.
 Due to its favourable properties, it has been recently applied for efficient discretization
 of high order finite element methods for numerical weather forecasting in <cit.>, <cit.>.
 Its appealing features  also make it   an interesting candidate for multirate approaches targeted at  large scale stiff problems, resulting either from  the simulation industrial systems or from the spatial discretization of partial differential equations.
 In particular, the  cubic Hermite interpolant can be employed to maintain higher accuracy also for stiff
 problems in the framework of a multirate approach.
  The proposed self-adjusting multirate TR-BDF2 method   is  equipped here with a partitioning 
  and time step selection criterion based on the technique proposed in  <cit.>.  
  
  The stability of the method has been analyzed in the framework of the classical linear model problem.
  Even though no complete theoretical results can be achieved, numerical computation of the stability function norm
  for a range of relevant model problems shows that the method has  good stability properties.
  This is true   not only for contractive problems with strictly dissipative behaviour, but also for problems
  with purely  imaginary eigenvalues, which suggests that the method could be advantageous also for
  the application to hyperbolic PDEs and structural mechanics problems.
  The numerical results obtained on several benchmarks  show that the application of the proposed method leads
  to significant efficiency gains, that are analogous to those achieved by other self-adjusting multirate approaches
  with respect to their single rate counterparts. In particular, the method appears to perform well when applied to the time discretization of nonlinear, hyperbolic partial differential equations, allowing to achieve automatic detection
  of complex localized phenomena such as shock waves and significant reductions in computational cost.
 

 In section <ref>, the self-adjusting implicit multirate approach is described.
    The time step refinement criterion and the resulting
   partitioning principle between active and latent components are presented  in section <ref>. 
 The TR-BDF2 method is reviewed in section <ref>.
 The interpolation procedures we have considered are described in section
 <ref>. A linear stability analysis is then presented in section <ref>.
 The approach proposed for such analysis is rather general and allows to study generic multirate methods based on
 any one step solver and interpolation procedure. 
 Numerical simulations are presented in section <ref>, while some conclusions
 and perspectives for future work are presented in <ref>.
 

  
  

§ A SELF ADJUSTING IMPLICIT MULTIRATE APPROACH

  
 We focus on multirate methods for the solution of the Cauchy initial value problem

    𝐲^'=𝐟(t,𝐲)       𝐲(0)=𝐲_0∈𝐑^m,      t∈[0,T].

We   consider time discretizations associated to discrete time levels
t_n,  n=0,…,N such that h_n=t_n+1-t_n and
we will denote by
𝐮_n the numerical approximation of 𝐲(t^n). 
 We will also denote by u^n+1= S ( u^n,h_n) the
implicitly defined operator S :𝐑^m→𝐑^m
 whose application is equivalent to the computation of  one step of size h_n  of
a given single step method. While here only implicit methods will be considered,
the whole framework can be extended to explicit and IMEX methods.
Notice that, if 𝐏 is the projector onto a linear subspace V⊂𝐑^m with dimension p < m,
the operator S^ V :𝐑^p×𝐑^m-p→𝐑^p
that represents the solution of the subsystem obtained freezing the components
of the unknown belonging to  V^⊥  to the value z∈𝐑^m-p can be defined by
y =    S^ V ( x, z, h_n )= 𝐏 S (  x⊕𝐳, h_n) .
Furthermore, we will denote by 𝐐 ( u^n+1, u^n,ζ) 
an interpolation operator, that provides an approximation of the numerical solution at
intermediate time levels t_n+ζ, where ζ∈[0,h_n]. Linear interpolation is often employed, but, for
 multistage methods, knowledge of the intermediate stages  also allows the application of more accurate interpolation procedures  without 
substantially increasing the computational cost.

We propose a self-adjusting  multirate algorithm that
  is a generalization of the method introduced in <cit.>, 
in which the choice of the local time step is left completely to the time step selection criterion.  
The multirate algorithm can be described as follows.

	
  1)  Perform  a tentative global (or macro) time step of size h_n with the standard single rate method and compute
	 û^n+1=  S ( u^n,h_n). The specific way in which error control mechanisms are
	  applied to the choice of the global and local time steps
	are crucial for the efficiency of the algorithm and will be discussed in detail in section
	<ref>.  
	 
  2)  Compute the error estimator 
	 and  partition the state space into  active and  latent variables, based on the value of the error estimate.
	  The projection onto the subspace V_0
	of the active variables is denoted by 𝐏_n=𝐏_n^(0),  
	while the projection onto the complementary
	subspace will be denoted by 𝐏̅_n^(0).
	Define 𝐏̅_n^(0) u^n+1 =    𝐏̅_n^(0)û^n+1
	as well as u^n,0=  u^n and
	  t_n,0=t_n. 
	
  3)  For k≥ 1,
	choose a local (or micro) time step  h_n^(k-1) for the active variables, based on the value of the error estimator.  Set  t_n,k=min{t_n,k-1+h_n^(k-1),  t_n+1}. 
	
		
  3.1)      Update the latent variables by interpolation
	
    𝐏̅_n^(k-1) u^n,k=
    	𝐐  ( 𝐏̅_n^(k-1) u^n+1,  𝐏̅_n^(k-1) u^n,k-1, h_n^(k-1) ).

	
  3.2)  Update the active variables  by computing
		
    𝐏^(k-1)_n u^n,k=
    	   S^ V_k-1 (𝐏^(k-1)_n u^n,k-1,  𝐏̅_n^(k-1) u^n,k, h_n^(k-1)).

	
  3.3) Compute the error estimator  for the active variables only and   partition again  V_k-1
	into latent and active variables. Denote by V_k⊂ V_k-1 the new
	subspace of active variables and  by  𝐏_n^(k)  the corresponding projection.

	
	
  3.4) Repeat  3.1) - 3.3) until t_n,k=t_n+1.





















		
	

 
Clearly, the effectiveness of the above procedure depends in a crucial way on the accuracy and stability of the basic ODE 
solver S, as well as on the time step
refinement and partitioning criterion.  Furthermore, as well known in multirate methods,
unconditional stability of the solver S does not necessarily entail that the same property holds for
the derived multirate solver.
The refinement and partitioning criterion will be described in detail in section
<ref>. Here, we only remark that, as it will be shown by the numerical experiments
in section <ref>, the approach outlined above is able to reduce significantly the computational
cost with respect to the equivalent single rate methods, without a major reduction in stability
for most of the envisaged applications.
 

 
   

§ THE TIME STEP REFINEMENT AND PARTITIONING CRITERION

  
 We now describe in detail the time step refinement and partitioning strategy that we have used  in the multirate algorithm described in section <ref>.
Our approach is based on the strategy proposed for an explicit Runge Kutta multirate method in
 <cit.>, where the time steps for refinement are obtained from the error estimates of the global step. The user specified tolerance plays an important role in the partitioning of the system. In <cit.>, a simple absolute error tolerance was considered. However, in most engineering system simulations, whenever the typical values of different components can vary greatly, the tolerance used is in general
 a combination  of absolute and relative error tolerances, see e.g. <cit.>.
  We have thus extended the strategy  proposed in
 <cit.> in order to employ such a combination of absolute and relative error tolerances.
 
More specifically, denote by   τ_r, τ_a the user defined error tolerances for relative and absolute errors, respectively.
Furthermore, assume that the tentative global step û^n+1= S ( u^n,h_n^(0))
has been computed and that   an   error estimator ϵ^n+1,0 is available.
The first task of the time step selection criterion is to asses whether the global time step h_n^(0) has been properly chosen.
Denoting by ϵ^n+1,0_i, û^n+1_i, i=1, …,m the i-th components of the error estimator and of the
tentative global step numerical solution, respectively, 
we define a   vector η^n+1,0 of normalized errors with components
   
    η^n+1,0_i= ϵ^n+1,0_i/τ_r|   û^n+1_i| + τ_a.

   Since clearly the condition max_i  η^n+1,0_i ≤ 1
has to be enforced, before proceeding to the partitioning into
active and latent variables, this condition is checked and the global step
is repeated with a smaller value of h_n^(0) whenever it is violated.
Numerical experiments have shown that, while an increase of efficiency with respect
to the single rate version of the method is always achieved, independently of the choice of the tentative
step,  the greatest improvements are only possible if the global time step  does not have
to  be repeated too often.  


   Once the condition  max_i  η^n+1,0_i ≤ 1 is satisfied,
   the set of indices of the components flagged for the first level of time step refinement is identified by 


    A^n+1,0 = { i :  η^n+1,0_i > δη^n+1,0_∞},

where δ∈(0,1) is a user defined coefficient. The smaller the value of δ,
the larger is the fraction of  components marked for refinement. Notice that,
if δ is set to unity,  the algorithm effectively operates in  single rate mode.
For each iteration k=1,…,k_max of the algorithm described in section <ref>,
active variable sets A^n+1,k are then defined analogously.
In each iteration k=0,…,k_max, the time step h_n^(k) is chosen according
to the standard criterion (see e.g. <cit.>, <cit.>) 


    h_n^(k) = νmin_j∈ A^n+1,k( τ_r|û^n+1|_j +τ_aϵ^n+1,k_i )^1/p+1,

where p is the convergence order of the solver S and ν is an user defined 
safety parameter.
 
 

§ THE SINGLE RATE TR-BDF2 METHOD

  
 The TR-BDF2 method was originally proposed in <cit.>. It is  a composite one step, two stage method, consisting
of one stage of the trapezoidal method followed by another of the BDF2 method. The stages are so adjusted that both the trapezoidal and the BDF-2 stages use the same Jacobian matrix. This composite method  has been reinterpreted 
in <cit.> 
as a one step Diagonally Implicit Runge Kutta (DIRK) method with two internal stages  The TR-BDF2 method is also in some sense the optimal method among all 3-stage DIRK methods,  owing to the following properties:

	
  1) it is first same as last (FSAL), so that only two implicit stages need to be evaluated; 
	
  2) the Jacobian matrix for both the stages is the same;
	
  3) it has an embedded third order companion that allows for a cheap error estimator;
	
  4) the method is strongly S-Stable;
	
  5) all the stages are evaluated within the time step;
	
  6) it is endowed with a cubic Hermite interpolation algorithm for dense output that yields globally C^1 continuous
	trajectories.

Features 3), 5) and 6) will play an important role in the multirate method proposed here.
Furthermore, as shown in  <cit.>, the method   has an explicit second order
 companion with which it can be combined to form a second order implicit-explicit additive Runge Kutta method.
 Due to its favourable properties, it has been recently applied for efficient discretization
 of high order finite element methods for numerical weather forecasting in <cit.>, <cit.>.

The TR-BDF2 method, considered as a composite method consisting of a step with the trapezoidal method followed by a step of the BDF2 method, can be written as


    u^n+γ   = u^n + γ h_n/2 ( f(t, u^n) +  f(t, u^n+γ) )
     u^n+1   = 1/γ(2-γ) u^n+γ -(1-γ)^2/γ(2-γ) u^n + 1-γ/2-γh_n
    	 f(t, u^n+γ)



The TR-BDF2 method viewed as a DIRK method has the following Butcher tableau:








0      0      0      0

γ     d      d      0

1     w      w      d
      w      w      d
 *    (1-w)/3     (3w+1)/3    d /3
	
	



where γ=2-√(2), d=γ/2 and w=√(2)/4 and the row * corresponds to the embedded
 third order method that can be used to build a convenient  error estimator.
 Here, the value  γ=2-√(2) is chosen for the two stages to have the same Jacobian matrix, 
 as well as in order to achieve L-stability, as it will be shown shortly.
It can be seen that the method has the First Same As Last (FSAL) property, i.e., the first stage of any step is the same as the last stage of the previous step. Thus, in any step, the first explicit stage need not be computed. 

 The implementation of the two implicit stages is done as suggested in <cit.>.
  The two stages according to the Butcher tableau are given by equation
(<ref>).
Instead of iterating on the variable u^n+1, 
we define another variable z  =h f(t, u) and solve for this variable by iteration. Thus, for the first implicit stage, we take  and z^n+γ,k is computed by Newton iterations as


    (𝐈-dh𝐉^n)Δ^k   =h_n f(t_n+γ,_n+γ^k) -  z^n+γ,k
     z^n+γ,k+1   =  z^n+γ + Δ^k,

where  J^n=  J(t_n,  u_n) denotes the Jacobian of f.  
Similarly, for the second implicit stage we take u_n+1^k =  u_n + w z^n + w  z^n+γ + d z^n+1,k and the Newton iteration is given by

    (𝐈-dh𝐉^n)Δ^k   =h f(t_n+γ,_n+γ^k) -  z^n+γ,k
     z^n+1,k+1   =  z^n+1 + Δ^k.


The reason for doing so is that if the problem is stiff, a function evaluation will amplify the numerical error in the stiff components. Iterating on z, however, ensures that it is computed to within the specified tolerance. For a more detailed discussion the interested reader is referred to <cit.>.
An important point to be noted is that, although the TR-BDF2 method is L-stable, its third order companion formula is not. Therefore, the error estimate at time level n+1,   given by

    ϵ^n+1,*=(b_1^* - b_1) z^n +(b_2^* - b_2) z_n+γ +(b_3^* - b_3) z_n+1

cannot be expected to be accurate for stiff problems. To overcome this problem, it is suggested
in  <cit.>   
to modify the error estimate by considering as error estimator the quantity ϵ_n+1 defined as the solution
of the linear system
 

    (𝐈-dh𝐉^n)ϵ^n+1=ϵ^n+1,*.

This modification of the error estimate  allows to improve it for stiff components, while preserving its accuracy in the limit of small time steps. Notice that the stability function of the TR-BDF2 method is given by

    R(z) =[1 + (1 - γ)^2] z + 2(2-γ)/z^2(1-γ)γ + z(γ^2-2) +2(2-γ),

whence it can be seen that the method is L-stable for γ = 2-√(2).
 
























   

§ THE INTERPOLATION PROCEDURES

  
 An essential component of any multirate algorithm is the  procedure employed
 to reconstruct the values of the latent variables at those intermediate time levels
 for which only the active variables are computed. Self-adjusting multirate procedures
 based on implicit methods, such as the one proposed in <cit.> and that
 presented in this paper, can avoid the use of extrapolation, thus increasing their overall stability.
 The simplest and most commonly used procedure is linear interpolation, that is defined for ζ∈[0,h_n] as
 
    𝐐 ( u^n+1, u^n,ζ) = ζ/h_n u^n+1 + (h_n-ζ)/h_n u^n.

 In case a multi-stage solver is employed, higher order interpolation can also be employed at no extra cost.
 Assuming that a sufficiently accurate approximation of the solution  u^n+λ is available at time t_n+λ∈ (t_n,t_n+1)
 and setting h_λ= t_n+λ-t_n,
 quadratic Lagrange interpolation can then be written as 
 
    𝐐 ( u^n+1, u^n,ζ)    =   (ζ-h_λ)ζ/h_n(h_n-h_λ) u^n+1  
    + ζ (ζ -h_n) /(h_λ-h_n)h_λ u^n+λ
       +   (h_λ-ζ)(h_n-ζ) /h_λh_n u^n.

 
 An interesting feature of the TR-BDF2 method in the multirate framework is that the method, as explained in <cit.>,
is endowed with a cubic Hermite, globally C^1 interpolation algorithm for dense output. 
Using the notation of section   <ref>, the Hermite cubic interpolant can be defined as

    𝐐 ( u^n+1, u^n,ζ)    =    (α_3 - 2α_2)β^3(ζ) 
       +    (3α_2 - α_3)β^2(ζ)   
      + α_1β(ζ) +α_0,

	where, for 0≤ζ≤γ h_n one has
	
    α_0   =_n, α_1=γ𝐳_n, α_2=_n+γ-_n -γ z_n, 
    α_3   =γ( z_n+γ -  z_n), β(ζ)= ζ/γ h_n

and  for γ h_n≤ζ≤ h_n one has instead

    α_0   =_n+γ, α_1=(1-γ) 𝐳_n+γ, α_2=_n+1 - _n+γ - α_1,  
    α_3   =(1-γ)( 𝐳_n+1 - 𝐳_n+γ) β(ζ)=ζ-γ h_n/(1-γ)h_n.

This interpolant allows to have accurate
approximations of the latent variables without extra computational cost or memory storage requirements, since  
it employs the intermediate stages of the TR-BDF2 method in order to achieve higher order accuracy. This contributes
to making  the TR-BDF2 method an attractive choice as the basis for a multirate approach.

 
  

§ LINEAR STABILITY ANALYSIS

  
 In this section,  we will study the linear stability of the multirate algorithm outlined in the
 previous sections in the case of a generic  one step  solver.
Following the approach employed in <cit.> for the θ-method,
 we will consider the special case in which the
macro time step is constant and denoted by h and a single refinement is carried out, with the
refined time step used for the active components given by h/2.
As a result, a single subspace of active variables is introduced,
V⊂𝐑^m, with dimension p < m. 
As in section <ref>, 𝐏 will denote the projector onto  this subspace, while
𝐏^⊥ will denote the corresponding projector onto the complementary
subspace V^⊥  of the latent variables. We will also denote by 
𝐄, 𝐄^⊥ the corresponding natural embedding operators
of V,   V^⊥ into 𝐑^m. As a result,
the operators  Π=𝐄𝐏 and Π^⊥=𝐄^⊥𝐏^⊥ 
represent the zeroing of the components of a given vector belonging to  V,   V^⊥
respectively.  It is to be remarked that this special setting is a major simplification with respect to the generality of the algorithm
described in section <ref>, so that the results of the analysis may not fully explain the behaviour of the algorithm
in practical applications. In particular, large values of the stability function norm in this context do not necessarily imply that
the method is unable to achieve the desired accuracy when applied in its most general form.

We  now consider the linear problem 	'=(t,)=𝐀  and we set h𝐀=. 
We assume that the amplification matrix of  the basic one step solver can be written  as

    𝐑= 𝐑( 𝐙)= 𝐃( 𝐙)^-1𝐍( 𝐙)=𝐃^-1𝐍,

where 𝐃, 𝐍 are polynomials in 𝐙. Notice that this includes most standard one step
solvers, either explicit or implicit. We will also use the notations

    𝐃_1/2  =  𝐃( 𝐙/2 )     𝐍_1/2=𝐍 ( 𝐙/2).


Notice that, in the special case considered here, the interpolation operator can be represented by the application to the data
of a matrix _1/2 whose entries only depend on h  and on the amplification matrix of the single rate method employed.
For example, the linear interpolant (<ref>) is represented in this case by the matrix

    _1/2 = 1/2 (  𝐈+𝐑 ),
 
while for the cubic Hermite interpolant   (<ref>) one has

    𝐐_1/2 =   𝐈 +  βγ𝐙+   β^2𝐅  
    + β^3  𝐆,

where β=1/(2γ) and

    𝐅=3 (𝐑_γ  - 𝐈 -γ𝐙)  -γ𝐙 (𝐑_γ  - 𝐈)
          𝐆=γ𝐙 (𝐑_γ  - 𝐈)
     -2(𝐑_γ  - 𝐈 -γ𝐙),


    𝐑_γ=  (  𝐈 -γ/2𝐙 )^-1 ( 𝐈 +γ/2𝐙 ).

 
 Given these definitions, one step of 
 the self-adjusting multirate algorithm can be rewritten in this special setting as


    𝐏^⊥ u^n+1= 𝐏^⊥𝐑 u^n=𝐏^⊥ u^n,2    𝐏^⊥ u^n,1=𝐏^⊥_1/2 u^n 	
          𝐏𝐃_1/2 (Π u^n,1+ Π^⊥_1/2 u^n
    	 )=  𝐏𝐍_1/2 u^n 	
          𝐏𝐃_1/2 (Π u^n,2+Π^⊥𝐑 u^n ) =  𝐏𝐍_1/2(Π u^n,1+Π^⊥_1/2 u^n ).

Recalling the previous definitions, the equation for 	𝐏 u^n,1 can be rewritten as

    𝐏𝐃_1/2𝐄𝐏 u^n,1= 𝐏𝐍_1/2 u^n
    -𝐏𝐃_1/2Π^⊥_1/2 u^n,

which implies	

    𝐏 u^n,1= (𝐏𝐃_1/2𝐄)^-1 [𝐏𝐍_1/2 
    -𝐏𝐃_1/2Π^⊥_1/2]  u^n.

Substituting this expression in the last equation of (<ref>) one obtains

    𝐏𝐃_1/2𝐄𝐏 u^n,2
    	 =    - 𝐏𝐃_1/2Π^⊥𝐑 u^n
            +𝐏𝐍_1/2𝐄   
    	(𝐏𝐃_1/2𝐄)^-1 [𝐏𝐍_1/2 -𝐏𝐃_1/2Π^⊥_1/2]  u^n 
          +𝐏𝐍_1/2Π^⊥_1/2 u^n .
	
The stability matrix of the complete multirate method	can then be obtained	
deriving from (<ref>) an explicit expression for
𝐏 u^n,2 and substituting in

    u^n+1=	Π u^n,2 +Π^⊥𝐑 u^n.

Introducing  for compactness the notation

    𝐀^∥,∥= 𝐏𝐀𝐄,   𝐀^∥,⊥= 𝐏𝐀𝐄^⊥

one can express the stability matrix as follows 
  
    𝐑_mr   =   𝐄{ (𝐃_1/2^∥,∥)^-1𝐍_1/2^∥,∥ (𝐃_1/2^∥,∥)^-1 [   𝐏𝐍_1/2
    -  𝐃_1/2^∥,⊥𝐏^⊥_1/2] . 
       +   .  (𝐃_1/2^∥,∥)^-1𝐍_1/2^∥,⊥𝐏^⊥_1/2 
    -  (𝐃_1/2^∥,∥)^-1𝐃_1/2^∥,⊥𝐏^⊥𝐑} +Π^⊥𝐑.
		
 
 We will not  attempt a complete theoretical analysis based on the previously derived expression of
 𝐑_mr. However, the same expression can be employed to compute numerically the stability matrix norm and
 spectral radius for relevant model problems. We have done so for three classes of linear ODE systems '=𝐀. 
 Notice that, for simplicity, all the systems have been rearranged so that the latent variables are the first
 of the vector of unknowns.
 
 Firstly, we consider a 2× 2 system  with the same characteristics as those 
 employed in <cit.>, defined by 

    A=[
      [    -1     1; -1000 -1000;       ]].

In particular, in the notation of <cit.>, these coefficient values correspond to the case of κ=1000,  β=-1.
The amplification matrix and its spectral and l_1, l_2, l_∞ matrix norms have been computed for
different values of the time step, up to a maximum value h_max
such that h_max / max|λ( 𝐀)|=100.
Therefore, the maximum time step considered is approximately 100 times larger than that of the customary
stability restriction  for explicit schemes. 
The values of the matrix norms are reported in figures <ref>, <ref>, for the spectral
norm and the other matrix norms, respectively.
It can be observed that no stability problems arise and that the multirate algorithm appears to be insensitive to
the choice of the interpolation method.

 

 
  
  We have then considered the 4× 4 dynamical system defined

    A=[
      [        0        1        0        0; -k_1/m_1     -γ_1        0        0;        0        0        0        1;  k_2/m_2        0 -k_2/m_2     -γ_2 ]].

  The system represents two masses m_1,m_2 such that the first is tied to a wall
  by a spring of elastic constant k_1 and the second is tied to the first by a spring of elastic constant k_2,
  while γ_1, γ_2 represent friction parameters. Similar systems
  have been often used to analyze empirically the stability of  numerical methods for structural mechanics,
  see e.g. <cit.>. We have considered as an example the case
 m_1=m_2=1, k_1=1, k_2=10^6, γ_1=0, γ_2=100.
 The values of the matrix norms are reported in figures <ref>, <ref>, for the spectral
norm and the other matrix norms, respectively.
We observe again that, even though the stability function norm can achieve large
values, the method is still stable in the spectral norm sense. 
Furthermore, it does not appear   to
be sensitive to the choice of the interpolation operator. In figure <ref>, the same quantities
are plotted for the same system in the case γ_2=0. The good stability behaviour displayed in these 
tests seems to justify further work on the application
 of this multirate approach to structural mechanics problems with multiple time scales, such as those
 considered e.g. in <cit.>. As a comparison,
 we also plot in figure  <ref> the corresponding norms for the single rate TR-BDF2 algorithm,
 that display an entirely analogous behaviour.

 
  

 
  
  
   

   
  
 Finally, we have considered a 40× 40 linear ODE system resulting from
 the discretization of a linear heat equation by second finite differences,
 with a diffusivity parameter 10^6 times larger for the last 20 variables than
 for the first 20. All matrix norms are essentially identical to one over the whole range
 of time step values. Analogous results are obtained considering the corresponding
 system for linear advection-diffusion equation, taking in this case centered finite differences
 for the advection discretization and assuming the velocity and diffusivity parameters
 to be 10^4 times larger for the last 20 variables than
 for the first 20.  The values of the matrix norms are reported in figure <ref> for the cubic interpolation case.
 In figure <ref>, the analogous quantities are plotted instead in the case of a pure advection system
 with the same velocity values as in the previous case.
 The good stability properties that result in this case seem to justify further work on the application
 of this multirate approach to hyperbolic equations with multiple time scales, such as those
 considered e.g. in <cit.>.
 
 
  
   
 
  
 
  
 
 
  

§ NUMERICAL EXPERIMENTS

  
 Several benchmark problems already considered in the literature have been employed to assess the performance of
 the multirate algorithm outlined in section <ref>.  The self adjusting multirate approach based on the TR-BDF2 solver has been
 compared to its single rate counterpart, considering several performance indicators, such as
 the CPU time required by the codes,   the number of scalar function calls and
 the number of components being computed at any given time step,  as an estimate of the computational workload for each time step.
 In all the numerical experiments, the TR-BDF2 multirate algorithm employed the cubic Hermite interpolant and a nonlinear solver
 based on the Newton method with approximate computation of the Jacobian to solve the nonlinear systems involved.
	


 §.§ The Inverter Chain Problem


For the first test we consider the inverter chain problem which is an important test problem from the field of electrical circuits. It has also been considered e.g.  in <cit.>, <cit.>.  
The system of equations is given by


    y_1'(t)   =U_op - y_1(t) -Γ g(u(t),y_1(t)),
    
    	y_j'(t)   =U_op -y_j(t) - Γ g(y_j-1(t),y_j(t)),  j=2,3,… ,m,

where y_j, j=1,…,m  represent a chain of m inverters,
 U_op is the operating voltage corresponding to the logical value 1 and Γ serves as a stiffness parameter. 
The function g is defined by 


    g(y,z) = ( max( y-U_τ,0))^2 - ( max( y-z-U_τ,0))^2.

The   parameter values used were m=500, Γ=100, U_τ=1 and U_op=5.
The initial condition was defined as


    y_j(0)={
           6.247× 10^-3for j even
    
           5 for j odd,
    .

while the input signal was given by


    u(t) = {   t-5     for  5≤ t≤10
       5     for  10≤ t ≤15
       5/2(17-t)     for  15≤ t≤17
       0    otherwise.

In this test case, the values of all the components vary in the range (0,5), so that only the absolute tolerance was used. 


Figure <ref> shows the solution of the inverter chain problem described above.   The input signal can be seen to be propagating across the chain of inverters. It is seen that the input signal is also smoothed by the inverter chain. The reference solution was computed by the  MATLAB solver with stringent tolerance and small maximum time step. The same reference solution
was employed to assess the accuracy of the multirate approach with respect to its single rate counterpart.  
 The maximum norm  errors at the final time of the simulation are shown in figure <ref>.
 Even though the   multirate  approach introduces a larger error than that of the
 single rate counterpart,  the errors committed are in general close to the specified tolerances. 







The performance of the algorithms with respect to the different metrics was analyzed for different values of the tolerance. 
Figure <ref> shows the plot of CPU time taken by the multirate and single rate   TR-BDF2 algorithms.  
 At a tolerance of 10^-5,  a speed-up  of a factor 3 was  achieved. Figure <ref> shows  the total number of scalar function evaluations  at different tolerances. Again it can be seen that the performance of the multirate algorithm is better than that of
  the corresponding single rate algorithm. At a tolerance of 10^-5 the number of function evaluations decreased by a factor of 3 for the TR-BDF2.   For a tolerance of 10^-3, the multirate TR-BDF2 made a larger number of function evaluations than it did at tolerance 10^-4. 
  Indeed, if the tolerance is not too stringent, the time step adaptation mechanism causes the algorithm to try larger steps, 
  which makes it harder for the Newton solver to converge.  
  








In figure <ref>, a portion of the so-called space time diagram of the multirate algorithm is shown,
denoting all the active components at each time step. For clarity,  only a zoomed in view  of 200 inverters for the first 500 time steps
is shown. SInce the simulation is started with a very conservative value of the time step,
 the first few time steps involve all the degrees of freedom. The time step is then
 gradually increased by the   adaptation mechanism, so that at a later stage there are only a few components being refined at each time step. 
 The set of refined components moves forward across the inverters as the simulation goes on. This is as expected, because the signal is propagating across the inverter chain and the  active components   are being refined, while the others are being integrated with larger time-steps.  The number of space-time points in the   figure <ref> can be regarded as a measure of the computational workload of the algorithm.  Obviously, for a single rate algorithm this would simply be equal to the number of components times the number of time-steps. Figure <ref> shows this estimate of computational workload for the algorithms for the inverter chain problem. With regard to this metric, using a tolerance of 10^-5 the multirate algorithm leads to a gain in efficiency of 3.4 times with respect to the single rate algorithm.
  
  
  










 


 §.§ A reaction diffusion problem


We then consider a test in which the ODE problem to be solved results from the space discretization
of a partial differential equation. In particular, we consider the reaction diffusion problem also
discussed in <cit.>, whose associated PDE can be written as


    ∂ y/∂ t=ϵ ∂ y^2/∂ x^2 + γ y^2(1-y),

in the domain 0<x<L, 0<t<T. The initial and boundary conditions are given by

    y(0,x)= 1/1+exp(λ(x-1))∂ y(t,0)/∂ x =   ∂ u(t,L)/∂ x  =0.

The values used for parameters were λ=12√(2γϵ), γ=1ϵ=100, L=5 and T=3. 
The equation was discretized in space by simple second order finite differences on a mesh of constant spacing
Δ x=L/N, with different values of N. Even though no analytic solution is available,  the solution is well known
to consist in a   wave moving in the positive direction and connecting the two stable states
of the nonlinear reaction potential. Figure <ref> shows the solution as a function of space at different time instants.









As in the previous section, we have compared the accuracy of the proposed algorithm with respect to the reference solution computed by the MATLAB  solver with tight tolerance  and small maximum time step values.
Figure <ref> shows the plot of the errors at the final simulation time T=3  s for different values of user specified tolerance.
 The performance of the algorithms with respect to the different metrics was analyzed for different values of the tolerance. 
Figure <ref> shows the plot of CPU time and the number of function evaluations required by the multirate and single rate   TR-BDF2 algorithms, respectively.  
 At a tolerance of 10^-5, the speed-up of the multirate algorithm with respect to the single rate algorithm was     2.7. 
 
 



 
 
 

 §.§ Linear Advection Equation

 
 As first example of hyperbolic equations we consider the linear problem


    ∂ u∂ t + ∂ u∂ x = 0,      x∈[-20,20]

with a Gaussian initial datum.
To discretize in space we used a first order upwind method 
with periodic boundary conditions.
The interval time is [0,3] with a number of cells equal to 400.  The relative error tolerance  was
set to 10^-6, while the absolute error tolerance was set to
10^-8. The initial size of the time step was taken to be 10^-2.
The upwind discretization entails a relevant amount of numerical diffusion,  
which is responsible for the spreading of the initial profile. This effect also entails  that the number of computed components increases as time progresses, see Figure <ref>.

 

 














 
 
 
 





























































The CPU times and total number of computed time steps are reported in Table <ref> for
the multirate and single rate TR-BDF2 methods, respectively. For the multirate method, cubic interpolation was found to be necessary in order to achieve an improvement over the performance of the single rate method, while
linear interpolation lead to an excessive number of rejections of macro steps. 
In order to assess also the accuracy of the multirate approach, the relative errors obtained in this test case
are reported in Table <ref>. More specifically, errors obtained by the multirate TR-BDF2  method with cubic Hermite interpolator, the  multirate TR-BDF2  method with linear interpolator  and the  single-rate TR-BDF2 method are compared. It can be seen that all three approaches yield comparable errors.
Since the first order upwind space discretization is responsible for the largest part of the error,
we also compute the error with respect to the reference solution obtained by discretizing in time  with the MATLAB ode45  solver set to a suitably stringent tolerance. The results, reported in
Table <ref>, show that also in this case the errors obtained with the three variants are
entirely analogous, albeit smaller by several orders of magnitude due the factoring out of the spatial discretization error.  

  
 

 


 §.§ Burgers equation

Finally, as a nonlinear test for hyperbolic equations, we consider the Riemann problem for the inviscid Burgers equation 


    ∂ u∂ t + ∂∂ x(12u^2) = 0,      x∈[-1,3],   t>0,

which amounts to assigning a piecewise constant  initial datum
u_0(x) = u_l      x<0,
u_0(x) = u_r      x>0.
We first consider the case  u_l > u_r, for which a shock wave solution is obtained. 
The space discretization is given by   finite volume method with Rusanov numerical flux,
see e.g. <cit.>.
The interval time is [0,1] with a number of cells equal to 400.
The relative error tolerance  was
set to 10^-4, while the absolute error tolerance was set to
10^-6 and   the tolerance for the Newton solver required by the implicit TR-BDF2 method
was taken to be 10^-8.
We took the initial size of the time step equal to 10^-2. For the Riemann problem we used in this first case u_l=1, while u_r=0. 




























































Figure <ref> shows how the shock wave
solution evolves in space as time progresses.
It can be seen how the automatically selected set of
 active components moves along with the shock wave and does not increase in size.
 The same result is apparent from the space time diagram in Figure <ref>.
 The CPU times and total number of computed time steps are reported in Table <ref> for
the multirate and single rate TR-BDF2 methods, respectively. For the multirate method, cubic interpolation was found to be necessary in order to achieve an improvement over the performance of the single rate method, while
linear interpolation lead to an excessive number of rejections of macro steps. 
In order to assess also the accuracy of the multirate approach, the relative errors obtained in this test case
are reported in Table <ref>. More specifically, errors obtained by the multirate TR-BDF2  method with cubic Hermite interpolator, the  multirate TR-BDF2  method with linear interpolator  and the  single-rate TR-BDF2 method are compared. It can be seen that all three approaches yield comparable errors.
Since the first order upwind space discretization is responsible for the largest part of the error,
we also compute the error with respect to the reference solution obtained by discretizing in time  with the MATLAB ode45  solver set to a suitably stringent tolerance. The results, reported in
Table <ref>, show that also in this case the errors obtained with the three variants are
entirely analogous, albeit smaller by several orders of magnitude due the factoring out of the spatial discretization error.  

 
 


In order to highlight the capability of the adaptive strategy to select automatically a time step that
is compatible with the underlying dynamics, we 
  have also reported the maximum Courant numbers at each time step, computed as

    max_i=1,...,N_x|f'(c_i^n)|h_nΔ x

where N_x is the total number of cells. In Figure <ref> a),
the Courant number values associated to the global time steps are represented by red circles,
while those associated to local, refined time steps are represented by blue symbols.
It can be seen that, while   large Courant numbers are used for the latent components,
much smaller values are employed for the active components that correspond to the locations
actually crossed by the shock wave.  In Figure <ref> b) the corresponding fractions of active
components are displayed, highlighting the significant reduction in computational cost in this case.
 






















 We then consider the case u_l > u_r, for which a rarefaction wave solution is obtained.
More specifically, we consider u_l=0 and u_r=1. The other parameters are the same as in the previous case.


























































Figure <ref> shows how the rarefaction wave
solution evolves in space as time progresses.
It can be seen how the automatically selected set of
 active components moves along with the   wave.
 The same result is apparent from the space time diagram in Figure <ref>.
 Since in this case the solution  undergoes significant changes in a  larger number of cells,
 the number of active components is correspondingly larger than in the shock wave case. 
 
 
 The CPU times and total number of computed time steps are reported in Table <ref> for
the multirate and single rate TR-BDF2 methods, respectively. For the multirate method, cubic interpolation was found to be necessary in order to achieve an improvement over the performance of the single rate method, while
linear interpolation lead to an excessive number of rejections of macro steps. 
In order to assess also the accuracy of the multirate approach, the relative errors obtained in this test case
are reported in Table <ref>. More specifically, errors obtained by the multirate TR-BDF2  method with cubic Hermite interpolator, the  multirate TR-BDF2  method with linear interpolator  and the  single-rate TR-BDF2 method are compared. It can be seen that all three approaches yield comparable errors.
Since the first order upwind space discretization is responsible for the largest part of the error,
we also compute the error with respect to the reference solution obtained by discretizing in time  with the MATLAB ode45  solver set to a suitably stringent tolerance. The results, reported in
Table <ref>, show that also in this case the errors obtained with the three variants are
entirely analogous, albeit smaller by several orders of magnitude due the factoring out of the spatial discretization error.  

 

 



In Figure <ref> a)  we have reported the Courant number values associated to each time step. Comparing to the shock wave case, we can see that larger Courant numbers
 are allowed by the error estimator also in the case of refined time steps, as a consequence of the different dynamics of this kind of solution.  In Figure <ref> b) the corresponding fractions of active
components are displayed, highlighting the lesser reduction in computational cost in this case with respect
to the shock wave solution.















  
 

§ CONCLUSIONS

  

 We have  proposed a self-adjusting multirate method  
that relies on the TR-BDF2 method as fundamental single rate solver. 
 The TR-BDF2 method has  a number  of interesting properties, such as cheap error estimation
 via an embedded third order method and  continuous output via a cubic Hermite interpolant, that can
 be exploited to increase efficiency and accuracy of a multirate approach. 
 
 Our  self-adjusting multirate TR-BDF2 method  has been coupled to a
    partitioning and time step selection criterion based on the technique proposed in  <cit.>.  
  The stability of the method has been analyzed in the framework of the classical linear model problem.
  Even though the results achieved were only based on the numerical computation of the stability function norm,
  we show that  for a range of relevant  stiff model problems  the method has remarkable stability properties.
  This is true not only for contractive problems with strictly dissipative behaviour, but also for problems
  with purely  imaginary eigenvalues, which suggests that the method could be advantageous also for
  the application to hyperbolic PDEs and structural mechanics problems.
  The numerical results obtained on several standard benchmarks show that application of the proposed method leads  to significant efficiency gains, that are analogous to those achieved by other self-adjusting approaches when compared   to their corresponding single rate variants. In particular,  the has been applied to the time discretization of nonlinear, hyperbolic partial differential equations, allowing to achieve automatic detection
  of complex localized phenomena such as shock waves and significant reductions in computational cost.
  
  In a companion work, we will focus on the extensive application of the method described in this paper
  to nonlinear hyperbolic equations and on the derivation of mass conservative versions
  of  this multirate approach.
  
  



§ ACKNOWLEDGEMENTS
 This paper contains developments of the results 
presented in the PhD thesis by A.R. <cit.>  and in the Master thesis by L.D. <cit.>, both at Politecnico di Milano.
 A.R.'s PhD grant was supported by the Erasmus Mundus Heritage project, while he is presently supported by the Science Foundation Ireland grant 13/RC/2094.
L.B. was partially supported by the INDAM - GNCS 2015 project Metodi numerici semi-impliciti e semi-Lagrangiani  per sistemi iperbolici di leggi di bilancio.
 Useful discussions with Luca Formaggia are gratefully acknowledged, as well the comments by
  Nicola Guglielmi and Claus Führer on the first draft of the thesis by Akshay Ranade.
   

 
plain




