## Derivation of System Equations

The electrical network of example 1 is shown above. The diode is modeled as a current source whose value is a nonlinear function of its terminal voltage ($$v_D$$) in parallel with two voltage dependent capacitors and a series resistance. Refer to the diode section of the device models for a detailed description. The elaborated network consists of $$l = 7$$ branches and $$k = 4$$ nodes.

The proper tree and the basic cuts of the graph corresponding to the elaborate network are shown above. Since capacitors $$C_E$$ and $$C_J$$ are members of the proper tree, the state vector may be written as:

$\mathbf{x} = \begin{bmatrix} v_{C_E} \\ v_{DCJ} \\ \end{bmatrix}$

The differential equation for $$v_{DCJ}$$ is obtained from the basic cut involving capacitor $$C_J$$. Adopting the convention that current flowing out of the surface containing node 2 is positive, then the equation may be written as:

$$i_{D} + i_{DCJ} + i_{DCT} - i_{RD} = i_{D} + C_J\frac{dv_{DCJ}}{dt} + i_{DCT} - i_{RD} = 0 \\ \label{eq:ex1_kcl1}$$

Equation \eqref{eq:ex1_kcl1} may be further simplified by noting that, since capacitors $$C_J$$ and $$C_T$$ are in parallel, $$v_{DCJ} = v_{DCT}$$.

$$i_{D} + C_J\frac{dv_{DCJ}}{dt} + C_T\frac{dv_{DCJ}}{dt} - i_{R_D} = 0 \\ \therefore (C_J + C_T)\frac{dv_{DCJ}}{dt} = i_{R_D} - i_{D} \\ \label{eq:ex1_kcl2}$$

The differential equation for $$v_{C_E}$$ is obtained from the basic cut involving capacitor $$C_E$$. Adopting the convention that current flowing out of the surface containing nodes 0 and 1 is positive, then the equation may be written as:

$$i_{C_E} + i_{R} + i_{R_D} = C_E\frac{dv_{C_E}}{dt} + i_{R} - i_{R_D} = 0 \\ \label{eq:ex1_kcl3}$$

Differential equations \eqref{eq:ex1_kcl2} and \eqref{eq:ex1_kcl3} contain the terms $$i_{R_D}$$ and $$i_{R}$$ which are neither state variables or system inputs. Resistors $$R_D$$ and $$R$$ are both members of the proper tree and therefore equations relating these terms to the state variables and system inputs may be derived by writing the KVL equations relating to these branches:

$$v_{R} - v_{C_E} = 0 \\ \label{eq:ex1_kvl1}$$ $$v_{R_D} - v_{in} + v_{C_E} + v_{DCJ} = 0 \\ \label{eq:ex1_kvl2}$$

Equations \eqref{eq:ex1_kcl2}, \eqref{eq:ex1_kcl3}, \eqref{eq:ex1_kvl1} and \eqref{eq:ex1_kvl2} may be compactly written in matrix notation:

$\begin{bmatrix} 1 & 0\\ 0 & 1 \\ \end{bmatrix} \begin{bmatrix} v_{R_D} \\ v_{R} \\ \end{bmatrix} = \begin{bmatrix} v_{in} - v_{C_E} - v_{DCJ} \\ v_{C_E} \\ \end{bmatrix} \label{eq:ex1_sys1}$ $$i_{R_D} = \frac{v_{R_D}}{R_D} \\ i_{R} = \frac{v_{R}}{R} \\ \label{eq:ex1_sys2}$$ $\begin{bmatrix} (C_J + C_T) & 0\\ 0 & C_E \\ \end{bmatrix} \begin{bmatrix} \frac{dv_{DCJ}}{dt} \\ \frac{dv_{C_E}}{dt} \\ \end{bmatrix} = \begin{bmatrix} i_{R_D} - i_{D} \\ i_{R_D} - i_{R} \\ \end{bmatrix} \label{eq:ex1_sys3}$

The diode current $$i_{D}$$ is a function of the diode terminal voltage $$v_{D}$$. From the basic loop involving this Voltage Controlled Current Source (VCCS), it can easily be shown that $$v_{D} = v_{DCJ}$$ and therefore $$v_{DCJ}$$ may be substituted for all occurrences of $$v_{D}$$ in the device characteristic equations. Equations \eqref{eq:ex1_sys1}, \eqref{eq:ex1_sys2} and \eqref{eq:ex1_sys3} therefore represent a complete system description.

## Solving the System of Equations

Octave's function lsode can be used to solve Ordinary Differential Equations (ODEs) of the form:

$$\frac {d\mathbf{x}}{dt} = f(\mathbf{x}, t) \\ \label{eq:oct1}$$

The ODE solver requires a function to compute the vector of derivatives of equation \eqref{eq:oct1}. For this example, the derivatives are evaluated by the function dypc_example1.

	  function dy = dypc_example1(y,t)

dy = zeros(2,1);

% external capacitance and reistance values.
R = 10e3;
CE = 0.01e-6;

% diode model parameters 1N4148
Is = 2.52e-9;
N = 1.752;
BV = 100.1;
Ibv = 10;
Gmin = 1e-12;
Cjo = 4e-12;
Vj = 0.3905;
m = 0.4;
tt = 20e-9;
FC = 0.5;
Rs = 0.568;

% diode p-n junction temperature
temp = 27;

% calculate the input voltage at time t.
A=0.85;
u = 0.3;
fc = 1000e3;
fm = 5000;
V_VIN = A*(1+ u*cos(2*pi*fm*t))*cos(2*pi*fc*t);

% Variable y holds the current state vector for time t.
% y(1) = VDCJ
% y(2) = VCE
V_DCJ = y(1);
V_CE = y(2);

% declare the Bm and Vm matrices and solve for Zm.
Bm = [1 0; 0 1];

Vm = [V_VIN-V_DCJ-V_CE; V_CE];

Zm = Bm \ Vm;
V_RD = Zm(1);
V_R = Zm(2);

% calculate the diode current, series resistance and junction capacitance.
d = diode(V_DCJ, Is, N, BV, Ibv, Gmin, Cjo, Vj, m, tt, FC, Rs, temp);
RD = d(1);
CJ = d(2);
CT = d(3);
i_ID = d(4);

% calculate resistor currents.
i_RD = V_RD / RD;
i_R = V_R / R;

% declare the Qm and Tm matrices and solve for Rm.
Qm = [CJ+CT 0; 0 CE];

Tm = [i_RD-i_ID; i_RD-i_R];

Rm = Qm \ Tm;

% assign the function outputs.
dy(1) = Rm(1);
dy(2) = Rm(2);


The function lsode uses a numerical integration algorithm to solve the system of equations. It must be linked to the dypc_example1 function and invoked with suitable arguments. For this example, the following code sequence is used to solve the system of equations for time $$t=0$$ to $$t=600 \mu s$$.

	  clear all;

abTol = 1e-6;

t = linspace (0, 600e-6, 10e3);
x0 = [0 0];

lsode_options("absolute tolerance",[abTol abTol]);
lsode_options("relative tolerance", 1e-6);
lsode_options("integration method", "stiff");

tic
Y = lsode("dypc_example1", x0, t);
toc