1

Solving nonlinear circuits with pulsed excitation by multirate
partial differential equations
Andreas Pels1 , Johan Gyselinck2 , Ruth V. Sabariego3 , and Sebastian Schöps1
School of Computational Engineering and Institut für Theorie Elektromagnetischer Felder,
Technische Universität Darmstadt, Germany
2 BEAMS Department, Université libre de Bruxelles, Belgium
3 Department of Electrical Engineering, EnergyVille, KU Leuven, Belgium

In this paper the concept of Multirate Partial Differential Equations (MPDEs) is applied to obtain an efficient solution
for nonlinear low-frequency electrical circuits with pulsed excitation. The MPDEs are solved by a Galerkin approach and a
conventional time discretization. Nonlinearities are efficiently accounted for by neglecting the high-frequency components
(ripples) of the state variables and using only their envelope for the evaluation. It is shown that the impact of this
approximation on the solution becomes increasingly negligible for rising frequency and leads to significant performance
gains.
Index Terms—Finite element analysis, Nonlinear circuits, Numerical simulation, Multirate partial differential equations.

I. I NTRODUCTION
ULTISCALE and multirate problems occur naturally in many applications from electrical engineering. Classical discretization schemes are often inefficient in these cases and it is preferable to address
the dynamics of each scale separately. In time domain
the multirate phenomenon is often characterized by the
fact that some solution components are active while the
majority is latent (e.g. behind a low-pass filter [1]) or by
problems with oscillatory solutions that are composed of
multiple frequencies. An example is depicted in Fig. 1, it
consists of a fast periodically varying ripple and a slowly
varying envelope.
Multirate Partial Differential Equations (MPDEs) have
been successfully applied in nonlinear high-frequency
applications with largely separated time scales [2], [3].
Various methods have been proposed for the numerical
solution of the MPDEs, e.g. harmonic balance [2], shooting methods, classical time stepping [4] or a combination
of both [5]. Instead of the Fourier basis functions, one can
also use classical nodal basis functions or more sophisticated problem-specific basis functions as for example the
pulse width modulation (PWM) basis functions [6], [7].
This contribution focuses on improving the efficiency
of the numerical solution of nonlinear MPDEs with
a high-frequent pulsed excitation. In contrast to prior
works, we propose neglecting the influence of highfrequency components on the nonlinearity. Consequently,
the nonlinear relation is evaluated using only the envelope
of the state variables.

M

Corresponding author: A. Pels (e-mail: pels@gsc.tu-darmstadt.de).

voltage (V) / current (A)

arXiv:1710.06278v1 [] 17 Oct 2017

1 Graduate

100

50
Current iL,ref
Voltage vC,ref
Excitation vi

0
0

2

4

6

8

10

time t (ms)
Fig. 1. Solution of the buck converter shown in Fig. 2 at fs = 1000 Hz.
It consists of a slowly varying envelope and fast ripples. The switching
cycle of the pulsed input voltage is Ts = 1 ms, the duty cycle D = 0.7.

The paper is structured as follows: after this introduction, Section II presents the multirate formulation. Section III discusses the Galerkin approach in time domain
and the extraction of the envelope to evaluate nonlinear
behavior. Section IV gathers numerical results based on
the buck converter benchmark example and discusses the
accuracy and efficiency of the modified method. Finally,
conclusions are drawn.
II. M ULTIRATE FORMULATION
Spatial discretization of low-frequency field formulations as, e.g. electro- or magneto-quasi-statics [8], or
network models of power converter circuits as, e.g. the
buck converter in Fig. 2 lead to (nonlinear) systems of
ordinary or differential algebraic equations of the form
A x(t)

d

x(t) + B x(t) x(t) = c(t)
dt

(1)

c 2017 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media,
including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution
to servers or lists, or reuse of any copyrighted component of this work in other works.

2

L

RL

iL

vi

vC

C R

vC (V)

80
Fig. 2. Simplified buck converter.

with an initial condition x(t0 ) = x0 , where x(t) ∈ RNs is
the vector of Ns state variables, A(x), B(x) ∈ RNs ×Ns
are matrices that may depend on the solution and c(t) ∈
RNs is the excitation vector.
The system of equations (1) is hereafter rewritten as
MPDEs [2], [3], [7] in terms of the two time scales t1
and t2


∂b
x
∂b
x
b(t1 , t2 ) = b
+
+ B(b
x) x
c(t1 , t2 ) , (2)
A(b
x)
∂t1
∂t2
b(t1 , t2 ) and b
where x
c(t1 , t2 ) are the multivariate forms
of x(t) and c(t). If b
c(t1 , t2 ) fulfills the relation b
c(t, t) =
c(t), the solution of the original problem can be exb(t1 , t2 ) of the MPDEs by
tracted from the solution x
b(t, t) [2].
x(t) = x
Let t1 denote the slow time scale and t2 the fast
time scale, which is furthermore assumed to be periodic.
Without limiting the generality of the approach, we
choose b
c(t1 , t2 ) := c(t2 ) such that the right-hand-side
only depends on the fast time scale [7].
III. N UMERICAL METHOD
System (2) can be numerically solved either using a
Galerkin framework, shooting methods, classical time
stepping or a combination. Here, we propose to use a
variational setting, i.e. a Galerkin approach, for the fast
time scale and a conventional time stepping for the slowly
varying envelope.
A. Galerkin in time domain
We represent the solution by an expansion of Np +
1 suitable basis functions pk (τ (t2 )) and coefficients
bj (t1 , t2 ) can
wj,k (t1 ). The approximated state variables x
be written as the series
bj (t1 , t2 ) =
x

Np
X

pk (τ (t2 ))wj,k (t1 ),

(3)

k=0

with τ (t2 ) = Tt2s mod 1, where Ts is the switching cycle
of the excitation and mod denotes the modulo operation.
Since we deal with pulsed right-hand-sides, we choose
the PWM basis functions of [7], although constructed for

8

40
6
0
0

4
2

4

2

6

8

t2 (ms)

t1 (ms)

0

Fig. 3. Illustration of the multirate solution for fs = 500 Hz. The
solution of the original system of differential equations is denoted in
black.

linear problems [6]. The zero-th basis function is constant
p0 (τ ) = 1 and the corresponding coefficient wj,0 defines
the envelope of the j-th solution component. The first
basis function is defined by
( √
−D
if 0 ≤ τ ≤ D
3 2τD
p1 (τ ) =
,
√ 1+D−2τ
3 1−D
if D ≤ τ ≤ 1
where D is a free parameter, which can be chosen
according to the duty cycle of the PWM, i.e. between
0 and 1. The basis functions of higher degree pk (τ ), 2 ≤
k ≤ Np are obtained recursively by integrating pk−1 (τ )
and orthonormalizing with respect to the L2 (0, 1) scalar
product, [6], [7].
Finally, a Galerkin approach is applied on the interval
of one period, which yields
ZTs
A(b
x)
0

 ∂b
x
∂t1

+

∂b
x
∂t2

(4)


b−b
+B(b
x) x
c pl (τ (t2 )) dt2 = 0,
for all l = 0, . . . , Np . The integration with respect to t2
leads to a system of differential (algebraic) equations in
t1 , which can be solved by conventional time integration.
B. Treatment of nonlinearity
During time integration of (4), the integrals have to
be evaluated every time step due to their nonlinear
dependency on the solution. This may lead to an unnecessary increase of the computational effort as the
ripple components are often small in comparison with the
magnitude of the envelope. Consequently, one can save
computational time by ignoring the ripple components
of the solution and only using its envelope w̄ for the
evaluation of the nonlinearity.

As mentioned before, the envelope is stored in the
vector of coefficients w(t1 ) given by
w = [w1,0 , . . . , w1,Np , w2,0 , . . . , wNs ,Np ]> .

(5)

Let us abstractly define a function f , which extracts the
envelope w̄(t1 ) from w(t1 ), i.e.

>

w̄ = [w1,0 , w2,0 , . . . , wNs ,0 ] .

(7)

Therefore the matrices A, B only depend on t1 and
the evaluation of (4) simplifies significantly. Equation (4)
becomes
A(w̄)

0

∂t1

+

(8)

dw
+ B(w(t1 )) w(t1 ) = C(t1 ) ,
dt1

(11)

where the matrices are given by
A(w) = A(f (w)) ⊗ I,

(12)

B(w) = B(f (w)) ⊗ I + A(f (w)) ⊗ Q,

(13)

and ⊗ denotes the Kronecker product. The right-handside vector is given by

C(t1 ) =
0


b
c1 (t1 , t2 ) p(τ (t2 ))


..

 dt2 .
.
b
cNs (t1 , t2 ) p(τ (t2 ))

100

150

current iL (A)
Fig. 4. Inductance of nonlinear coil versus current iL through the coil.
original
simplified

10−2
10−4
10−6
103

104

105

Fig. 5. Error  of the original and simplified MPDE approach (Np = 4)
versus frequency fs .

IV. N UMERICAL R ESULTS

in equation (8), we get

Ts

50

frequency fs (Hz)

for all l = 0, . . . , Np , where the matrices A and B are
independent of t2 . Introducing
Z 1
I = Ts
p(τ ) p>(τ ) dτ ,
(9)
0
Z 1
dp >
p (τ ) dτ ,
(10)
Q =−
0 dτ

Z

0

∂b
x

∂t2

b−b
+B(w̄) x
c pl (τ (t2 )) dt2 = 0,

A(w(t1 ))

0.5

0

error 

which is in the case of PWM basis functions the zero-th
components, i.e.

 ∂b
x

1

(6)

w̄(t1 ) = f (w(t1 )),

ZTs

inductance (mH)

3

(14)

Eventually, system (11) can be more efficiently solved by
conventional time discretization in the sense that larger
time steps can be used than for the original problem (1).
However, drawbacks are the approximation of the nonlinearity and the increased size of the matrices, i.e. A(w),
B(w) ∈ RNs (Np +1)×Ns (Np +1) and C(t) ∈ RNs (Np +1) ,
which is a similar tradeoff as in harmonic balance [2].

The numerical tests are performed on the simplified
buck converter model [6] using a nonlinear coil, whose
characteristic is shown in Fig. 4. The code is implemented
in GNU Octave [9]; for time integration the high-order
implicit Runge-Kutta method Radau5 from odepkg is
used, [10], [11]. As basis functions pk (τ ) we choose the
problem-specific PWM basis functions introduced earlier.
The reference solution to which all results are compared is calculated directly by solving (1) with a very
accurate time discretization (tol = 10−12 ). The buck
converter is operated in the range of frequencies from
500 Hz to 100 kHz. To determine the accuracy of the
MPDE approach the relative L2 -error of the buck converter voltage
=

kvC,ref − vC kL2 (Ω)
kvC,ref kL2 (Ω)

(15)

with respect to t ∈ Ω = [0, 10] ms is approximated by
numerical quadrature. This calculation is performed for
each frequency fs = T1s .
The MPDE solution is expanded using Np = 4
basis functions and (11) is solved using a tolerance of
tol = 10−6 . The result is exemplary depicted in Fig. 3.
Fig. 5 shows the error of the approach without (denoted
as “original approach”) and with the simplified evaluation
of the nonlinearity (8) (denoted as “simplified approach”)
with respect to the frequency. Without the simplification,
the integrals in (4) are evaluated using Gauss-Kronrod
quadrature and lead to an accuracy of  < 10−4 for
frequencies fs > 10 kHz. As expected, the higher the
frequency, the higher the accuracy of the method since

102
101

original
simplified

100
103

104

105

frequency fs (Hz)
Fig. 6. Computational time for the original MPDE formulation and the
one with simplified evaluation of nonlinearity (Np = 4).

the magnitude of the ripples in relation to the envelope
decreases. The simplified approach introduces an additional error due to the approximation of the integrals.
However, the accuracy for fs > 10 kHz, i.e.,  < 10−3 ,
is still sufficient for most applications. Fig. 7 shows
the solutions of reference, simplified and original MPDE
approach versus time at a low frequency fs = 1000 Hz.
Here, the error committed in the simplified approach is
clearly distinguishable.
Table I shows the speedup (in terms of time for solving
the equation systems) of the simplified MPDE approach
compared to a conventional time discretization at the
same accuracy. For higher frequency, the conventional
time discretization of (1) becomes more and more inefficient as a higher number of ripples has to be resolved.
The MPDE approach on the contrary resolves the ripples
with the Galerkin approach so that the time discretization
resolves only the envelope. This leads to a solution time
almost independent of the frequency, i.e., approximately
1 s for the simplified and 100 s for the original approach,
see Fig. 6. The higher solution time of the original
approach is a result of the evaluation of the integrals in
each time step. Thus the speedup of the original approach
is much lower compared to the simplified approach.
V. C ONCLUSION
The MPDE approach is applied to a nonlinear lowfrequency example with pulsed excitation. The solution is
obtained by a Galerkin approach and time discretization.
To evaluate the nonlinearity the ripple components due to
the pulsed excitation are neglected and only the envelope
is used. The accuracy of the proposed method rises with
increasing excitation frequency and the method offers
a considerable speedup compared to conventional time
discretization with the same accuracy.
ACKNOWLEDGEMENT
This work is supported by the ‘Excellence Initiative’ of
German Federal and State Governments and the Graduate
School CE at TU Darmstadt and in part by the Walloon Region of Belgium (WBGreen FEDO, grant RW1217703).

output voltage vC (V)

solution time (s)

4

80
60
40
original MPDE, Np = 4
simplified MPDE, Np = 4
reference solution

20
0
0

2

4

6

8

10

time (ms)
Fig. 7. Comparison of solutions for fs = 1000 Hz.
TABLE I
S PEEDUP OF MPDE APPROACH (Np = 4) COMPARED TO
CONVENTIONAL TIME DISCRETIZATION FOR DIFFERENT
FREQUENCIES .

fs (kHz)
10
50
100

approx. speedup
60
400
1000

approx. error
8 · 10−4
3 · 10−5
7 · 10−6

R EFERENCES
[1] S. Schöps, H. De Gersem, and A. Bartel, “A cosimulation
framework for multirate time-integration of field/circuit coupled
problems,” IEEE Trans. Magn., vol. 46, no. 8, pp. 3233–3236,
Jul. 2010.
[2] H. G. Brachtendorf, G. Welsch, R. Laur, and A. Bunse-Gerstner,
“Numerical steady state analysis of electronic circuits driven by
multi-tone signals,” Electr. Eng., vol. 79, no. 2, pp. 103–112, 1996.
[3] J. Roychowdhury, “Analyzing circuits with widely separated time
scales using numerical PDE methods,” IEEE Trans. Circ. Syst.
Fund. Theor. Appl., vol. 48, no. 5, pp. 578–594, May 2001.
[4] T. Mei, J. Roychowdhury, T. Coffey, S. Hutchinson, and D. Day,
“Robust, stable time-domain methods for solving MPDEs of
fast/slow systems,” IEEE Trans. Comput. Aided. Des. Integrated
Circ. Syst., vol. 24, no. 2, pp. 226–239, Feb. 2005.
[5] K. Bittner and H. G. Brachtendorf, “Adaptive multi-rate wavelet
method for circuit simulation,” Radioengineering, vol. 23, no. 1,
Apr. 2014.
[6] J. Gyselinck, C. Martis, and R. V. Sabariego, “Using dedicated
time-domain basis functions for the simulation of pulse-widthmodulation controlled devices – application to the steady-state
regime of a buck converter,” in Electromotion 2013, Cluj-Napoca,
Romania, Oct. 2013.
[7] A. Pels, J. Gyselinck, R. V. Sabariego, and S. Schöps, “Multirate partial differential equations for the efficient simulation of
low-frequency problems with pulsed excitations,” 2017, arXiv
1707.01947.
[8] H. A. Haus and J. R. Melcher, Electromagnetic Fields and Energy.
Prentice-Hall, 1989.
[9] J. W. Eaton, D. Bateman, S. Hauberg, and R. Wehbring,
The GNU Octave 4.0 Reference Manual 1/2: Free Your
Numbers. Samurai Media Limited, Oct. 2015. [Online].
Available: http://www.gnu.org/software/octave/doc/interpreter
[10] E. Hairer and G. Wanner, “Stiff differential equations solved by
radau methods,” J. Comput. Appl. Math., vol. 111, no. 1–2, pp.
93–111, 1999.
[11] T. Treichl and J. Corno, ODEpkg – A package for solving
ordinary differential equations and more., 0th ed., GNU Octave,
2015. [Online]. Available: https://octave.sourceforge.io/odepkg/

