Frequency Response and Gap Tuning for Nonlinear Electrical Oscillator Networks

We study nonlinear electrical oscillator networks, the smallest example of which consists of a voltage-dependent capacitor, an inductor, and a resistor driven by a pure tone source. By allowing the network topology to be that of any connected graph, such circuits generalize spatially discrete nonlinear transmission lines/lattices that have proven useful in high-frequency analog devices. For such networks, we develop two algorithms to compute the steady-state response when a subset of nodes are driven at the same fixed frequency. The algorithms we devise are orders of magnitude more accurate and efficient than stepping towards the steady-state using a standard numerical integrator. We seek to enhance a given network's nonlinear behavior by altering the eigenvalues of the graph Laplacian, i.e., the resonances of the linearized system. We develop a Newton-type method that solves for the network inductances such that the graph Laplacian achieves a desired set of eigenvalues; this method enables one to move the eigenvalues while keeping the network topology fixed. Running numerical experiments using three different random graph models, we show that shrinking the gap between the graph Laplacian's first two eigenvalues dramatically improves a network's ability to (i) transfer energy to higher harmonics, and (ii) generate large-amplitude signals. Our results shed light on the relationship between a network's structure, encoded by the graph Laplacian, and its function, defined in this case by the presence of strongly nonlinear effects in the frequency response.


Introduction
Networks of nonlinear electrical oscillators have found recent application in several microwave frequency analog devices [1][2][3][4][5][6]. The fundamental unit in these networks is a nonlinear oscillator wired as in Figure 1; this oscillator consists of one inductor, one voltage-dependent capacitor, one source, and one sink (a resistor). While many nonlinear oscillatory circuits have been studied for their chaotic behavior, the particular oscillator in Figure 1 does not exhibit sensitive dependence on initial conditions in the regime of operation that we consider. Instead, assuming the source is of the form A cos (vtzw), the oscillator reaches a steady-state consisting of a sum of harmonics with fundamental frequency v [7].
When networks of these oscillators have been studied, the network topology has either been a one-dimensional linear chain, in which case the circuit is called a nonlinear transmission line [8][9][10][11][12][13]-see Figure 2, or a two-dimensional rectangular lattice [14][15][16][17][18]-see Figure 3. Even if each individual block in the chain/ lattice is weakly nonlinear, the overall circuit can exhibit strongly nonlinear behavior. It is this property that is exploited for microwave device applications, enabling low-frequency, lowpower inputs to be transformed into high-frequency, high-power outputs.
The first objective of this work is to develop numerical algorithms to compute the frequency response of a nonlinear electrical network with topology given by an arbitrary connected graph. Here we are motivated by the successful application of computational techniques in the design of the high-frequency analog devices referenced above. As we show, to compute steadystate solutions with comparable accuracy, both the perturbative and iterative algorithms developed in this paper require orders of magnitude less computational time than standard numerical integration. While the perturbative algorithm generalizes derivations given in prior work [7,18], the iterative algorithm has not been previously applied to nonlinear electrical networks. Both new algorithms show exponential convergence in the number of iterations, and for a test problem on a network with N~400 nodes, less than 20 iterations are required to achieve machine precision errors.
The second objective of this work is to relate structural properties of the network to the dynamics of the nonlinear oscillator system. The derivation of the perturbative algorithm indicates that nonlinearity in the electrical network manifests itself through energy transfer from the fundamental forcing frequency to higher harmonics. This helps us understand why properties such as amplitude boosting [7,18] and frequency upconversion [1], observed in nonlinear electrical networks with regular lattice topologies, can be expected when the topology is that of a random, disordered network. Additionally, we observe that an inductanceweighted graph Laplacian matrix features prominently in both algorithms for computing the steady-state solution. This graph Laplacian matrix encodes the structure of the network, and its eigenvalues are the squares of the resonant frequencies for the undamped, linear version of the circuit. Driving the damped, linearized circuit at one of these resonances results in large amplitude outputs. It is reasonable to hypothesize that the locations of these resonances play a large role in the dynamics of the nonlinear network.
This motivates the following question: how do the eigenvalues of the graph Laplacian influence the nonlinear network's properties of frequency upconversion and amplitude boosting? While it is possible to alter the spectrum of the graph Laplacian by changing the node-edge relationships in the graph, we can also change its spectrum by keeping the topology fixed and manipulating the network's inductances. We formulate and solve the inverse problem of finding the inductances such that the graph Laplacian achieves a prescribed spectrum. The solution proceeds via a Newton-type algorithm that takes the desired spectrum as input and iteratively alters the inductances until a convergence criterion is met.
For three types of random graphs, we find that the Newton-type method effectively finds circuit inductances that close the gap between the first two eigenvalues of the graph Laplacian. We conduct a series of numerical experiments to examine the effect of closing this eigenvalue gap on a given circuit's ability (i) to transfer energy from the fundamental driving frequency to higher harmonics, and (ii) to generate high-amplitude output signals.
The results indicate that the two metrics (i-ii) can be improved dramatically by closing the gap between the graph Laplacian's first two eigenvalues. Table 1 shows results we obtained for graphs with N~175 nodes. Though this a small portion of the results we describe later, this table already illustrates the effect of gap tuning on network performance. Note that each pre and post circuit have the same graph topology, differing only in their edge inductances.
Note that we have made available open-source Python implementations of all algorithms described in this work. The Python code, together with R code used for plotting, has been posted on a public repository. This enables the reader to reproduce all results in this paper. Instructions on how to download this code is given below.

Connections to Other Systems
We can make several connections between the problem studied in this paper and other problems of interest: N Random elastic networks. Using a mechanical analogy between inductorscapacitors and massessprings, the nonlinear electronic network can be transformed into a mathematically equivalent network of masses and anharmonic springs [19,Appendix I]. Such random elastic networks have been of recent interest as models of amorphous solids [20][21][22]. For such networks, quartic spring potential energies have been considered [23]. Nonlinear random elastic networks have also been used to model molecular machines; in this context, tuning the gap between the first two eigenvalues of the linearized system enables the construction of networks with properties similar to those of real proteins [24]. Despite this activity, algorithms for computing and manipulating the frequency response of nonlinear elastic networks have not been developed. Our work addresses this issue directly.
N Nonlinear electromagnetic media. The circuit we analyze, for particular values of the circuit parameters, arises naturally as a finite volume discretization of Maxwell's equations for TE/TM modes in a nonlinear medium [25,26]. The arbitrary connected graph topology of the circuit corresponds to a finite volume discretization on an arbitrary unstructured mesh. The algorithms developed here can be used to compute and optimize the frequency response of nonlinear electromagnetic media.
N Coupled phase oscillator networks. There has been intense interest in nonlinear phase oscillator networks, primarily due to the ability of such networks to model biophysical systems featuring synchronization. Though synchronization is not of primary interest in our system, we may still draw parallels. The effect of network topology on the properties of coupled phase oscillators has been studied extensively [27][28][29][30]. Manipulating eigenvalues of the Laplacian matrix enables one to enhance a network's synchronization properties [31]. More recently, several authors have developed algorithms for optimizing the synchronization of phase oscillator networks [32][33][34][35][36][37]. The questions considered in this subset of the coupled phase oscillator literature are related to the issues addressed in the present work.

Problem Formulation
Let H(N,e) be a connected, simple graph with N nodes and e edges. Each edge corresponds to an inductor that physically connects two nodes. Each node corresponds to a capacitor and resistor, wired in parallel, that physically connect the node to a common ground. Let f ƒN be the number of nodes that are  Frequency Response of Electrical Networks PLOS ONE | www.plosone.org driven by prescribed sources. Since the voltage at the prescribed source is known, we do not model it using a node. The connection between the source and the node that it drives is modeled by a half-edge, also known as a dangling edge since one end is connected to a driven node and the other end does not connect to any node. We let H(N,e,f ) denote the graph together with the f half-edges.
The capacitance and conductance (inverse resistance) at node j are C j and G j , respectively. We let V j (t) denote the voltage from node j to ground at time t. The inductance of edge k is L k , while the current through edge k at time t is I k (t). The exact dimensions for each component of H, along with the currents and voltages, are tabulated in Table 2.
In order to write down Kirchhoff's laws, we must choose an orientation of the edges. The orientation of an edge records the direction of positive current flow through the edge. If we solve the problem with opposite orientations, the only difference we will notice is that the currents will pick up a factor of {1. Consequently, the orientation we choose does not affect the solution in any material way. In what follows, we will choose a random orientation of the edges.
In Figure 4 we show an example graph corresponding to H (6,9,2). The edges are oriented randomly. The inputs are connected at nodes 1 and 6 through two inductors. These input nodes correspond to half-edges in H. On the right we view node 3 in detail. Each of the two edges connected to this node correspond  to an inductor. A capacitor with capacitance C 3 and a resistor with conductance G 3 connect node 3 to ground.
To arrange Kirchhoff's laws compactly, we use the N|(ezf ) incidence matrix of H(N,e,f ), denoted by B. Let j be an edge connecting the nodes i 0 and i. If j is oriented such that positive current starts at node i 0 and flows to node i, we write j~(i 0 ,i). If j is a half-edge attached to node i, we write j~(1,i), leaving the first slot empty and orienting the half-edge so it always points toward the forced node. The entries of the incidence matrix B are i fj~(i,i 0 ) for some node i 0 0 otherwise: This paper will only consider single frequency time-harmonic forcing of the form ae ivt zae {ivt where a[ N . Let P be an N|(ezf ) matrix with entries P i,j~1 if node i is connected to an input edge j and 0 otherwise. Using the projection matrix P we define the forcing Using the notation summarized in Table 2, Kirchhoff's laws for the nonlinear circuit on the graph H(N,e,f ) can now be written compactly as Here L(dI=dt), C(dV =dt), and GV are examples of component-wise multiplication of vectors. For a,b[ m , we define c~ab[ m by c j~aj b j for 1ƒjƒm. Note that in this case, we can also write c~diag(a)b. Here diag(a) is the m|m matrix that contains the vector a along its diagonal (½diag(a) ii~a i ) and is zero elsewhere.
The formulation (2-3) generalizes previous formulations [25,38] in which the capacitors were constant and the systems considered were linear.
By differentiating (3) and inserting it into (2), we obtain a second-order system for the voltages: Here Note that D is the weighted Laplacian for the network with edge weights given by reciprocal inductance.
We assume that the capacitance at node i depends on the voltage at node i: where C 0 [R is a constant. Note that this choice of capacitance function means that (4) features a quadratic nonlinearity. We can then formulate the frequency response problem for the nonlinear electrical network: given the amplitude vector a and frequency v for the forcing function (1), determine the steady-state solution V (t) of (4).

Perturbative Algorithm
We first solve the frequency response problem using a perturbative expansion in powers of ". We use dots to denote differentiation with respect to time. Substituting the capacitance function (7) in (4) and rearranging, we obtain  Each node corresponds to a voltage-dependent capacitor to ground, wired in parallel with a resistor to ground, as depicted in the zoomed-in schematic for node 3.
Each edge corresponds to an inductor. Each half-edge connects one prescribed voltage source to one given node. In this paper, all methods that are developed are valid for connected graphs with at least one half-edge. Note that the circuits in Figures 1-3  We expand Inserting (9) into (8), we obtain equations for each order of ". At zeroth order, we obtain For k §1, the k-th order equation is We now solve (10)(11). Let us introduce the Fourier transform in time,ŷ with inverse Fourier transform Note that with these definitions, _ y y _ y y(a)~iaŷ y(a): This implies that the Fourier transforms of both sides of (10-11) can be summarized by writing and By (5) and (1), we see that where d is the Dirac delta. Then the k~0 branch of (14) yieldŝ Using the inverse Fourier transform, we have where ''c.c.'' stands for the complex conjugate of the previous terms. Here we have used the property that L({a)~L(a).
Once we have computed V 0 (t), we can insert it into (16) to compute w 1 (t). We will find that w 1 (t) is a linear combination of e {2ivt , e 0ivt , and e 2ivt . Using this fact in the k~1 branch of (14), we can solve forV V 1 (a) and then apply the inverse Fourier transform to compute V 1 (t). We will find that V 1 (t) contains the same modes as w 1 (t).
The above shows how we get the perturbative solution algorithm started. Now let us move to the more general case where we seek V k (t) for any k §1. Assume that we have already computed V j (t) for 0ƒjƒk{1, and that the solution takes the following form: In words, V 2m contains odd modes 1,3, . . . ,2mz1, and V 2mz1 contains even modes 0,2, . . . ,2mz2. Here we assume that 0ƒ2mv2mz1ƒk{1, and that the a i,j [ N coefficients are known.
In order to solve for V k (t), we use the k §1 branch of (14), which requires us to compute (16). We have two cases, when k is odd and when k is even. In both cases, it is a simple (if tedious) algebraic exercise to show that w k (t) yields: In both cases, it is clear that using (14) to solve forV V k (a) results in a sum of Dirac delta's. Applying the inverse Fourier transform yields V k (t), which will be a sum of Fourier modes. One can check that V k (t) will have precisely the form (19a) or (19b) depending on whether k is even or odd, respectively.
The algorithm is then clear. Starting with (19), we apply component-wise multiplication to particular pairs of the a i,j vectors in order to compute the coefficients of the Fourier modes of w k (t) defined in (16). Next, we combine the step of solving forV V k (a) using the k §1 branch of (14) together with the step of computing the inverse Fourier transform. After component-wise multiplication of the Fourier coefficients of w k by {(a 2 =2)C 0 , we multiply each coefficient on the left by ½L(a) {1 with a set to match the frequency of the corresponding Fourier mode. Dividing these coefficients by 2p yields the Fourier coefficients of V k (t), as desired.
While we have presented the algorithm in an intuitive way, the statements made above can be made rigorous, and a convergence theory for the perturbative expansion (9) can be established. This is the subject of ongoing work.
There are a few brief remarks to make about the algorithm presented above: N As described above, we consider only those networks that contain resistance at all nodes, i.e., G i w0 for all nodes i. Such an assumption is not only physically realistic; it also guarantees that for all a[R, the matrix L(a) is invertible. The invertibility for the a~0 case is a consequence of Corollary 1 proved below.
N In this work, we are interested in the weakly nonlinear regime where "EaE is sufficiently small such that the perturbative method converges. As the nondimensional constant "EaE is increased beyond the breakdown point of the perturbative method, direct numerical solutions of the equations of motion reveal subharmonic oscillations, and eventually, chaotic oscillations.
N The fact that the Fourier transform yields the steady-state solution has been explained in our earlier work [7]. By fixing an arbitrary set of initial conditions and using the Laplace transform to derive the full solution, one can show that after the decay of transients, the part of the solution that remains is precisely what we obtain using the Fourier transform. This also explains why it was not necessary for us to specify initial conditions for (4) in our derivation above-the initial conditions only influence the decaying transient part of the solution.

Iterative Algorithm
The perturbative method developed above shows us that the solution V (t) is a sum of harmonics where the fundamental frequency is given by the input frequency v. This implies that the steady-state solution V (t) is periodic with period T~2p=v. This observation leads us to ask whether it is possible to directly solve for the Fourier coefficients of V (t) without first expanding in powers of E. In this section, we develop a fixed point iteration scheme that accomplishes this task.
First, we integrate both sides of (8) from t~0 to t~T to derive We show below that as long as the network contains at least one half-edge, D is invertible. Hence (20) implies This means there is no zero/DC mode present in V (t), motivating the Fourier series expansion a k e ikvt zc:c: ð22Þ In order to compute the solution, we truncate at k~M, leading to an approximation V &V : a k e ikvt zc:c: ð23Þ Using orthogonality we derive Using the T-periodicity of V and integration by parts, we have To simplify notation, we combine (1) and (5) and write V in~w e ivt zc:c: where Now let d m,n denote the Kronecker delta function which equals 1 if m~n, and 0 otherwise. We multiply both sides of (8) by e {ikvt , integrate from t~0 to t~T, and finally divide by T to obtain where L was defined in (15) and Because the form of the nonlinearity is simple, we can insert (23) into (26) and derive with the understanding that a 0~0 , a {j~aj for jw0, and a j~0 for jjjwM. We insert (27) into (25) and obtain We convert this into an iterative scheme in a natural way. Let a (j) k denote the j-th iterate, and assume that a k terms appearing on the left-hand side are at iteration jz1, while those appearing on the right-hand side are at iteration j. Let A (j) denote the N|M complex matrix whose k-th column is a (j) k . Then the scheme is where the k-th column of the matrix F M (A (j) ) is Here we assume 1ƒkƒM, which is also why we have deleted the second Kronecker delta from the right-hand side.
Starting at A (0) , we iterate forward using (28), stopping the computation when EF(A (j) ){A (j) E is below a specified tolerance. Note that in our implementation of F, we precompute and store the LU factorization for the M matrices fL(kv)g M k~1 , since this part of the computation of the right-hand side of (29) does not change from one iteration to the next.
Again, we have derived the algorithm but have not proven its convergence. Instead, we will demonstrate empirically that the algorithm converges using several numerical tests.

Inverse Problem
In this section, we consider the inverse problem of finding a set of inductances such that D, the Laplacian defined by (6), achieves a desired spectrum. Before describing an algorithm to solve this inverse problem, we review basic spectral properties of D.
Lemma 1. Assume all inductances are positive. Then D as defined in (6) is symmetric positive semidefinite, and all its eigenvalues must be nonnegative.
Proof. Let ½diag(L) {1=2 be the diagonal matrix whose Let fl i g N i~1 denote the spectrum of D, with eigenvalues arranged in nondecreasing order: l 1 ƒl 2 ƒ Á Á Á ƒl N . The above argument shows that l 1 §0. We can be more precise about this: if there are no half-edges, then l 1~0 , while the presence of at least one half-edge causes l 1 w0.
Lemma 2. Let H~H(N,e) be a connected graph with N nodes, e edges, and zero half-edges. For a particular orientation of the graph, let B denote the signed incidence matrix. Then rank(B)~N{1.
Proof. Let r be any integer from 1 to N{1. Consider any subset S of r vertices of the graph. Take the sum of the rows of the incidence matrix corresponding to the elements of S. This sum cannot be zero; if it were, there would be no path connecting S to the complement S c and the resulting graph would not be connected. Hence the sum of these rows must contain a nonzero entry. As the same would be true if we considered linear combinations of the rows corresponding to S, we conclude that any subset of at most N{1 rows must be linearly independent. At the same time, if we take the sum of all the rows we get a zero row, because each column contains precisely one z1 and one {1. Lemma 3. Let H 0~H (N,e,f ) be a connected graph with N nodes, e edges, and f w0 half-edges. For a particular orientation of the graph, let B 0 denote the signed incidence matrix. Then rank(B 0 )~N.
Proof. Without loss of generality, we can assume that the N|(ezf ) incidence matrix B 0 is organized such that the first e columns correspond to full edges, while columns ez1, . . . ,ezf correspond to half-edges. Now choose any j such that 1ƒjƒf , and examine column ezj of B. Let k be the unique row in which this column contains +1. Since row k of B is the only row that contains an entry in column ezj, row k is linearly independent from the other N{1 rows of B. By Lemma 2, the submatrix of B consisting of all rows other than row k has rank N{1. Including row k increases the rank by one, yielding a rank N matrix. Proof. Combine Lemmas 1 and 4. We now describe an algorithm that quantifies how we must change the vector of inductances L in order to make D have a desired set of eigenvalues. In what follows, we assume we work with a system that satisfies the hypotheses of Corollary 1.
For nƒN, let l Ã~( l Ã 1 , . . . ,l Ã n ) T denote a vector of desired eigenvalues satisfying We treat the vector of inductances L as a variable, and let l(L) denote the sorted vector of eigenvalues of the graph Laplacian D defined in (6). Since D is symmetric, it possesses an orthonormal basis of eigenvectors. We assume that v j (L) is the normalized eigenvector corresponding to l j (L). Now let F : R ezf ?R n be the function We now apply a version of Newton's method to find a zero of this function. To use Newton's method we will need to compute the Jacobian J(F (L)). Let primes denote differentiation with respect to L k . To form the Jacobian we need to find l 0 j :~L LL k l j (L): We proceed by implicit differentiation, starting from the eigenvector equation Differentiating both sides with respect to L k , and omitting the dependence on L, we obtain Multiplying (34) on the left by v T j and using (35) together with where Using l 0 j we can compute the entries of the Jacobian matrix and the corresponding Newton's method with pseudoinverse becomes where { denotes the Moore-Penrose pseudoinverse. Using (34) as shown might produce inductances such that the ratio of the largest to smallest inductance is too large. In order to avoid these large variations, we constrain E {1 ƒL i ƒE. We incorporate these constraints using an active set approach, replacing F by the function G i : R ezf zm ?R nzm , where i denotes the iteration number and m denotes the number of constraints violated by L (i) . Let Q + denote the functions For every constraint p violated from below, we set . For every constraint q violated from above, we set G i q (L)~Q z (L q {E). Since the Q + functions are continuously differentiable, it is easy to compute the Jacobian J(G i (L)) and then apply the algorithm Algorithm (36) can be used to alter all the eigenvalues of the system if n~N and l Ã [R N . Alternatively, one can set n~2 and only request the two smallest eigenvalues to be changed to l Ã 1 and l Ã 2 , respectively.
In the next section we show that altering the lowest eigenvalue l 1 is enough to cause higher energy transfer to the higher modes. To show, we will use (39) to change l 1 to some desired value, keeping l 2 constant. We note that since we do not constrain l 3 , Á Á Á ,l N , they can change as a result of altering L, but l 2 ƒl j for j[f3, Á Á Á ,Ng will be maintained.  For all applications of this inverse problem algorithm described in the next section, we use (36) with the initial conditions L (0)~½ 1,1, . . . ,1 T and the constraint violation parameter E~10 3 .

Gap Tuning: Methodology
How does the steady-state voltage in the nonlinear circuit change as a function of the gap between the first two eigenvalues of the graph Laplacian D? In this section, we address this question by combining the perturbative/iterative algorithms with the inverse problem algorithm. We describe numerical experiments designed to test the effect of closing the graph Laplacian's first eigenvalue gap on the circuit's ability to (a) transfer more energy to higher harmonics, and (b) generate higher-amplitude output signals.
We conduct our numerical experiments on three types of random graphs, all generated using the NetworkX package [39]: N Barabási-Albert (BA), a preferential attachment model with one parameter, m, the number of edges to draw between each new node and existing nodes [40]. We set m~3 in our experiments. N Watts-Strogatz (WS), a small world model with two parameters, k, the number of nearest neighbor nodes to which each node is initially connected, and p, the probability of rewiring each edge [41]. In our experiments, we set k~5 and p~0:3. When we produce realizations of any of these graphs, we accept only those graphs that are connected. Suppose we have used one of these three models to generate a connected, random graph with N nodes. To make this a concrete circuit problem, we set C 0 i~1 for all nodes i, and L j~1 for all edges j. We fix the nonlinearity parameter E~0:5. We select N=10 nodes uniformly at random, and attach half-edges to these nodes with inductance L j~1 . For each node i, we set the conductance G i~0 :15 for the BA and WS graphs, and G i~0 :5 for the ER graphs. This selection will be explained in more detail below.
With these parameters set, we have enough information to compute the graph Laplacian D defined by (6). As we did before, let l 1 , . . . ,l N denote the eigenvalues of D sorted in increasing order. We set the forcing frequency v~ffi ffiffiffi ffi l 2 p . Since this value is a resonant frequency of the linear, undamped system we expect it lies close to a resonance for the nonlinear, damped system. The type of forcing we consider is A sin vt, a special case of (1) with a~A=(2i).
With this setup, we use both the perturbative method and the iterative method to compute the steady-state solution V (t). For the perturbative method, we solve up to order 9, and for the iterative method, we solve using 20 modes. This means that the iterative scheme captures ten modes-11v through 20v-that are not captured by the perturbative scheme. We compare the two solutions as a check for whether the number of modes we have considered is sufficient. In all tests, we find that there is no significant difference between the solutions, implying that the first 10 harmonics-v through 10v-are sufficient to resolve the solution.
Since E~0:5, the capacitance model (7) is valid only for V i v2. For all computed solutions, we check that the maximum voltage across all nodes in one cycle satisfies this constraint.
One quantity of interest in our simulations is the amount of energy in the higher harmonics. Let Y be an N|M complex matrix such that the j-th column of Y contains the Fourier coefficients of the zjv mode over all N nodes. Here j goes from 1 to M, the maximum mode to which the solution is computed. We then define the fraction of energy in modes z2v and higher, averaged over all nodes. We also compute the maximum magnitude voltage produced anywhere in the circuit during one period. For both k and V, the subscript ''pre'' denotes that these quantities have been computed before we change L to manipulate the eigenvalues of D.
Having computed k pre , we now study how this fraction changes when we reduce the gap between the first two eigenvalues of D. For a fixed d[(0,1), we set l Ã 1~d l 2 and l Ã 2~l 2 , and then apply the inverse problem algorithm.
Using (36), we solve for a vector of inductances L Ã such that the first two eigenvalues of D are given by l Ã 1 and l Ã 2 . When we iterate forward using (36), if we find that jjG i (L)jj 2 §10 {12 after 200 iterations, we generate a new random graph and restart the experiment.
We recompute the graph Laplacian D using the new vector inductances L Ã , and again apply the perturbative and iterative algorithms to solve for the steady-state solution V (t). Using this solution, we compute the energy in the higher harmonics using the right-hand side of (37), now labeling this average fraction as k post . We also compute the right-hand side of (38) and label this quantity as V post .
Let us now describe how we choose the particular values of the conductance G i and the eigenvalue fraction d. In Table 3, we tabulate l 2 {l 1 , the gap between the second and first eigenvalue for each of the three types of random graphs described above. The eigenvalue gaps we present are averaged over 200 simulations for each of four graph sizes: N[f25,75,125,175g.
We observe that the eigenvalue gaps for the BA and WS graphs do not change appreciably as a function of N, while for ER graphs, the gaps grow quickly as a function of N. When we solve for the steady-state voltages on these three types of graphs, we also notice a difference. For ER graphs, the maximum voltage grows quickly as a function of N, while for BA and WS graphs, the same phenomenon does not occur. To counteract the large growth of maximum voltages for large graph sizes, we set the conductance G i to the larger value of 0:5 for ER graphs, causing more energy to dissipate through resistors. For BA and WS graphs, we set G i to 0:15.

Comparison of Steady-State Algorithms
In this section, we compare steady-state solutions computed by numerical integration against the solutions computed using the perturbative and iterative methods derived earlier.
For the tests described in this section, the domain is a 20|20 square lattice with N~400 nodes. Nodes along the left and bottom boundaries of the lattice are driven by input forcing. The input provided is 0:03 sin (vt) with v~1. For the capacitance model (7), we set C 0~1 and E~0:5. For each edge j, we set L j~1 . The conductance G i is set to 0:01 at all points except for the topright corner, where it is set to 1:0.
To compare the results of the perturbative and iterative methods against the numerical integrator, we will need to obtain the steady-state solution using the numerical integrator. To do this, we start at t~0 and numerically integrate the first-order system (2-3) forward in time for 1500 cycles. The ODE solver uses the Dormand-Prince (dopri5) method with relative and absolute tolerances equal to 10 {10 and 10 {12 , respectively. For the parameters given above, this number of cycles is sufficient so that, from one cycle to the next, the change in the solution is on the order of the relative tolerance of the numerical integrator. Hence we take the solution over the last cycle to be the steady-state solution.
As a preliminary check, we directly compare the three steadystate solutions. We treat the solution obtained from numerical time integration as a reference solution z ref (t). Let z (i) (t) denote either the perturbative or iterative solution after i iterations-for the perturbative method, the iteration count is defined as the largest mode number present in the solution. Let T be the period of the steady-state solution, and for an integer tw0, consider the equispaced discretization of the interval ½0,T given by ft k~k T=tg k~t k~1 . For each iteration i, we evaluate both the the perturbative/iterative and reference solution on this equispaced grid with t~64 points, and we compute the error In Figure 5 we have plotted log 10 E(i) as a function of the iteration i. While both methods initially tend towards the reference solution, we see from Figure 5 that the error does not drop below 10 {9 . In the following subsections, we provide evidence that the reference solution is less accurate than the solutions computed using the perturbative/iterative methods. This explains why the error in Figure 5 does not go to zero, i.e., why the perturbative/ iterative methods will not converge to the solution produced by time integration.
Our first tests concern the Fourier coefficients of the computed solutions. In what follows, we use A to denote the vector of Fourier series coefficients associated with a steady-state solution computed using any of the three methods discussed above.
Fixed Point Error. Suppose that V (t) is an exact T-periodic steady-state solution of (4). If we were to expand this solution in a Fourier series as in (22), the resulting (infinite) coefficient vector A would satisfy A~F ? k (A) for all k, with F as in (29).  Table 5. Comparison of the three solutions' preservation of the energy balance (44).
Scheme In both the perturbative and iterative methods, what we seek is a finite-mode truncation of the exact solution. For the iterative method we fix M~20 so that the highest mode has frequency 20v. In the perturbative method we solve, we solve up to order 19, which implies that the highest mode in the solution has frequency 20v.
Combining the ideas of the last two paragraphs, it is natural to use as an error metric for the M-mode truncation of the exact solution. In Table 4, we record (40) for solutions computed using the perturbative, iterative, and numerical integration methods. Note that the iterative and perturbative methods directly provide us with the Fourier coefficients necessary for this calculation. We compute the Fourier coefficients of the numerical integrator's solution using the FFT. Table 4 shows that the perturbative and iterative solutions are about five orders of magnitude closer to being fixed points of F M k than the solution obtained from numerical integration.
For the perturbative and iterative methods, let us examine how the fixed point error (40) decreases as a function of iteration count. In Figure 6, we plot log E M (A (i) ) versus the iteration number i. Here A (i) is the vector of Fourier coefficients for the solution computed after only i iterations. The plot shows that, for both the perturbative and iterative methods, approximately 10 iterations are required to match the fixed point error of the solution computed using time integration. The error of this latter solution, taken from Table 4, is represented on Figure 6 by a horizontal black line. Figure 6 also shows that the perturbative and iterative methods converge exponentially in the number of iterations. From iteration 1 until iteration 16, fitting lines of best fit to the perturbative and iterative error curves results in slopes of {2:1407 and {2:2326, respectively. For both methods, this can be approximated by  Table 4. For the perturbative/iterative methods, Time I records the time required to compute 10 iterations, resulting in a fixed point error comparable to that of the numerical method-see the crossing of the curves with the black horizontal line in Figure 6. For the perturbative/iterative methods, we also record under Time II the time required to achieve the O(10 {16 ) errors as in Table 4. All times are averaged over 10 runs. doi:10.1371/journal.pone.0078009.t006 Figure 8. Barabá si-Albert random graph results. From left to right, we present results for Barabási-Albert random graphs with N~25, 75, 125, and 175 nodes. For each graph, we use algorithm (36) to modify the inductances L such that the ratio of the smallest to the second smallest eigenvalue is d. We use pre and post to denote, respectively, the graphs before and after algorithm (36) is applied. By shrinking the gap between the first two eigenvalues, the energy transferred to higher harmonics (37) can be increased from approximately k pre &0:5% to k post &5% (for all graph sizes), and the maximum voltage (38) can be increased from V pre v0:05 volts to V post [½0:3,0:5 volts (depending on the graph size). We also note that for larger graphs, choosing d~1 (i.e., no gap between the first two eigenvalues) does not yield optimal behavior. doi:10.1371/journal.pone.0078009.g008 For each graph, we use algorithm (36) to modify the inductances L such that the ratio of the smallest to the second smallest eigenvalue is d. We use pre and post to denote, respectively, the graphs before and after algorithm (36) is applied. By shrinking the gap between the first two eigenvalues, the energy transferred to higher harmonics (37) can be increased from k pre &1% to k post §8:75% (for all graph sizes), and the maximum voltage (38) can be increased from V pre &0:05 volts to V post &0:5 volts (for all graph sizes). The values of k post for Watts-Strogatz graphs are about twice as large as the values of k post for Barabási-Albert graphs in Figure 8. doi:10.1371/journal.pone.0078009.g009 Figure 10. Erdö s-Ré nyi random graph results. From left to right, we present results for Erdö s-Rényi random graphs with N~25, 75, 125, and 175 nodes. For each graph, we use algorithm (36) to modify the inductances L such that the ratio of the smallest to the second smallest eigenvalue is d. We use pre and post to denote, respectively, the graphs before and after algorithm (36) is applied. By shrinking the gap between the first two eigenvalues, the energy transferred to higher harmonics (37) can be increased to k post [½1,8% (depending on the graph size), and the maximum voltage (38) can be increased to V post [½0:1,0:8 (depending on the graph size). The results for Erdö s-Rényi graphs are much more strongly dependent on the number of nodes N than the results shown in Figures 8 or 9. Note that the peak voltages for the N~175 graphs with forcing amplitude 0:01 are the largest voltages for any graphs considered in this paper. We can increase the peak voltages for smaller graphs by choosing a smaller value of the conductance than G i~0 :5 (for all nodes i), the value used to compute the results in this figure. doi:10.1371/journal.pone.0078009.g010 writing E M (A (i) )*e {2i . After 16 iterations, the error has approached machine epsilon, and both curves level off before reaching the final values shown in Table 4. Decay Rate. Suppose we write the first-order system (2)(3) in the form _ z z~F(z,t), with z~(I,V ) T . Then on the open set D~fI[R ezf ,V [R N ,t[R j jV i jvE {1 , 1ƒiƒNg, the vector field F : D?R ezf zN is C ? , i.e., F is j times continuously differentiable for any integer j §1. By the standard existence/uniqueness theorem for ordinary differential equations, it follows that wherever the solution z(t)~(I(t),V (t)) T exists, it must also be C ? in t.
The above observation enables us to test the decay of the Fourier coefficients of all three solutions. For if the steady-state solution V (t) is C ? in t, then the Fourier series coefficients of V must satisfy the following decay property: To examine the decay of the Fourier coefficients for the three computed solutions, we plot log Ea k E versus k in Figure 7. For the perturbative and iterative solutions, the curves on the plot coincide and are nearly linear with slope {2:8004. This implies that Ea k E*e {2:8k , which is sufficient to satisfy the theoretical decay rate given by (41).
The Fourier coefficients obtained from the numerical integrator's solution, on the other hand, do not decay at all beyond mode 10. This violates the theoretical decay rate (41) even for '~1.
Energy Conservation. Next we test the energy conservation properties of the three computed solutions. We proceed to derive an energy balance equation. Because our capacitors are voltagedependent, the charge Q and voltage V are related via dQ~C(V )dV , which implies dQ dt~C (V ) dV dt : Using this in (3) together with (2), we obtain Let M(t) be the total energy stored in the magnetic fields of all inductors at time t. Then the first term on the left-hand side of (42). Let E(t) be the total energy stored in the electric fields of all capacitors at time t. Then the second term on the left hand side of (42). Hence (42) reads: If the system has reached a T-periodic steady state, then I(t), V (t), M(t), and E(t) will all be T-periodic. Therefore, integrating (43) in t from t~0 to t~T, we find that the left-hand side vanishes. The remaining terms yield the following energy balance equation: The left-hand side is the energy pumped into the system over one cycle, while the right-hand side denotes the energy dissipated through resistors, again over one cycle. Table 5 shows the absolute difference between the left-and right-hand sides of (44), computed for each of the three methods. We find that for the perturbative and iterative methods' energy balance errors are below machine epsilon. The numerical integrator yields an error approximately five orders of magnitude larger than that of the two other methods.
Computational Time. The results presented thus far indicate that whether we measure error using the fixed point error (40) or the violation of the energy balance (44), the solution obtained via numerical integration has errors that are approximately five orders of magnitude larger than that of the perturbative/iterative methods. The actual values of the errors committed by the numerical integrator in Tables 4 and 5, as well as the final error values for the curves in Figure 5, are close to the numerical integration relative and absolute tolerances of 10 {10 and 10 {12 , respectively. We hypothesize that, if computational time were not an issue, we could run the numerical integrator with smaller tolerances and obtain steady-state solutions that more closely match, in the same error metrics described above, the perturbative and iterative solutions.
As we now proceed to show, computational time is a major issue for the time integration method. In Table 6, we record the time required to compute steady-state solutions using the three methods. We see from the Time I column that to achieve the error of &2|10 {11 in Table 4, the numerical integrator requires over 483 seconds. We know from Figure 6 that the perturbative/ iterative methods require 10 iterations to achieve approximately the same error as the time integrator; the remaining entries in the Time I column show that both the perturbative and iterative algorithms compute such a solution hundreds of times faster than the time integrator.
The Time II column in Table 6 records how long it takes the perturbative/iterative algorithms to achieve the errors recorded in Table 4. Observe that even if we run the perturbative/iterative algorithms all the way to full convergence, they are much faster than time integration. In this case, the time integrator is 542 (respectively, 66) times slower than the iterative (respectively, perturbative) algorithm.
Note that the perturbative and iterative algorithms were implemented in Python using the Numpy/Scipy packages. The dopri5 implementation used for numerical time integration is the implementation provided by the scipy.integrate.ode module. All times reported are average times across 10 runs on the same machine.

Gap Tuning
For each N[f25,75,125,175g, and for each of 10 values d chosen in an equispaced fashion from the intervals given above, we compute 100 runs of the complete procedure described abovesee Gap Tuning: Methodology. For each such run, we compute pre/post values of k and V for three values of the input forcing amplitude, which we take to be the same at all input nodes k: A k [f0:001,0:005,0:01g. These results for k and V, averaged over the 100 runs, are plotted in Figures 8, 9, and 10. Figure 8 shows the results for Barabási-Albert (BA) graphs. By shrinking the gap between the first two eigenvalues, the percentage of energy transferred to higher harmonics (37) can be increased by approximately one order of magnitude, for all graph sizes, while the maximum magnitude voltage (38) can be increased by a factor of 6 to 20, depending on the graph size. Note that for larger graphs, choosing d~1, i.e., forcing the first two eigenvalues to coincide, does not yield optimal energy transfer to higher harmonics.
The results in Figure 9 for Watts-Strogatz (WS) graphs are similar to those for BA graphs. We again find that by shrinking the gap between the first two eigenvalues, the energy transferred to higher harmonics (37) can be increased. However, the values of k post for Watts-Strogatz graphs are about twice as large as the values of k post for Barabási-Albert graphs in Figure 8. For all graph sizes, tuning the eigenvalue gap can increase the percentage of energy transferred to higher harmonics by a factor of up to 8, while the maximum magnitude voltage can be increased by approximately one order of magnitude.
In Figure 10, we present the results for Erdös-Rényi graphs. The results again support the finding that by shrinking the gap between the first two eigenvalues, the circuit can transfer more energy to higher harmonics and boost the peak magnitude of output signals. Specifically, we see that the energy transferred to higher harmonics (37) can be increased to k post [½1,8%, and the maximum voltage (38) can be increased to V post [½0:1,0:8.
The results for Erdös-Rényi graphs are much more strongly dependent on the number of nodes N than the results shown in Figures 8 or 9. Note that the peak voltages for the N~175 graphs with forcing amplitude 0:01 are the largest voltages for any graphs considered in this paper. We can increase the peak voltages for smaller graphs by choosing a smaller value of the conductance than G i~0 :5 (for all nodes i), the value used to compute the results in Figure 10.
For all three types of graphs, both pre and post values of k and V increase as a function of the input forcing amplitude.

Code
All code necessary to reproduce the above results have been posted as a public repository on GitHub, accessible at the following URL: https://github.com/GarnetVaz/Nonlinear-electrical-oscillators We use Python together with the numpy, scipy, matplotlib, and networkx modules for all numerical computing. The code that generates Figures 8, 9, and 10 is set to utilize 10 processors using the open-source multiprocessing module. For plotting, we use R together with the ggplot2, plyr, and reshape packages. All languages, packages and modules used are open source.
Assuming all packages and modules have been correctly installed, one can reproduce all results by running the Python codes numerical_comparison.py and graphmulti.py. The latter code may require several hours to run. The Python codes will generate figures using R; the R codes we provide need not be run independently.
Further details on how to run the codes, including the specific versions of required packages and modules, are given in the README.md file at the URL given above.
The code that we provide can easily be modified to run simulations not described here. For example, one can compare the perturbative/iterative algorithms against numerical integration using graphs other than the 20|20 grid graph used above. One can also explore gap tuning results for random graphs with different parameters than the ones we have chosen.

Conclusion
For nonlinear electrical circuits on arbitrary connected graphs, we have developed two numerical methods to compute the steadystate voltage. Using both absolute metrics and relative comparisons with a solution obtained via direct numerical integration, we validated the new algorithms. The results show that for the same error tolerance, both the perturbative and iterative methods are orders of magnitude faster than the solution obtained by time stepping. Moreover, these methods are able to capture the behavior in high Fourier modes and converge to machine precision in a fixed point error metric.
In future work, we plan to apply the steady-state algorithms developed above to solve Maxwell's equations in nonlinear electromagnetic media [26]. This application makes use of the correspondence between the nonlinear electrical network and a finite volume discretization of Maxwell's equations on an unstructured mesh.
In order to enhance the nonlinearity-driven features of these circuits, we developed a Newton-like algorithm that alters the eigenvalues of a network's graph Laplacian. The algorithm leaves the topology of the network untouched, changing only the inductances, i.e., the edge weights. By applying the Newton-like algorithm to three types of random graphs, we showed that reducing the gap between the graph Laplacian's first two eigenvalues leads to enhanced nonlinear behavior. Comparing pre-and post-optimized circuits, it is evident that optimizing the eigenvalue gap significantly increases (i) energy transfer from the fundamental driving frequency to higher harmonics, and (ii) the maximum magnitude output voltage.
In both the perturbative and iterative algorithm, the only way in which the network's structure influences its frequency response is through the graph Laplacian matrix. In our experiments, we have tuned this graph Laplacian's first eigenvalue gap by holding the topology of the graph constant and altering the edge weights, i.e., inductances. What if we had instead tuned the gap by holding the edge weights constant and altering the topology of the graph? This would amount to altering the incidence matrix B instead of the inductance vector L. So long as both types of alterations result in the same graph Laplacian D, our results indicate that the nonlinear electrical network's functionality should be enhanced significantly. This, of course, leads to the question of whether it is possible to algorithmically alter the network topology to achieve a particular graph Laplacian matrix, an interesting avenue for further work.