Analog computing of partial differential equations

We show in this paper that a class of Periodic Networks of Resistances (called PNR) can be used to solve partial differential equations. More precisely, each cell of this periodic network computes locally a finite difference approximation of the partial differential equation. The PNR circuits are composed of current source, voltage-controlled voltage source and positive and negative resistances. An example of application of PNR could be the real time control of distributed MEMS.


INTRODUCTION
Actuation, sensoring and control of very large sized arrays of MEMS, says 1000×1000 MEMS or cells on the same chip, are out of reach of digital computing (microprocessor). The revival of analog computing might be the recourse for filling the gap of our computing power. Control laws are formulated in terms of ordinary or partial differential equations. The electronic analog computers used in the seventies were limited to solve ordinary differential equations. They were composed of operational amplifiers and capacitances. Unfortunatly, the use of numerous capacitances renders them unsuitable for high density integration. Today, advent of nanotechnology could make conceivable molecular or quantum computing [1].
We propose in this paper a much simpler solution to bring back to life analog computing. Our method is based on the design of electronic circuits which "computes" locally the finite difference approximation of a partial differential equation. These unit cell circuits must have the properties of being possibly assembled then periodically recopied. The whole periodic network is then composed of the numerous repeated identical cells, except at the boundaries where the cell takes into account the boundary conditions. So we can say that this Periodic Network or Resistances (PNR) "solves" or "is described" by a PDE.
The unit cell involves only a current source, VCVSs (Voltage-Controlled Voltage Sources) and positive and negative resistances. The pratical realization of a PNR is not an issue of this paper, but we can just say that the designs of negative resistances can be based on variety of techniques, for example negative input impedance [2][3] [4], negative impedance converters and inverters [5] or translinear circuits [6]. The last reference is itself electronically tunable.
We now present how to construct an electronic circuit of a linear PDE up to the fourth order. The discrete solution at each voltage node of the center of each cell composing this circuit converges towards the solution of the PDE, with its boundary conditions, for a given input signal (the right hand side of the PDE) when the number of periodic cells increases. The choice of the finite difference scheme used to represent the PDE is a mathematical problem and is not relevant here. It is clear that if this finite difference scheme is stable and accurate, its electronic implementation will have a good behaviour too.
The 14+1 partial differential operators up to order 4 are summarized in table I. It is sufficient to study 9 of them (at the left side of the table), all the other ones can be deduced by symmetry.
The following sections propose a study of eight unit circuits which "code" a finite difference (FD) approximation of these 8 partial differential operators. The section names "Molecule αβ" refers to the centered FD approximation of the derivative ∂ α+β ∂x α ∂y β with an error proportional to the step size of order 2. The calculation of the harmonic derivative (Molecule 20) is well known, all the other proposed circuits are originals. The study of each molecule is presented through an increasing difficulty and by gathering the similar solutions.
Next, we present how to mix these elementary cells to construct the circuit of a linear partial differential equation, and finally we show how to implement the boundary conditions in a circuit form.
The uppercase letter V i,j is used to denote the finite difference approximation of the function v (ih, jh). The two-dimensional grid used in the finite difference approximation is represented on Fig. 1 in the case N = 4. The points of the domain are labeled by a filled circle, those on the boundary are labeled by a filled square, and finally the point labeled by a cross are used to compute Neumann boundary conditions.
The easiest computational molecule to implement under the form of an electronic circuit is the one associated with the second derivative. The molecule is symmetrical.
We recall that a computational molecule is a graphical depiction of an approximate partial derivative formula. The previous one stands for: The elementary cell which "calculates" (1) is represented on Fig. 2 with its too adjacent cells. We call ρ 1 the resistances along the paths located by the index 1 in the scheme and linking two adjacent potentials (i.e. ρ 1 is the resistance between V i and V i−1 or V i and V i+1 ). By applying the Kirchhoff's Current Law (KCL) at node V i , rearranging some terms, the equation of the cell i can be written under the form (2). 1 If we choose resistances which scaled down in h 2 with the cell size, i.e. ρ 1 = −h 2 ρ 0 , the potential of a cell i estimates the (second) derivative of the function v(x) at a single point. The whole electronic circuit composed of N cells computes a finite difference approximation of the differential equation (3). The particular choice to take ρ 1 = −h 2 ρ 0 instead of ρ 1 = h 2 ρ 0 is very important, it allows to superpose the electrical network of this molecule to other one (all the circuit equations must have the same left hand side to be superposed).
The values of the resistances inside a cell depend only on the circuit topology and are easily expressed as a function of ρ 1 , here r 1 = r 2 = ρ 1 /2.
The computational molecule of the fourth derivative is given by (4). The molecule is symmetrical and expands over two steps at each side.  The elementary cell which "calculates" (4) is represented on Fig. 3 with its four adjacent cells. The resistances ρ 1 connect the potentials V i and V i±2 and the resistances ρ 2 connect the potentials V i and V i±1 . Applying KCL at node V i and doing some algebraic manipulations leads to (5).
If we choose ρ 1 = −h 4 ρ 0 and ρ 2 = h 4 ρ 0 /4, then the potential at node V i estimate the fourth derivative of the function v(x) at a given point (Eq. (6)).
The molecule of the first derivative is represented by (7). This molecule is not symmetrical and thus cannot be represented by a symmetrical network.
To implement the molecule, we must construct an asymmetrical circuit. The proposed one uses the potential −V i,j . It is obtained by a linear VCVS e 0 of value −1. (VCVS stands for Voltage-Controlled Voltage Source). The circuit is drawn on Fig. 4.
The resistances ρ 1 connect the potentials V i and V i±1 and the resistance ρ 2 connects the potentials V i and −V i+1 . Applying KCL at node V i and after some algebraic manipulations, one can obtain (8).
The molecule of the third derivative is represented by (10). Like every odd order derivative, the third molecule is not symmetrical and thus, it cannot be represented by a symmetrical network. Moreover, it expands over two steps at each side.
The proposed circuit is represented on Fig. 5. The circuit is a combination of the method used in molecule 40 to obtain the value of the potential, and those used in molecule 10 to get an asymmetrical value in the molecule. The paths 1 and 2 code the outer values of the molecule {−1, 1}, while the paths 3 and 4 code its inner values {2, −2}. Note that the minus sign is inverted between these outer and inner values, this is translated into a direction change in the network (compare paths 2 and 4).
The resistances ρ 1 connect the potentials V i and V i±2 and the resistance ρ 2 connects the potentials V i and −V i+2 . In a similar way, the resistances ρ 3 connect the potentials V i and V i±1 and the resistance ρ 4 connects the potentials V i and −V i−1 . Applying KCL at node V i and after some algebraic manipulations, one can obtain (11).
and r 0 = −2h 3 ρ 0 /12, then the potential at node V i estimates the third derivative of the function v(x) at a given point (Eq. (12)).
The molecule of the 11 derivative is represented by (13). This molecule is the tensor product of two asymmetrical molecules (the molecules 10 and 01) and is therefore symmetrical.
The proposed circuit is given on Fig 6. Following the same convention and notation as before, ρ 1 corresponds to the resistance values between the potentials V i,j and V i±1,j±1 , and ρ 2 are the resistance values between the potentials V i,j and V i∓1,j±1 . The electronic network which evaluates the molecule 11 is drawn on the center cell of Fig. 6, its eight adjacent cells are also represented. To facilitate the understanding and the analysis of the network, only the center cell is completely drawn. On all the other cells, only the resistances along the paths which link there cell center (V i−1,j+1 , V i,j+1 , etc) to V i,j have been represented. Fig. 6. Network of molecule 11.
Applying KCL at node V i and after some algebraic manipulations, one can obtain (14).
The actual resistances are represented by thick vertical and horizontal segments. The resistances represented by long segments can be taken as r {long} = ρ 1 /2, while the ones represented by short segments can be taken as r {short} = ρ 2 /2. For example, along the path from V i,j to V i+1,j+1 , there are 4 resistances (of the same value) in series which are in parallel with 4 others, so the equivalent resistance is The network has been designed in such a way that one gets a maximum number of resistances having the same value, up to the sign. It seems that such resistances are easier to integrate in practical circuits. It is not difficult to derive an equivalent network, from the proposed one, with a much lower number of resistances.
The molecule 22 derivative is represented by (16). This molecule is the tensor product of two symmetrical molecules (the molecules 20 and 02) and is therefore symmetrical.
The proposed circuit is given on Fig 8. As shown in the network, ρ 1 correspond to the resistance values between the potentials V i,j and V i±1,j and between V i,j and V i,j±1 , and ρ 2 are the resistance values between the potentials V i,j and V i±1,j±1 and between V i,j and V i∓1,j±1 . The electronic network which evaluates the molecule 22 is drawn on the center cell of Fig. 8, its eight adjacent cells are also represented. To facilitate the understanding and the analysis of the network, only the center cell is completely drawn. Applying KCL at node V i and after some algebraic manip-ulations, one can obtain (17).
The actual resistances are represented by thick vertical and horizontal segments. The eight resistances represented by very long segments can be taken as r {very long} = ρ 1 /2, while all the other ones (i.e. resistances represented by long and short segments) can be taken as r {long,short} = ρ 2 /2. The explanation of this choice is identical to that of the molecule 11.

VII. COMPUTATION OF ∂ x 2 y (MOLECULE 21)
The molecule 21 derivative is represented by (19). This molecule is the tensor product of one symmetrical molecule and one asymmetrical molecule (the molecules 20 and 01) and is therefore asymmetrical.
To define its electronic circuit, the molecule is splitted into 3 smaller molecules crossing the point (i, j), the vertical one (−2, 0, 2), and the molecules along the two diagonals (−1, 0, 1). Each submolecule can then be implemented by the same method as the one used for the molecule 10.
A careful drawing of the location of resistances to obtain a periodicity of the cells leads to the network of Fig. 9. Only the center cell is completely drawn to facilitate the reading of the network. The resistance between one of the 9 paths linking the potentials of adjacent cell to V i,j are called ρ α . The bold index α on the scheme indicates their corresponding path.
If we choose r 0 = −2h 3 /12 and the sets of ρ α so that, the potential at node V i estimates the 2-1 partial derivative of the function v(x) at a given point (Eq. (21)).
The actual resistances are represented by segments between two dots. The 6 resistances along the paths 1 and 2 can be taken as r 1 , 2 = ρ 1 /2, while the resistances along the paths 3, 4, 5 and 6 can be taken as r 1 , 2 , 3 , 4 , 5 , 6 = ρ 2 /3. The value of the actual resistances along the paths define, by periodicity, the resistances in the elementary cell.

VIII. COMPUTATION OF ∂ x 3 y (MOLECULE 31)
The molecule 21 derivative is represented by (22). This molecule is the tensor product of two asymmetrical molecules (the molecules 30 and 01) and is therefore symmetrical.
The molecule is separated into two submolecules arround the center point (i, j). The first submolecule is the inner one (−2, 2, −2, 2), it can be defined under a circuit form from the circuit of the molecule 11 by inverting and scaling the resistance ρ α . The circuit of the second molecule, the outer one (1, −1, 1, −1) is drawn on Fig 10. As indicated in the network, ρ 1 correspond to the resistances between the potentials V i,j and V i±2,j±1 and ρ 2 are the resistances between the potentials V i∓2,j∓1 . The circuit of the molecule 31 is the "schematic sum" of the network of Fig. 10 and the network of molecule 11 (with resistances renamed ρ ′ 1 and ρ ′ 2 ). A "schematic sum" means that the two circuits must be superposed with separated dots in the cell border, and that these two circuits have the same current source i 0 . Applying KCL at node V i and after some algebraic manipulations, one can obtain (23).
The actual resistances of Fig. 10 are represented by segments between two dots. The 4 resistances along the paths 1 and 2 can be taken as r 1 , 2 = ρ 1 /4. The value of the actual resistances along the paths define, by periodicity, the resistances in the elementary cell.

IX. COMPUTATION OF A PDE
A linear PDE up to the order four can be written under the general form (25). The coefficient ρ ij in front of each derivative is the factor ρ 0 defined in the study of the circuit of the corresponding molecule. The coefficient a ij is simply a multiplicative term. For the PDEs coming from physics (Poisson's equation, deflection of plates, etc... ) the ratio a ij /ρ ij is usually an integer.
a 00 ρ 00 V + a 10 ρ 10 ∂V ∂x + a 01 ρ 01 ∂V ∂y + · · · + a 04 ρ 04 The elementary cell of a linear PDE is obtained by superposing the corresponding cell of each of its terms. This superposition is carried out in the following way: • The nodes on the boundaries of each collected cell must remain disjoint. • The voltage nodes V i,j , as the current sources i 0 , of all the collected cell must be merged. • The resistances of each assembled cell must be scaled with the corresponding factor a ij . In the case of nonlinear PDE, the parts under the form of a product of derivative, could be implemented by multiplying the voltage nodes V i,j of each term and then by connecting the result to the current source i 0 .

X. BOUNDARY CONDITIONS
Consider for example the deflection of a square clamped plate (26-28) to illustrate the implementation of the boundary conditions. The Dirichlet boundary conditions (27) impose the value of the function along the boundary. The Neumann boundary conditions (28) impose the normal derivative along the boundary.  The cell which implements the Neumann boundary conditions is built by simply replacing the current source i 0 of each elementary cell by a voltage-controlled voltage source v −1 , in the case of condition (28), the VCVS has the value of the voltage node V 1 at the inner side of the boundary. This is the PNR-circuit counterpart of how the finite difference method imposes a null derivative of a function along the boundaries. Other Neumann conditions are built in the same way. The cell of this Neumann boundary conditions for the fourth derivative term is given in Fig. 12. These cells are located arround a cross on Fig. 1. Notice that for the cell at the left part of the boundary the resistances r 1 , r 2 and r 3 are not connected and could be removed.

XI. SIMULATION AND COMMENTS
We have developed a program which generates the Spice network corresponding to the chosen finite difference discretisation scheme from a (linear) partial differential equation. The generated file contained the N or N 2 cells of the domain, and the cells at the boundaries to take account of the boundary conditions. All the circuits defined in the previous sections have been independently checked by simulating the PDE: with α+β 4 on the domain Ω = [0, 1]×[0, 1], with the right boundary conditions. The Spice solution for the potential V i,j has been compared with the PDE solutions v (x, y) = x α y β for different mesh sizes. On these academic examples, there is no significant difference between the solutions.
The case α = β = 0 corresponds to the Molecule 00, which is the trivial molecule ?>=< 89:; 1 . The unit circuit which "codes" the equation V i = ρ 0 i 0 is given in Fig 13. This case is important, it allows to manage the terms of PDE which involves the function itself (or derivative of order 0 of a function to justify the terminology Molecule 00). Fig. 13. Network of molecule 00.
We can illustrate the method on a classical partial differential equation. The analog computing of the Poisson's equation A comparison of the solutions computed by the direct circuit simulation with those obtained by a numerical simulation has been done for the Poisson's equation circuit with M × N = 10 × 10 cells (see Figure 15). In the following simulation results, the comparison focus only on node voltages magnitude v.
The direct simulation of the periodic circuit has been made with Spice. The maximum amplitude is 1.4723V . The numerical computation of v has been done using the Finite Difference Method (FDM) on a regular mesh of 40 × 40 squares. We must emphasize that the mesh size was chosen to obtain an accurate numerical solution of the PDE and is not related at all to the number of the cells of the circuit. The maximum amplitude is 1.4742 V . In Figure 15 the continuous FDM solution v (x, y) is represented by the mesh while all voltage nodes v i,j computed by Spice are located by bullets.

XII. CONCLUSION
We have proposed a concept of periodic network of resistances (PNR) to define electronic circuits which are able to solve PDE. We adopt an electronician point of view, i.e. designing the electronic circuits and simulating them. The matrix equations of our periodic circuit are the same as the corresponding set of finite difference equations to resolve a given PDE. Hovewer, the accuracy and the stability of the numerical schemes have to be carrefully studied.
By doing a "schematic sum" of some elementary networks defined in the previous sections, it is possible to build the 2D electronics circuit which is a solution of a given linear partial differential equation, for example the plate equation: Moreover, generalization to non constant coefficient PDE, non-elliptic PDE or PDE with an arbitrary function in the right hand side is faisable by the same method. By using multiplior, even nonlinear PDE could be manageable. All these generalizations suppose implicitly that the finite difference equations relative to the periodic network are mathematicaly valid. An example of application of PNR could be the real time control of a system, with a PDE control law implemented under the form of an analog integrated circuit instead of a numerical one. It is difficult to quatify the gain in speed, that may be of a factor 1000.
Of course, a practical implementation of the method described in this paper under the form of integrated circuits is challenging, some of the numerous problems are the integration of the many NIC (Negative Impedance Convertor) to simulate the negative resistances: How to choose the magnitude of the elements ? Should the scaling operate on the resistances or on the current sources ? What is the effect of the scattering of the component values (some finite difference schemes are very sensitive to it, but finite difference schemes deriving from physical equations have "generally" a good behavior).