Non-Iterative Simulation Methods for Virtual Analog Modelling

The simulation of nonlinear components is central to virtual analog simulation. In audio effects, circuits often include devices such as diodes and transistors, mostly operating in a strongly nonlinear regime. Mathematical models are in the form of systems of nonlinear ordinary differential equations (ODEs), and traditional integrators, such as the trapezoid and midpoint methods, can be employed as solvers. These methods are fully implicit and require the solution of a nonlinear algebraic system at each time step, introducing further complications regarding the existence and uniqueness of the solution, as well as the choice of halting conditions for the iterative root finder. On the other hand, fast explicit methods such as Forward Euler, are prone to unstable behaviour at standard audio sample rates. For these reasons, in this work, a family of linearly-implicit schemes is presented. These schemes take the form of a perturbation expansion, making the construction of higher-order schemes possible. Compared with classic implicit designs, the proposed methods have the advantage of efficiency, since the update is computed in a single iteration. Furthermore, the existence and uniqueness of the update are proven by simple inspection of the update matrix. Compared to classic explicit designs, the proposed schemes display stable behaviour at standard audio sample rates. In the case of a single scalar ODE, sufficient conditions for numerical stability can be derived, imposing constraints on the choice of the sampling rate. Several theoretical results are provided, as well as numerical examples for typical stiff equations used in virtual analog modelling.


Non-Iterative Simulation Methods for Virtual Analog Modelling Michele Ducceschi and Stefan Bilbao , Senior Member, IEEE
Abstract-The simulation of nonlinear components is central to virtual analog simulation.In audio effects, circuits often include devices such as diodes and transistors, mostly operating in a strongly nonlinear regime.Mathematical models are in the form of systems of nonlinear ordinary differential equations (ODEs), and traditional integrators, such as the trapezoid and midpoint methods, can be employed as solvers.These methods are fully implicit and require the solution of a nonlinear algebraic system at each time step, introducing further complications regarding the existence and uniqueness of the solution, as well as the choice of halting conditions for the iterative root finder.On the other hand, fast explicit methods such as Forward Euler, are prone to unstable behaviour at standard audio sample rates.For these reasons, in this work, a family of linearly-implicit schemes is presented.These schemes take the form of a perturbation expansion, making the construction of higherorder schemes possible.Compared with classic implicit designs, the proposed methods have the advantage of efficiency, since the update is computed in a single iteration.Furthermore, the existence and uniqueness of the update are proven by simple inspection of the update matrix.Compared to classic explicit designs, the proposed schemes display stable behaviour at standard audio sample rates.In the case of a single scalar ODE, sufficient conditions for numerical stability can be derived, imposing constraints on the choice of the sampling rate.Several theoretical results are provided, as well as numerical examples for typical stiff equations used in virtual analog modelling.
Index Terms-Virtual analog modeling, finite difference schemes, ordinary differential equations.
In an SSM, the system is represented directly in terms of physical "Kirchhoff" quantities, such as voltages and currents, and ultimately as a first order system of nonlinear ordinary differential equations (ODEs), or more generally differentialalgebraic equations (DAEs).Discretisation may be performed using a number of well-known numerical integration methods, often of implicit character, such as the trapezoid and midpoint methods [9].A PHS representation extends this direct representation in terms of currents and voltages, but using a Hamiltonian structure.In this framework, energetically passive behaviour may be translated directly to the numerical setting when an implicit discrete gradient method is employed for time discretisation [10].A numerically-stable simulation algorithm results from this approach.WDFs, on the other hand, are a scattering formulation employing wave variables, rather than currents and voltages directly.In the linear case, WDFs have the property of unconditional stability, through a passivity property encoded into power-preserving elements and scattering junctions (and thus are similar to PHS in this respect).Until recent years, WDF-based VA methods were limited to systems with a single one-port nonlinearity, though this limitation has since been overcome [11], [12].The use of the trapezoid rule (equivalent to a bilinear transformation in the frequency domain in the case of a linear system) is assumed in almost all WDF models, though the use of general linear multistep methods has been proposed [13].
Though all such frameworks have proven successful in handling many virtual analog systems, they rely mostly on fullyimplicit numerical methods, generally requiring iterative root finders such as Newton-Raphson [14].An exception is the case of circuits including a single scalar nonlinear element, which may be resolved explicitly in the wave domain by placing the nonlinearity at the "root" of the WDF network [6], [15].For typical VA systems, the iterative search for zeros is a computational bottleneck and furthermore is serial in nature, though in some cases the number of Newton-Raphson iterations may be kept low through a judicious selection of the free parameters in a nonlinear WDF [16].Beyond this basic computational cost, numerous theoretical and computational difficulties emerge: iterative methods typically operate at a variable cost at each time step (depending on the number of iterations required in order to achieve a predefined error); the user is faced with new choices regarding appropriate tolerance thresholds and the maximum number of iterations; and existence and uniqueness of numerically-computed solutions is generally not ensured [17].
This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/Though standard explicit methods avoiding iterative solvers are available, their behaviour is prone to numerical instability in the case of the strong nonlinearities that occur in VA modeling, and may require high oversampling factors to yield reasonable results [18].
This context motivates the current study: it will be shown that the typical form of virtual analog nonlinearities lends itself naturally to linearly implicit discretisations, in that at each time step, the unknown update vector is expressed as the solution of a linear system of equations [19], [20].These new schemes, obtained in the form of a perturbation expansion of increasing order of accuracy, are thus non-iterative.These schemes may be viewed as an application of modified equation techniques, employed typically in the context of linear partial differential equations [21], [22], [23], [24].Compared to fully-implicit methods, the proposed schemes have the advantage of efficiency, and of constant operation count per time step.The results presented here expnd considerably on previous work by the same authors [25].Stability results are available in the case of a VA model described by a single ODE (such as the diode clipper), and partial results are available in the case of a vector nonlinearity.A general framework to obtain the expansion coefficients for higher-accurate schemes is presented, and the theoretical results are validated by numerous numerical examples.
The article is organised as follows.Section II presents the test case of a scalar nonlinear equation, and standard discretisations of explicit (forward Euler and Runge-Kutta 4) and implicit (trapezoid and midpoint) type.The ability of these schemes to handle stiff problems is tested in the case of the diode clipper.Section III introduces the proposed non-iterative schemes, in the form of a general perturbation expansion.Stability conditions are given in the scalar case for functions that are commonly used to model nonlinear audio components.It is shown that, unlike fully explicit designs, the proposed schemes remain stable even under large input amplitudes, as demonstrated in a number of numerical experiments.Section IV generalises the scalar case to the vector case, and typical virtual analog systems (the ring modulator, and a nonlinear filter) are framed in this context.Section V assesses the properties of the numerical designs in a number of experiments.Partial results regarding the stability of non-iterative schemes in the vector case are given in Appendix.

II. BASIC NUMERICAL METHODS
Consider the basic scalar test problem, defined as x = x(t) is the unknown depending on time t ≥ 0, and initialized with x(0) = x 0 .u = u(t) is a known driving function, and f = f (x) is a nonlinear function.Equation (1), when generalized to the vector case, serves as the basis for many virtual analog models written in a state space form [1], [2], [26].The zero-input case is useful for analysis purposes: In this article, we assume that so that the nonlinearity is sector-bounded to [0, ∞], and thus corresponds to a passive loss mechanism [27].In the zero-input case (2), it immediately follows that d dt and the solution magnitude decreases monotonically, so that There are many families of numerical methods that can be used in order to solve (1) or its vector generalisations approximately; these differ widely in terms of computational cost and their stability properties.In this article, we will only consider numerical methods operating with a fixed time step T , in s, where F s = 1/T is the sample rate, as is standard in audio applications.In this case, the unknown x(t) will be approximated by a time series x n , for integer n ≥ 0, where t = nT ; similarly, u n is an input sequence (possibly drawn from samples of a continuous waveform u(t) at times t = nT ).All such schemes will be assumed initialised with x 0 = x 0 .

A. Explicit Methods
As a first step, define the forward difference operator δ + , through its application to a general time series z n , by The most basic explicit numerical method for solving (1) is the Forward Euler (FE) method, which may be written as where f n f (x n ).x n+1 may be computed directly from previously computed values of x and the input signal u.FE is first order accurate in the time step T .The many generalizations include the linear multistep and Runge-Kutta families of numerical methods (and RK4 in particular), which trade increased computational cost for increased order of accuracy [9].All such methods must observe a stability condition, usually derived in the linearized case, and framed in terms of an upper bound on the time step T .

B. Implicit Methods
First define the averaging operator μ + , through its application to a time series z n , by Two important and widely used methods may be defined, with reference to the model problem (1), by: Trapezoid Rule : Both are implicit schemes, requiring the solution of a nonlinear algebraic equation in the unknown value x n+1 at every time step.
For the midpoint rule, one must solve (10) and for the trapezoid rule, In either case, an iterative method (such as Newton Raphson [14]) will be necessary in general.Beyond greatly increased computational cost, new issues emerge, including choices of stopping criteria, initialisation, variable computational load (generally undesirable for real time applications) as well as existence/uniqueness of solutions.When f is linear, both the midpoint and trapezoid rules coalesce into the equivalent of the application of the bilinear transform to approximate the time derivative in the frequency domain (and thus are equivalent to wave digital filtering [28] in this case).Both the midpoint rule and trapezoid rule are formally second order accurate.
Both methods are unconditionally stable, allowing for very robust behaviour suitable for the modeling of strong nonlinearities that occur in virtual analog circuits.In the zero-input case, for the midpoint rule, it is true that, for a passive nonlinearity, and thus For the trapezoid rule, one may arrive at the energy balance where and again, the solution may be bounded in terms of the initial conditions, but behaviour is no longer necessarily monotonic.

C. Example: Diode Clipper
As a practical example of the behaviour of these methods in virtual analog modeling, consider the case of the diode clipper, as described in [29].Here, f (x) and u(t) are defined by Following Yeh [18], we use parameters R = 2.2 kΩ, C = 10 nF, I s = 2.52 nA, and V T = 45.3 mV.As a comparison, illustrating some of the features of the basic methods just described, consider operation at an oversampled rate of 4 × 48 kHz, using a sinusoidal input voltage signal v(t) = v 0 sin(2πf 0 t).In Fig. 1, numerical results are shown for FE, RK4, and the midpoint and trapezoid rules, for a peak input voltage of v 0 = 1.3 V, and a frequency of f 0 = 1 kHz.Even for this relatively low input voltage, severe numerical distortion is visible in the case of the FE and RK4 methods; beyond v 0 = 1.3 V, the explicit methods are unstable.This can be alleviated by increasing the sample rate, but, as pointed out by Yeh, upwards of 30× oversampling may be required for good  behaviour at higher input voltages (up to 4.5 V).In contrast, both the midpoint and trapezoid rule produce good results at an audio rate, but now at increased computational cost per sample, due to the need for an iterative root-finding method (here Newton-Raphson).Using a relative error stopping threshold (in this case, such that the relative difference 1 − x k+1 /x k of the k th iteration falls below 10 −15 ), the average number of iterations per sample for both the midpoint rule and trapezoid rule are as shown in Fig. 2. For this case of a low input voltage, it is approximately 4 for both the midpoint and trapezoid rules, but this number grows with the amplitude of the input signal and the driving frequency, or if the sample rate is reduced.

III. A FAMILY OF NON-ITERATIVE SCHEMES
In Section II, we saw the key features of two families of integrators.Explicit methods are relatively cheap, computationally, but have poor stability properties, requiring operation at an oversampled rate.In contrast, implicit methods are much better behaved numerically, but incur additional computational expense per sample due to the need for iterative solvers.
In this section, a new family of methods is proposed, following from preliminary work presented in [25], which allows good numerical behaviour without the need for iterative techniques.These methods exploit the properties of typical nonlinearities that occur in virtual analog modeling, and are thus less general than the methods described in Section II.A key new requirement is that f (x) satisfy the property: More simply, f (x) must approach 0 smoothly in the limit of small x.This condition is reasonable for virtually all nonlinearities that appear in the virtual analog modeling literature.The sector-boundedness condition ( 3) is assumed to continue to hold and implies, furthermore, that Note that these conditions require that f approaches zero for small x: some nonlinear circuits elements do not meet this requirement, such as triodes, where the grid current does not approach zero as the anode voltage becomes small [30].

A. Zero-Input Case
Consider first the zero-input case for the model problem, as in (2).Under condition (17), this may be rewritten as An explicit family of schemes may be written, using the operator definitions from ( 6) and ( 8), as Here, the discrete-time sequence g n is defined as The factor σ (P ) (x n ) will take the form of a general perturbation expansion with scale T , the time step: Here, ζ (0) = 1, and the functions ζ (p) (z), p = 1, . . ., P , can be set to achieve desired formal orders of accuracy in the scheme above (up to order P + 1), and settings will be given in the next subsection; because the leading term in σ (P ) (x n ) is of order 1, scheme ( 20) is formally consistent with the model system (19).The construction of higher-order, one-step schemes of this kind falls within the larger context of modified equation methods, typically employed for the numerical integration of partial differential equations [21], [22].To the knowledge of the authors, the construction of higher-order schemes for model problem (1), under conditions ( 17), ( 18) was presented for the first time in a previous publication by the authors [25], and the results therein are extended here.Scheme (20) allows explicit updating, when rewritten as: In this form, a basic stability condition may be arrived at: when κ n > 0, or when we have, immediately, and the numerical solution is monotonically decreasing.Condition ( 23) may hold over a restricted range of values for the time step T , and thus stability is conditional in this case.
Under the stronger condition that we have unconditional stability-the scheme will produce bounded solutions regardless of the choice of time step T .The conditions above are sufficient, but not necessary for stability.That is, there may be cases where these conditions are violated, and the scheme still behaves in a stable manner.Such conditions may be arrived at via a numerical study of the function σ (P ) (z), as shown in the examples of Section III-C.

B. Order of Accuracy
The general order of accuracy of scheme (20), despite being a one-step explicit method, may be set arbitrarily, though new order-dependent conditions on the nonlinearity f intervene.
To analyse this, we would like to show that the time series x n generated by scheme (20) are samples of a continuous function x(t) satisfying a modified equation that coincides with the model problem (2) up to a term in T P +1 : The operators δ + and μ + , when applied to x, and expanded in Taylor series about the point t = nT may be written as Employing the form of σ (P ) from ( 21), and the expansions above gives the following form for scheme (20): ) Arranging this in a series in terms of ascending powers of T yields the form where G (0) = H and G (p) , 1 ≤ p ≤ P , have the form: At this point, one may arrive at settings such that Up to p = 3, this yields the settings: where f , f and f refer to the first to third derivatives of the function f in terms of its argument, evaluated at x. (Note that by consistency, the condition (31) for G (0) holds by default.)  (1, ζ (2) AND ζ (3)   Under these conditions, and from ( 29) and (31), which implies (26) or accuracy to (P + 1)th order.

C. Examples
In this section, examples of nonlinear functions f (x) that occur in virtual analog modelling will be examined with regard to the family of schemes (20).All satisfy the properties (3) and (17).The examples of the cubic, tanh, sinh and exponential nonlinearities, all parameterized by a constant a > 0 are given in Table I, alongside expressions for ζ (1) , ζ (2) and ζ (3) .The results in this section are derived by studying numerically the behaviour of the function σ (P ) (z), as pointed out at the end of Section III-A.
The cubic and sinh nonlinearities share most features.For second order accuracy, σ (1) ≥ 1 > 0, so we have unconditional stability in this case.We cannot obtain such a bound for the third order accurate scheme, as σ (2) is negatively unbounded.But in either case, we have σ (3) ≥ 1 > 0, and thus unconditional stability for the fourth order accurate scheme.
The hyperbolic tangent nonlinearity appears frequently as a soft-clipping mechanism in virtual analog, including in the Moog 4-pole ladder filter [31], [32], [33].Here, the situation is more complex, as the functions ζ (p) generally change sign.Conditional stability conditions, in terms of an upper bound on the time step T may be arrived at: For the exponential nonlinearity, conditions for accuracy to second and fourth order are: A stability condition for the third order accurate scheme is not available in this case.
As examples of the behaviour of the integrators described above, consider the very basic case of the four nonlinearities described here, all initialized with x(0) = 1, and with a = 1, and run for a fixed duration of 1 s.In Fig. 3, the relative error between the numerical and exact solution at the end of the simulation interval is shown, as a function of sample rate, and for orders of accuracy (P + 1) up to P = 3.In all cases, error decreases as Fig. 3. Relative error between exact and numerically-computed solutions, for the four nonlinearities described in this section, and using schemes of order of accuracy as indicated.
expected, as the (P + 1)th power of the sample rate (visible as a slope of −(P + 1) in this log/log plot).

D. Sources and the Diode Clipper
Returning now to the diode clipper, presented in Section II-C, consider now scheme (20) including a source term: where f and u are as defined in (16), and g = f/x.Consider operation at a high input voltage, using scheme (36) to second order accuracy at a moderately oversampled rate of 4 × 48 kHz, and using a sinusoidal input with peak voltage amplitude v 0 and frequency f 0 .For low-to medium-range frequencies, the behaviour of the scheme is essentially equivalent to that produced by either the midpoint or trapezoid methods, as illustrated in Fig. 4, at left; recall that neither FE nor RK4 is numerically stable under these conditions, even at a 4× oversampled rate.Under extreme conditions, spurious oscillations appear in both the midpoint method and the non-iterative method relative to the trapezoid rule, but there is no explosive instability.See Fig. 4, right column.The operation count and run time of the non-iterative scheme are now a small fraction of the cost of running an iterative method (indeed, the operation count for the non-iterative method is approximately equivalent to one iteration of Newton Raphson).As can be seen from Fig. 2, at f 0 = 5 kHz, Fig. 4. Output voltage for the diode clipper, using a sinusoidal input of peak voltage v 0 = 4.5 V, and at a moderate frequency f 0 (1 kHz, left column) and a high frequency (5 kHz, right column).Results are shown using the midpoint method (top row), the trapezoid rule (middle row) and a second-order accurate non-iterative scheme (36), all running at 4 × 48 kHz.and at a sample rate of 4 × 48 kHz, both the midpoint and trapezoid rules require between 5 and 6 iterations per sample.

IV. VECTOR SYSTEMS
The extension to the vector case is direct, and allows noniterative emulation of a variety of VA systems.In continuous time, and after scaling, system (1) is generalized to: Here, x = x(t) ∈ R M is an M -element state vector.The nonlinearity is expressed through the vector f ∈ R M : where B ∈ R M ×M and F ∈ R M ×N are constant coefficient matrices and B s = (B + B T )/2 is assumed positive-semidefinite.q = q(η) ∈ R N is a vector whose components q k , k = 1, . . ., N are nonlinear functions of the scalar arguments η k , i.e. q(η) = [q 1 (η 1 ), . .., q N (η N )] .It is assumed that each q k (η k ), k = 1 . . ., N, satisfies the sector-boundedness property (3) and the limit condition (17) with respect to its argument η k .Furthermore, η ∈ R N is defined as Here, u(t) ∈ R M and c(t) ∈ R N are time-dependent sources.
For ease of notation, in (38) a linear term has been extracted within f ; the vector q contains the residual nonlinear effects.System ( 37) is a nonlinear state-space model, obtained from the general form found in e.g.[2], under restrictions on the forms of B and η.System (37) possesses a natural Lyapunov function x T x/2 in the zero-input case, leading to the following dissipative behaviour: At this point, we turn our attention to a pair of examples: the ring modulator, and the Korg35 filter, both of great practical interest in VA modeling.

A. Ring Modulator
From [34], a model for the ring modulator can be written as in (37).with M = 5 and N = 4.Given the matrices one has Here, C, C p are capacitances, L is an inductance, and R m , R a , R i are resistances; all are constant, and assumed positive.See Section V-A for numerical values for these constants.The state vector is See Fig. 5.The source vectors u(t) and c(t) are defined in terms of scalar modulator and carrier source signals u m (t) and u c (t) respectively by The nonlinear functions are given as where here, I s , V T are the saturation current and thermal voltage, respectively.The nonlinear functions (43) clearly individually satisfy the sector-boundedness (3) and limit (17) conditions.It is remarked that B 0 (and thus B), is not symmetric, but has a positive semidefinite symmetric part.

B. Korg35 Filter
An analysis of the Korg35 filter is given in [35].In this model, one has M = 2, N = 1 and, under a linear transformation, the model can be written in form (37), with and c = 0. Here, ω is the cutoff frequency of the filter, and α ≥ 0 is a resonance parameter, to be specified shortly.u(t) is the input source voltage.In this system, the state x is a two-element vector whose components are where v k (t) is the voltage across the k th capacitor in the physical circuit, and V T is the thermal voltage.A circuit diagram is available in [35].The nonlinear function is given in this case by and it is straightforward to show that the sector-boundedness and limit conditions are satisfied.Here, W (z) is the Lambert W function of z.See also Fig. 6.By inspection of the form of B in (44), it is immediate to verify that its symmetric part is diagonal, with eigenvalues 0, ω(2 − α), and is thus positive-semidefinite under the condition that While this condition certainly suffices for absolute stability (i.e., monotonic decay of the Lyapunov function (40)), it is however not necessary, as it depends only on the properties of the linearized system, through B. A more refined analysis, taking into account the combination of B and the nonlinearity is tractable in this case, and a sufficient and necessary condition for monotonic decay is obtained when Since β is small, the sufficient condition (47) is in fact not very restrictive.This system is able to self-oscillate, since q(η) η is bounded for large |η|.When (i.e. for 2 < α < 8), the filter is partly dissipative, and limit cycles are observed [35].When α ≥ 8, the system exhibits unbounded growth.

V. NUMERICAL METHODS AND EXAMPLES
System (37) may be integrated using vector extensions of the trapezoid, midpoint, and the non-iterative schemes.In the following, the zero input case is shown.The numerical schemes are given as follows: Midpoint: Non-iterative: Here, in analogy with (20), and suppressing the time index n, where D = diag(q i /η i ).For the non-iterative scheme, a perturbation expansion analogous to (21) may be obtained in the vector case to a desired order P + 1.The first two terms are given here explicitly as σ (1) where I denotes the k × k identity matrix, and Λ = diag(dq i /dη i ).
When sources are included, these can be approximated using the averaging operator, as (53) The use of the operator μ + in front of the sources serves two purposes: the local truncation error is maintained to O(T 2 ) and the input signal is low-passed.

A. Example: Ring Modulator
Consider the ring modulator, as described in Section IV-A, with circuit parameters I s = 40.63nA; V t = 56.3mV; C = C p = 10 nF; L = 0.8 H; R a = 600 Ω; R i = 50 Ω; R m = 80 Ω.Consider sinusoidal inputs u m (t) = 1.2 sin(800πt) and u c (t) = u c,max sin(3780πt), for different peak carrier amplitudes u c,max .This is a very stringent test of numerical method performance in VA.Even for relatively small carrier amplitudes u c,max , standard explicit methods fail entirely.(For, example, if u c,max = 0.5 V, then FE is unstable at sample rates below 3.3 MHz.) In Fig. 7, numerical results are shown under the choices u c,max = 0.5 V and u c,max = 2 V.In both cases, errors are smallest in the case of the trapezoid rule, and largest for the non-iterative method.However, at the relative tolerance threshold of 10 −10 employed for Newton Raphson, at u c,max = 2 V, the trapezoid and midpoint rules require 6.2 and 26.8 iterations per sample, respectively.At higher carrier input amplitudes (u c,max = 3 V), however, the non-iterative method begins to exhibit spurious oscillatory behaviour (see Fig. 8), and at even higher amplitudes will become unstable, signalling an as-yet undetermined stability condition.Such non-iterative methods allow operation at rates far below those of standard explicit methods such as FE.

B. Example: Korg35
In Figs. 9 and 10, the Korg35 filter is tested, again under large input amplitudes.Here, the output of the schemes is relatively well behaved, since the nonlinearity is quite soft compared to that of the diode clipper or ring modulator.Fig. 11 shows the average number of iterations needed for convergence of the iterative schemes up to a tolerance threshold of 10 −10 .While the trapezoid method requires significant fewer iterations than the midpoint method, around 4 such iterations are required for higher input frequencies, compared with a single iteration of the non-iterative scheme.One problem here concerns the evaluation of the Lambert W function.While in these examples Matlab's lambertw was used, this generally requires some kind of iterative solver, such as Fritsch's iteration, see e.g.[36], [37].A separate iterative method is required in order to evaluate the Lambert function, in addition to any iterative solver required in order to perform an update, and remains as a fixed cost, and is independent of the choice of numerical method (i.e.iterative or non-iterative).

VI. CONCLUSION
Analog circuits often include nonlinear components to process audio signals.Virtual models, particularly when geared toward real-time applications, must ensure the stability of the resulting time stepping schemes.It has been shown here that the stability of typical virtual analog state-space models can be analysed by means of Lyapunov functions.Lyapunov stability carries over to the discrete case, when discretisation is performed using the trapezoid and the midpoint methods.While stable, such methods are fully implicit, and require the solution of a nonlinear algebraic system of equations at each time step.Further complications arise in the choice of tolerance thresholds, and halting conditions.For these reasons, in this work a family of higher-order linearly implicit schemes was derived in the form of a perturbation expansion.Since the update is expressed as a linear system of equations, no iterative methods are required.Compared with standard implicit schemes, these schemes have the advantage of efficiency, and avoid entirely the problem of tolerance thresholds.Stability of the schemes in the absence of sources can be inferred using Lyapunov arguments, and sufficient stability conditions can be given in the case of models described by a single scalar ODE.The extensions to vector systems is possible, and numerical examples show that the non-iterative schemes can run very efficiently, compared with fully implicit methods.Stability results are more difficult to obtain in the vector case; some partial results appear in Appendix.Stability results for systems including a source, as well as for vector systems, are worthy of future investigations.

APPENDIX NUMERICAL STABILITY
The non-iterative schemes presented here outperform standard expicit methods (such as FE or RK4).That is, they are able to produce output comparable to that of a standard iterative solver (e.g.trapezoid), at high input levels where standard explicit methods, such as FE or RK4 fail, or require very high oversampling factors in order to operate.In the case of a scalar system (with M = N = 1), stability results have been shown in Section III.In the case of vector systems, the situation is less clear, but the desirable attribute of robust performance with a low oversampling factor, and at high input amplitudes remains.Some partial results are presented here.
To this end, consider system (37)- (38), with an arbitrary state size M , but with a single nonlinearity, so that N = 1.In this case, F = v is a constant M × 1 vector, but B remains a constant M × M matrix, with positive definite symmetric part.Using the non-iterative scheme (50c), with the second-order choice of σ (1)  Noting that the matrix left-multiplying δ + x is of the form of a rank 1 perturbation of the identity, the Sherman-Morrison identity [38] may be used to arrive at: (56) This leads, finally, to where and z = μ + x.
In order for the natural choice of 1 2 x T x to be a Lyapunov function for this scheme, we require Q ≥ 0 for all z.Conditions on λ and d ensuring this are not immediately forthcoming.It is easily seen that Q ≥ 0 under the edge case of a linear system (when λ = d = 0), in which case In the case where B = 0, the expression for Q reduces to in which case, given that d ≥ 0 from the sector-boundedness condition, if λ ≥ d, then again Q ≥ 0. This is indeed the case for the Korg35, as illustrated in Fig. 6.If λ ≤ d, we arrive at conditional stability in the same manner as in Section III.Instability has not been observed using this scheme in the case of the Korg35, which has N = 1 and is of the form above.This could imply that the natural choice of Lyapunov function is not suitable in this case, but this remains unproven.Given the good behaviour of the scheme, it may be that a different form is more appropriate.
The more general case of a system with multiple nonlinearities (i.e.N > 1) is more difficult, and in this case the scheme presented here does exhibit numerical instability (as in the case of the ring modulator), though at much higher driving amplitudes than standard explicit methods such as FE or RK4.

Fig. 1 .
Fig. 1.Simulation results for the diode clipper, at a sample rate of 4 × 48 kHz, and using v 0 = 1.3 V, f 0 = 1 kHz, for simulation methods as indicated.

Fig. 2 .
Fig. 2. Diode clipper: Average number of iterations/sample as a function of driving frequency f 0 , for the midpoint rule (left) and trapezoid rule (right) at different input voltages v 0 as indicated, and at a sample rate of 4 × 48 kHz.
40)The last inequality follows immediately from positive semidefiniteness of B s , and sector-boundedness of q.

Fig. 7 .
Fig. 7. Ring modulator simulation results.with parameters and input signals as given in Section V-A, and for u c,max = 0.5 V (left column) and 2 V (right column).Reference solutions (top row) are computed using the trapezoid rule at a sample rate f s = 200 • 48 kHz.The errors are computed as the difference between the reference solution and the output of the three methods using a sample rate f s = 4 • 48 kHz.Newton-Raphson is run with a tolerance threshold τ = 10 −10 , and the maximum number of iterations is limited to 100.The output is y = v 2 .

Fig. 11 .
Fig. 11.Korg35.Average Number of iterations for the trapezoid and midpoint rules, using an input square wave of input frequency f 0 , and peak amplitude u 0 .The sample rate used for all simulations is 2 • 48 kHz.Here, Newton-Raphson is run with a tolerance threshold τ = 10 −10 .

TABLE I NONLINEAR
FUNCTIONS f (x) AND ASSOCIATED FUNCTIONS ζ