Analog simulator of integro-differential equations with classical memristors

An analog computer makes use of continuously changeable quantities of a system, such as its electrical, mechanical, or hydraulic properties, to solve a given problem. While these devices are usually computationally more powerful than their digital counterparts, they suffer from analog noise which does not allow for error control. We will focus on analog computers based on active electrical networks comprised of resistors, capacitors, and operational amplifiers which are capable of simulating any linear ordinary differential equation. However, the class of nonlinear dynamics they can solve is limited. In this work, by adding memristors to the electrical network, we show that the analog computer can simulate a large variety of linear and nonlinear integro-differential equations by carefully choosing the conductance and the dynamics of the memristor state variable. We study the performance of these analog computers by simulating integro-differential models related to fluid dynamics, nonlinear Volterra equations for population growth, and quantum models describing non-Markovian memory effects, among others. Finally, we perform stability tests by considering imperfect analog components, obtaining robust solutions with up to 13% relative error for relevant timescales.


Introduction
Analog computers employ continuously tunable quantities of a system, such as its electrical, mechanical, or hydraulic properties, to codify and solve a given problem. Usually, the dynamics of an analog computer perfectly matches the dynamics of the simulated system. The advantage of analog computers lies on the computational power provided by real-time operation and complete parallelism, so that they require less resources than a digital counterpart for the simulation. On the other hand, the accuracy of an analog computer is limited by its computing elements and by the quality of its analog components. Despite the astonishing evolution of digital computers during the last decades, the interest in analog computers has re-emerged in recent years 1 . For example, posible applications as math co-processors in VLSI architectures 2,3 , and wave-based analog computation in metamaterials have recently been proposed [4][5][6] .
An electrical analog computer is an active network composed of electrical elements, namely, resistors, capacitors and operational amplifiers which, connected together, are capable of simulating any set of linear ordinary differential equations 7 . In this case, the solutions of these equations are encoded into the time evolution of the voltage waveform produced by the analog computer. Indeed, electrical analog devices can only work with a single independent variable when codified in time 8,9 . Nonetheless, by making use of finite-difference methods, it is also possible to solve partial differential equations 7,10 .
Resistive memory effects, in which the resistance depends on the history of charges crossing through a material, can be described by Kubo's response theory 11 , and have been observed as resistance switching phenomena in thin films 12,13 . The connection between resistance switching devices and memristive devices 14,15 was first made in 2008 by HP Labs 16 . Afterwards, various other physical implementations of memristors have been reported [17][18][19][20] . The memory and non-volatility properties of the memristor 21 have sparked great interest in its applications in the fabrication of digital memories 22 , performance of computational tasks [23][24][25][26][27] and recently, in the design of a perceptron 28 , which motivates the possibility of building memristor-based hardware for a physical neural network and proved their universality as approximators. While most of the interest in memris-tors lies in applications in digital computation, there have also been proposals of utilization of memristors in programmable analog components of electrical networks 29 and in the design of highly efficient operational amplifiers 30 . Regarding the solution of differential equations, memristor arrays have been employed as accelerators for digital solvers 31,32 and, recently, a digital approach based on memristor crossbar arrays have been employed to solve partial differential equations 33 with high precision. Nonetheless, a natural question is whether the characteristic properties of memristors could extend the applications of electrical analog computers.
In this work, by adding classical memristors to the network of an electrical analog computer, we show that it can simulate a large variety of linear and nonlinear integro-differential equations with interest in mathematics and physics. This is possible by appropriately choosing the conductance of the memristors and the dynamics of its state variable. We study the performance of these analog computers by simulating integro-differential models of fluid dynamics, nonlinear Volterra equations for population growth, and quantum models describing non-Markovian memory effects, among others. Finally, by considering 10% error in the analog components, we perform stability tests of the dynamics, showing the robustness of the simulations with up to 13% error for relevant timescales. It is noteworthy to mention that the operation conditions of an off-the-shelf memristor fit those of usual circuit toolboxes.

Electric Analog Computer
To simulate a linear ordinary differential equation, the analog computer only requires the following operations: (i) summation, (ii) sign inversion, (iii) integration and (iv) multiplication by a constant.
We will describe briefly the implementation of the aforementioned operations for the analog computer. First, let us consider the summation and sign inversion operations, which are performed by an adder circuit whose transfer characteristics for n inputs are given by where V out is the output voltage and K i = R f /R i . The symbolic representation and equivalent circuit of the adder is shown in Fig. 1(a). Sign inversion can be implemented by an adder with a single input with K 1 = 1. Next, integration, is meant to be the most important operation in an analog computer. It can be implemented by an integrator circuit, whose transfer characteristics for n inputs are given by where the constant IC stands for the initial condition. The symbolic representation and equivalent circuit of the integrator is shown in Fig. 1(b). For the adder and integrator circuits, we consider ideal operational amplifiers. The last linear operation, multiplication by a constant, can be implemented with a potentiometer, which is described by where α < 1 is a constant that characterizes the device. In addition, analog computers can also simulate nonlinear dynamics to some extent, by considering signal generators and function multipliers, which enable the mathematical multiplication of two signals.
We will consider a simple example to illustrate the configuration of an analog computer for the simulation of a linear ordinary differential equation. Let us consider the following second-order linear differential equation This equation is written in this form to remark thatÿ can be obtained as the sum of −ẏ and y, which are obtained from integratingÿ once and twice, respectively. The computer diagram that implements Eq. (4) is shown in Fig. 1(c). Therefore, we see that integrator circuits enable the solution of ordinary differential equations by making use of feedback connections, which make the input and output voltage in Eq. (2) to be equal. In addition, several integrators connected in series stand for subsequent integrations of the desired variable. In this way, the simulation of a linear differential equation of order n would require n integrators. This method can be extended to a system of m linear ordinary differential equations of order n, which can be written as The simulation of such a system requires nm integrators, and since the summation can be absorbed into the integrators, it requires at most nm 2 sign inverters. It is important to note that, in order to solve a differential equation with an analog computer, it is not necessary to know the voltage waveform at any point of the electrical network, provided that the circuit was properly constructed, the desired solution is obtained by measuring the voltage signal at the output of the circuit. We remark that analog computers are best suited for solving systems of ordinary differential equations, and aditionally, it is possible to use them to solve partial differential equations by making use of finite differences approximations 7 . Since analog computers are better suited to solve time-dependent problems, the non-temporal variable is discretized to obtain a system of nonlinear ordinary differential equations for the time variable.

Memristive analog computation
A voltage-controlled memristor can be described by the following current-voltage relation 14 where i(t) and v(t) are the current and voltage across the device, respectively. The quantity ω is an internal state variable specified by the physical implementation of the memristor. For ideal voltage-controlled memristors the state variable corresponds to the magnetic flux φ , which means f (v) = v. The quantity g(ω) is the so-called memductance which has units of conductance and is a function of the state variable. Since the memristor is a passive circuit element, the memductance g(ω) must always be positive. The notion of memristive behavior can be extended to a broader class of system whose characteristics resemble those of the memristor, and are called memristive systems 15 . These are defined by Where in general, ω represents a set of n state variables, and f and g are continuous functions. Now, to study the inclusion of memristive devices into an analog computer we will consider the case of an integrator circuit with a single input, where we substitute the resistor for a memristor, as is shown in Fig. 2. We will consider that the memristor has a single state variable ω. In the figure, the current across the memristor, i M , is given by Eq. (8). Then, by Kirchhoff's law at node A, we have Integrating Eq. (9) allows to write Eq. (11) aṡ which in general will be a nonlinear integro-differential equation that will be specified by g and f , which describe the memductance and the memristor state variable dynamics, respectively. It must be noticed that these functions are specified by memristor engineering. We now show the type of equations that can be solved using this circuit. Let us consider the equation that Volterra introduced for single populations in the study of growth of biological populations 34 , in which N(t) is the population size at time t, a and b are positive rate constants, and k(t) is the "hereditary" influence. Now, in order to simulate this equation, we consider a memristive system as described by Eq. (8) and Eq. (9). By assuming where we have chosen C = 1 and k(t, τ) = k 1 (t)k 2 (τ). This is restricted, however, to kernels that can be separated in this way, otherwise it is only possible to simulate single-variable kernels. This equation can be implemented with a single integrator and memristor as shown in Fig. 2.
By connecting integrator circuits in series, we can generate nonlinear integro-differential equations of higher-order. For example, if we consider n integrators in series, where only the first integrator circuit contains a memristor, then we can generate the nth-order integro-differential equation On the other hand, if each of the n integrators contains a memristor, then the nth-order integro-differential equation will involve the composition of n memductance functions g 1 , . . . , g n corresponding to each memristor. The resulting equation can be written as follows This allows for the simulation of a wide class of nonlinear integro-differential equation. It is also possible to solve first-order linear integro-differential equations with a suitable change of variables. If we consider Eq. (12) and define the memristor by Then, with u = ln(v), we havė Notice that we can avoid the sign inversion of the integrand by normalizing the voltage such that 0 ≤ v ≤ 1. In this way, it is possible to simulate several linear integro-differential equations coming from quantum models. An example of this is a Volterra equation describing non-Markovian quantum memory effects 35 , given by Since L is a linear superoperator, the integral kernel will involve linear terms in the elements of ρ. This leads to a system of integro-differential equations for the elements of ρ which can be implemented by appropriately choosing the dynamics of the internal state variable of the memristor.
Aditionally, linear integro-differential equations that appear in the model of turbulent diffusion 36 , have the forṁ Equations of this type can be implemented with a memristor defined by g(ω, v,t) = α(t)k 1 (t)ω andω = k 2 (t) α(t) 2 ln(v) 2 , such that K(t, s) = k 1 (t)k 2 (s), and where α(t) = exp t 0 p(s)ds . By using z(t) = α(t)u(t), we can configure Eq. (20) into the analog computer as the following equation, with k ′ (s) = k(s) α(s) 2 . Next, connecting the output z(t) to a signal multiplier with the signal 1 α(t) , we recover the solution u(t) of Eq. (20).

Equivalent memristors for enhancing simulations
We have seen that the integro-differential equation that can be implemented with the integrator circuit of Fig. 2 depends on the memductance, g, and internal variable dynamics, f , in Eqs. (8) and (9), respectively. However, the choice of the memductance function g is restricted by currently available memristors. A possible way for broadening the class of available memductance functions is to consider a composite memristor circuit, connected in series or in parallel 37,38 .
The behavior of composite memristor circuits is nontrivial, it involves transient phenomena which depend on the initial states of each memristor and is also affected by their polarities. Nonetheless, it has been shown that, for ideal voltagedependent memristors, a composite memristor circuit involving either in parallel or in series connection can obey the same rules as with standard resistors 37,38 . This behavior depends on the memristor internal relaxation time and on the state variable dynamics. Assuming that the timescale of memristor relaxation is much shorter than the timescale associated with the state variable dynamics, then we can neglect the former and ignore transient phenomena.

6/11
We can expect that, under the aforementioned assumptions on timescales, a composite memristor circuit involving memristive systems will obey a dynamical model similar to Eqs. (8) and (9). Then, we look to describe the composite memristor circuit as an equivalent memristor with an equivalent memductance. For a composite memristor circuit connected in series, as shown in Fig. 3 (a), the equivalent memductance is given by whereas for memristors connected in parallel, it holds with v = v 1 + v 2 , and ω = ω 1 + ω 2 . Then, it is possible to approximate a memductance function, g(ω, v,t), by considering its Taylor expansion, g(ω, v,t) = c 1 + c 2 v + 1 2! c 3 v 2 + ..., and construct a composite memristor circuit connected in parallel such that each memristor in the circuit provides one of the terms of the desired Taylor series. Then, the equivalent memductance, given by is an approximation of the desired memductance function, g(ω, v,t), provided that g( 2! c 3 v 2 and so on. In this way, the desired memductance function, g(ω, v,t), can be approximated with a controllable accuracy which depends on the amount of terms of the Taylor series implemented. In addition, when memristors are set in proximity, they may mutually interact 39,40 , affecting the dynamics of their internal variables, and consequently, their memductances. This leads to a different class of integro-differental equations that can be implemented in the analog computer by including coupled memristors into the network. Let us consider an ideal voltagecontrolled memristor system described by where we consider g 1 (φ ) = αφ + β , and α, β are constants. Then, we have We can consider again two coupling cases, memristors connected in series and in parallel. For two memristors connected in series, as shown in Fig. 4(a), we define v 12 = v 1 + v 2 , and φ 12 = φ 1 + φ 2 . By applying Kirchhoff's voltage law, together with the constitutive relations of the memristors, one obtains a set of coupled differential equations for the internal variables of the memristors 39 . In the special case of α 1 = α 2 = α, β 1 = β 2 = β , κ 1 = κ 2 = κ, the equivalent memductance, defined by i(t) = g 12 (φ 12 )v 12 , is consequently given by (see Ref. 39 ) This result may allow us to implement an additional kind of integro-differential equation, related to non-Markovian dynamics, given by the following expression: On the other hand, for coupled memristors connected in parallel, as shown in Fig. 4(b), we have i = i 1 + i 2 , and φ 12 = φ 1 + φ 2 .
In this case the total memducante of the dual coupled memristors in parallel connection is given by the expression (see Ref. 39 ) This is, the memristors coupled in parallel operate as a new memristor, and the equivalent memductance is obtained as the sum of the individual memductances. We see that, in this case, coupled memristors connected in parallel behave very similar

7/11
to the uncoupled case. The main difference is the appearance of additional parameters. We want to mention that to consider a large memristor array in our case could benefit the proposal, specially in obtaining a good approximation of any desired memductance function by constructing its Taylor expansion with an equivalent memristor array connected in parallel. The main limitation of using a large number of memristors lies in the error accumulation due to imperfections in the devices, as will be seen later. In this sense, the flexibility with large arrays increases as the device fabrication improves.

Adaptability with other classes of memdevices
Another possibility for broadening the class of integro-differential equations that can be implemented by the analog simulator is to consider other classes of memdevices, such as memcapacitive and meminductive devices 41 . Here we will consider the substitution of the capacitor in the integrator circuit of Fig. 2 for a voltage-controlled memcapacitive device which can be specified as follows: where q, v, x and C, are the charge, voltage, internal variable and memcapacitance in the memcapacitive system, respectively. Then, performing the same calculation in node A of Fig. 2 with Kirchhoff's law, as in Eq.(11), we obtain his leads to the following integro-differential equatioṅ This naturally adds a new parameter, the memcapacitance C(x, v,t), that increases the flexibility in the class of integrodifferential equations that can be implemented with the analog simulator. For instance, we can see that by considering C(x, v,t) = x 41 we can include a linear term in the voltage that will depend on the dynamics of the internal variable of the memcapacitive system. This can ease the requirements on the memductance which for some applications can become a very intricate function of the internal variable, voltage and time.

Numerical robustness
In order to provide an example of the previous procedures and to test the stability of the solutions, let us consider the following integro-differential equatioṅ This equation can be implemented by taking C = 1 and chosing g(ω, v,t) = −2 + 0.001v + e −t ω andω = e s s 1+s v in Eq. (12). In Fig. 5(a), we show the exact solution of this equation compared against the case with up to 10% error in the coefficients of Eq. (38), which stands for imperfections in the analog components of the computer. We can see that the average relative error stabilizes at ∼ 10% indicating that the solution is robust against imperfections in the analog components. As a second example we consider the following integro-differential equatioṅ This equation can be taken into the form of Eq. (21) by considering α(t) = exp(−1/16(e −2t−1 )), k 1 (t) = 1/2e −t , k 2 (s) = e −s andω = k 2 (s) ln(v) 2 . In Fig. 5(b), we show the solution of this equation compared against the case with up to 10% error in the coefficients of Eq. (39), as well. We observe that, in this case, the error increases continuously, as it is usually the case in analog computers, reaching up to 13% within the timescale considered.

Discussion
We have studied the inclusion of memristors into the network of electric analog computers. We have found that this addition enables analog computers to simulate linear and nonlinear integro-differential equations by appropriately choosing the memductance and the dynamics of the memristor state variable. This broadens the applicability of analog computers which, otherwise, could only solve systems of ordinary differential equations. Additionally, we have numerically studied the performance of these analog computers by simulating integro-differential models related to fluid dynamics, nonlinear Volterra equations for population growth, and quantum models describing non-Markovian memory effects. Moreover, we have tested and evaluated the stability of the solutions provided against imperfections of analog components by introducing up to 10% error in their corresponding coefficients, finding that the relative error reaches up to 13% for relevant timescales. From these results, we have concluded the robustness and resilience of the results of the simulation provided by the methods introduced in this work. Moreover, it is noteworthy to mention that our results have been obtained with minimal architecture variations, leaving the possibility of more intrincated arrangements, whose complexity grows in complexity extremely fast, open for future studies. The analog simulator can benefit from a large array of memristors, coupled in series or parallel, specially if we consider the Taylor expansion for any memductance function, which is directly related to the class of integro-differential equations that can be implemented by the circuit. The main limitation in using a large number of memristors is the error accumulation due to imperfect fabrication of the devices. In this sense, the flexibility of the analog simulator with the number of memristors increases as the device fabrication and design improves. As a further scope, it would be interesting to consider the proposed scenarios in the context of quantum memristors [42][43][44][45] . These devices are the quantized version of the memristors employed in this work, i.e. open non-Markovian controllable open quantum systems which are able to process quantum information and can be used for simulations and for implementation of quantum neural networks 46,47 . The intuition is that quantum entanglement could be useful in the search of quantum speed-ups for analog computers, in a similar manner as it already happens with quantum simulations.