- Email: [email protected]
Accurate dense output formula for exponential integrators using the scaling and squaring method
Applied Mathematics Letters 43 (2015) 101–107
Contents lists available at ScienceDirect
Applied Mathematics Letters journal homepage: www.elsevier.com/locate/aml
Accurate dense output formula for exponential integrators using the scaling and squaring method Chengshan Wang a , Xiaopeng Fu a , Peng Li a,∗ , Jianzhong Wu b a
Key Laboratory of Smart Grid of Ministry of Education, Tianjin University, Tianjin, China
b
Institute of Energy, School of Engineering, Cardiff University, Cardiff, UK
article
abstract
info
Article history: Received 9 October 2014 Received in revised form 7 December 2014 Accepted 8 December 2014 Available online 17 December 2014 Keywords: Exponential integrator Dense output formula Matrix exponential
This paper proposes an accurate dense output formula for exponential integrators. The computation of matrix exponential function is a vital step in implementing exponential integrators. By scrutinizing the computational process of matrix exponentials using the scaling and squaring method, valuable intermediate results in this process are identified and then used to establish a dense output formula. Efficient computation of dense outputs by the proposed formula enables time integration methods to set their simulation step sizes more flexibly. The efficacy of the proposed formula is verified through numerical examples from the power engineering field. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction Exponential integrators are a class of time integration methods for differential equations and are especially suitable for stiff equations [1]. For the initial value problems u′ (t ) = F (u (t )) ,
u (t0 ) = u0 ,
(1)
these methods perform a linearization procedure on the original problem, u′ (t ) = F (u (t )) = Au (t ) + g (u (t )) ,
(2)
to extract a matrix A that retains the stiffness property. The associated matrix exponential and ϕ functions are embedded into integration schemes, and the nonlinear remainder g (u (t )) is then treated approximately. The ϕ functions are defined as follows [1,2]:
ϕ0 (z ) = ez , ϕk (z ) = z ϕk+1 (z ) + 1/k!.
(3) (4)
For the linear differential equations with polynomial inhomogeneity (as a truncated Taylor expansion of (2)), u′ (t ) = Au (t ) +
p
gk /k! (t − t0 )k ,
u (t0 ) = u0 ,
k=0
∗
Corresponding author. E-mail addresses: [email protected] (C. Wang), [email protected] (X. Fu), [email protected] (P. Li), [email protected] (J. Wu).
http://dx.doi.org/10.1016/j.aml.2014.12.008 0893-9659/© 2014 Elsevier Ltd. All rights reserved.
(5)
102
C. Wang et al. / Applied Mathematics Letters 43 (2015) 101–107
their solution can be expressed exactly using the matrix exponential and the ϕ functions [2], as follows: u (t ) = e(t −t0 )A u0 +
p
(t − t0 )k+1 ϕk+1 ((t − t0 ) A) gk .
(6)
k=0
The introduction of the matrix exponential and ϕ functions into the integration schemes allows for the analogous linear problem (5) to be analytically solved, which thus provides the exponential integrators with better numerical performance over that of traditional methods. However, the computational burden per step also grows. In spite of breakthroughs that have been achieved in recent years [3–5], matrix exponential-related functions are more computationally intense compared to the solution of a linear system that has the same dimension. For applications in the engineering fields, certain supplements to the algorithmic implementation are required to fit the intrinsic properties of exponential integrators to the actual demands of engineering analysis. The development of the dense output formula is one possible approach to facilitate engineering applications of exponential integrators. In engineering applications, the choice of the time integration step size involves consideration of both the numerical requirements, i.e., the deviation of the solution at selected time discretization nodes, and the non-numerical constraints, such as visualization requirements for the simulation results. Since simulation is typically not the final stage of the entire analysis process, the simulation result itself is an input to subsequent analysis stages (‘‘advanced applications’’), such as in a frequency spectrum scan, which often demands a certain density in time discretization. In this case, for methods that offer high accuracy, such as exponential integrators, if the step size is chosen solely from a numerical perspective, an insufficient sampling frequency is provided; if a smaller step size is adopted due to external constraints, then the computational cost is high and it is difficult to compete with low-order methods. Taking power system engineering as an example, the dominant time integration method in this field is the implicit trapezoidal rule [6]. More advanced and recent numerical methods are seldom used. The dense output formula decouples the output step size from the computational step size, which therefore provides flexibility in addressing the above-mentioned dilemma. For an interval [tn , t n +1 ] that is formed by consecutive time discretization nodes (which is now chosen only for the numerical aspect), a dense output formula provides cheap outputs on {t n , k |k = 0, 1, . . . , σ , s.t. t n ,0 = tn , tn,σ = t n +1 , tn,k+1 − tn,k = (t n +1 − tn )/σ }. Going one step further, if the accuracy of these dense outputs does not deteriorate the overall simulation, it can be seen that the simulation step size shrinks to 1/σ , and the per-step computational cost is spread to these σ small steps. This action equivalently boosts the efficiency of the exponential integrators and enables them to gain an edge in computational speed over low-order methods. Dense output is also important in addressing other practical issues, such as event location and the treatment of discontinuities in the simulation. In this paper, a dense output formula that is embedded into the matrix exponential computational process is proposed. Unlike other generic approaches, e.g., interpolation, our approach is tailored for exponential integrators and fully utilizes the peculiar properties of the matrix exponential function. Our approach has a clear physical meaning and an efficient implementation, and it can accurately recover the oscillatory dynamics between time discretization nodes. One of the most widely used exponential integrators, the Exponential Rosenbrock–Euler (ERE) method [7,8], was chosen to demonstrate the proposed approach. 2. Core steps of exponential integrators: ERE method as an example The ERE method is a Rosenbrock-type exponential integrator with a stiff order of 2 [7,8]. For the general nonlinear system (1), the ERE method performs continuous linearization along the numerical solution at each step, i.e.,
∂F ∂F u (t ) = (un ) u (t ) + F (u (t )) − (un ) u (t ) := An u (t ) + gn (u (t )) , ∂u ∂u ′
(7)
where un ∈RN is the numerical approximation of the accurate u (tn ). Using the variation-of-constant formula [1,2], the following matrix exponential function is introduced: (t −tn )An
u (t ) = e
u (tn ) +
t
e(t −τ )An gn (u (τ )) dτ .
(8)
tn
For simplicity, gn (u (τ )) is approximated with gn (un ) over the interval [tn , t n +1 ]. Let hn = t n +1 − tn , the ERE integration scheme is derived as follows: un+1 = ehn An un + hn ϕ1 (hn An ) gn (un ) .
(9)
In practice, a separate computation for e and ϕ1 (hn An ) is not required. Using a technique from [9], the computation of ϕ1 (hn An ) is absorbed into a larger matrix exponential computation. The final integration formula is as follows: hn An
u (tn+1 ) ≈ un+1 = IN
Note that matrix IN
0N ×1 e
hn
An 01×N
gn (un ) un 0 := IN 1
˜
0N ×1 ehn An u ˜ n.
0N ×1 is a vector slicing operator that does not imply real matrix–vector multiplication.
(10)
C. Wang et al. / Applied Mathematics Letters 43 (2015) 101–107
103
˜
For large-scale problems, the matrix exponential ehn An itself is a high-dimensional dense matrix that cannot be easily ˜n hn A
˜ n is computed as a whole using the Krylov subspace approximacomputed. In this case, the matrix–vector product e u tion [9,10]. For arbitrary A∈RN ×N and v ∈RN , the m dimensional Krylov subspace Km is defined as Km = span v , Av , . . . , Am−1 v .
(11)
The standard Arnoldi algorithm returns an orthonormal basis matrix Vm ∈RN ×m and a reduced order matrix Hm = VmT AVm for Km . The following approximation can thus be made: eA v ≈ ∥v ∥ Vm eHm e1 ,
(12)
where e1 is the standard basis vector [1, 0, 0, . . . , 0] , which does not imply real matrix–vector multiplication. The matrix exponential computation eHm on the right-hand side of (12) now has a lower dimension (m ≪ N ). T
˜
˜ n in (10), we again let Hm be the reduced order Applying formula (12) to the matrix exponential-vector product ehn An u ˜ n and let Vm be the basis vector matrix, and we utilize the scaling-invariant property of the Krylov subspace: matrix of A Hm = VmT AVm → α Hm = VmT (α A) Vm ,
(13)
the integration scheme with the Krylov subspace approximation is obtained as follows: 0N ×1 u ˜ n Vm ehn Hm e1 .
un+1 = IN
(14)
The computation of the matrix exponential appears in both the original scheme (10) and the approximated scheme (14), for which the scaling and squaring method is one of the most popular methods [4,5]. This method uses the property eA = (eA/σ )σ . Let σ = 2s , and use scaling to make ∥eA/σ ∥ sufficiently small; eA/σ is then approximated as eA/σ ≈ rpq (A/σ ), where rpq is the [p/q] Padé approximant of the matrix exponential function [4,11]. Function rpq is in form of a rational function, with p and q being the polynomial degree of numerator and denominator, respectively. The matrix exponential eA can then be recovered from the successive squaring of rpq (A/σ ) s times, as follows: eA = eA/σ
σ
2s A ≈ rpq . σ
(15)
3. Dense output formula using the scaling and squaring method 3.1. Dense output formula for low-dimensional problems ˜
Applying the scaling and squaring formula (15) to the matrix exponential ehn An in the integration scheme (10), we obtain
˜
hn
˜n A
ehn An = e σ
σ
2s hn ˜n ≈ rpq A . s
(16)
2
This computational process has s intermediate results, i.e.,
Mi =
hn
˜n A 2s
rpq
2i
hn 2i hn ˜ ˜ A ≈ e 2s An = e 2s−i n ,
i = 0, 1, . . . , s − 1.
(17)
These Mi matrices contain valuable state transition information. For illustration, use a small step size hn /2s−i in (10):
u tn +
hn
2s−i
≈ IN
hn
0N ×1 e 2s−i
˜n A
˜ n ≈ IN u
0N ×1 Mi u ˜ n,
(18)
i.e., the intermediate results {Mi } lead naturally to dense outputs at nodes {tn + hn /2s−i |i = 0, 1, . . . , s − 1 }, with only one ˜ n required per output. Because these dense outputs amount to simulation results of repeatedly multiplication operation Mi u halved step sizes, their accuracy is at least higher than that of the original un+1 . Dense outputs from formula (18) are distributed on the time axis in a 2’s power pattern. To obtain equidistant results, dense outputs un,k at {tn + khn /2s |k ∈ Z, 20 6 k ≤ 2s and k ̸= 2i , i = 0, 1, . . . , s − 1 } must be filled. By formula (10), this goal k
˜ h A
˜ n,k = e 2s n n u˜ n . Each un,k is then obtained by discarding the last element of the corresponding requires computing vectors u ˜ n,k . u Index k can be represented in the binary form as follows: k=
s−1 j =0
kj 2j ,
kj ∈ {0, 1} .
(19)
104
C. Wang et al. / Applied Mathematics Letters 43 (2015) 101–107 Table 1 Summary of computation time. ODE methods
Dense output formula
Comp. step size/ output step size (µs)
Comp. time (s)
TRAP ERE ERE ERE ERE
None None Proposed formula Linear interpolation Hermite interpolation
50/50 50/50 400/50 400/50 400/50
0.304 0.490 0.105 0.091 0.137
On the other hand, because the matrix exponential satisfies the property e(α+β)A = eα A eβ A , we have
e
s −1
k ˜ h A 2s n n
=e
kj ˜ h A s−j n n j=0 2
=
s−1
kj
e 2s−j
˜n hn A
.
(20)
j =0
For every un,k to be computed, there exists l ∈ Z, s.t. 2l−1 < k < 2l ; thus, kl−1 = 1, and for every j ≥ l, kj = 0. From (20), we have k
˜ n,k = e 2s u
˜n hn A
˜n = u
s−1
kj
e 2s−j
˜n hn A
hn
˜ n = e 2s−l+1 u
˜n A
l −2
kj
e 2s−j
˜n hn A
˜ n ≈ Ml−1 u˜ n,k−2l−1 . u
(21)
j =0
j =0
˜ n,k is computed in ascending order of k, then the u˜ n,k−2l−1 vector in (21) has already been computed. Therefore, only one If u ˜ n,k . matrix–vector multiplication is required to obtain the current u In summary, for matrix exponential computation with s successive squaring stages, the proposed dense output formula provides 2s − 1 additional outputs; each consumes one matrix–vector multiplication. 3.2. Dense output formula for large-scale problems When Krylov subspace approximation is used for large-scale problems, the dense output formula described in Section 3.1 requires some adaptation. Again, assume that s successive squaring stages are available in computing ehn Hm in (14). The dense outputs un,k at nodes {tn + khn /2s |k ∈ Z, 20 ≤ k < 2s } can be expressed as follows:
un,k = IN
k
0N ×1 u ˜ n Vm e 2s hn Hm e1 := IN
0N ×1 u ˜ n Vm v˜ n,k .
k h H 2s n m
(22) k ˜ h A 2s n n
˜ n ,k = e ˜ n in Section 3.1, following the The term v˜ n,k = e e1 in (22) can be handled in a similar way as the term u u same deduction process. Briefly speaking, each v˜ n,k requires only a low dimensional (m ≪ N ) matrix–vector multiplication, ˜ n Vm v˜ n,k is required. In summary, the which can be neglected. To obtain un,k , one additional multiplication operation u computational expense in high dimensional case remains practically the same. 4. Numerical examples This section presents two numerical case studies from the electric power system field. The fixed-step trapezoidal method (TRAP) is a widely used time integration method for power system simulation, and is used in this section for comparison. As a benchmark for error analysis, variable step methods with very small error tolerance were used to generate reference results for both cases. 4.1. AC electric motor dynamic simulation The starting dynamics of an AC electric motor were simulated, where numerical solutions at 50 µs spaced grids are required. The model is as follows: x′1 = −x1 /τs1 + x3 /τs2 − ωx2 + (α sin 2ωt + β sin (2ωt + 2π/3) + γ sin (2ωt − 2π/3)) /3 x′2 = −x2 /τs1 + x4 /τs2 + ωx1 + 2 α sin2 ωt + β sin2 (ωt − 2π/3) + γ sin2 (ωt + 2π/3) /3
x′3 = −x3 /τr1 + x1 /τr2 − ωx4 + x4 x5 x′4 = −x4 /τr1 + x2 /τr2 + ωx3 − x3 x5 x′5 = −x5 /τf + µx1 x4 − µx2 x3 y1 = x5 y2 = (x1 cos ωt + x2 sin ωt ) / (Rs τs1 ) − (x3 cos ωt + x4 sin ωt ) / (Rs τs2 ) .
(23)
Simulations have been done using both the TRAP and the ERE methods, with different step sizes and dense output formulas. The computation time and the error curves are summarized in Table 1 and Figs. 1 and 2.
C. Wang et al. / Applied Mathematics Letters 43 (2015) 101–107
105
Fig. 1. Error comparison of the TRAP and the ERE methods, with different simulation step sizes: (a) error comparison for output variable y1 (rotation speed), (b) error comparison for output variable y2 (winding current).
a
b
Fig. 2. Error comparison of three dense output formulas, including the proposed dense output formula, linear interpolation-based (interp1-based) formula, and Hermite interpolation-based (spline-based) formula.
Fig. 1 compares the error curves of the TRAP method and the ERE method. Under the required step size of 50 µs, the ERE method achieves higher accuracy than the TRAP method, as the black curve is two magnitudes lower than the blue curve; however, the consumed computation time is also longer, as shown in Table 1. For engineering analysis where speed is often in higher priority, the original ERE method can hardly compete with the traditional TRAP method. With help of the proposed dense output formula, the ERE method can take larger step size and produce multiple outputs per step to match the required fine output grid. In this case, 400 µs step size is used and 8 outputs are computed per step. As shown in Fig. 1 and Table 1, this enhanced ERE method achieves better accuracy and faster computation speed than the TRAP method at the same time. The proposed dense output formula is further compared with two interpolation-based dense output formulas. The linear interpolation-based formula (interp1-based) is the simplest and fastest dense output formula one can expect, but the generated outputs have low accuracy, as shown in Fig. 2(a). The Hermite interpolation-based formula (spline-based) utilizes the first derivative and achieves good accuracy, but the computation is slower, as it involves linear system solves. The proposed dense output formula achieves the best accuracy among the three, and has high efficiency that is very close to the linear interpolation-based formula, as shown in Fig. 2(b) and Table 1. 4.2. Electric power distribution network simulation The large-scale power distribution network presented in [12] was simulated. The modeled components include distribution transformers, overhead lines, underground cables, compensation devices and residential loads. Eigen analysis revealed high-frequency oscillation modes, suggesting stiffness; thus, a small simulation step size is required to fully capture the dynamic system process. The ERE method with Krylov subspace approximation (14) along with its accompanied dense output formula (22) were used in this case. Fig. 3 displays a segment of simulated high-frequency dynamics. A step size of 5 µs is adopted, and an increasing number of dense outputs are set from subplots (a) to (d). For comparison, the TRAP method with the same step size and the 4th-order Runge–Kutta method (RK4) with a 0.1 µs step size (for numerical stability) are also applied. Fig. 3(a) shows that although the ERE method returns results with sufficient accuracy on discretization nodes without the dense output formula, the 5 µs step size fails to cover the frequency range of the studied dynamics. When dense output functionality is added, finer scale dynamics are captured. With 64 outputs per step (Fig. 3(d)), the exponential integrator’s result nearly coincides with the reference curve. In comparison, the TRAP and the RK4 methods show larger deviations in their results and miss the fine-scale dynamics. To demonstrate the efficiency of the proposed dense output formula, Table 2 lists the computation time for different dense output settings. It is clear that the additional computational burden due to the dense outputs has a low ratio.
106
C. Wang et al. / Applied Mathematics Letters 43 (2015) 101–107
Fig. 3. Distribution network simulation results under different dense output settings: (a) no dense output; (b) dense output setting: 4 outputs/step; (c) dense output setting: 16 outputs/step; (d) dense output setting: 64 outputs/step. Table 2 Computation time for different dense output settings. ODE methods
Dense output setting
Output step size (µs)
Comp. time (s)
ERE ERE (dense output) ERE (dense output) ERE (dense output)
None 4 outputs/step 16 outputs/step 64 outputs/step
5 5/4 5/16 5/64
76.113 76.852 77.672 79.189
Table A.1 Parameter Value
τs1
τs2
τr1
τr2
τf
7.792746e−03
7.993539e−03
8.402578e−03
8.619083e−03
8.019246e+00
Parameter Value
α
β
γ
ω
µ
1.05×469.4855
1.00×469.4855
0.95×469.4855
2π × 60
1.828072e+04
Rs 2.053000e+00
5. Conclusions This paper introduces a dense output formula of exponential integrators by exploiting valuable intermediate results of the scaling and squaring method. Numerical examples are shown that the proposed formula helps the exponential integrators to increase the efficiency and reveal fine-scale dynamic processes when applied to engineering problems. The proposed method is a valuable supplement to exponential integrators in practical applications, especially when numerical solutions are required on find output grids. Acknowledgment This work was supported by the National High Technology Research and Development Program of China (863 Program) (2011AA05A114). Appendix. Parameters for AC electric motor simulation case See Table A.1.
C. Wang et al. / Applied Mathematics Letters 43 (2015) 101–107
107
References [1] M. Hochbruck, A. Ostermann, Exponential integrators, Acta Numer. 19 (2010) 209–286. [2] A.H. Al-Mohy, N.J. Higham, Computing the action of the matrix exponential, with an application to exponential integrators, SIAM J. Sci. Comput. 33 (2011) 488–511. [3] C. Moler, C. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM Rev. 45 (1) (2003) 3–49. [4] N.J. Higham, The scaling and squaring method for the matrix exponential revisited, SIAM Rev. 51 (4) (2009) 747–764. [5] A.H. Al-Mohy, N.J. Higham, A new scaling and squaring algorithm for the matrix exponential, SIAM J. Matrix Anal. Appl. 31 (2009) 970–989. [6] H. Dommel, Digital computer solution of electromagnetic transient in single and multiphase networks, IEEE Trans. Power Appar. Syst. 88 (4) (1969) 388–399. [7] M. Hochbruck, A. Ostermann, J. Schweitzer, Exponential Rosenbrock-type methods, SIAM J. Numer. Anal. 47 (2009) 786–803. [8] M. Caliari, A. Ostermann, Implementation of exponential Rosenbrock-type integrators, Appl. Numer. Math. 59 (2009) 568–581. [9] Y. Saad, Analysis of some Krylov subspace approximations to the matrix exponential operator, SIAM J. Numer. Anal. 29 (1992) 209–228. [10] M. Hochbruck, C. Lubich, On Krylov subspace approximations to the matrix exponential operator, SIAM J. Numer. Anal. 34 (1997) 1911–1925. [11] C. Brezinski, U. Ieea, J. Van Iseghem, A taste of Padé approximation, Acta Numer. 4 (1995) 53–103. [12] W.H. Kersting, Radial distribution test feeders, IEEE Trans. Power Syst. 6 (3) (1991) 975–985.
Contents lists available at ScienceDirect
Applied Mathematics Letters journal homepage: www.elsevier.com/locate/aml
Accurate dense output formula for exponential integrators using the scaling and squaring method Chengshan Wang a , Xiaopeng Fu a , Peng Li a,∗ , Jianzhong Wu b a
Key Laboratory of Smart Grid of Ministry of Education, Tianjin University, Tianjin, China
b
Institute of Energy, School of Engineering, Cardiff University, Cardiff, UK
article
abstract
info
Article history: Received 9 October 2014 Received in revised form 7 December 2014 Accepted 8 December 2014 Available online 17 December 2014 Keywords: Exponential integrator Dense output formula Matrix exponential
This paper proposes an accurate dense output formula for exponential integrators. The computation of matrix exponential function is a vital step in implementing exponential integrators. By scrutinizing the computational process of matrix exponentials using the scaling and squaring method, valuable intermediate results in this process are identified and then used to establish a dense output formula. Efficient computation of dense outputs by the proposed formula enables time integration methods to set their simulation step sizes more flexibly. The efficacy of the proposed formula is verified through numerical examples from the power engineering field. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction Exponential integrators are a class of time integration methods for differential equations and are especially suitable for stiff equations [1]. For the initial value problems u′ (t ) = F (u (t )) ,
u (t0 ) = u0 ,
(1)
these methods perform a linearization procedure on the original problem, u′ (t ) = F (u (t )) = Au (t ) + g (u (t )) ,
(2)
to extract a matrix A that retains the stiffness property. The associated matrix exponential and ϕ functions are embedded into integration schemes, and the nonlinear remainder g (u (t )) is then treated approximately. The ϕ functions are defined as follows [1,2]:
ϕ0 (z ) = ez , ϕk (z ) = z ϕk+1 (z ) + 1/k!.
(3) (4)
For the linear differential equations with polynomial inhomogeneity (as a truncated Taylor expansion of (2)), u′ (t ) = Au (t ) +
p
gk /k! (t − t0 )k ,
u (t0 ) = u0 ,
k=0
∗
Corresponding author. E-mail addresses: [email protected] (C. Wang), [email protected] (X. Fu), [email protected] (P. Li), [email protected] (J. Wu).
http://dx.doi.org/10.1016/j.aml.2014.12.008 0893-9659/© 2014 Elsevier Ltd. All rights reserved.
(5)
102
C. Wang et al. / Applied Mathematics Letters 43 (2015) 101–107
their solution can be expressed exactly using the matrix exponential and the ϕ functions [2], as follows: u (t ) = e(t −t0 )A u0 +
p
(t − t0 )k+1 ϕk+1 ((t − t0 ) A) gk .
(6)
k=0
The introduction of the matrix exponential and ϕ functions into the integration schemes allows for the analogous linear problem (5) to be analytically solved, which thus provides the exponential integrators with better numerical performance over that of traditional methods. However, the computational burden per step also grows. In spite of breakthroughs that have been achieved in recent years [3–5], matrix exponential-related functions are more computationally intense compared to the solution of a linear system that has the same dimension. For applications in the engineering fields, certain supplements to the algorithmic implementation are required to fit the intrinsic properties of exponential integrators to the actual demands of engineering analysis. The development of the dense output formula is one possible approach to facilitate engineering applications of exponential integrators. In engineering applications, the choice of the time integration step size involves consideration of both the numerical requirements, i.e., the deviation of the solution at selected time discretization nodes, and the non-numerical constraints, such as visualization requirements for the simulation results. Since simulation is typically not the final stage of the entire analysis process, the simulation result itself is an input to subsequent analysis stages (‘‘advanced applications’’), such as in a frequency spectrum scan, which often demands a certain density in time discretization. In this case, for methods that offer high accuracy, such as exponential integrators, if the step size is chosen solely from a numerical perspective, an insufficient sampling frequency is provided; if a smaller step size is adopted due to external constraints, then the computational cost is high and it is difficult to compete with low-order methods. Taking power system engineering as an example, the dominant time integration method in this field is the implicit trapezoidal rule [6]. More advanced and recent numerical methods are seldom used. The dense output formula decouples the output step size from the computational step size, which therefore provides flexibility in addressing the above-mentioned dilemma. For an interval [tn , t n +1 ] that is formed by consecutive time discretization nodes (which is now chosen only for the numerical aspect), a dense output formula provides cheap outputs on {t n , k |k = 0, 1, . . . , σ , s.t. t n ,0 = tn , tn,σ = t n +1 , tn,k+1 − tn,k = (t n +1 − tn )/σ }. Going one step further, if the accuracy of these dense outputs does not deteriorate the overall simulation, it can be seen that the simulation step size shrinks to 1/σ , and the per-step computational cost is spread to these σ small steps. This action equivalently boosts the efficiency of the exponential integrators and enables them to gain an edge in computational speed over low-order methods. Dense output is also important in addressing other practical issues, such as event location and the treatment of discontinuities in the simulation. In this paper, a dense output formula that is embedded into the matrix exponential computational process is proposed. Unlike other generic approaches, e.g., interpolation, our approach is tailored for exponential integrators and fully utilizes the peculiar properties of the matrix exponential function. Our approach has a clear physical meaning and an efficient implementation, and it can accurately recover the oscillatory dynamics between time discretization nodes. One of the most widely used exponential integrators, the Exponential Rosenbrock–Euler (ERE) method [7,8], was chosen to demonstrate the proposed approach. 2. Core steps of exponential integrators: ERE method as an example The ERE method is a Rosenbrock-type exponential integrator with a stiff order of 2 [7,8]. For the general nonlinear system (1), the ERE method performs continuous linearization along the numerical solution at each step, i.e.,
∂F ∂F u (t ) = (un ) u (t ) + F (u (t )) − (un ) u (t ) := An u (t ) + gn (u (t )) , ∂u ∂u ′
(7)
where un ∈RN is the numerical approximation of the accurate u (tn ). Using the variation-of-constant formula [1,2], the following matrix exponential function is introduced: (t −tn )An
u (t ) = e
u (tn ) +
t
e(t −τ )An gn (u (τ )) dτ .
(8)
tn
For simplicity, gn (u (τ )) is approximated with gn (un ) over the interval [tn , t n +1 ]. Let hn = t n +1 − tn , the ERE integration scheme is derived as follows: un+1 = ehn An un + hn ϕ1 (hn An ) gn (un ) .
(9)
In practice, a separate computation for e and ϕ1 (hn An ) is not required. Using a technique from [9], the computation of ϕ1 (hn An ) is absorbed into a larger matrix exponential computation. The final integration formula is as follows: hn An
u (tn+1 ) ≈ un+1 = IN
Note that matrix IN
0N ×1 e
hn
An 01×N
gn (un ) un 0 := IN 1
˜
0N ×1 ehn An u ˜ n.
0N ×1 is a vector slicing operator that does not imply real matrix–vector multiplication.
(10)
C. Wang et al. / Applied Mathematics Letters 43 (2015) 101–107
103
˜
For large-scale problems, the matrix exponential ehn An itself is a high-dimensional dense matrix that cannot be easily ˜n hn A
˜ n is computed as a whole using the Krylov subspace approximacomputed. In this case, the matrix–vector product e u tion [9,10]. For arbitrary A∈RN ×N and v ∈RN , the m dimensional Krylov subspace Km is defined as Km = span v , Av , . . . , Am−1 v .
(11)
The standard Arnoldi algorithm returns an orthonormal basis matrix Vm ∈RN ×m and a reduced order matrix Hm = VmT AVm for Km . The following approximation can thus be made: eA v ≈ ∥v ∥ Vm eHm e1 ,
(12)
where e1 is the standard basis vector [1, 0, 0, . . . , 0] , which does not imply real matrix–vector multiplication. The matrix exponential computation eHm on the right-hand side of (12) now has a lower dimension (m ≪ N ). T
˜
˜ n in (10), we again let Hm be the reduced order Applying formula (12) to the matrix exponential-vector product ehn An u ˜ n and let Vm be the basis vector matrix, and we utilize the scaling-invariant property of the Krylov subspace: matrix of A Hm = VmT AVm → α Hm = VmT (α A) Vm ,
(13)
the integration scheme with the Krylov subspace approximation is obtained as follows: 0N ×1 u ˜ n Vm ehn Hm e1 .
un+1 = IN
(14)
The computation of the matrix exponential appears in both the original scheme (10) and the approximated scheme (14), for which the scaling and squaring method is one of the most popular methods [4,5]. This method uses the property eA = (eA/σ )σ . Let σ = 2s , and use scaling to make ∥eA/σ ∥ sufficiently small; eA/σ is then approximated as eA/σ ≈ rpq (A/σ ), where rpq is the [p/q] Padé approximant of the matrix exponential function [4,11]. Function rpq is in form of a rational function, with p and q being the polynomial degree of numerator and denominator, respectively. The matrix exponential eA can then be recovered from the successive squaring of rpq (A/σ ) s times, as follows: eA = eA/σ
σ
2s A ≈ rpq . σ
(15)
3. Dense output formula using the scaling and squaring method 3.1. Dense output formula for low-dimensional problems ˜
Applying the scaling and squaring formula (15) to the matrix exponential ehn An in the integration scheme (10), we obtain
˜
hn
˜n A
ehn An = e σ
σ
2s hn ˜n ≈ rpq A . s
(16)
2
This computational process has s intermediate results, i.e.,
Mi =
hn
˜n A 2s
rpq
2i
hn 2i hn ˜ ˜ A ≈ e 2s An = e 2s−i n ,
i = 0, 1, . . . , s − 1.
(17)
These Mi matrices contain valuable state transition information. For illustration, use a small step size hn /2s−i in (10):
u tn +
hn
2s−i
≈ IN
hn
0N ×1 e 2s−i
˜n A
˜ n ≈ IN u
0N ×1 Mi u ˜ n,
(18)
i.e., the intermediate results {Mi } lead naturally to dense outputs at nodes {tn + hn /2s−i |i = 0, 1, . . . , s − 1 }, with only one ˜ n required per output. Because these dense outputs amount to simulation results of repeatedly multiplication operation Mi u halved step sizes, their accuracy is at least higher than that of the original un+1 . Dense outputs from formula (18) are distributed on the time axis in a 2’s power pattern. To obtain equidistant results, dense outputs un,k at {tn + khn /2s |k ∈ Z, 20 6 k ≤ 2s and k ̸= 2i , i = 0, 1, . . . , s − 1 } must be filled. By formula (10), this goal k
˜ h A
˜ n,k = e 2s n n u˜ n . Each un,k is then obtained by discarding the last element of the corresponding requires computing vectors u ˜ n,k . u Index k can be represented in the binary form as follows: k=
s−1 j =0
kj 2j ,
kj ∈ {0, 1} .
(19)
104
C. Wang et al. / Applied Mathematics Letters 43 (2015) 101–107 Table 1 Summary of computation time. ODE methods
Dense output formula
Comp. step size/ output step size (µs)
Comp. time (s)
TRAP ERE ERE ERE ERE
None None Proposed formula Linear interpolation Hermite interpolation
50/50 50/50 400/50 400/50 400/50
0.304 0.490 0.105 0.091 0.137
On the other hand, because the matrix exponential satisfies the property e(α+β)A = eα A eβ A , we have
e
s −1
k ˜ h A 2s n n
=e
kj ˜ h A s−j n n j=0 2
=
s−1
kj
e 2s−j
˜n hn A
.
(20)
j =0
For every un,k to be computed, there exists l ∈ Z, s.t. 2l−1 < k < 2l ; thus, kl−1 = 1, and for every j ≥ l, kj = 0. From (20), we have k
˜ n,k = e 2s u
˜n hn A
˜n = u
s−1
kj
e 2s−j
˜n hn A
hn
˜ n = e 2s−l+1 u
˜n A
l −2
kj
e 2s−j
˜n hn A
˜ n ≈ Ml−1 u˜ n,k−2l−1 . u
(21)
j =0
j =0
˜ n,k is computed in ascending order of k, then the u˜ n,k−2l−1 vector in (21) has already been computed. Therefore, only one If u ˜ n,k . matrix–vector multiplication is required to obtain the current u In summary, for matrix exponential computation with s successive squaring stages, the proposed dense output formula provides 2s − 1 additional outputs; each consumes one matrix–vector multiplication. 3.2. Dense output formula for large-scale problems When Krylov subspace approximation is used for large-scale problems, the dense output formula described in Section 3.1 requires some adaptation. Again, assume that s successive squaring stages are available in computing ehn Hm in (14). The dense outputs un,k at nodes {tn + khn /2s |k ∈ Z, 20 ≤ k < 2s } can be expressed as follows:
un,k = IN
k
0N ×1 u ˜ n Vm e 2s hn Hm e1 := IN
0N ×1 u ˜ n Vm v˜ n,k .
k h H 2s n m
(22) k ˜ h A 2s n n
˜ n ,k = e ˜ n in Section 3.1, following the The term v˜ n,k = e e1 in (22) can be handled in a similar way as the term u u same deduction process. Briefly speaking, each v˜ n,k requires only a low dimensional (m ≪ N ) matrix–vector multiplication, ˜ n Vm v˜ n,k is required. In summary, the which can be neglected. To obtain un,k , one additional multiplication operation u computational expense in high dimensional case remains practically the same. 4. Numerical examples This section presents two numerical case studies from the electric power system field. The fixed-step trapezoidal method (TRAP) is a widely used time integration method for power system simulation, and is used in this section for comparison. As a benchmark for error analysis, variable step methods with very small error tolerance were used to generate reference results for both cases. 4.1. AC electric motor dynamic simulation The starting dynamics of an AC electric motor were simulated, where numerical solutions at 50 µs spaced grids are required. The model is as follows: x′1 = −x1 /τs1 + x3 /τs2 − ωx2 + (α sin 2ωt + β sin (2ωt + 2π/3) + γ sin (2ωt − 2π/3)) /3 x′2 = −x2 /τs1 + x4 /τs2 + ωx1 + 2 α sin2 ωt + β sin2 (ωt − 2π/3) + γ sin2 (ωt + 2π/3) /3
x′3 = −x3 /τr1 + x1 /τr2 − ωx4 + x4 x5 x′4 = −x4 /τr1 + x2 /τr2 + ωx3 − x3 x5 x′5 = −x5 /τf + µx1 x4 − µx2 x3 y1 = x5 y2 = (x1 cos ωt + x2 sin ωt ) / (Rs τs1 ) − (x3 cos ωt + x4 sin ωt ) / (Rs τs2 ) .
(23)
Simulations have been done using both the TRAP and the ERE methods, with different step sizes and dense output formulas. The computation time and the error curves are summarized in Table 1 and Figs. 1 and 2.
C. Wang et al. / Applied Mathematics Letters 43 (2015) 101–107
105
Fig. 1. Error comparison of the TRAP and the ERE methods, with different simulation step sizes: (a) error comparison for output variable y1 (rotation speed), (b) error comparison for output variable y2 (winding current).
a
b
Fig. 2. Error comparison of three dense output formulas, including the proposed dense output formula, linear interpolation-based (interp1-based) formula, and Hermite interpolation-based (spline-based) formula.
Fig. 1 compares the error curves of the TRAP method and the ERE method. Under the required step size of 50 µs, the ERE method achieves higher accuracy than the TRAP method, as the black curve is two magnitudes lower than the blue curve; however, the consumed computation time is also longer, as shown in Table 1. For engineering analysis where speed is often in higher priority, the original ERE method can hardly compete with the traditional TRAP method. With help of the proposed dense output formula, the ERE method can take larger step size and produce multiple outputs per step to match the required fine output grid. In this case, 400 µs step size is used and 8 outputs are computed per step. As shown in Fig. 1 and Table 1, this enhanced ERE method achieves better accuracy and faster computation speed than the TRAP method at the same time. The proposed dense output formula is further compared with two interpolation-based dense output formulas. The linear interpolation-based formula (interp1-based) is the simplest and fastest dense output formula one can expect, but the generated outputs have low accuracy, as shown in Fig. 2(a). The Hermite interpolation-based formula (spline-based) utilizes the first derivative and achieves good accuracy, but the computation is slower, as it involves linear system solves. The proposed dense output formula achieves the best accuracy among the three, and has high efficiency that is very close to the linear interpolation-based formula, as shown in Fig. 2(b) and Table 1. 4.2. Electric power distribution network simulation The large-scale power distribution network presented in [12] was simulated. The modeled components include distribution transformers, overhead lines, underground cables, compensation devices and residential loads. Eigen analysis revealed high-frequency oscillation modes, suggesting stiffness; thus, a small simulation step size is required to fully capture the dynamic system process. The ERE method with Krylov subspace approximation (14) along with its accompanied dense output formula (22) were used in this case. Fig. 3 displays a segment of simulated high-frequency dynamics. A step size of 5 µs is adopted, and an increasing number of dense outputs are set from subplots (a) to (d). For comparison, the TRAP method with the same step size and the 4th-order Runge–Kutta method (RK4) with a 0.1 µs step size (for numerical stability) are also applied. Fig. 3(a) shows that although the ERE method returns results with sufficient accuracy on discretization nodes without the dense output formula, the 5 µs step size fails to cover the frequency range of the studied dynamics. When dense output functionality is added, finer scale dynamics are captured. With 64 outputs per step (Fig. 3(d)), the exponential integrator’s result nearly coincides with the reference curve. In comparison, the TRAP and the RK4 methods show larger deviations in their results and miss the fine-scale dynamics. To demonstrate the efficiency of the proposed dense output formula, Table 2 lists the computation time for different dense output settings. It is clear that the additional computational burden due to the dense outputs has a low ratio.
106
C. Wang et al. / Applied Mathematics Letters 43 (2015) 101–107
Fig. 3. Distribution network simulation results under different dense output settings: (a) no dense output; (b) dense output setting: 4 outputs/step; (c) dense output setting: 16 outputs/step; (d) dense output setting: 64 outputs/step. Table 2 Computation time for different dense output settings. ODE methods
Dense output setting
Output step size (µs)
Comp. time (s)
ERE ERE (dense output) ERE (dense output) ERE (dense output)
None 4 outputs/step 16 outputs/step 64 outputs/step
5 5/4 5/16 5/64
76.113 76.852 77.672 79.189
Table A.1 Parameter Value
τs1
τs2
τr1
τr2
τf
7.792746e−03
7.993539e−03
8.402578e−03
8.619083e−03
8.019246e+00
Parameter Value
α
β
γ
ω
µ
1.05×469.4855
1.00×469.4855
0.95×469.4855
2π × 60
1.828072e+04
Rs 2.053000e+00
5. Conclusions This paper introduces a dense output formula of exponential integrators by exploiting valuable intermediate results of the scaling and squaring method. Numerical examples are shown that the proposed formula helps the exponential integrators to increase the efficiency and reveal fine-scale dynamic processes when applied to engineering problems. The proposed method is a valuable supplement to exponential integrators in practical applications, especially when numerical solutions are required on find output grids. Acknowledgment This work was supported by the National High Technology Research and Development Program of China (863 Program) (2011AA05A114). Appendix. Parameters for AC electric motor simulation case See Table A.1.
C. Wang et al. / Applied Mathematics Letters 43 (2015) 101–107
107
References [1] M. Hochbruck, A. Ostermann, Exponential integrators, Acta Numer. 19 (2010) 209–286. [2] A.H. Al-Mohy, N.J. Higham, Computing the action of the matrix exponential, with an application to exponential integrators, SIAM J. Sci. Comput. 33 (2011) 488–511. [3] C. Moler, C. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM Rev. 45 (1) (2003) 3–49. [4] N.J. Higham, The scaling and squaring method for the matrix exponential revisited, SIAM Rev. 51 (4) (2009) 747–764. [5] A.H. Al-Mohy, N.J. Higham, A new scaling and squaring algorithm for the matrix exponential, SIAM J. Matrix Anal. Appl. 31 (2009) 970–989. [6] H. Dommel, Digital computer solution of electromagnetic transient in single and multiphase networks, IEEE Trans. Power Appar. Syst. 88 (4) (1969) 388–399. [7] M. Hochbruck, A. Ostermann, J. Schweitzer, Exponential Rosenbrock-type methods, SIAM J. Numer. Anal. 47 (2009) 786–803. [8] M. Caliari, A. Ostermann, Implementation of exponential Rosenbrock-type integrators, Appl. Numer. Math. 59 (2009) 568–581. [9] Y. Saad, Analysis of some Krylov subspace approximations to the matrix exponential operator, SIAM J. Numer. Anal. 29 (1992) 209–228. [10] M. Hochbruck, C. Lubich, On Krylov subspace approximations to the matrix exponential operator, SIAM J. Numer. Anal. 34 (1997) 1911–1925. [11] C. Brezinski, U. Ieea, J. Van Iseghem, A taste of Padé approximation, Acta Numer. 4 (1995) 53–103. [12] W.H. Kersting, Radial distribution test feeders, IEEE Trans. Power Syst. 6 (3) (1991) 975–985.