Analysis of switched-capacitor circuits using driving-point signal-flow graphs

This paper extends the driving-point signal-flow graphs to switched-capacitor (SC) circuits by introducing a new theoretical element: an auxiliary voltage source that transfers no charge. In contrast to existing SFG methods, our method has no restrictions as to what types of SC circuits can be analysed, it requires no equivalent circuits or tables, and it works with two-phase as well as multi-phase SC circuits of any complexity. Compared to charge-equation matrix methods, it requires more effort, but is better suited for hand analysis because it makes causal relationships visible. Three illustrative examples are given to show the efficiency of the method and present a few application hints: a voltage doubler, the standard SC integrator, and a four-phase circuit simulating an inductor.


Introduction
Switched-capacitor (SC) circuits are analog discrete-time circuits that consist of switches, capacitors, and amplifiers. They are as diverse and powerful as they are hard to analyse. This is mainly so because SC circuits reconfigure between phases, so they will not have the same topology, nor even the same number of nodes, in the different phases.
Signal-flow graphs (SFGs) are often used for hand analysis of circuits (and for teaching), because they make causal relationships visible and provide insight into a circuit even before transfer functions are calculated [1][2][3][4]. While the SFG analysis of arbitrary continuous-time circuits is solved, there is no SFG method yet that can be applied to all SC circuits without restrictions on their structure or number of phases. In this paper we propose such an SFG method, compare it to previous SFG method to highlight the improvement, and compare it to the established charge-equation method to show that it describes the same equation systems in a different way.
The DPSFG was introduced in 1998 [5] when Ochoa combined driving-point impedance techniques with Mason signal-flow graphs (SFGs) [6,7]. Ochoa used auxiliary voltage generators (which we call auxiliary voltage sources or aux sources) and explained the SFG derivation by first splitting schematics into sub-circuits and then coupling the sub-circuits with voltage-controlled current sources.
We cast Ochoa's method into a different cognitive framework in [8], making it unnecessary to resort to intermediate circuit representations and simplifying the application of the method to the point where no written material or tables at all are required to use it. This way of explaining the method followed closely Mason's original idea to use SFGs because they ''offer a visual structure, a universal graph language, a common ground upon which causal relationships among a number of variables may be laid out and compared [7]''. The one-to-one correspondence between circuit and graph also made it straightforward to derive transposed circuits, as shown in [8].
Concerning SC networks, two groups devised and published general analysis methods in parallel. The first publication used the concepts of nodal analysis on SC circuits and also extended classical two-port theory to SC circuits [9]. In the second publication, the derivation of the matrices is simpler, because it uses a switch matrix that is straightforward to derive, but since the author does not use the z transform, but only discrete-time equations with time indices, the derivations do not result in z transfer functions [10].
The nodal analysis was then further simplified by using the indefinite admittance matrix (IAM) in [11] and was extended to the analysis of multiphase SC networks [12]. Since there is a very close correspondence between signalflow graphs and IAMs, SFGs started to be used at that point in time because, by making causal relationships and signal flow visible, authors gained more insight into the circuits.
The first comprehensive paper about SFG analysis of SC circuits was published in 1984 [13], it bases on the theory introduced in [11] and makes it possible to easily assemble signal-flow graphs of SC filters by first tabulating basic building blocks and the corresponding graph parts. This method was extended to multiphase SC networks [14].
Both [13] and [14] succeed in making the signal flow visible, but the graphs are not one-to-one maps of the circuits. Furthermore, the methods only work for sourcesink-node (SSN) networks [14], which means that all nodes are either at a voltage source or at an amplifier input. Many circuits, e.g., those doing correlated double sampling and voltage multiplication, are therefore not covered.
These methods assume from the start that there are no parasitic capacitors and that the gain of the op-amps used in the SC circuits is infinite. They cannot possibly analyse the effect of finite opamp gain on the transfer functions, which is quite important, as was, e.g., discussed in [15] for SC integrators.
This was remedied in Chichocki et al. [16] who made it possible to analyse any SC circuit, with non-ideal op-amps and also parasitic capacitors. They used a small number of equivalent z-domain circuits and the Coates flow-graph technique. In their own words, they ''require a smaller number of elements than the equivalent circuits previously proposed''. But they still require equivalent circuits containing some additional elements.
In a previous paper we applied DPSFGs to the noise analysis in SC circuits, [17]. In this paper, we introduce a DPSFG analysis method for SC circuits that requires no equivalent circuits whatsoever and imposes no restrictions on amplifier gains, parasitic capacitors, or the number of clock phases. What results is a general method for visualbased (i.e., intuitive) analysis of switched-capacitor networks of any complexity, with any number of phases.
As in all cited papers, the theory is simple (here it is the introduction of a new type of auxiliary voltage source). We want to show that this one intuitive understanding-the understanding of that new aux source-makes it possible to analyse continuous-time and switched-capacitor circuits with the same SFG technique. The purpose of this paper is to extend and unify SFG analysis and enable SFG users to tackle all SC circuits as easily as they can treat continuoustime circuits.
Therefore, we proceed as follow: In Sect. 2, we give a short introduction to SFGs and DPSFGs of continuous-time circuits. Our explanation complements [4,5,8]; it is structured such that the introduction of the new theoretical element is as simple as possible in Sect. 3, where we also give a minimal example. Section 4 contains the main example in which we show how to efficiently apply our method. Section 5 extends this to multi-phase networks and shows how to calculate a discrete-time input impedance. Finally, in Sect. 6 we show what to do when graphical evaluation becomes too complicated.

Driving-point signal-flow graphs (DPSFG) of continuous-time circuits
This section presents the DPSFG method introduced in [5], brought to its present form in [8], and provided as a sequence of video tutorials in [4].

Signal-flow graphs and Mason's gain rule
Mason graphs are just graphical representations of linear equation systems. The so-called branches are multipliers, and the so-called nodes are variables. All branches going into a node are added. Therefore, each node that has branches going into it describes a linear equation, and the example shown in Fig. 1 represents the following linear equation system: This can also be expressed in matrix notation: where we have simply replaced all numbers 0 by Á such that it becomes immediately apparent that this graph corresponds to a sparse-matrix equation system. It would be equally simple to take any matrix equation, e.g., an IAM representation of a circuit, and draw it as a Mason graph, but that would not necessarily be helpful as it would then not show the causal relationships.
Equation (2) could now be solved by matrix manipulations, but Mason showed a better way in [7]: A transfer function from a node u to a node z can be calculated as The graph determinant D can be calculated as where r 1 is the sum of all loops, r 2 is the sum of products of 2 loops that do not have nodes in common, and r 3 is the sum of products of 3 loops neither of which have nodes in common, and so on. The sum (4) can, in principle, go to r N in a system with N loops, but most feedback systems derive their special properties from interacting loops, so the sum often ends soon in practice. In more technical terms, whenever a matrix representation of the equation system would be sparse, then Mason's formula is efficient. In Fig. 1, we have three loops: L 1 and L 2 have nodes in common, but L 3 touches neither of the other loops, so (4) becomes The P i in the numerator of (3) are the forward paths, the possible ways to get from u to z. The D i are their respective sub-determinants, which are calculated just like D, but only using loops that have no nodes in common with the forward path in question. In Fig. 1, we have two forward paths: Path P 1 touches all loops, so when calculating the subdeterminant, all terms containing loops in (6) disappear and D 1 ¼ 1. Path P 2 only touches loop L 3 , so all terms containing L 3 in (6) have to be omitted to calculate Inserting everything into (3) results in Note that if a different transfer function has to be calculated (e.g., T y ¼ y=u), then only the numerator can change, because D does not depend on the forward paths.
Mathematically, an SFG is just a graphical representation of an equation system, which can even be manipulated and simplified graphically [6]. But as Mason stated there, it can also be used to visualise causal relationships-if it is drawn in the right way.

The driving-point signal-flow graph (DPSFG)
The most important aspect of a signal-flow graph that is to be used for hand analysis is that it maintains and represents causal relationships. The common tool used in literature is source superposition, as it is also done when circuit equations are derived from circuits.
Source superposition will give individual contributions from all sources to all circuit nodes for the condition that all other sources are set to zero. So unless there is a voltage source in every node of a circuit, it is necessary to tabulate all possible networks that could reasonably appear between sources and sinks, which invariably means restricting methods to certain classes of circuits, or to use another analysis method to derive the individual contributions.
As explained in [8], the DPSFG method solves this problem by making sure there is a voltage source at every circuit node. In order to achieve this, the first step in circuit analysis is to introduce an auxiliary source (aux source) at every node that does not yet have a voltage source. This aux source is defined as follows: It produces the voltage that is already present at a node, and therefore conducts no current.
To give an example, Fig. 2(a) shows the standard continuous-time integrator with a resistor having conductance G and a capacitor having capacitance C, and an opamp having gain A. Node 1 is connected to the input source, node 3 to the opamp output, which is also a voltage source, but there is no voltage source on node 2. Therefore an aux source is connected to node 2, resulting in Fig. 2 To draw the DPSFG, we can first draw all nodes. Observe that we also draw the current flowing into the aux source, I 2 , in Fig. 3(a).
Then we simply do source superposition, voltage by voltage. V in controls V 1 , and we have Fig. 3(b). V 1 lets a current I 2 ¼ G V 1 flow, which gives Fig. 3(c). V 2 controls the opamp output as V 3 ¼ ÀA V 2 , giving Fig. 3(d). V 3 lets a current I 2 ¼ sC V 3 flow, and it also is the same as V out , and we have Fig. 3(e). This is almost all, but now we have a problem: I 2 is not zero! The aux source was defined to conduct no current. However, I 2 is not the true current into the source, but the current flowing into the aux source for V 2 ¼ 0: the shortcircuit current of the aux source. So one part of the source superposition is missing. In order to have zero current under superposition, we now need to calculate the necessary voltage V 2 ¼ Z 2 I 2 to have I 2 flow out of the aux source when all other sources are set to zero. This is Fig. 3 This Z 2 is the driving-point impedance (DPI) which gives the DPSFG part of its name. Since we have made sure that every node now has a voltage source attached, all of which are set to zero during superposition, this DPI is always the inverse of the sum of all conductances connected to the node. Here, This is all, and it becomes apparent that the result is an SFG that shows the causality as well as describing the equation system of the circuit. There is only one loop, and only one path, which even touches the loop. Therefore: where the last approximation is for A ! 1.
Additional things can now be calculated with comparatively small effort. For example, what happens if a parasitic capacitance C p is connected to node 2? Then only Z 2 changes to What if we need to know the input admittance? The input current is I in ¼ ðV 1 À V 2 Þ G. Drawn as SFG branches, this gives Fig. 3(g). The only thing that changes are the forward paths, and we get: Analog Integrated Circuits and Signal Processing For A ! 1, this gives Y in ¼ G, as expected, because then node 2 is virtual ground. And this method always works, for every linear circuit. The only thing that happens for larger circuits is that there are more aux sources, and of course Mason's rule becomes more tedious to evaluate, but the derivation of the branches is never different from, or more difficult than, what we just showed. Note, however, that from a didactic point of view, it can be better to fill in the DPSFG in a different order (c.f. [4,Chap. 7]).
This method works because the aux sources make Kirchhoff current-law (KCL) equations appear in the SFG. We illustrate this by looking at the equations for the auxsource-related nodes I 2 and V 2 : We insert (22) into (21), take all sum terms to the same side, and collect expressions with G and sC: And we see that this is the KCL for node 2.
In summary, the connection of an aux source causes the KCL equation for that node to appear in the graph. This is basically why the DPSFG method must work for any linear circuit.

The charge-mode auxiliary source
At this point we are ready to extend the method to SC circuits. After the lengthy introduction it is clear that (a) SC circuits reconfigure when the phases change, and (b) if the proper auxiliary source is chosen, then source superposition becomes simple and drawing an SFG becomes straightforward.

Theorem 1
The DPSFG method as introduced in [5], and as presented without the need to split a circuit into subcircuits in Sect. 2, works for all SC circuits with an arbitrary number of phases if (a) one circuit diagram per phase is drawn, with switches shorted or opened, and (b) in each of those circuit diagrams, auxiliary sources are attached to all nodes that are not connected to a voltage source. Such an aux source is defined as follows: It produces the voltage that is already present at a node, and therefore conducts no charge.
The validity of this theorem will now be demonstrated by giving examples. Figure 4 shows an SC voltage doubler, a simple SC circuit, which nevertheless is not a source-sink network. The figure shows the circuit with the switches, and the two circuit configurations in the odd phase (phase 1) and the even phase (phase 2). We have given all nodes individual numbers.
Node 1 in the odd phase is at a voltage source terminal, but nodes 2 and 3 in the even phase are not, so there we attach aux sources. This time the sources have a voltage V e 2;3 and a charge Q e 2;3 that is the charge flowing into the source while the even-phase circuit settles to its final state. Now we can construct the DPSFG just as in Sect. 2.2 but by using charges and voltages instead of currents and voltages. First we write down all nodes of the SFG, in two rows for the two phases, giving Fig. 5(a). Now source by source: V o in is V o 1 : Fig. 5(b). If all sources other than V o 1 are set to zero, then C 1 and C 2 are discharged into the aux sources V e 2 and V e 3 . Observing the respective plates of the capacitors discharged into the nodes, and since this is another phase, Fig. 5(c). If V e 2 is switched on and all others are zero, a charge Q e 3 ¼ C 2 V e 2 flows: Fig. 5(d). Similarly, if V e 3 is switched on and all others are zero, a charge Q e 2 ¼ C 2 V e 3 flows, and this also is V e out : Fig. 5(e). As in Sect. 2.2, we now have calculated the short-circuit charges of the aux sources and need a last step answering the question: what should the voltage of the aux sources be such that, under superposition, the aux source conducts no charge? This gives Fig. 5(f). Before, we got driving-point impedances in this step, and now we get S e 2;3 , which are two driving-point elastances, an elastance being an inverse capacitance (unit F À1 ), and obviously for the exact same reason why driving-point impedances have this structure. Mason's rule is simple since there is only one loop and two forward paths that touch the loop: This is very straightforward to extend, particularly with an analysis of parasitic capacitances, which eluded most of the hitherto presented SFG methods. For example, if a C p is attached to node 2, the only thing that changes is S e 2 in (24): Then (30) can be calculated again: As with the continuous-time DPSFG in Sect. 2, every problem could be solved by step-wise source superposition, but changing the order of how the DPSFG is drawn may make an analysis quicker and less error-prone.

Main example: the SC integrator
This is already all of the theory, but the strength of the method lies in the systematic application, so it is best exemplified with a simple (but sufficiently complex) SC circuit: the stray-sensitive SC integrator in Fig. 6.
In the previous section we just numbered all nodes, but now we introduce a technique that makes the node numbers in all phases unique: the nodes are numbered with powers of 2, i.e., 1; 2; 4; 8; 16; . . ., in the schematic with switches, Fig. 6 (top). During the derivation of the odd and even schematics in Fig. 6, nodes of the original circuit are connected by closed switches, and the number of such a new node is the sum of the numbers of the nodes connected.
Since the original numbers were powers of 2, the new nodes are guaranteed to have unique numbers and contain local topological information, with the effect that all elastances S / i with the same i are the same, irrespective of their phase /.

Calculation with charge equations
The circuit in Fig. 6 can now be analysed using charge equations. The capacitor charges are in the odd phase and in the even phase, where Q 1;2 is the charge on C 1;2 with the positive plate of the capacitors indicated in Fig. 6. The amplifier with gain A works in both phases, giving two equations: Now it is straightforward to write down the charge differences on the capacitors during the phase transitions: Then there is charge conservation. In each phase, there are nodes without voltage sources on them where no charge can appear or disappear. In the odd phase, this is only node 4, giving the equation In the even phase, only node 6 has no voltage source attached, giving one more equation: Observe how both DQ e 1 and DQ e 2 are on the same side of the equation; this is because both capacitor plates connected to node 6 in the even phase have the same sign on them.
The final step is to specify the input and the output of the circuit: The equation system (33)-(40) can now be solved for V o out : Considering that our z À2 means we have one delay from an odd phase to the next odd phase (see ''Appendix''), this is the well-known transfer function of the SC integrator (e.g., [15]).

Calculation with a driving-point signal-flow graph
The same calculation can now be done with signal-flow graphs. In this section we show how to do it systematically. First, auxiliary sources have to be attached to all circuit nodes that do not already have voltage sources on them. These are the nodes 4 and 6 in Fig. 7. Second, a signal-flow graph can be prepared, as shown in Fig. 8(a). This graph contains the input and output paths, the amplifying paths, and, in every node where there is an Fig. 6 Schematics of the simple, stray-sensitive SC integrator. The top schematic has all the switches; the centre schematic shows the odd phase 1, and the bottom schematic shows the even phase 2 Third, all paths due to all voltage sources can be drawn, leading from Fig. 8(a, b). We will just explain two of these paths. First, let all voltage sources be zero, and then let V o 8 be switched to its value. Then the charge C 2 V o 8 will flow into the aux source at node 4, giving a contribution to Q o 4 . This is the path with weight C 2 from V o 8 to Q o 4 . However, this is not the only effect V o 8 will have, because in this step of calculating source superposition, all other voltages are zero, particularly all voltages in the even phase. The capacitor C 2 has now been charged to ÀC 2 V o 8 in the odd phase (observe the þ sign on the left plate of C 2 to understand why the charge is negative), and in the even phase it is completely discharged because the voltage is zero on both sides. Therefore the charge ÀC 2 V o 8 flows into the aux source at node 6. Since this is after a phase change, we get a path Àz À1 C 2 from V o 8 to Q e 6 . All other missing paths can then be found by applying the same method to V o 3 , V o 4 , V e 6 , and V e 8 , and then the DPSFG is complete.
Finally, it would already be possible to apply Mason's gain rule, and if our only goal were to get a transfer function as quickly as possible, then we would stop here. However, using one of the graphical simplification methods mentioned in Sect. 2, Fig. 8(b) can readily be re-drawn as Fig. 9. Then a term C 2 ðA þ 1Þ appears, and it becomes visible in the graph already that C 2 is subject to the Miller effect.
The graph in Fig. 9 has only three loops: According to Mason's gain rule, the graph determinant becomes Note that the term L 1 L 3 does not appear, because Loop 1 and Loop 3 have nodes in common, and the same is true for L 2 L 3 . There is one path from V o in to V o out , and since this path touches all three loops, its sub-determinant D 1 ¼ 1.
With this information, the transfer function of the circuit can be calculated as which, of course, gives the same result obtained in (41). As before, calculating the effects of parasitic capacitors would be simple with this method as they would only change the elastance equations (42).

Comparison of the two methods
This DPSFG method gives a direct, graphical way to draw a graph from a circuit and solve it using Mason's gain rule. What remains to be shown now is how the two methods relate.
In the classical method, we wrote down charge conservation equations at all nodes where there were no voltage sources. In the DPSFG method, we introduced aux sources to all nodes where there were no voltage sources. So it   Fig. 8(b) by graphical simplification. In Fig. 8(b), there were two ways to go from V o 4 to Q e 6 , one had the weight z À1 C 2 , the other the weight A z À1 C 2 . Their sum is z À1 C 2 ðA þ 1Þ, which is the path from V o 4 to Q e 6 in this SFG. The same was done for the path from V e 6 to Q o stands to reason that the SFG branches pertaining to the aux-source nodes should be charge conservation equations, just as they were KCL equations in the continuous-time case.
Proof for (39): First, substitute (36) into (39): Now substitute (33) and (34) into (48): This can be re-ordered as follows: Now let us look at Fig. 8(b), at the DPSFG branch from Q e 6 to V e 6 . It has the weight S e 6 , which corresponds to the equation Substituting (42) into this and solving for Q e 6 gives: This can now be used together with (50) to calculate Q e 6 : and it is immediately apparent that this is the equation obtained from adding all branches going into the node Q e 6 in the DPSFG. In summary, the connection of an aux source causes the charge conservation equation for this node to appear in the graph. This is why the method will work for any SC circuit operating with voltage settling and charge transfer.

Second example: multi-phase SC networks and discrete-time input impedance
The real power of this method lies in the fact that it only uses very basic principles, and therefore can be applied to any switched-capacitor circuit that works with settling and charge transfer. We will now demonstrate this using the simulated inductor in Fig. 10. This circuit does not only have four phases, but it requires the calculation of the input charge difference DQ in that the input source V in has to provide. The first step is to draw the circuit in its four phases; we call them a, b, c and d instead of e and o. This is shown in Fig. 11.
Drawing the DPSFG can be done as before: first, we draw the graph in Fig. 12(a) by observing the source superposition theorem in all four phases individually, and then we obtain Fig. 12(b) by inserting the paths showing charge transfer between the phases.
The elastances on the branches of Fig. 12 are: and as explained in Sect. 4, all elastances with the same index number are the same, even if they are in different phases. Two more things should be noticed when going from Fig. 12(a, b): first, we have applied the simplification demonstrated in Fig. 9, which is why the integrating capacitor C 2 gives rise to four paths with the weight ð1 þ AÞz À1 C 2 . Second, the input charge DQ in has two contributions. At the beginning of phase a, the charge of C 0 is C 0 V d 18 . At the end it is C 0 V a 2 . The charge DQ in that has to be supplied by the voltage source is therefore This gives the two branches from V d 18 and from V a 2 to DQ in , and the signal-flow graph is complete.
At this point we can evaluate Mason's gain rule: There is one loop each in each of the four phases, and one big loop going through everything: The graph determinant is: There are three forward paths, with the respective sub-determinants Now it is possible to compute the discrete-time input admittance, which has the unit F, as follows:  Fig. 12 The DPSFGs derived from Fig. 11. Graph a shows the portions of each graph derived from the individual phases, b adds the paths for charge transfer between phases and for calculating the input charge according to (58) the latter line comes from the fact that D 1 ¼ D.
Substituting (54) to (70) into (72) gives a rather complicated expression, which we will not reproduce here. In [16], the whole calculation is made for A ! 1; if we do this as well, we get which is the same as Eq. (15) in [16]. This acts like a discrete-time simulated inductor for a specific choice of capacitor ratios: Since this is a four-phase system which only samples the input at phase 1, this corresponds to relative to phase one (see ''Appendix''). As expected for an inductance, the admittance has one pole at z ¼ 1 (which corresponds to f ¼ 0).
6 When Mason is too difficult...
The example in the previous section resulted in a very complicated and tedious evaluation of Mason's gain rule. Clearly, this example is at the border of what is reasonably calculatable by hand. Note, however, that this only concerns Mason's rule. We can still use the DPSFG method to obtain a correct graph that corresponds to the circuit diagram (which makes it easily reviewable), but then use a computer algebra tool (Mathematica, Python/NumPy/SciPy,...) to actually evaluate the equation system. Such a computation entails listing the dependent variables, listing all equations (one per node), and then solving the equations and substituting the elastances.
This can easily be done for the graph in Fig. 12, but we will show the script for the graph in Fig. 9 such that it is more instructive to the reader. In Mathematica 1 notation, this would be: A brief explanation: In Mathematica, = is an assignment, == an equality, -[ a substitution rule, and the operator /. performs substitutions. The solution sols consists of one substitution rule for every dependent variable. Therefore, the last line calculates the transfer function voo/vio, substitutes the solution for voo and then the elastances elast, takes the lim A!1 , and finally simplifies. The result of this code is: This is evidently the same as (41).

Conclusions
In this paper, we have adapted the driving-point signal-flow graphs (SFGs) to switched-capacitor circuits by one main modification: we have exchanged the zero-current auxiliary source by a zero-charge aux source. This let the drivingpoint impedances become driving-point elastances.
We have demonstrated the application of the method by examples, in which we also showed a few tricks helping simplify the analysis even further, and have shown by argument that this method is applicable to any SC circuit relying on charge transfer just as the original method is applicable to any linear continuous-time circuit. Therefore, we have unified the SFG analysis methodology for continuous-time and SC circuits. Our method only requires that the users remember Mason's gain rule and understand source superposition and the aux source. Written material like tables of partial transfer functions, equivalent circuits, etc., are unnecessary. In this respect our method is unique among graph-based SC analysis methods.
Of course, the same analysis can be done by deriving and solving charge equations, which would give the same symbolic results with fewer steps and more quickly using a computer algebra tool. The graph method, however, maintains causal relationships and makes them visible, and when studying and comparing different SC structures, this insight is often worth the additional effort.
While the last example gave a very complicated graph determinant and is at the limit of what can be done by hand (in fact, we, the authors, evaluated Mason's rule using Mathematica to be sure we do not make stupid mistakes), the one-to-one correspondence makes it straightforward to review a drawn SFG versus its circuit. It is uncommon with this method that an equation is omitted, and near impossible to introduce redundant equations, which is something that everyone who has analyzed circuits by hand using equations knows all too well.
Using Mason's gain rule to derive transfer functions can give more insight into a circuit, but if evaluating Mason's rule becomes too complicated, then calculating transfer functions from a signal-flow graph can also be done directly by entering the euqations into a computer algebra tool.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creative commons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.