Electrical Circuits as Dynamical Systems

Abstract

An electrical circuit containing at least one dynamic circuit element (inductor or capacitor) is an example of a dynamic system. The behavior of inductors and capacitors is described using differential equations in terms of voltages and currents. The resulting set of differential equations can be rewritten as state equations in normal form. The eigenvalues of the state matrix can be used to verify the stability of the circuit. The most fitted numerical methods to integrate electrical circuit differential equations are the Euler Method (Forward and Backward), the Trapezoidal Rule, and the Gear Method of second to sixth degree, for circuits having stiff equations. These methods are implemented, with adjustable time-step integration, in the majority of circuit simulation software, such as SPICE. The analytical solution can also be computed, for small-size circuits, applying the Laplace Transform. It is interesting to compare the graphical presentation of numerically and analytically obtained solutions. While the numerical methods can be used for both linear and nonlinear circuits, the Laplace Transform is mostly used for linear circuits. A method of using it for nonlinear circuits is also presented.

1. Introduction

The existence and uniqueness of a dynamic circuit solution are strongly tied to the existence of the state equation in normal form. If the equation $\dot{x}=f\left(x,t\right)$ exists, then it can be proved, as it can be found in the mathematic literature, that if $f$ is Lipschitzian (for any $x_1$ and $x_2$ and any $t$, $‖f\left(x_1,t\right)-f\left(x_2,t\right)‖\le k\bullet ‖x_1-x_2‖$ where $k>0$ and $‖‖$ is the Euclidian norm) and if the function $f\left(0,t\right)$ is uniformly continuous and bounded, then the state equation has an unique solution for any initial state $x_1=x\left(t\right)$. The existence of the normal form of the state equation is related to the existence and uniqueness of the resistive multiport solution for any values of the source parameters (which replace the dynamic elements) connected to the ports, and the existence of nonzero dynamic capacitances and inductances for any values of the control parameters of these elements.

It can be seen as in the case of linear circuits $f\left(x,t\right)$ is Lipschitzian: $‖f\left(x_1,t\right)-f\left(x_2,t\right)‖=‖A\bullet \left(x_1-x_2\right)‖$ because there is always a constant $k$ such that $‖A\bullet \left(x_1-x_2\right)‖\le k\bullet ‖x_1-x_2‖$ [1, 2, 3].

In the case of nonlinear circuits without excess state quantities (the circuit does not contain any loop consisting only of independent or controlled voltage sources and/or capacitors and does not contain any cut-set consisting only of independent or controlled current sources and/or inductors), if the characteristics of the dynamic elements are strictly increasing and derivable, then they have a nonzero dynamic parameter at any point of operation. If the resistors’ characteristics are strictly increasing and the resistive multiport does not have loops formed only from voltage sources and cut-sets formed only from current sources, then this resistive multiport has a unique solution. So, there are state equations in the normal form $\dot{x}=f\left(x,t\right)$. For the dynamic circuit to have a unique solution for any initial state $x\left(t_0\right)$, it is enough that $f$ is Lipschitzian [1, 2, 3]. When we have excess state quantities, the problem is treated similarly.

2. State equations for dynamic circuits

The simplest dynamic circuit elements are the linear capacitor and the linear inductor. The operating equation of the linear capacitor is $i_c\left(t\right)=C\bullet \frac{d{v}_c\left(t\right)}{\mathit{dt}}$ where $v_c\left(t\right)$ is the voltage at the capacitor terminals, $i_c\left(t\right)$ is the current through the capacitor, and C is a constant called the capacitor capacity. The operating equation of the linear inductor is $v_L\left(t\right)=L\bullet \frac{d{i}_L\left(t\right)}{\mathit{dt}}$ where $v_L$ is the voltage at the inductor terminals, $i_L\left(t\right)$ is the current through the inductor, and L is a constant called the inductance of the inductor. The ideal dynamic elements are, unlike resistors, lossless elements, i.e., they do not dissipate energy but accumulate it. The energy accumulated at a given moment by such an element can be subsequently transferred to the circuit in which the respective element is connected.

A circuit that contains at least one dynamic element is called a dynamic circuit. The behavior of dynamic circuits, consisting of independent sources, inductors, capacitors, and resistors, is described by a system of differential equations.

2.1 First-order dynamic circuits

A first-order linear circuit contains only one dynamic element (an inductor or a capacitor), other linear circuit elements (resistors, linear controlled sources), and independent sources. The resistive two poles have an equivalent voltage generator (Thevenin) in Figure 1 [4] and/or an equivalent current generator (Norton) in Figure 2 [5] at the input of the dynamic element. With this consideration in mind, it is sufficient to consider only the next two first-order linear circuits.

Equations of Figure 1’s circuits	Equations of Figure 2’s circuits
$R\bullet i_c+v_c=e\left(t\right)$	$\frac{v_L}{R}+i_L=i_s\left(t\right)$
$i_c=C\bullet \frac{d{v}_c}{\mathit{dt}}$	$v_L=L\bullet \frac{d{i}_L}{\mathit{dt}}$
$R\bullet C\bullet \frac{d{v}_c}{\mathit{dt}}+v_c=e\left(t\right)$	$\frac{L}{R}\bullet \frac{d{i}_L}{\mathit{dt}}+i_L=i_s\left(t\right)$
$\frac{d{u}_c}{\mathit{dt}}+\frac{v_c}{R\bullet C}=\frac{e\left(t\right)}{R\bullet C}$	$\frac{d{i}_L}{\mathit{dt}}+\frac{R}{L}\bullet i_L=\frac{R}{L}\bullet i_s\left(t\right)$
If we note:
$v_c=x$	$i_L=x$
$\frac{d{v}_c}{\mathit{dt}}=\dot{x}$	$\frac{d{i}_L}{\mathit{dt}}=\dot{x}$
$R\bullet C=\tau$	$\frac{L}{R}=\tau$
$e\left(t\right)=s\left(t\right)$	$i_s\left(t\right)=s\left(t\right)$
then
$\dot{x}+\frac{x}{\tau }=\frac{s\left(t\right)}{\tau }$	$\dot{x}+\frac{x}{\tau }=\frac{s\left(t\right)}{\tau }$

Figure 1.

A Thevenin first-order circuit.

Figure 2.

A Norton first-order circuit.

So, a first-order linear circuit satisfies the equation:

$$\dot{x}+\frac{x}{\tau }=\frac{s\left(t\right)}{\tau }$$E1

where $x$ is the state variable of the circuit, $\tau$ is the time constant of the circuit, and $s\left(t\right)$ is the parameter of the equivalent independent source.

2.2 Dynamic circuits of second order or greater

The aim is to write the state equations of state in normal form $\dot{x}=f\left(x,t\right)$. Bringing the equations to this form has two advantages, the qualitative properties of the circuit can be studied more easily, and certain numerical methods for circuit analysis can be applied, formulated to solve a system of differential equations written in this form.

There are two cases. In the first case, the circuit does not contain any loop consisting only of independent or controlled voltage sources and/or capacitors and does not contain any cut-set consisting only of independent or controlled current sources and/or inductors. In this case, the state variables are independent of each other, and we say that the circuit has no excess state quantities. In the second case, the circuit does not satisfy the above restrictions, the state variables are not independent of each other, and we say that the circuit has excess state quantities.

2.2.1 Circuits without excess state quantities

Any linear dynamic circuit without excess state quantities can be considered as a linear resistive multiport (containing linear resistors and independent sources) with dynamic elements connected at the ports. The dynamic circuit of order $n>2$ that contains $p$ capacitors and $n-p$ inductors can be represented as follows, in Figure 3.

Figure 3.

A $p$ capacitors and $n-p$ inductors dynamic circuit.

Capacitors are replaced with independent voltage sources and inductors with independent current sources, as in Figure 4. If the following notations are made:

Figure 4.

Capacitors/inductors replaced with independent voltage/current sources.

$x={\left[v_1,\cdots ,v_p,i_{p+1},\cdots ,i_n\right]}^t$—vector of state variables,

$y={\left[\begin{array}{cccccc}{i}_1,& \cdots ,& i_p,& v_{p+1},& \cdots ,& v_n\end{array}\right]}^t$—the vector of the output quantities and

$u={\left[\begin{array}{cccccc}{e}_1,& \cdots ,& e_{\alpha },& i_{\alpha +1},& \cdots ,& i_{\mu }\end{array}\right]}^t$—the vector of the input quantities (independent sources), according to the superposition theorem [6], the circuit equations become $y=k_0\bullet x+k_1\bullet u$, where $k_0$ and $k_1$ are matrices with constant elements.

Using the notation: $\mathrm{\Delta }=diag\left(\begin{array}{cccccc}{C}_1,& \cdots ,& C_p,& {\mathit{L}}_{\mathit{p}+1},& \cdots & {\mathit{L}}_{\mathit{n}}\end{array}\right)$ we obtain $\dot{x}={\mathrm{\Delta }}^{-1}\bullet y$ and $y=\mathrm{\Delta }\bullet \dot{x}$ and the state equations are $\dot{x}={\mathrm{\Delta }}^{-1}\bullet \left(k_0\bullet x+k_1\bullet u\right)$ or:

$$\dot{x}=A\bullet x+B\bullet u$$E2

where $A$ is called the circuit state matrix.

As an example, for the circuit without excess state quantities in the Figure 5 [3, 7], the following state equations are obtained.

Figure 5.

Circuit without excess state quantities.

2.2.2 Circuits with excess state quantities

In this case, there is at least one loop consisting only of capacitors and voltage sources or a cut-set consisting only of inductors and current sources. This means that there will be at least one capacitor that cannot be placed in the tree or an inductor that cannot be placed in the co-tree. The normal tree therefore contains all voltage sources and as many capacitors as possible whose voltages are independent state quantities. The other capacitors will be contained in the normal co-tree, and their voltages are excess state quantities. The normal co-tree contains all current sources and as many inductors as possible whose currents are independent state quantities. The other inductors are contained in the normal tree and their currents are excess state quantities. We consider the circuit as a resistive multiport with dynamic elements connected to the ports. According to the substitution theorem, the capacitors and inductors in the tree are replaced with voltage sources and the capacitors and inductors in the co-tree with current sources. The result is the circuit in Figure 6, in which for simplicity, only one source was represented for each category of element (tree capacitors, tree inductors, co-tree inductors, co-tree capacitors). The second size index represents tree (t) or co-tree (c).

$$\text{Using the notation}:x={\left[\begin{array}{cc}{v}_{\mathit{Ct}},& i_{\mathit{Lc}}\end{array}\right]}^t,x^{\ast }={\left[\begin{array}{cc}{v}_{\mathit{Cc}},& i_{\mathit{Lt}}\end{array}\right]}^t,y={\left[\begin{array}{cc}{i}_{\mathit{Ct}},& u_{\mathit{Lc}}\end{array}\right]}^t,y^{\ast }={\left[\begin{array}{cc}{i}_{\mathit{Cc}},& v_{\mathit{Lt}}\end{array}\right]}^t.$$

Figure 6.

The capacitors and inductors in the tree replaced with voltage sources, and the capacitors and inductors in the co-tree with current sources.

The resistive circuit being linear, according to the superposition theorem we obtain:

$$y=k_0\bullet x+k_1\bullet {\mu }_S+k_2\bullet y^{\ast }$$E3

where $k_0$, $k_1$, and $k_2$ are matrices with constant parameters. The operating equations of the dynamic elements can be written:

$$y=-\mathrm{\Delta }\bullet \dot{x}\text{where}\mathrm{\Delta }=diag\left(\begin{array}{cc}{C}_t,& L_c\end{array}\right),$$

$$y^{\ast }=-{\mathrm{\Delta }}^{\ast }\bullet {\dot{x}}^{\ast }\text{where}{\mathrm{\Delta }}^{\ast }=diag\left(\begin{array}{cc}{C}_c,& L_t\end{array}\right).$$

The excess state quantities can be expressed, with the help of the Kirchhoff's laws, depending on the independent state quantities and the parameters of some independent sources: $x^{\ast }=k_3\bullet x+k_4\bullet {\mu }_S$ and we obtain:

$$y^{\ast }=-{\mathrm{\Delta }}^{\ast }\bullet k_3\bullet \dot{x}-{\mathrm{\Delta }}^{\ast }\bullet k_4\bullet {\dot{\mu }}_S$$E4

But:

$$\dot{x}=-{\mathrm{\Delta }}^{-1}\bullet y=-{\mathrm{\Delta }}^{-1}\bullet \left[k_0\bullet x+k_1\bullet {\mu }_S-k_2\bullet \left({\mathrm{\Delta }}^{\ast }\bullet k_3\bullet \dot{x}+{\mathrm{\Delta }}^{\ast }\bullet k_4\bullet {\dot{\mu }}_S\right)\right]$$E5

so, the last relationship can be written as:

$$\dot{x}=A\bullet x+B\bullet {\mu }_S+B^{\ast }\bullet {\dot{\mu }}_S$$E6

where $A$ is the state matrix of the circuit.

As an example, for the circuit with excess state quantities (capacitor loop and inductor cut-set circuit) and without sources for simplicity, in Figure 7 [8, 9], the state equations are obtained:

Figure 7.

Circuit with excess state quantities.

2.3 Dynamic circuits with resistive nonlinearities

The methods presented above for writing state equations can also be used and in the case of dynamic circuits with resistive nonlinearities, with only a few changes. A first approach is to approximate the nonlinear characteristic of the resistor with piece-wise linear characteristic. The state equations for the nonlinear circuit can be written separately for each linear portion of the piece-wise linear characteristic. For example, for the rectifier with a diode and the $\mathit{RC}$ load in Figure 8, we have two states. One in which the diode conducts and is equivalent to a very low value resistor, ideal $0\mathrm{\Omega }$ as in Figure 9 (state 1) and the other in which the diode does not conduct and is equivalent to a very high value resistor, ideal $\infty \mathrm{\Omega }$ as in Figure 10 (state 2).

Figure 8.

A diode rectifier with $\mathit{RC}$ load.

Figure 9.

Diode conducts (state 1).

Figure 10.

Diode does not conduct (state 2).

From the equations obtained separately for the two states, the following state equation can be written:

$$\frac{d{v}_c\left(t\right)}{\mathit{dt}}=\left(-\frac{1-D}{R_E·C}-\frac{1}{R·C}\right)·v_c\left(t\right)-\frac{1-D}{R_E·C}·e\left(t\right)$$E7

Where $e\left(t\right)=5\bullet \mathit{sin}\left(2\mathit{\pi f}\bullet \left(t+t_0\right)\right)\left[V\right]$ and $f=10\mathit{kHz}$.

When $D=0$ the equation for state 1 is obtained, and when $D=1$ the equation for state 2 is obtained. $D$ acts as a closed-open switch over a well-defined time interval depending on the state of the diode. If the rectifier diode conducts, $D$ acts as an open switch, and if the diode does not conduct, $D$ acts as a closed switch [10].

Another approach to writing the state equations for dynamic circuits with resistive nonlinearities is to replace in the state equations the nonlinear resistance with the function that describes the nonlinearity (Figure 11). In the example of a rectifier with a diode and $\mathit{RC}$ load, the diode can be modeled as a nonlinear piece-wise linear function: $R_D=0.1\mathrm{\Omega }$ if the voltage at its terminals is positive $v>0$, and $R_D=10M\mathrm{\Omega }$ if the voltage at its terminals is negative $v<0$, Figure 12.

Figure 11.

A rectifier with $\mathit{RC}$ load; diode replaced with a nonlinear resistance.

Figure 12.

The nonlinear resistance function.

The following state equation is obtained:

$$\frac{d{v}_c\left(t\right)}{\mathit{dt}}=-\frac{R_D+R_E+R}{\left(R_D+R_E\right)·R·C}·v_c\left(t\right)+\frac{1}{\left(R_D+R_E\right)·C}·e\left(t\right)$$E8

As the handling of functions is more difficult than that of symbols, the second approach is suitable for small circuits and the first for larger circuits.

3. Qualitative behavior of dynamic circuits

3.1 First-order dynamic circuits

In the previous chapter, it has been shown that a linear circuit of the first order satisfies the equation:

$$\dot{x}+\frac{x}{\tau }=\frac{s\left(t\right)}{\tau }$$E9

where $x$ is the state variable of the circuit, $\tau$ is the time constant of the circuit, and $s\left(t\right)$ is the parameter of the equivalent independent source.

The solution of a first-order linear circuit with independent direct current sources consists of two terms: the solution of the homogeneous equation $x_v$ and a particular solution $x_p$:

$$x=x_v+x_p$$E10

The homogeneous equation is ${\dot{x}}_{\nu }+\frac{x_{\nu }}{\tau }=0$. This can be solved using the separation of variables method:

$${\dot{x}}_{\nu }=-\frac{x_{\nu }}{\tau }$$

$$\frac{d{x}_v}{\mathit{dt}}=-\frac{x_{\nu }}{\tau }$$

$$\frac{d{x}_v}{x_{\nu }}=-\frac{\mathit{dt}}{\tau }$$

$${\int }_{x_0}^{x_v}\frac{d{x}_v}{x_{\nu }}=-{\int }_{t_0}^t\frac{\mathit{dt}}{\tau }$$

$$\mathit{ln}\left(x_v\right)-\mathit{ln}\left(x_0\right)=-\frac{\left(t-t_0\right)}{\tau }$$

$$\mathit{ln}\left(x_v\right)=C_1-\frac{\left(t-t_0\right)}{\tau }$$

$$x_v=e^{C_1}\bullet e^{-\frac{\left(t-t_0\right)}{\tau }}=C_2\bullet e^{-\frac{\left(t-t_0\right)}{\tau }}$$

The particular solution is a constant, in the form of free term $x_p=C_3,$so:

$$x\left(t\right)=C_2\bullet e^{-\frac{\left(t-t_0\right)}{\tau }}+C_3$$

The constant $C_3$ is calculated by replacing $t=t_{\infty }$:

$$x\left(t_{\infty }\right)=C_2\bullet e^{-\infty }+C_3=0+C_3\mathrm{so}{C}_3=x\left(t_{\infty }\right)$$

The solution is determined only if the initial condition is known $x\left(t_0\right)$. Replacing $t=t_0$ we obtain:

$$x\left(t_0\right)=C_2\bullet e^0+x\left(t_{\infty }\right)=C_2+x\left(t_{\infty }\right)\mathrm{so}{C}_2=x\left(t_0\right)-x\left(t_{\infty }\right)\text{and results}:$$

$$x\left(t\right)=\left[x\left(t_0\right)-x\left(t_{\infty }\right)\right]·e^{-\frac{\left(t-t_0\right)}{\tau }}+x\left(t_{\infty }\right)$$E11

We distinguish two cases in which the behavior of the solution is different: $\tau >0$ and $\tau <0$. If $\tau <0$, when $t\to \infty$, $x\left(t\right)$ decreases exponentially over time and the solution tends to equilibrium (Figure 13). If $\tau <0$, when $t\to \infty$, $x\left(t\right)$ increases exponentially over time and the solution tends to an infinite value (Figure 14).

Figure 13.

The solution decreases exponentially.

Figure 14.

The solution increases exponentially.

3.2 Dynamic circuits of second order or greater

The qualitative behavior of the circuit is determined by the eigenvalues of the state matrix $A$. These can be calculated as roots of the equation:

$$\mathit{det}\left(A-\lambda \bullet I\right)=0$$E12

where $I$ is the unit matrix of the same order as the state matrix $A$.

First case: The state matrix $A$ has real and distinct eigenvalues. There are three possibilities:

1. the state matrix A has negative eigenvalues. In this case the solution moves towards the steady state when $t\to \infty$ (Figure 15a).

2. the state matrix A has null eigenvalues. In this case the solution is a constant, the steady state, when $t\to \infty$ as well as when $t\to -\infty$ (Figure 15b).

3. the state matrix A has positive eigenvalues. In this case the solution moves towards the steady state when $t\to -\infty$ (Figure 15c).

Figure 15.

Circuit response for real and distinct eigenvalues (a) stable response; (b) constant response; and (c) unstable response.

Second case: The state matrix $A$ has complex conjugated eigenvalues. In this case, there are also three possibilities:

1. state matrix A has eigenvalues with a negative real part. In this case, the solution contains damped harmonic components that move toward the steady state when $t\to \infty$ (Figure 16a).

2. The state matrix $A$ has eigenvalues with null real part. In this case, the solution contains undamped harmonic components when $t\to \infty$ as well as when $t\to -\infty$ (Figure 16b).

3. The state matrix $A$ has eigenvalues with a positive real. In this case, the solution contains damped harmonic components that go toward the steady state when $t\to -\infty$ (Figure 16c).

Figure 16.

Circuit response for complex conjugated eigenvalues (a) stable response; (b) constant response; and (c) unstable response.

It can be observed that circuits that have the state matrix $A$ with negative real or complex conjugated eigenvalues with the negative real part are stable. If at least one real eigenvalue is positive or a pair of complex conjugated eigenvalues has a positive real part, they are unstable.

For example, for the circuit in paragraph 2.2.1 (Figure 5), by replacing the numerical values for the circuit elements in the state matrix, the following eigenvalues are obtained:

$$\mathit{det}\left(A-\lambda \bullet I\right)=\mathit{det}\left(\begin{array}{ccc}-1-\lambda & 0& {10}^5\\ 0& -\lambda & {10}^{10}\\ -4\bullet {10}^7& -4\bullet {10}^7& -4\bullet {10}^6-\lambda \end{array}\right)=0$$

So:

$$-{\lambda }^3-4\bullet {10}^6\bullet {\lambda }^2-4\bullet {10}^{17}\bullet \lambda -4\bullet {10}^{17}=0\to \left\{\begin{array}{c}{\lambda }_1=-1\\ {\lambda }_2=-2\bullet {10}^6+6.3246\bullet {10}^8\bullet i\\ {\lambda }_3=-2\bullet {10}^6-6.3246\bullet {10}^8\bullet i\end{array}\right.$$

It is observed that all eigenvalues have a negative real part, so the circuit is stable.

4. Analytical solving of dynamic circuits equations

With the help of the Laplace transform, it is possible to construct a system of algebraic equations in the complex variable $s$ domain, which corresponds to a system of linear differential equations in the time domain. The use of algebraic equations instead of differential equations shows evident benefits in the study of electrical circuits.

Although section 2 shows how to write state equations in normal form $\dot{x}=f\left(x,t\right)$, this method is cumbersome for a human operator even for circuits with a small number of dynamic elements. For this reason, in the analysis of electrical circuits, it is preferred to apply the Laplace transform first to the equations that describe the operation of the circuit elements and then to write the circuit equations.

4.1 Circuit analysis with Laplace transform

The study of the time-varying regime can be performed if the initial conditions ($i_k\left(0_{-}\right)$ inductor currents and $v_k\left(0_{-}\right)$ capacitor voltages) at $t=0_{-}$ are known. It is considered that the independent sources are connected in the circuit at $t=0_{-}$. In this case, all voltages and all currents are original functions. For example, a direct current source (voltage or current) having $E=\mathit{ct}$ or $I_S=\mathit{ct}$, being connected at $t=0_{-}$ has $e\left(t\right)=E\bullet u\left(t\right)$ or $i_S\left(t\right)=I_S\bullet u\left(t\right).$ Thus, the quantities $e\left(t\right)$ and $i_S\left(t\right)$ become original functions due to the presence of the factor $u\left(t\right)$. For such a circuit there are Laplace images of voltages and currents: $I_k\left(s\right)=L\left\{i_k\left(t\right)\right\}$ and $V_k\left(s\right)=L\left\{v_k\left(t\right)\right\}$. Laplace image computation is called operational computation. Applying the linearity property of the Laplace transform, Kirchhoff's Laws ${\sum }_N{i}_k\left(t\right)=0$ and ${\sum }_B{v}_k\left(t\right)=0$ can be written in the complex variable $s$ domain: ${\sum }_N{I}_k\left(s\right)=0$ and ${\sum }_B{V}_k\left(s\right)=0$.

According to the linearity and derivation properties of the original function of the Laplace transform, the differential equations that express the connections between $v_k\left(t\right)$ and $i_k\left(t\right)$ for the circuit elements, correspond to algebraic equations that express the connections between $V_k\left(s\right)=L\left\{v_k\left(t\right)\right\}$ and $I_k\left(s\right)=L\left\{i_k\left(t\right)\right\}$ [11].

4.2 Equivalent operational circuits

A circuit in which currents and voltages are Laplace images is called the equivalent operational circuit. Below are the most common circuit elements and their Laplace images [11].

1. The ideal direct current voltage (or current) source

   

2. The ideal alternating current voltage (or current) source

   

3. The resistor

   The ideal resistor has the operation equation $v\left(t\right)=R·i\left(t\right)$ and if $V\left(s\right)=L\left\{v\left(t\right)\right\}$, $I\left(s\right)=L\left\{i\left(t\right)\right\}$ then $V\left(s\right)=R·I\left(s\right)$. The factor multiplied with $I\left(s\right)$ to get $V\left(s\right)$ is called the operational impedance $\left(Z\left(s\right)=\frac{V\left(s\right)}{I\left(s\right)}\right)$. For the ideal resistor $Z_R\left(s\right)=R$. The circuit corresponding to this relationship is shown below:

   

4. The inductor

   For the ideal inductor, according to the derivation property of the original function, the equation $v\left(t\right)=L·\frac{\mathit{di}\left(t\right)}{\mathit{dt}}$ transforms to:

   $$V\left(s\right)=L·\left[s·I\left(s\right)-i\left(0_{-}\right)\right]=s·L·I\left(s\right)-L·i\left(0_{-}\right)$$E13

   where $i\left(0_{-}\right)$ is the initial condition at the time $t=0_{-}$ for the current through the inductor. The equivalent operational circuit is:

   

5. The capacitor

   Similar to the inductor, according to the derivation property of the original function, for the ideal capacitor, the relationship $i\left(t\right)=C·\frac{\mathit{dv}\left(t\right)}{\mathit{dt}}$ transforms to:

   $$I\left(s\right)=C·\left[s·V\left(s\right)-v\left(0_{-}\right)\right]=\frac{V\left(s\right)}{\frac{1}{s·C}}-C·v\left(0_{-}\right)$$E14

   where $v\left(0_{-}\right)$ is the initial condition at the $t=0_{-}$ for the voltage drop on the capacitor. The equivalent operational circuit is:

   

4.3 Response computation of linear dynamic circuits

All theorems and analysis methods valid for direct or alternating current circuits are valid for the equivalent operational circuits: equivalent generator theorems, superposition theorem, the two-port representation theorems, Kirchhoff’s laws analysis, nodal analysis, mesh analysis, etc.

The Laplace transform analysis algorithm of a linear circuit in time-varying regime consists of:

1. determine the initial conditions for the inductors and capacitors in the circuit, $i_L\left(0_{-}\right)$ and $v_C\left(0_{-}\right)$;

2. build the circuit with operational sources and operational impedances using the operational equivalent circuits;

3. write the equations and solve for the unknowns $V_K\left(s\right)$ and $I_K\left(s\right)$;

4. determine the analytical solution, the original functions $v_k\left(t\right)$ and $i_k\left(t\right)$ as the inverse Laplace transform (using Heaviside's theorems).

If the circuit has more than four dynamic elements, determining the circuit response using analytical methods requires substantial effort for a human operator. In this case, a software product such as Mathematica [12] or Maple [13] capable of performing symbolic computations can be used.

To illustrate the above algorithm, consider the following example, the circuit in Figure 17 [14]. The current through the inductor $i_L\left(t\right)$ and the voltage at the terminals of the capacitor $v_C\left(t\right)$ in the transient mode that occurs after the switch K is closed at the time $t=0$ in the circuit in the figure below are required. We know $R=1\Omega$, $C=1/3F$, $L=3/2H$, $E=6V$, $i_L\left(0_{-}\right)=3A$ and $v_C\left(0_{-}\right)=3V$.

Figure 17.

Dynamic linear circuit in transient mode.

Using the equivalent operational circuits (Figure 18), we obtain:

Figure 18.

The equivalent operational circuit.

To make it easier to calculate the transient solution, we breakdown the expression into simple fractions:

$$I_L\left(s\right)=\frac{6}{s}+\frac{-4}{s+1}+\frac{1}{s+2}\text{and}{V}_C\left(s\right)=\frac{6}{s}+\frac{-6}{s+1}+\frac{3}{s+2}$$

So:

$$i_L\left(t\right)=L^{-1}\left\{I_L\left(s\right)\right\}=6-4·e^{-t}+e^{-2·t}\left[A\right]$$

and

$$v_C\left(t\right)=L^{-1}\left\{V_C\left(s\right)\right\}=6-6·e^{-t}+3·e^{-2·t}\left[V\right]$$

4.4 Response computation of the dynamic circuits with resistive nonlinearities

Solving systems of differential equations by a human operator is difficult even in the case of linear circuits and even more so in the case of nonlinear ones. To solve them, it is preferred to use programs capable of performing symbolic computations, such as Mathematica, Maple, etc. To solve the circuit state equation obtained by the second approach in paragraph 2.2:

$$\frac{d{v}_c\left(t\right)}{\mathit{dt}}=-\frac{R_D+R_E+R}{\left(R_D+R_E\right)·R·C}·v_c\left(t\right)+\frac{1}{\left(R_D+R_E\right)·C}·e\left(t\right)$$E15

where: $R_D=\left\{\begin{array}{c}0.1\mathrm{\Omega }\text{for}v>0\\ 10M\mathrm{\Omega }\text{for}v<0\end{array}\right.$, $e\left(t\right)=5\bullet \mathit{sin}\left(2\mathit{\pi f}\bullet \left(t+t_0\right)\right)\left[V\right]$ and $f=10\mathit{kHz}$, it is possible to proceed in this way. The above equation is solved separately for each linear region of the nonlinear characteristic for $R_D$. When $R_D$ value changes, the initial condition $v_c\left(t_0\right)=v_{c0}$ is considered as the value of the voltage $v_c\left(t\right)$ obtained at the end of the previous interval. These calculations were made in the Maple program for two periods of the source $e\left(t\right)$, with the following solution obtained by concatenating the results (Figure 19):

Figure 19.

The analytical solution (Maple) and the numerical solution (SPICE).

As a verification of this method, we can compare the results with those obtained by using the SPICE program [15], in the same time frame. It is observed that the analytical solution obtained with Maple is almost identical to numerical solution obtained with SPICE.

4.5 Transfer functions

Consider a circuit with constant parameters and a unique solution. We assume that we have only one independent source and initial conditions equal to zero. The transfer function from gate i to gate j is defined as the ratio between the output quantity on side j and the input quantity of side i (Figure 20). The output quantity can be a voltage or a current, and the input quantity can be the electromotive voltage of an independent voltage source or the current of an independent current source.

Figure 20.

Circuit with one input (i) and one output (j).

Generally, a circuit function also called a transfer function is:

$$H\left(s\right)=\frac{P\left(s\right)}{Q\left(s\right)}$$E16

The roots of $Q\left(s\right)$ are called the poles of$H\left(s\right),$ and the roots of $P\left(s\right)$ are called the zeros of$H\left(s\right)$. The dynamic behavior of the network depends upon the location of the poles and zeros on the network function curve [16]. In general, one can graphically deduce the magnitude and phase curve of any network function from the location of its poles and zeros. The poles of $H\left(s\right)$ are the circuit natural frequencies. Not all natural frequencies are the poles of any function of that circuit because certain common expressions that appear in both $Q\left(s\right)$ and $P\left(s\right)$ can disappear by simplification.

The poles of the transfer function determine the stability of the circuit, similarly to the eigenvalues of the state matrix presented in paragraph 3. If all the poles have a negative real part, the circuit is stable. If at least one pole has a positive real part, the circuit is unstable.

HSPICE uses the Muller method to calculate the roots of polynomials $P\left(s\right)$ and $Q\left(s\right)$. This method approximates the polynomial with a quadratic equation that fits through three points in the vicinity of a root. Successive iterations toward a particular root are obtained by finding the nearest root of a quadratic equation whose curve passes through the last three points [16]. Pole/zero analysis results are based on the circuit's DC operating point, so the operating point solution must be accurate [16].

For the circuit in paragraph 2.2.1, Figure 5, the poles in Table 1 are obtained using the HSPICE program.

Poles	Real	Imaginary
1	$-9.9999\bullet {10}^{-1}$	$0$
2	$-2\bullet {10}^6$	$+6.3246\bullet {10}^8$
3	$-2\bullet {10}^6$	$-6.3246\bullet {10}^8$

Table 1.

The poles obtained with HSPICE for the circuit in Figure 5.

It is observed that they are identical with the eigenvalues of the state matrix, calculated in paragraph 3.2.

5. Numerical approach of solving dynamic circuits equations

We aim to solve a system of state equations written in normal form $\dot{x}=f\left(x,t\right)$. Even for linear circuits whose equations have an analytical solution, the use of numerical methods is preferred because even for a second-order circuit, the analytical calculations are quite complicated to be performed efficiently by a human operator.

Automatic solving of analytical solutions requires considerable computational effort, as computers are designed to operate with numbers rather than symbols. For this purpose, specialized software such as Mathematica or Maple that performs analytical calculations can be used.

5.1 Numerical integration methods used in circuit simulation

Any numerical method starts from the initial condition $x\left(t_0\right)$ and determines successively:

$$x\left(t_0+h\right),x\left(t_0+2h\right),\dots ,$$E17

where h is the time step.

The simplest numerical integration methods used in circuit simulation programs are [17]:

5.1.1 The forward Euler method (FE)

Expanding $x\left(t_{k+1}\right)$ in Taylor series in the vicinity of the point $t_k$ we obtain:

$$x\left(t_{k+1}\right)=x\left(t_k\right)+{\left.\frac{\mathit{dx}}{\mathit{dt}}\right|}_{t_k}·\frac{\left(t_{k+1}-t_k\right)}{1!}+{\left.\frac{d^{2}x}{d{t}^2}\right|}_{t_k}·\frac{{\left(t_{k+1}-t_k\right)}^2}{2!}+\dots$$E18

Using the notation $x\left(t_k\right)=x_k$, neglecting the higher-order terms, and taking into account the fact that $\frac{\mathit{dx}}{\mathit{dt}}=\dot{x}=f\left(x,t\right)$, we obtain:

$$x_{k+1}=x_k+h·{\dot{x}}_k\text{or}{x}_{k+1}=x_k+h·f\left(x_k\right)$$E19

where $h=t_{k+1}-t_k$ (time step).

The graphical interpretation of this method of numerical integration is presented below in Figure 21.

Figure 21.

The forward Euler method (FE) graphical interpretation.

5.1.2 The backward Euler method (BE)

Similarly, expanding $x\left(t_k\right)$ in Taylor series in the vicinity of the point $t_{k+1}$ results:

$$x\left(t_k\right)=x\left(t_{k+1}\right)+{\left.\frac{\mathit{dx}}{\mathit{dt}}\right|}_{t_{k+1}}·\frac{\left(t_k-t_{k+1}\right)}{1!}+{\left.\frac{d^{2}x}{d{t}^2}\right|}_{t_{k+1}}·\frac{{\left(t_k-t_{k+1}\right)}^2}{2!}+\dots$$E20

With the above notations we obtain:

$$x_{k+1}=x_k+h·{\dot{x}}_{k+1}\text{or}{x}_{k+1}=x_k+h·f\left(x_{k+1}\right)$$E21

In Figure 22 is presented the graphical interpretation of this numerical integration method.

Figure 22.

The backward Euler method (BE) graphical interpretation.

5.1.3 The trapezoidal rule (TR)

From the graphic interpretation of the two numerical integration methods presented above, it is observed that for the same function, a method approximates the integral with a value smaller and the other with a bigger value. This drawback can be overcome by using the trapezoidal rule.

Adding and dividing by two the two formulas obtained for $x_{k+1}$, the one obtained with the forward Euler method and the one obtained with the backward Euler method, the formula for the trapezoidal rule is obtained:

$$\left.\begin{array}{c}{x}_{k+1}=x_k+h·{\dot{x}}_k\\ x_{k+1}=x_k+h·{\dot{x}}_{k+1}\end{array}\right\}\Rightarrow 2·x_{k+1}=2·x_k+h·\left({\dot{x}}_k+{\dot{x}}_{k+1}\right)$$EQ5

or

$$x_{k+1}=x_k+\frac{h}{2}·\left({\dot{x}}_k+{\dot{x}}_{k+1}\right)$$E22

The graphical interpretation of the trapezoidal rule is presented below, in Figure 23.

Figure 23.

The trapezoidal rule (TR) graphical interpretation.

It is easy to see that the forward Euler method, in which $x\left(t_{k+1}\right)$ is determined from $x\left(t_k\right)$, is an explicit method, while the backward Euler method and the trapezoidal rule are implicit methods.

The solutions obtained by numerical integration being approximate, the errors introduced at each step are calculated (local error and global error of the method). Local error at the time $t=t_{n+1}$ is ${\epsilon }_n=x_{\mathit{\text{exact}}}\left(t_{n+1}\right)-x_{\mathit{\text{approx}}}\left(t_{n+1}\right)$ where $x_{\mathit{\text{exact}}}\left(t_{n+1}\right)$ was calculated starting from the $x\left(t_n\right)$ approximate calculation. The total error at $t=t_{n+1}$ is ${\epsilon }_n=x_{\mathit{\text{exact}}}\left(t_{n+1}\right)-x_{\mathit{\text{approx}}}\left(t_{n+1}\right)$, where $x_{\mathit{\text{exact}}}\left(t_{n+1}\right)$ was calculated from the initial $x\left(0\right)$.

The numerical integration method for which the total error decreases as time passes is a stable method. A method that does not have this property is numerically unstable, even if the local error is small and decreases over time.

In a stable method, the size of the time step is limited only by the imposed local error, which depends on the studied problem. Working with very low values of h leads to a large number of time steps required to determine the solution that take up an unjustifiably high computation time.

5.1.4 The gear methods (G)

For circuits with eigenvalues that differ by a few orders of magnitude ("stiff") the numerical methods presented above do not give correct results. In this case, special gear methods are used. The relationships that define the second to sixth-order gear methods are:

Second-order gear: $x_{k+1}=\frac{3}{4}{x}_k-\frac{1}{3}{x}_{k-1}+h\frac{2}{3}·{\dot{x}}_{k+1}$

Third-order gear: $x_{k+1}=\frac{18}{11}{x}_k-\frac{9}{11}{x}_{k-1}+\frac{2}{11}·x_{k-2}+h\frac{6}{11}·{\dot{x}}_{k+1}$

Fourth-order gear 4: $x_{k+1}=\frac{48}{25}{x}_k-\frac{36}{25}{x}_{k-1}+\frac{16}{25}·x_{k-2}-\frac{3}{25}·x_{k-3}+h\frac{12}{25}·{\dot{x}}_{k+1}$

Fifth-order gear 5: $x_{k+1}=\frac{300}{137}{x}_k-\frac{300}{137}{x}_{k-1}+\frac{200}{137}·x_{k-2}-\frac{75}{137}·x_{k-3}+\frac{12}{137}·x_{k-4}+h\frac{60}{137}·{\dot{x}}_{k+1}$

Sixth-order gear 6: $x_{k+1}=\frac{360}{147}{x}_k-\frac{450}{147}{x}_{k-1}+\frac{400}{147}·x_{k-2}-\frac{225}{147}·x_{k-3}+\frac{72}{147}·x_{k-4}-\frac{10}{147}·x_{k-5}+h\frac{60}{137}·{\dot{x}}_{k+1}$

In the case of an explicit method, after choosing the time step, we start from the initial condition $x\left(t_0\right)$ and successively compute $x\left(t_0+h\right),x\left(t_0+2h\right),\cdots$ covering the entire time interval of interest. In the case of an implicit method, several iterations are made at each step. At the first iteration, an explicit method is used, and a predictor is obtained ${x_{k+1}}^{\left(0\right)}=x_k+h·f\left(x_k\right)$. The value obtained for the predictor is entered on the right-hand side of the equation of the backward method obtaining a new value for ${x_{k+1}}^{\left(1\right)}$ (in the left hand side), which at the next iteration is introduced again on the right-hand side and ${x_{k+1}}^{\left(2\right)}$ is obtained and so on until $\left|{x_{k+1}}^{\left(N-1\right)}-{x_{k+1}}^{\left(N\right)}\right|<\epsilon$ imposed. Values ${x_{k+1}}^{\left(1\right)},{x_{k+1}}^{\left(2\right)}\cdots$ are called correctors.

In current SPICE circuit simulation programs, the forward Euler, backward Euler, trapezoidal rule, and gear method of the second order are used. In circuit design, a small simulation time is very important and integration methods of degree greater than 2 are not used because it involves many more computations and the improvement of the solution is not considerable visible.

5.2 Resistive (companion) models for dynamic circuit elements

The numerical methods described in the previous paragraph require writing the state equations in normal form $\dot{x}=f\left(x,t\right)$. Writing the circuit equations in this form requires the elimination of some variables and the expanding of others starting from the circuit equations. This process involves the numerical solution of some systems of linear or nonlinear algebraic equations, operation that may be affected by significantly errors and involves a certain calculation effort. This process of integrating circuit equations is simplified by replacing the dynamic elements with so-called resistive (or companion) models. By doing so, the circuit response is determined by solving a linear or nonlinear resistive circuit at each time step.

Companion models for a linear capacitor or an inductor derive from the numerical integration method used in the previous paragraph.

Starting from the relations given by the most used numerical integration method and taking into account that the state variable for the capacitor is the voltage ($i=C·\frac{\mathit{dv}}{\mathit{dt}}$ or $i=C·\dot{v}$ so $\dot{v}=\frac{1}{C}·i$), we obtain the models in Table 2 for the capacitor.

BE	$i_{k+1}=\frac{C}{h}·v_{k+1}-\frac{C}{h}·v_k$	$R=\frac{h}{C}$	$I_S=\frac{C}{h}·v_k$
TR	$i_{k+1}=\frac{2\bullet C}{h}·v_{k+1}-\left(\frac{2\bullet C}{h}·v_k+i_k\right)$	$R=\frac{h}{2·C}$	$I_S=\frac{2\bullet C}{h}·v_k+i_k$
2^nd G	$i_{k+1}=\frac{3C}{2h}·v_{k+1}-\left(\frac{9C}{8h}·v_k-\frac{C}{2h}·v_{k-1}\right)$	$R=\frac{2h}{3C}$	$I_S=\frac{9C}{8h}·v_k-\frac{C}{2h}·v_{k-1}$

Table 2.

The BE, TR and second G companion models for the capacitor.

These equations are describing the following equivalent circuit, the capacitor companion model.

In a similar manner but considering that the state variable for the inductor is the current ($v=L·\frac{\mathit{di}}{\mathit{dt}}$ or $v=L·\dot{i}$ so $\dot{i}=\frac{1}{L}·v$), we obtain the models in Table 3 for the inductor.

BE	$i_{k+1}=\frac{h}{L}·v_{k+1}+i_k$	$R=\frac{L}{h}$	$I_S=i_k$
TR	$i_{k+1}=\frac{h}{2\bullet L}·v_{k+1}+\left(\frac{h}{2\bullet L}·v_k+i_k\right)$	$R=\frac{2\bullet L}{h}$	$I_S=\frac{h}{2\bullet L}·v_k+i_k$
2^nd G	$i_{k+1}=\frac{2h}{3L}·v_{k+1}+\left(\frac{3}{4}·i_k-\frac{1}{3}·i_{k-1}\right)$	$R=\frac{3L}{2h}$	$I_S=\frac{3}{4}·i_k-\frac{1}{3}·i_{k-1}$

Table 3.

The BE, TR, and second G companion models for the inductor.

These equations are describing the following equivalent circuit, the inductor companion model.

For example, using the companion models starting from the trapezoidal rule, the dynamic circuit in paragraph 2.2.1, Figure 5, becomes a resistive circuit but with different values at each time step for the resistors and sources in the dynamic element models (Figure 24).

Figure 24.

Dynamic elements replaced with companion models for circuit in Figure 5.

The advantage of this approach is that solving a resistive circuit is much easier than a dynamic circuit, for example, using the modified nodal analysis (MNA) method [18]:

$$\left[\begin{array}{ccc}\frac{1}{R_1}+\frac{1}{R_{L_1}}& -\frac{1}{R_{L_1}}& 0\\ -\frac{1}{R_{L_1}}& \frac{1}{R_{L_1}}+\frac{1}{R_{C_1}}& -\frac{1}{R_{C_1}}\\ 0& -\frac{1}{R_{C_1}}& \frac{1}{R}+\frac{1}{R_{C_1}}+\frac{1}{R_C}\end{array}\right]\bullet \left[\begin{array}{c}{V}_2\\ V_3\\ V_4\end{array}\right]=\left[\begin{array}{c}\frac{V_1}{R_1}-I_{L_1}\\ I_{L_1}+I_{C_1}\\ \frac{V_1}{R}+I_C-I_{C_1}\end{array}\right]$$

This method, using a variable integration time step, is implemented in all SPICE circuit simulators.

Similar to linear dynamic element models, companion models can be built for nonlinear dynamic elements. Next, we will determine the parameters for these models for the nonlinear capacitor and the nonlinear inductor for the trapezoidal rule. Models corresponding to other numerical integration methods can be determined in a similar way.

For a nonlinear capacitor where the electric charge is a polynomial dependence of the voltage, of the form, $q=C\bullet \left(v+c_1\bullet v^2+c_2\bullet v^3\right)$ the $i=\frac{\mathit{dq}}{\mathit{dt}}=\dot{q}$ the companion model can be obtained from the following.

In this case, the state variable for the capacitor is the electric charge:

$$q_{k+1}=q_k+\frac{h}{2}·\left({\dot{q}}_k+{\dot{q}}_{k+1}\right)=q_k+\frac{h}{2}·\left(i_k+i_{k+1}\right)$$

$$i_{k+1}=\frac{2}{h}·q_{k+1}-\frac{2}{h}·q_k-i_k$$

$$i_{k+1}=\frac{2}{h}\bullet C\bullet \left(v_{k+1}+c_1\bullet v_{k+1}^2+c_2\bullet v_{k+1}^3\right)-\left[\frac{2}{h}\bullet C\bullet \left(v_k+c_1\bullet v_k^2+c_2\bullet v_k^3\right)+i_k\right]$$E23

This equation describes the following equivalent circuit, the nonlinear capacitor companion model.

For a nonlinear inductor where the magnetic flux is a polynomial dependency of the current, of the form $\phi =L\bullet \left(i+c_1\bullet i^2+c_2\bullet i^3\right)$ and $v=\frac{\mathit{d\phi }}{\mathit{dt}}=\dot{\phi }$, the companion model can be obtained from the following.

In this case, the inductor state variable is the magnetic flux:

$${\phi }_{k+1}={\phi }_k+\frac{h}{2}·\left({\dot{\phi }}_k+{\dot{\phi }}_{k+1}\right)={\phi }_k+\frac{h}{2}·\left(v_k+v_{k+1}\right)$$

$$v_{k+1}=\frac{2}{h}·{\phi }_{k+1}-\frac{2}{h}·{\phi }_k-v_k$$

$$v_{k+1}=\frac{2}{h}\bullet L\bullet \left(i_{k+1}+c_1\bullet i_{k+1}^2+c_2\bullet i_{k+1}^3\right)-\left[\frac{2}{h}·L\bullet \left(i_k+c_1\bullet i_k^2+c_2\bullet i_k^3\right)+v_k\right]$$E24

This equation describes the following equivalent circuit, the nonlinear inductor companion model.

The resistive circuit that is solved at each time step is nonlinear in this case, even if the resistors in the circuit are linear. This circuit is solved by an iterative method, usually the Newton method.

6. Conclusions

All companion models contain a linear resistor with the resistance depending on the parameter of the dynamic element ($L$ or $C$) and the time step $h$, and an independent source whose parameter depends on the value of the state variable at the previous time.

Companion models of linear dynamic elements can also be built with voltage sources (the actual current source can be transformed into a voltage source). Current sources were preferred because they are suitable for the modified nodal analysis (MNA) method.

For a given time step $h$, starting from the given initial state of the dynamic elements, the circuit response is calculated at $t_0+h$ using a first-order numerical integration method. In this way, the analysis of a linear dynamic circuit can be done by solving a linear resistive circuit at each time step. Starting from $t_0+2\bullet h$, a second-order method can be used (such as the trapezoidal rule or the second-order Gear method). If the time step does not change, only the independent sources change in the resistive circuit, the values of the resistors remain the same. If the time step changes, the model resistance values must be recalculated.

The integration of the circuit equations is usually done with a variable time step $h$. As $h$ decreases, the resistance in the capacitor companion model decreases and the resistance in the inductor companion model increases. In the case of an $\mathit{LC}$ branch circuit, there are two resistors in parallel whose resistors are different by several orders of magnitude. Such a circuit cannot be solved correctly by using just simple precision. In some cases, even double-precision computations can lead to incorrect results.