## Final notes for ECE1254, Modelling of Multiphysics Systems

I’ve now finished my first grad course, Modelling of Multiphysics Systems, taught by Prof Piero Triverio.

I’ve posted notes for lectures and other material as I was taking the course, but now have an aggregated set of notes for the whole course posted.
This is now updated with all my notes from the lectures, solved problems, additional notes on auxillary topics I wanted to explore (like SVD), plus the notes from the Harmonic Balance report that Mike and I will be presenting in January.

This version of my notes also includes all the matlab figures regenerating using http://www.mathworks.com/matlabcentral/fileexchange/23629-export-fig, which allows a save-as pdf, which rescales much better than Matlab saveas() png’s when embedded in latex.  I’m not sure if that’s the best way to include Matlab figures in latex, but they are at least not fuzzy looking now.

All in all, I’m pretty pleased with my notes for this course.  They are a lot more readable than any of the ones I’ve done for the physics undergrad courses I was taking (http://peeterjoot.com/writing/).  While there was quite a lot covered in this course, the material really only requires an introductory circuits course and some basic math (linear algebra and intro calculus), so is pretty accessible.

This was a fun course.  I recall, back in ancient times when I was a first year student, being unsatisfied with all the ad-hoc strategies we used to solve circuits problems.  This finally answers the questions of how to tackle things more systematically.

Here’s the contents outline for these notes:

Preface
Lecture notes
1 nodal analysis
1.1 In slides
1.2 Mechanical structures example
1.3 Assembling system equations automatically. Node/branch method
1.4 Nodal Analysis
1.5 Modified nodal analysis (MNA)
2 solving large systems
2.1 Gaussian elimination
2.2 LU decomposition
2.3 Problems
3 numerical errors and conditioning
3.1 Strict diagonal dominance
3.2 Exploring uniqueness and existence
3.3 Perturbation and norms
3.4 Matrix norm
4 singular value decomposition, and conditioning number
4.1 Singular value decomposition
4.2 Conditioning number
5 sparse factorization
5.1 Fill ins
5.2 Markowitz product
5.3 Markowitz reordering
5.4 Graph representation
6.1 Summary of factorization costs
6.2 Iterative methods
6.4 Recap: Summary of Gradient method
6.6 Full Algorithm
6.7 Order analysis
6.9 Gershgorin circle theorem
6.10 Preconditioning
6.11 Symmetric preconditioning
6.13 Problems
7 solution of nonlinear systems
7.1 Nonlinear systems
7.2 Richardson and Linear Convergence
7.3 Newton’s method
7.4 Solution of N nonlinear equations in N unknowns
7.5 Multivariable Newton’s iteration
7.6 Automatic assembly of equations for nonlinear system
7.7 Damped Newton’s method
7.8 Continuation parameters
7.9 Singular Jacobians
7.10 Struts and Joints, Node branch formulation
7.11 Problems
8 time dependent systems
8.1 Assembling equations automatically for dynamical systems
8.2 Numerical solution of differential equations
8.3 Forward Euler method
8.4 Backward Euler method
8.5 Trapezoidal rule (TR)
8.6 Nonlinear differential equations
8.7 Analysis, accuracy and stability (Dt ! 0)
8.8 Residual for LMS methods
8.9 Global error estimate
8.10 Stability
8.11 Stability (continued)
8.12 Problems
9 model order reduction
9.1 Model order reduction
9.2 Moment matching
9.3 Model order reduction (cont).
9.4 Moment matching
9.5 Truncated Balanced Realization (1000 ft overview)
9.6 Problems
Final report
10 harmonic balance
10.1 Abstract
10.2 Introduction
10.2.1 Modifications to the netlist syntax
10.3 Background
10.3.1 Discrete Fourier Transform
10.3.2 Harmonic Balance equations
10.3.3 Frequency domain representation of MNA equations
10.3.4 Example. RC circuit with a diode.
10.3.5 Jacobian
10.3.6 Newton’s method solution
10.3.7 Alternative handling of the non-linear currents and Jacobians
10.4 Results
10.4.1 Low pass filter
10.4.2 Half wave rectifier
10.4.3 AC to DC conversion
10.4.4 Bridge rectifier
10.4.5 Cpu time and error vs N
10.4.6 Taylor series non-linearities
10.4.7 Stiff systems
10.5 Conclusion
10.6 Appendices
10.6.1 Discrete Fourier Transform inversion
Appendices
a singular value decomposition
b basic theorems and definitions
c norton equivalents
d stability of discretized linear differential equations
e laplace transform refresher
f discrete fourier transform
g harmonic balance, rough notes
g.1 Block matrix form, with physical parameter ordering
g.2 Block matrix form, with frequency ordering
g.3 Representing the linear sources
g.4 Representing non-linear sources
g.5 Newton’s method
g.6 A matrix formulation of Harmonic Balance non-linear currents
h matlab notebooks
i mathematica notebooks
Index
Bibliography

## Thoughts on the “new” C++ style cast operators.

I happen to maintain the DB2 coding standards, which are mostly concerned with portability, and not style. I’ve joked that I was given that job since I had “broken the build” more than anybody else, so was most qualified to let others know how not to do so.

In our coding standards we have a prohibition against the use of exceptions. This is a historical restriction because we’ve built with compilation flags like -qnoeh (no exception handling) on some platforms to get a bit of additional performance. These days the compilers do much better at not degrading performance when exception handling is allowed and not used, but since our performance folks will sell their kids for a 1% improvement, we’ve kept using flags like this and the associated restriction. Components that must (or want) to use exception handling must “firewall” any exceptions, not letting them get thrown to external code (and also explicitly enable exceptions for their code).

We had a note in the coding standards not to use RTTI (Run Time Type Identification), because exceptions are required. That was a confusing and incomplete statement to include in our standards. It was interpreted by one developer as meaning that none of:

dynamic_cast<>()
reinterpret_cast<>()
static_cast<>()
const_cast<>()


were allowed. However, only the dynamic_cast is a RTTI operation, and only the dynamic_cast will throw an exception when the cast doesn’t match the underlying type.

I’ve now fixed up our coding standards. It now references dynamic_cast instead of RTTI.

The DB2 code is still very C’ish, compiled with a C++ compiler. I’d say the bulk of the casts in our code are old style C casts, and that most of our developers (including myself) don’t even know when to use the “new” cast operations. Here’s some thoughts on these:

• I’ve seen a fair amount of const_cast use in our code, and I know personally how this can be very useful. An example:


volatile int x ;
void foo( int * y ) ;

foo( const_cast(&x) ) ; // compilation error
foo( const_cast(&x) ) ; // allowed.  Just strips the volatile (or const) off the pointer type.

• A second nice use of const_cast<>() is to enforce type checking in macros. If the macro parameter is supposed to be a pointer to type T, you can enforce that by using a const_cast<T*>. This assumes that you don’t actually care about the const-ness of the pointer, and will force a compile error if the macro is used with any other type.
• In general I don’t think it’s a bad idea to use the new cast operators, since they represent a hierarchy of weaker than C style casts. You can also search for cast operations by name in a given file if they are used, which is much harder to do with a C style cast.
• I’m not sure how much use of static_cast<> and reinterpret_cast<> we have in the code.

Declaring my own stupidity, without looking them up, I didn’t personally know how to use the “new style” static_cast and reinterpret_cast operations correctly. I use a C cast by habit unless I’m trying to strip const or volatile attributes (or enforce a type in a macro).

If DB2 coders start using these, it will likely confuse old guys like me for a while, but I figure I can learn about these.

As a step in this direction, I see some helpful looking info in the C++ reference page on explicit_cast

My take on reference page is roughly: If there is a place to use a C cast, then you can use:

1. const_cast if it compiles. If not you can use:
2. static_cast. If that doesn’t compile you can use:
3. reinterpret_cast. If that doesn’t compile you can use:
4. (In code where exceptions are allowed:) dynamic_cast.
5. If that doesn’t work or is not allowed you can use:
A C style cast.

## New version of phy450 (Relativistic Electrodynamics) notes now posted.

I’d taken Professor Poppitz’s “Relativistic Electrodynamics” course in 2011, and wrote up my notes and problem set solutions in latex at that point.

That was only the second course that I tried this in, and was much less structured than any of the subsequent class notes collections that I produced later.  I now put problems with the chapter material using the latex exercise environment, and have now retrofitted the use of that package into my old phy450 notes.

I’ve also moved some of the tutorial content into the problem sections, and the rest into the appendix.  The result is a much more streamlined layout than was the case when I took the course originally.  There are still no figures, nor is there any index.  Only my later class notes collections have those and it’s a fair amount of work to retrofit that.  The whole thing still needs a complete review (and probably a rewrite).

## A matrix formulation of the Harmonic Balance method non-linear currents

Because it was simple, a coordinate expansion of the Jacobian of the non-linear currents was good to get a feeling for the structure of the equations. However, a Jacobian of that form is impossibly slow to compute for larger $$N$$. It seems plausible that eliminating the coordinate expansion, expressing both the currrent and the Jacobian directly in terms of the Harmonic Balance unknowns vector $$\BV$$, would lead to a simpler set of equations that could be implemented in a computationally more effective way. To aid in this discovery, consider the simple RC load diode circuit of fig. 1. It’s not too hard to start from scratch with the time domain nodal equations for this circuit, which are

fig. 1. Simple diode and resistor circuit

1. $$0 = i_s – i_d$$
2. $$Z v^{(2)} + C dv^{(2)}/dt = i_d$$
3. $$i_d = I_0 \lr{ e^{(v^{(1)} – v^{(2)})/V_T} – 1}$$

To setup for matrix form, let

\label{eqn:diodeRLCSample:1240}
\Bv(t) =
\begin{bmatrix}
v^{(1)}(t) \\
v^{(2)}(t) \\
\end{bmatrix}

\label{eqn:diodeRLCSample:1140}
\BG =
\begin{bmatrix}
0 & 0 \\
0 & Z \\
\end{bmatrix}

\label{eqn:diodeRLCSample:1160}
\BC =
\begin{bmatrix}
0 & 0 \\
0 & C \\
\end{bmatrix}

\label{eqn:diodeRLCSample:1180}
\Bd =
\begin{bmatrix}
1 \\
-1
\end{bmatrix}

\label{eqn:diodeRLCSample:1200}
\Bb =
\begin{bmatrix}
1 \\
0
\end{bmatrix},

so that the time domain equations can be written as

\label{eqn:diodeRLCSample:1220}
\BG \Bv(t)
+ \BC \Bv'(t)
=
\Bb i_s(t)
+
I_0
\Bd
\lr{
e^{ (v^{(1)}(t) – v^{(2)}(t))/V_T} – 1
}
=
\begin{bmatrix}
\Bb & -I_0 \Bd
\end{bmatrix}
\begin{bmatrix}
i_s(t) \\
1
\end{bmatrix}
+
I_0 \Bd
e^{ (v^{(1)}(t) – v^{(2)}(t))/V_T}.

Harmonic Balance is essentially the assumption that the input and outputs are assumed to be a bandwidth limited periodic signal, and the non-linear components can be approximated by the same

\label{eqn:diodeRLCSample:1260}
i_s(t) = \sum_{n=-N}^N I^{(s)}_n e^{ j \omega_0 n t },

\label{eqn:diodeRLCSample:1280}
v^{(k)}(t) =
\sum_{n=-N}^N V^{(k)}_n e^{ j \omega_0 n t },

\label{eqn:diodeRLCSample:1300}
\epsilon(t) =
e^{ (v^{(1)}(t) – v^{(2)}(t))/V_T}
\simeq
\sum_{n=-N}^N E_n e^{ j \omega_0 n t },

The approximation in \ref{eqn:diodeRLCSample:1300} is an equality only at the Nykvist sampling times $$t_k = T k/(2 N + 1)$$. The Fourier series provides a periodic extension to other times that will approximate the underlying periodic non-linear relation.

With all the time dependence locked into the exponentials, the derivatives are really easy to calculate

\label{eqn:diodeRLCSample:1281}
\frac{d}{dt} v^{(k)}(t) =
\sum_{n=-N}^N j \omega_0 n V^{(k)}_n e^{ j \omega_0 n t }.

Inserting all of these into \ref{eqn:diodeRLCSample:1220} gives

\label{eqn:diodeRLCSample:1320}
\sum_{n=-N}^N e^{ j \omega_0 n t} \lr{ \BG + j \omega_0 n \BC }
\begin{bmatrix}
V^{(1)}_n \\
V^{(2)}_n \\
\end{bmatrix}
=
\sum_{n=-N}^N e^{ j \omega_0 n t}
\lr{
-I_0 \Bd \delta_{n 0}
+
\Bb I^{(s)}_n
+ I_0 \Bd E_n
}.

The periodic assumption requires equality for each $$e^{j \omega_0 n t}$$, or

\label{eqn:diodeRLCSample:1340}
\lr{ \BG + j \omega_0 n \BC }
\begin{bmatrix}
V^{(1)}_n \\
V^{(2)}_n \\
\end{bmatrix}
=
-I_0 \Bd \delta_{n 0}
+
\Bb I^{(s)}_n
+ I_0 \Bd E_n.

For illustration, consider the $$N = 1$$ case, where the block matrix form is

\label{eqn:diodeRLCSample:1360}
\begin{bmatrix}
\BG + j \omega_0 (-1) \BC & 0 & 0 \\
0 & \BG + j \omega_0 (0) \BC & 0 \\
0 & 0 & \BG + j \omega_0 (1) \BC
\end{bmatrix}
\begin{bmatrix}
\begin{bmatrix}
V^{(1)}_{-1} \\
V^{(2)}_{-1} \\
\end{bmatrix} \\
\begin{bmatrix}
V^{(1)}_{0} \\
V^{(2)}_{0} \\
\end{bmatrix} \\
\begin{bmatrix}
V^{(1)}_{1} \\
V^{(2)}_{1} \\
\end{bmatrix}
\end{bmatrix}
=
\begin{bmatrix}
\Bb I^{(s)}_{-1} \\
\Bb I^{(s)}_{0} – I_0 \Bd \\
\Bb I^{(s)}_{1} \\
\end{bmatrix}
+
I_0
\begin{bmatrix}
\Bd E_{-1} \\
\Bd E_{0} \\
\Bd E_{1} \\
\end{bmatrix}.

The structure of this equation is

\label{eqn:diodeRLCSample:1380}
\BY \BV = \BI + \mathcal{I}(\BV),

The non-linear current $$\mathcal{I}(\BV)$$ needs to be examined further. How much of this can be precomputed, and what is the simplest way to compute the Jacobian? With

\label{eqn:diodeRLCSample:1400}
\BE =
\begin{bmatrix}
E_{-1} \\
E_{0} \\
E_{1} \\
\Bepsilon =
\begin{bmatrix}
\epsilon_{-1} \\
\epsilon_{0} \\
\epsilon_{1} \\
\end{bmatrix},

the non-linear current is

\label{eqn:diodeRLCSample:1420}
\mathcal{I} =
I_0
\begin{bmatrix}
\Bd E_{-1} \\
\Bd E_{0} \\
\Bd E_{1} \\
\end{bmatrix}
=
I_0
\begin{bmatrix}
\Bd \begin{bmatrix} 1 & 0 & 0 \end{bmatrix} \BE \\
\Bd \begin{bmatrix} 0 & 1 & 0 \end{bmatrix} \BE \\
\Bd \begin{bmatrix} 0 & 0 & 1 \end{bmatrix} \BE
\end{bmatrix}
=
I_0
\begin{bmatrix}
\Bd & 0 & 0 \\
0 & \Bd & 0 \\
0 & 0 & \Bd
\end{bmatrix}
\BF^{-1} \Bepsilon

In the last step $$\BE = \BF^{-1} \Bepsilon$$ has been factored out (in its inverse Fourier form). With

\label{eqn:diodeRLCSample:1480}
\BD =
\begin{bmatrix}
\Bd & 0 & 0 \\
0 & \Bd & 0 \\
0 & 0 & \Bd \\
\end{bmatrix},

the current is

\label{eqn:diodeRLCSample:1540}
\boxed{
\mathcal{I}(\BV) =
I_0 \BD \BF^{-1} \Bepsilon(\BV).
}

The next step is finding an appropriate form for $$\Bepsilon$$

\label{eqn:diodeRLCSample:1440}
\begin{aligned}
\Bepsilon &=
\begin{bmatrix}
\epsilon(t_{-1}) \\
\epsilon(t_{0}) \\
\epsilon(t_{1}) \\
\end{bmatrix} \\
&=
\begin{bmatrix}
\exp\lr{ \lr{ v^{(1)}_{-1} – v^{(2)}_{-1} }/V_T } \\
\exp\lr{ \lr{ v^{(1)}_{0} – v^{(2)}_{0} }/V_T } \\
\exp\lr{ \lr{ v^{(1)}_{1} – v^{(2)}_{1} }/V_T }
\end{bmatrix} \\
&=
\begin{bmatrix}
\exp\lr{
\begin{bmatrix}
1 & 0 & 0
\end{bmatrix}
\lr{ \Bv^{(1)} – \Bv^{(2)} }/V_T } \\
\exp\lr{
\begin{bmatrix}
0 & 1 & 0
\end{bmatrix}
\lr{ \Bv^{(1)} – \Bv^{(2)} }/V_T } \\
\exp\lr{
\begin{bmatrix}
0 & 0 & 1
\end{bmatrix}
\lr{ \Bv^{(1)} – \Bv^{(2)} }/V_T } \\
\end{bmatrix} \\
&=
\begin{bmatrix}
\exp\lr{
\begin{bmatrix}
1 & 0 & 0
\end{bmatrix}
\BF
\lr{ \BV^{(1)} – \BV^{(2)} }/V_T } \\
\exp\lr{
\begin{bmatrix}
0 & 1 & 0
\end{bmatrix}
\BF
\lr{ \BV^{(1)} – \BV^{(2)} }/V_T } \\
\exp\lr{
\begin{bmatrix}
0 & 0 & 1
\end{bmatrix}
\BF
\lr{ \BV^{(1)} – \BV^{(2)} }/V_T } \\
\end{bmatrix}.
\end{aligned}

It would be nice to have the difference of frequency domain vectors expressed in terms of $$\BV$$, which can be done with a bit of rearrangement

\label{eqn:diodeRLCSample:1460}
\begin{aligned}
\BV^{(1)} – \BV^{(2)}
&=
\begin{bmatrix}
V^{(1)}_{-1} – V^{(2)}_{-1} \\
V^{(1)}_{0} – V^{(2)}_{0} \\
V^{(1)}_{1} – V^{(2)}_{1} \\
\end{bmatrix} \\
&=
\begin{bmatrix}
1 & -1 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & -1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & -1 \\
\end{bmatrix}
\begin{bmatrix}
V_{-1}^{(1)} \\
V_{-1}^{(2)} \\
V_{0}^{(1)} \\
V_{0}^{(2)} \\
V_{1}^{(1)} \\
V_{1}^{(2)} \\
\end{bmatrix} \\
&=
\begin{bmatrix}
\Bd^\T & 0 & 0 \\
0 & \Bd^\T & 0 \\
0 & 0 & \Bd^\T \\
\end{bmatrix}
\BV \\
&= \BD^\T \BV,
\end{aligned}

\label{eqn:diodeRLCSample:1520}
\BH
=
\BF \BD^\T /V_T
=
\begin{bmatrix}
\Bh_1^\T \\
\Bh_2^\T \\
\Bh_3^\T
\end{bmatrix},

which allows the non-linear current to can now be completely expressed in terms of $$\BV$$.

\label{eqn:diodeRLCSample:1560}
\boxed{
\Bepsilon(\BV)
=
\begin{bmatrix}
e^{\Bh_1^\T \BV} \\
e^{\Bh_2^\T \BV} \\
e^{\Bh_3^\T \BV} \\
\end{bmatrix}.
}

## Jacobian

With a compact matrix representation of the non-linear current, attention can now be turned to the Jacobian of the non-linear current. Let $$\BA = I_0 \BD \BF^{-1} = [ a_{ij} ]_{ij}$$, the current (with summation implied) is

\label{eqn:diodeRLCSample:1580}
\mathcal{I} =
\begin{bmatrix}
a_{ik} \epsilon_k,
\end{bmatrix}

with coordinates

\label{eqn:diodeRLCSample:1600}
\mathcal{I}_i = a_{ik} \epsilon_k = a_{ik} \exp\lr{ \Bh_k^\T \BV }.

so the Jacobian components are

\label{eqn:diodeRLCSample:1620}
[\BJ^{\mathcal{I}}]_{ij}
=
a_{ik} \epsilon_k = a_{ik}
\PD{V_j}{}
\exp\lr{ \Bh_k^\T \BV }
=
a_{ik}
h_{kj}
\exp\lr{ \Bh_k^\T \BV }.

Factoring out $$\BU = [h_{ij} \exp\lr{ \Bh_i^\T \BV }]_{ij}$$,

\label{eqn:diodeRLCSample:1640}
\BJ^{\mathcal{I}}
= \BA \BU
=
\BA
\begin{bmatrix}
\begin{bmatrix} h_{11} & h_{12} & \cdots h_{1, R(2 N + 1)}\end{bmatrix} \exp\lr{ \Bh_1^\T \BV } \\
\begin{bmatrix} h_{21} & h_{22} & \cdots h_{2, R(2 N + 1)}\end{bmatrix} \exp\lr{ \Bh_2^\T \BV } \\
\begin{bmatrix} h_{31} & h_{32} & \cdots h_{3, R(2 N + 1)}\end{bmatrix} \exp\lr{ \Bh_3^\T \BV } \\
\end{bmatrix}
=
\BA
\begin{bmatrix}
\Bh_1^\T \exp\lr{ \Bh_1^\T \BV } \\
\Bh_2^\T \exp\lr{ \Bh_2^\T \BV } \\
\Bh_3^\T \exp\lr{ \Bh_3^\T \BV } \\
\end{bmatrix}.

A quick sanity check of dimensions seems worthwhile, and shows that all is well

• $$\BA$$ : $$R(2 N + 1) \times 2 N + 1$$
• $$\BU$$ : $$2 N + 1 \times R(2 N + 1)$$
• $$\BJ^{\mathcal{I}}$$ : $$R(2 N + 1) \times R(2 N + 1)$$

The Jacobian of the non-linear current is now completely determined

\label{eqn:diodeRLCSample:1660}
\boxed{
\BJ^{\mathcal{I}}( \BV ) =
I_0 \BD \BF^{-1}
\begin{bmatrix}
\Bh_1^\T \exp\lr{ \Bh_1^\T \BV } \\
\Bh_2^\T \exp\lr{ \Bh_2^\T \BV } \\
\Bh_3^\T \exp\lr{ \Bh_3^\T \BV } \\
\end{bmatrix}.
}

## Newton’s method solution

All the pieces required for a Newton’s method solution are now in place. The goal is to find a value of $$\BV$$ that provides the zero

\label{eqn:diodeRLCSample:1680}
f(\BV) = \BY \BV – \BI – \mathcal{I}(\BV).

Expansion to first order around an initial guess $$\BV^0$$, gives

\label{eqn:diodeRLCSample:1700}
f( \BV^0 + \Delta \BV ) = f(\BV^0) + \BJ(\BV^0) \Delta \BV \approx 0,

where the full Jacobian of $$f(\BV)$$ is

\label{eqn:diodeRLCSample:1720}
\BJ(\BV) = \BY – \BJ^{\mathcal{I}}(\BV).

The Newton’s method refinement of the initial guess follows by inversion

\label{eqn:diodeRLCSample:1740}
\Delta \BV = -\lr{ \BY – \BJ^{\mathcal{I}}(\BV^0) }^{-1} f(\BV^0).

## Development of Harmonic balance equations (considering a sample diode RLC circuit)

Previously, the time domain MNA equations and first steps at producing the Harmonic Balance equations were performed. That was a frequency domain analysis with an assumed Fourier solution associated with discrete time sampling.

The next goal is to put this in block matrix form. First introducing discrete time sampling vectors

\label{eqn:diodeRLCSample:100}
\Bv_a =
\begin{bmatrix}
v_a(t_{-N}) \\
v_a(t_{1-N}) \\
\vdots \\
v_a(t_{N-1}) \\
v_a(t_{N}) \\
\Bu_a =
\begin{bmatrix}
u_b(t_{-N}) \\
u_b(t_{1-N}) \\
\vdots \\
u_b(t_{N-1}) \\
u_b(t_{N}) \\
\Bw_a =
\begin{bmatrix}
w_c(t_{-N}) \\
w_c(t_{1-N}) \\
\vdots \\
w_c(t_{N-1}) \\
w_c(t_{N}) \\
\end{bmatrix},

and Fourier component vectors

\label{eqn:diodeRLCSample:120}
\BV_a =
\begin{bmatrix}
V^{(a)}_{-N} \\
V^{(a)}_{1-N} \\
\vdots \\
V^{(a)}_{N-1} \\
V^{(a)}_{N} \\
\BU_b =
\begin{bmatrix}
U^{(b)}_{-N} \\
U^{(b)}_{1-N} \\
\vdots \\
U^{(b)}_{N-1} \\
U^{(b)}_{N} \\
\BW_c =
\begin{bmatrix}
W^{(c)}_{-N} \\
W^{(c)}_{1-N} \\
\vdots \\
W^{(c)}_{N-1} \\
W^{(c)}_{N} \\
\end{bmatrix}.

With $$\alpha = e^{ 2 \pi j /(2 N + 1) }$$, and

\label{eqn:diodeRLCSample:140}
\BF =
\begin{bmatrix}
\alpha^{ N N } & \alpha^{ \lr{N-1} N } & \cdots & 1 & \cdots & \alpha^{ -\lr{N-1} N } & \alpha^{ -N N } \\
\alpha^{ N \lr{N-1} } & \alpha^{ \lr{N-1} \lr{N-1} } & \cdots & 1 & \cdots & \alpha^{ -\lr{N-1} \lr{N-1} } & \alpha^{ -N \lr{N-1} } \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
1 & 1 & 1 & 1 & 1 & 1 & 1 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
\alpha^{ -N \lr{N-1} } & \alpha^{ -\lr{N-1} \lr{N-1} } & \cdots & 1 & \cdots & \alpha^{ {N-1} \lr{N-1} } & \alpha^{ N \lr{N-1} } \\
\alpha^{ -N N } & \alpha^{ -N N } & \cdots & 1 & \cdots & \alpha^{ \lr{N-1} N } & \alpha^{ N N } \\
\end{bmatrix},

in each case the time domain sampling vectors are related to the Fourier components by relations of the form

\label{eqn:diodeRLCSample:160}
\Bx = \BF \BX.

## Block matrix form, with physical parameter ordering

To understand how to put \ref{eqn:diodeRLCSample:240} in block matrix form, it is helpful to consider a specific example. Consider again the specific example of the RLC circuit above, which has the form

\label{eqn:diodeRLCSample:260}
\begin{bmatrix}
0 & 0 & 0 & 0 \\
0 & Z & -Z & 1 \\
0 & -Z & Z + j \omega_0 n C & 0 \\
0 & -1 & 0 & + j \omega_0 n L \\
\end{bmatrix}
\begin{bmatrix}
V_n^{(1)} \\
V_n^{(2)} \\
V_n^{(3)} \\
I_n^{(L)} \\
\end{bmatrix}
=
\begin{bmatrix}
I_n^{(1)} \\
I_n^{(2)} \\
I_n^{(3)} \\
I_n^{(4)} \\
\end{bmatrix}

Here the $$I^{(i)}$$ terms are the DFT representations of both the linear and non-linear sources.

Suppose for example that $$N = 1$$. One way to write \ref{eqn:diodeRLCSample:260} would be

\label{eqn:diodeRLCSample:320}
\begin{aligned}
&
\begin{bmatrix}
I_{-1}^{(1)} \\
I_0^{(1)} \\
I_{1}^{(1)} \\
I_{-1}^{(2)} \\
I_0^{(2)} \\
I_{1}^{(2)} \\
I_{-1}^{(3)} \\
I_0^{(3)} \\
I_{1}^{(3)} \\
I_{-1}^{(4)} \\
I_0^{(4)} \\
I_{1}^{(4)} \\
\end{bmatrix}
=
\left[
\begin{array}{c|c|c|c}
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} \\ \hline
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
Z & 0 & 0 \\
0 & Z & 0 \\
0 & 0 & Z \\
\end{matrix} &
\begin{matrix}
-Z & 0 & 0 \\
0 & -Z & 0 \\
0 & 0 & -Z \\
\end{matrix} &
\begin{matrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 \\
\end{matrix} \\ \hline
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
-Z & 0 & 0 \\
0 & -Z & 0 \\
0 & 0 & -Z \\
\end{matrix} &
\begin{matrix}
Z & 0 & 0 \\
0 & Z & 0 \\
0 & 0 & Z \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} \\ \hline
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
-1 & 0 & 0 \\
0 & -1 & 0 \\
0 & 0 & -1 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} \\
\end{array}
\right]
\begin{bmatrix}
V_{-1}^{(1)} \\
V_{0}^{(1)} \\
V_{1}^{(1)} \\
V_{-1}^{(2)} \\
V_{0}^{(2)} \\
V_{1}^{(2)} \\
V_{-1}^{(3)} \\
V_{0}^{(3)} \\
V_{1}^{(3)} \\
I_{-1}^{(L)} \\
I_{0}^{(L)} \\
I_{1}^{(L)} \\
\end{bmatrix} \\
&+
\left[
\begin{array}{c|c|c|c}
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} \\ \hline
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} \\ \hline
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
j \omega_0 (-1) C & 0 & 0 \\
0 & j \omega_0 (0) C & 0 \\
0 & 0 & j \omega_0 (1) C \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} \\ \hline
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{matrix} &
\begin{matrix}
j \omega_0 (-1) L & 0 & 0 \\
0 & j \omega_0 (0) L & 0 \\
0 & 0 & j \omega_0 (1) L \\
\end{matrix} \\
\end{array}
\right]
\begin{bmatrix}
V_{-1}^{(1)} \\
V_{0}^{(1)} \\
V_{1}^{(1)} \\
V_{-1}^{(2)} \\
V_{0}^{(2)} \\
V_{1}^{(2)} \\
V_{-1}^{(3)} \\
V_{0}^{(3)} \\
V_{1}^{(3)} \\
I_{-1}^{(L)} \\
I_{0}^{(L)} \\
I_{1}^{(L)} \\
\end{bmatrix} \\
\end{aligned}

With a vector of fourier coeffient vectors

\label{eqn:diodeRLCSample:280}
\BV =
\begin{bmatrix}
\BV^{(1)} \\
\BV^{(2)} \\
\BV^{(3)} \\
\BI^{(L)}
\BI =
\begin{bmatrix}
\BI^{(1)} \\
\BI^{(2)} \\
\BI^{(3)} \\
\BI^{(4)}
\end{bmatrix}.

and a $$(2 N + 1) \times (2 N + 1)$$ matrix of indexes

\label{eqn:diodeRLCSample:220}
\BN =
\begin{bmatrix}
-N & & & & \\
& 1-N & & & \\
& & \ddots & & \\
& & & N-1 & \\
& & & & N \\
\end{bmatrix},

the complete block diagonalization is

\label{eqn:diodeRLCSample:300}
\boxed{
{\begin{bmatrix}
g_{rs} \BI_{2 N + 1} +
j \omega_0 c_{rs} \BN
\end{bmatrix}
}_{rs}
\BV
=
\BI.
}

## Block matrix form, with frequency ordering

It turns out that a better way of ordering the vector of Fourier components is using a frequency ordering that interleaves the physical parameters. With such an ordering the DFT MNA equations are

\label{eqn:diodeRLCSample:340}
\begin{aligned}
\BI =
&\begin{bmatrix}
I_{-1}^{(1)} \\
I_{-1}^{(2)} \\
I_{-1}^{(3)} \\
I_{-1}^{(4)} \\
I_0^{(1)} \\
I_0^{(2)} \\
I_0^{(3)} \\
I_0^{(4)} \\
I_{1}^{(1)} \\
I_{1}^{(2)} \\
I_{1}^{(3)} \\
I_{1}^{(4)} \\
\end{bmatrix}
+
\left[
\begin{array}{c|c|c}
\begin{matrix}
0 & 0 & 0 & 0 \\
0 & Z & -Z & 1 \\
0 & -Z & Z & 0 \\
0 & -1 & 0 & 0 \\
\end{matrix} & 0 & 0 \\ \hline
0 &
\begin{matrix}
0 & 0 & 0 & 0 \\
0 & Z & -Z & 1 \\
0 & -Z & Z & 0 \\
0 & -1 & 0 & 0 \\
\end{matrix} & 0 \\ \hline
0 & 0 &
\begin{matrix}
0 & 0 & 0 & 0 \\
0 & Z & -Z & 1 \\
0 & -Z & Z & 0 \\
0 & -1 & 0 & 0 \\
\end{matrix} \\
\end{array}
\right]
\begin{bmatrix}
V_{-1}^{(1)} \\
V_{-1}^{(2)} \\
V_{-1}^{(3)} \\
I_{-1}^{(L)} \\
V_{0}^{(1)} \\
V_{0}^{(2)} \\
V_{0}^{(3)} \\
I_{0}^{(L)} \\
V_{1}^{(1)} \\
V_{1}^{(2)} \\
V_{1}^{(3)} \\
I_{1}^{(L)} \\
\end{bmatrix} \\
&+
j \omega_0
\left[
\begin{array}{c|c|c}
(-1)
\begin{bmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & C & 0 \\
0 & 0 & 0 & L \\
\end{bmatrix} & 0 & 0 \\ \hline
0 &
(0)
\begin{bmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & C & 0 \\
0 & 0 & 0 & L \\
\end{bmatrix} & 0 \\ \hline
0 & 0 &
(1)
\begin{bmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & C & 0 \\
0 & 0 & 0 & L \\
\end{bmatrix} \\
\end{array}
\right]
\begin{bmatrix}
V_{-1}^{(1)} \\
V_{-1}^{(2)} \\
V_{-1}^{(3)} \\
I_{-1}^{(L)} \\
V_{0}^{(1)} \\
V_{0}^{(2)} \\
V_{0}^{(3)} \\
I_{0}^{(L)} \\
V_{1}^{(1)} \\
V_{1}^{(2)} \\
V_{1}^{(3)} \\
I_{1}^{(L)} \\
\end{bmatrix} \\
\end{aligned}

This ordering matches that of [1].

## Representing the linear sources

Assuming real sources with frequencies that are only multiples of the fundamental harmonic, a reasonable way to represent them in storage is with a pair of matrices

\label{eqn:diodeRLCSample:360}
\begin{bmatrix}
\BI \sim \BB \Bomega
\end{bmatrix}.

If $$R$$ is the dimension of $$\BG$$ and $$\BC$$, then $$\BB$$ is a $$R \times S$$ dimension matrix, where $$S$$ is the sum of

• 1, if there are any DC sources, plus
• 2 times the number of unique frequency sources.

For example, if there is a DC source and one AC source with frequency $$\nu$$, then for column vectors $$\Bb_i$$ this pair is of the form

\label{eqn:diodeRLCSample:380}
\BU \Bomega =
\begin{bmatrix}
\Bb_{-1} & \Bb_0 & \Bb_1
\end{bmatrix}
\begin{bmatrix}
– 2 \pi \nu \\
0 \\
2 \pi \nu
\end{bmatrix}.

This representation produces the time domain representation exactly when there are only DC sources, and can be used to construct the Fourier coefficients by inspection when there are AC sources. For example, for $$N = 1$$ in the example above, the Fourier coefficent vector is

\label{eqn:diodeRLCSample:400}
\BI
=
\begin{bmatrix}
\Bb_{-1} \\
\Bb_{0} \\
\Bb_{1} \\
\end{bmatrix}.

If $$N = 2$$ was used, then we would have instead

\label{eqn:diodeRLCSample:420}
\BI
=
\begin{bmatrix}
\Bzero \\
\Bb_{-1} \\
\Bb_{0} \\
\Bb_{1} \\
\Bzero \\
\end{bmatrix}.

## Representing non-linear sources

The time domain MNA fig. 1.

fig. 1. Simple diode circuit

With $$Z = 1/R, Z_g = 1/R_g$$, the KCL equations are

1. $$\lr{ v_1 – v_2 } Z_s = i_s – i_d$$
2. $$\lr{ v_2 – v_1 } Z_s + v_2 Z_g = -i_s + i_d$$

Using the model $$i_d = I^{(0)} \lr{ e^{ (v_1 – v_2)/V_T } – 1 }$$, with source $$i_s = I^{(s)} \cos( \omega_0 t )$$,
this has the block matrix form

\label{eqn:diodeRLCSample:580}
\BG =
\begin{bmatrix}
Z_s & -Z_s \\
-Z_s & Z_s + Z_g \\
\Bx =
\begin{bmatrix}
v_1(t) \\
v_2(t) \\
\end{bmatrix}

\label{eqn:diodeRLCSample:600}
\BB =
\begin{bmatrix}
I^{(s)}/2 & -I^{(0)} & I^{(s)}/2 \\
-I^{(s)}/2 & I^{(0)} & -I^{(s)}/2
\Bu(t) =
\begin{bmatrix}
e^{-j \omega_0 t} \\
1 \\
e^{j \omega_0 t}
\end{bmatrix}

\label{eqn:diodeRLCSample:620}
\BD
=
\begin{bmatrix}
I^{(0)} \\
-I^{(0)}
\Bw(t) = e^{(v_1(t) – v_2(t))/V_T}.

If $$E_n$$ is the nth DFT coefficient for $$e(t) = e^{(v_1(t) – v_2(t))/V_T}$$, then the DFT equations for the $$N = 1$$ DFT is

\label{eqn:diodeRLCSample:640}
\begin{aligned}
\lr{ V_{-1}^{(1)} – V_{-1}^{(2)} } Z_s &= I^{(s)}/2 – I^{(0)} E_{-1} \\
\lr{ V_{-1}^{(2)} – V_{-1}^{(1)} } Z_s + V_{-1}^{(2)} Z_g &= -I^{(s)}/2 + I^{(0)} E_{-1} \\
\lr{ V_{0}^{(1)} – V_{0}^{(2)} } Z_s &= I^{(0)} – I^{(0)} E_{0} \\
\lr{ V_{0}^{(2)} – V_{0}^{(1)} } Z_s + V_{0}^{(2)} Z_g &= -I^{(0)} + I^{(0)} E_{0} \\
\lr{ V_{1}^{(1)} – V_{1}^{(2)} } Z_s &= I^{(s)}/2 – I^{(0)} E_{1} \\
\lr{ V_{1}^{(2)} – V_{1}^{(1)} } Z_s + V_{1}^{(2)} Z_g &= -I^{(s)}/2 + I^{(0)} E_{1}
\end{aligned}

Let $$\Bb = [ \Bb_{-1}\, \Bb_0\, \Bb_1 ]$$, and $$\BD = [ \Bd_1 ]$$. The block matrix equivalent form, by inspection, is

\label{eqn:diodeRLCSample:660}
\begin{bmatrix}
\BG & 0 & 0 \\
0 & \BG & 0 \\
0 & 0 & \BG
\end{bmatrix}
\begin{bmatrix}
V_{-1}^{(1)} \\
V_{-1}^{(2)} \\
V_{0}^{(1)} \\
V_{0}^{(2)} \\
V_{1}^{(1)} \\
V_{1}^{(2)}
\end{bmatrix}
=
\begin{bmatrix}
\Bb_{-1} \\
\Bb_{0} \\
\Bb_{1} \\
\end{bmatrix}
+
\begin{bmatrix}
\Bd_1 E_{-1} \\
\Bd_1 E_{0} \\
\Bd_1 E_{1}
\end{bmatrix}.

This shows the stamping pattern required to form the non-linear portion of the Harmonic balance equations. The general pattern can be written as

\label{eqn:diodeRLCSample:820}
\boxed{
\mathcal{G} \BV = \BI + \mathcal{I}(\BV),
}

Here $$\mathcal{G}$$ is block diagonal, and in general has blocks of $$\BG + j \omega_0 n \BC$$. The matrix $$\BI$$ was generated from the Fourier coeffients of all the linear sources, and $$\mathcal{I}(\BV)$$ encodes all the non-linear contributions to the system.

### More general non-linear structure, for multiple diodes

For the diode exponentials, these non-linear term will be of the form

\label{eqn:diodeRLCSample:680}
\BD \Bw(t)
=
\begin{bmatrix}
\Bd_1 & \Bd_2 & \cdots & \Bd_S
\end{bmatrix}
\begin{bmatrix}
w_1(t) \\
w_2(t) \\
\vdots \\
w_S(t) \\
\end{bmatrix},

where $$w_i(t) = \exp( (v_{i,1}(t) – v_{i,2}(t))/V_{T,i} )$$. If the DFT coordinates of these functions are $$E_n^{(i)}$$, then the frequency domain representation is

\label{eqn:diodeRLCSample:700}
\mathcal{I}(\BV)
=
\sum_{i = 1}^S
\begin{bmatrix}
\Bd_i E_{-N}^{(i)} \\
\Bd_i E_{1-N}^{(i)} \\
\vdots \\
\Bd_i E_{N-1}^{(i)} \\
\Bd_i E_{N}^{(i)} \\
\end{bmatrix}.

This is a $$R (2 N + 1) \times 1$$ matrix, as expected.

The computation of these DFT coordinates is a bit messy since they are time dependent, and thus dependent on the (unknown) values of $$V_n^{(1)}$$. Consider the above circuit as an example where we have

\label{eqn:diodeRLCSample:720}
w(t) = \exp\lr{ (v_1(t) – v_2(t))/V_T }.

The DFT inverse is

\label{eqn:diodeRLCSample:740}
E_n = \sum_{k=-N}^N \exp\lr{ (v_1(t_k) – v_2(t_k))/V_T } e^{-j \omega_0 n t_k }.

With $$\BE = ( E_{-N}, E_{1-N}, \cdots, E_{N-1}, E_N )$$, this is

\label{eqn:diodeRLCSample:760}
\BE =
\inv{2 N + 1} \overline{{\BF}}
\begin{bmatrix}
\exp\lr{ (v_{-N}^{(1)} – v_{-N}^{(2)})/v_T } \\
\exp\lr{ (v_{1-N}^{(1)} – v_{1-N}^{(2)})/v_T } \\
\vdots \\
\exp\lr{ (v_{N-1}^{(1)} – v_{N-1}^{(2)})/v_T } \\
\exp\lr{ (v_{N}^{(1)} – v_{N}^{(2)})/v_T } \\
\end{bmatrix}
=
\inv{2 N + 1} \overline{{\BF}}
\begin{bmatrix}
\exp\lr{ {[\BF (\BV^{(1)} – \BV^{(2)})/v_T]}_{-N} } \\
\exp\lr{ {[\BF (\BV^{(1)} – \BV^{(2)})/v_T]}_{1-N} } \\
\vdots \\
\exp\lr{ {[\BF (\BV^{(1)} – \BV^{(2)})/v_T]}_{N-1} } \\
\exp\lr{ {[\BF (\BV^{(1)} – \BV^{(2)})/v_T]}_{N} }
\end{bmatrix}.

With the introduction a term by term exponentiation operator

\label{eqn:diodeRLCSample:780}
\exp[ \Bx ] =
\begin{bmatrix}
\exp(x_1) \\
\exp(x_2) \\
\vdots \\
\end{bmatrix},

this one Fourier coefficient vector is

\label{eqn:diodeRLCSample:800}
\BE =
\inv{2 N + 1} \overline{{\BF}}
\exp[ \BF (\BV^{(1)} – \BV^{(2)})/v_T ].

This is now completely expressed in terms of the unknown Fourier component vectors, each a subset of the aggreggated “voltage”, range selectable with the Matlab operation $$\BV^{(i)} = \BV(i:R:end)$$.

## Newton’s method

In order to solve the system, Newton’s method on the Fourier coeffients is required. Solutions to $$\mathcal{F}(\BV) = 0$$ are sought, where

\label{eqn:diodeRLCSample:840}
\mathcal{F}(\BV) = \mathcal{G} \BV – \BI – \mathcal{I}(\BV).

Here the sources current vector DFT coordinates have been split into the linear contributions $$\BI$$ and nonlinear contributions $$\mathcal{I}$$ defined by \ref{eqn:diodeRLCSample:700}.

Working with ones-indexed coordinates of $$\BV = [V_k]_k$$, where $$k \in [1, R(2 N + 1) ]$$, the Jacobian is

\label{eqn:diodeRLCSample:860}
\BJ = \mathcal{G} – {[ \PDi{V_s}{\mathcal{I}_r} ]}_{rs}.

To calculate these partials we need the partials of the coordinates of $$\BE$$ of
\ref{eqn:diodeRLCSample:800}. The kth coordinate of $$\BV^{(1)}, \BV^{(2)}$$ in terms of the coordinates of the $$R(2 N + 1)$$ vector of unknowns $$\BV$$ are

\label{eqn:diodeRLCSample:880}
\begin{aligned}
[\BV^{(1)}]_k &= V_{1 + (k-1) R} \\
[\BV^{(2)}]_k &= V_{2 + (k-1) R}
\end{aligned}

Using summation convention, with sums over any repeated indexes implied, those coordinates are

\label{eqn:diodeRLCSample:900}
E_r =
\inv{2 N + 1} \overline{{F}}_{r a}
\exp\lr{ F_{ab}
(
V_{1 + (b-1) R} –
V_{2 + (b-1) R}
)/v_T }.

\label{eqn:diodeRLCSample:920}
\PD{V_s}{E_r} =
\inv{2 N + 1} \inv{v_T} \overline{{F}}_{r a} F_{ab}
\exp\lr{ F_{ab}
(
V_{1 + (b-1) R} –
V_{2 + (b-1) R}
)/v_T }
\lr{
\delta_{s,1 + (b-1) R} –
\delta_{s,2 + (b-1) R}
}.

### Generalization

To generalize this, suppose that the diode exponential was associated with voltages spanning nodes $$m, n$$ so that

\label{eqn:diodeRLCSample:980}
\BE =
\inv{2 N + 1} \overline{{\BF}}
\exp[ \BF (\BV^{(m)} – \BV^{(n)})/v_T ].

In this case, the coordinates of the physical “voltages” are

\label{eqn:diodeRLCSample:1020}
\begin{aligned}
[\BV^{(m)}]_k &= V_{m + (k-1) R} \\
[\BV^{(n)}]_k &= V_{n + (k-1) R}
\end{aligned},

so

\label{eqn:diodeRLCSample:1040}
E_r =
\inv{2 N + 1} \overline{{F}}_{r a}
\exp\lr{ F_{ab}
(
V_{m + (b-1) R} –
V_{n + (b-1) R}
)/v_T }.

The Jacobian partials are

\label{eqn:diodeRLCSample:1060}
\PD{V_s}{E_r} =
\inv{2 N + 1} \inv{v_T} \overline{{F}}_{r a} F_{ab}
\exp\lr{ F_{ab}
(
V_{m + (b-1) R} –
V_{n + (b-1) R}
)/v_T }
\lr{
\delta_{s,m + (b-1) R} –
\delta_{s,n + (b-1) R}
}.

Note that this Jacobian

\label{eqn:diodeRLCSample:1080}
\BJ_\BE =
{
\begin{bmatrix}
\PD{V_s}{E_r}
\end{bmatrix}}_{rs}

is a $$(2 N + 1) \times R(2N + 1)$$ matrix.

The full Jacobian of $$\mathcal{I}(\BV)$$ is

\label{eqn:diodeRLCSample:1120}
\BJ_{\mathcal{I}}(\BV)
=
\sum_{i = 1}^S
\begin{bmatrix}
\Bd_i \PD{\BV}{E_{-N}^{(i)}} \\
\Bd_i \PD{\BV}{E_{1-N}^{(i)}} \\
\vdots \\
\Bd_i \PD{\BV}{E_{N-1}^{(i)}} \\
\Bd_i \PD{\BV}{E_{N}^{(i)}} \\
\end{bmatrix}.

Where $$\PDi{\BV}{E_{k}^{(i)}}$$ is the kth row of $$\BJ_{\BE^{(i)}}$$. The complete Jacobian is

\label{eqn:diodeRLCSample:1100}
\BJ = \mathcal{G} – \BJ_{\mathcal{I}}(\BV).

# References

Giannini and Giorgio Leuzzi. Nonlinear Microwave Circuit Design. Wiley Online Library, 2004.

## A sample diode RLC circuit

To get a feel for how to generate the MLN equations for a circuit that has both RLC and non-linear components, consider the circuit of fig. 1.

fig. 1. An RLC circuit with a diode.

The KCL equations for this circuit are

1. $$0 = i_s – i_d$$
2. $$i_L + \frac{v_2 – v_3}{R} = i_d$$
3. $$\frac{v_3 – v_2}{R} + C \frac{dv_3}{dt} = 0$$
4. $$-v_2 + L \frac{d i_L}{dt} = 0$$
5. $$i_d = I_0 \lr{ e^{(v_1 – v_2)/v_T} – 1}$$

FIXME: for the diode, is that the right voltage sign with respect to the current direction?

With $$Z = 1/R$$, these can be put into the standard MLN matrix form as

\label{eqn:diodeRLCSample:20}
\begin{bmatrix}
0 & 0 & 0 & 0 \\
0 & Z & -Z & 1 \\
0 & -Z & Z & 0 \\
0 & -1 & 0 & 0 \\
\end{bmatrix}
\begin{bmatrix}
v_1 \\
v_2 \\
v_3 \\
i_L \\
\end{bmatrix}
+
\begin{bmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & C & 0 \\
0 & 0 & 0 & L \\
\end{bmatrix}
{\begin{bmatrix}
v_1 \\
v_2 \\
v_3 \\
i_L \\
\end{bmatrix}}’
=
\begin{bmatrix}
I_0 & 1 \\
-I_0 & 0 \\
0 & 0 \\
0 & 0 \\
\end{bmatrix}
\begin{bmatrix}
1 \\
i_s(t) \\
\end{bmatrix}
+
\begin{bmatrix}
-I_0 \\
I_0 \\
0 \\
0 \\
\end{bmatrix}
\begin{bmatrix}
e^{(v_2 – v_3)/v_T}
\end{bmatrix}

Let’s write this as

\label{eqn:diodeRLCSample:40}
\BG \BX(t) + \BC \dot{\BX}(t) = \BB \Bu(t) + \BD \Bw(t).

Here $$\Bu(t)$$ collects up all the unique signature sources (for example sources with each different frequency in the system), and $$\Bw(t)$$ is a vector of all the unique non-linear (time dependent) terms.

Assuming a bandwidth limited periodic source we know how to express any of the time dependent variables $$v_1, …$$ in terms of their (discrete) Fourier transforms. Suppose that

the Fourier coefficients for $$v_a(t), u_b(t), w_c(t)$$ are given by

\label{eqn:diodeRLCSample:60}
\begin{aligned}
v_a(t) &= \sum_{n = -N}^N V_n^{(a)} e^{j \omega_0 n t} \\
u_b(t) &= \sum_{n = -N}^N U_n^{(b)} e^{j \omega_0 n t} \\
w_c(t) &= \sum_{n = -N}^N W_n^{(c)} e^{j \omega_0 n t}.
\end{aligned}

For example, in this circuit, if the source was zero phase signal at the fundamental frequency of our Fourier basis ($$i_s(t) = e^{j \omega_0 t}$$), the only non-zero Fourier components $$U_n^{(a)}$$ would be $$U_0^{(1)} = 1, U_1^{(2)} = 1$$.

Equation \ref{eqn:diodeRLCSample:40} then becomes

\label{eqn:diodeRLCSample:80}
0 = \sum_{n=-N}^N
e^{j n \omega_0 t}
\lr{
\lr{
\BG + j \omega_0 n \BC
}
\begin{bmatrix}
V_n^{(1)} \\
V_n^{(2)} \\
\vdots
\end{bmatrix}
– \BB
\begin{bmatrix}
U_n^{(1)} \\
U_n^{(2)} \\
\end{bmatrix}
– \BD
\begin{bmatrix}
W_n^{(1)} \\
\end{bmatrix}
}

The time dependence in the linear terms is nicely taken of by this transformation to the frequency domain. However, we have a fairly messy structure with sums of Fourier components instead of the nice Fourier component vectors that we see in \S A.4 of [1]. That reference does consider multivariable problems like this one, so it looks like fully digesting that methodology is the next step.

# References

Giannini and Giorgio Leuzzi. Nonlinear Microwave Circuit Design. Wiley Online Library, 2004.

## Matrix form for discrete time Fourier transform

December 4, 2014 ece1254 No comments , ,

## Matrix form

The discrete time Fourier transform has been seen to have the form

\label{eqn:discreteFourierMatrixForm:200}
x_k = \sum_{n = -N}^N X_n e^{ 2 \pi j n k /(2 N + 1)}

\label{eqn:discreteFourierMatrixForm:220}
X_n = \inv{2 N + 1} \sum_{k = -N}^N x_k e^{- 2 \pi j n k /(2 N + 1)}.

A matrix representation of this form is desired. Let

\label{eqndiscreteFourierMatrixForm:240}
\Bx =
\begin{bmatrix}
x_N \\
\vdots \\
x_0 \\
\vdots \\
x_{-N}
\end{bmatrix}

\label{eqndiscreteFourierMatrixForm:260}
\BX =
\begin{bmatrix}
X_N \\
\vdots \\
X_0 \\
\vdots \\
X_{-N}
\end{bmatrix}

\ref{eqn:discreteFourierMatrixForm:200} written out in full is
\label{eqn:discreteFourierMatrixForm:280}
\begin{aligned}
x_k
&= X_N e^{ 2 \pi j N k /(2 N + 1)} \\
&+ X_{N-1} e^{ 2 \pi j \lr{N-1} k /(2 N + 1)} \\
&+ \cdots \\
&+ X_0 \\
&+ \cdots \\
&+ X_{1-N} e^{ -2 \pi j \lr{N-1} k /(2 N + 1)} \\
&+ X_{-N} e^{ -2 \pi j N k /(2 N + 1)}.
\end{aligned}

With $$\alpha = e^{ 2 \pi j /(2 N + 1) }$$ the matrix form is

\label{eqn:discreteFourierMatrixForm:300}
\Bx =
\begin{bmatrix}
\alpha^{ N N } & \alpha^{ \lr{N-1} N } & \cdots & 1 & \cdots & \alpha^{ -\lr{N-1} N } & \alpha^{ -N N } \\
\alpha^{ N \lr{N-1} } & \alpha^{ \lr{N-1} \lr{N-1} } & \cdots & 1 & \cdots & \alpha^{ -\lr{N-1} \lr{N-1} } & \alpha^{ -N \lr{N-1} } \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
1 & 1 & 1 & 1 & 1 & 1 & 1 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
\alpha^{ -N \lr{N-1} } & \alpha^{ -\lr{N-1} \lr{N-1} } & \cdots & 1 & \cdots & \alpha^{ {N-1} \lr{N-1} } & \alpha^{ N \lr{N-1} } \\
\alpha^{ -N N } & \alpha^{ -N N } & \cdots & 1 & \cdots & \alpha^{ \lr{N-1} N } & \alpha^{ N N } \\
\end{bmatrix}
\BX

Similarily, from \ref{eqn:discreteFourierMatrixForm:220}, the inverse relation expands out to

\label{eqn:discreteFourierMatrixForm:320}
\begin{aligned}
( 2 N + 1 ) X_n
&= x_N e^{- 2 \pi j n N /(2 N + 1)} \\
&+ x_{N-1} e^{- 2 \pi j n \lr{N-1}/(2 N + 1)} \\
&\cdots \\
&+ x_0 \\
&\cdots \\
&+ x_{1-N} e^{ 2 \pi j n \lr{ N-1 }/(2 N + 1)} \\
&+ x_{-N} e^{ 2 \pi j n N/(2 N + 1)},
\end{aligned}

with a matrix form of

\label{eqn:discreteFourierMatrixForm:340}
( 2 N + 1 ) \BX =
\begin{bmatrix}
\alpha^{- N N } & \alpha^{- N \lr{N-1}} &\cdots & 1 &\cdots & \alpha^{ N \lr{ N-1 }} & \alpha^{ N N} \\
\alpha^{- \lr{N-1} N } & \alpha^{- \lr{N-1} \lr{N-1}} &\cdots & 1 &\cdots & \alpha^{ \lr{N-1} \lr{ N-1 }} & \alpha^{ \lr{N-1} N} \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
1 & 1 & 1 & 1 & 1 & 1 & 1 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
\alpha^{\lr{N-1} N } & \alpha^{ \lr{N-1} \lr{N-1}} &\cdots & 1 &\cdots & \alpha^{ -\lr{N-1} \lr{ N-1 }} & \alpha^{ -\lr{N-1} N} \\
\alpha^{ N N } & \alpha^{ N \lr{N-1}} &\cdots & 1 &\cdots & \alpha^{ -N \lr{ N-1 }} & \alpha^{ -N N} \\
\end{bmatrix}

Lettting

\label{eqn:discreteFourierMatrixForm:360}
\BF =
\begin{bmatrix}
\alpha^{ N N } & \alpha^{ \lr{N-1} N } & \cdots & 1 & \cdots & \alpha^{ -\lr{N-1} N } & \alpha^{ -N N } \\
\alpha^{ N \lr{N-1} } & \alpha^{ \lr{N-1} \lr{N-1} } & \cdots & 1 & \cdots & \alpha^{ -\lr{N-1} \lr{N-1} } & \alpha^{ -N \lr{N-1} } \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
1 & 1 & 1 & 1 & 1 & 1 & 1 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
\alpha^{ -N \lr{N-1} } & \alpha^{ -\lr{N-1} \lr{N-1} } & \cdots & 1 & \cdots & \alpha^{ {N-1} \lr{N-1} } & \alpha^{ N \lr{N-1} } \\
\alpha^{ -N N } & \alpha^{ -N N } & \cdots & 1 & \cdots & \alpha^{ \lr{N-1} N } & \alpha^{ N N } \\
\end{bmatrix},

the discrete transform pair has the following compactly matrix representation

\label{eqn:discreteFourierMatrixForm:380}
\Bx = \BF \BX

\label{eqn:discreteFourierMatrixForm:400}
\BX = \inv{2 N + 1} \overline{{\BF}} \Bx,

where $$\overline{{\BF}}$$ is the complex conjugate of $$\BF$$.

# References

## Discrete Fourier Transform

In [1] a verification of the discrete Fourier transform pairs was performed. A much different looking discrete Fourier transform pair is given in [2] \$ A.4. This transform pair samples the points at what are called the Nykvist time instants given by

\label{eqn:discreteFourier:20}
t_k = \frac{T k}{2 N + 1}, \qquad k \in [-N, \cdots N]

Note that the endpoints of these sampling points are not $$\pm T/2$$, but are instead at

\label{eqn:discreteFourier:40}
\pm \frac{T}{2} \inv{1 + 1/N},

which are slightly within the interior of the $$[-T/2, T/2]$$ range of interest. The reason for this slightly odd seeming selection of sampling times becomes clear if one calculate the inversion relations.

Given a periodic ($$\omega_0 T = 2 \pi$$) bandwith limited signal evaluated only at the Nykvist times $$t_k$$,

\label{eqn:discreteFourier:60}
\boxed{
x(t_k) = \sum_{n = -N}^N X_n e^{ j n \omega_0 t_k},
}

assume that an inversion relation can be found. To find $$X_n$$ evaluate the sum

\label{eqn:discreteFourier:80}
\begin{aligned}
&\sum_{k = -N}^N x(t_k) e^{-j m \omega_0 t_k} \\
\sum_{k = -N}^N
\lr{
\sum_{n = -N}^N X_n e^{ j n \omega_0 t_k}
}
e^{-j m \omega_0 t_k} \\
\sum_{n = -N}^N X_n
\sum_{k = -N}^N
e^{ j (n -m )\omega_0 t_k}
\end{aligned}

This interior sum has the value $$2 N + 1$$ when $$n = m$$. For $$n \ne m$$, and
$$a = e^{j (n -m ) \frac{2 \pi}{2 N + 1}}$$, this is

\label{eqn:discreteFourier:100}
\begin{aligned}
\sum_{k = -N}^N
e^{ j (n -m )\omega_0 t_k}
&=
\sum_{k = -N}^N
e^{ j (n -m )\omega_0 \frac{T k}{2 N + 1}} \\
&=
\sum_{k = -N}^N a^k \\
&=
a^{-N} \sum_{k = -N}^N a^{k+ N} \\
&=
a^{-N} \sum_{r = 0}^{2 N} a^{r} \\
&=
a^{-N} \frac{a^{2 N + 1} – 1}{a – 1}.
\end{aligned}

Since $$a^{2 N + 1} = e^{2 \pi j (n – m)} = 1$$, this sum is zero when $$n \ne m$$. This means that

\label{eqn:discreteFourier:120}
\sum_{k = -N}^N
e^{ j (n -m )\omega_0 t_k} = (2 N + 1) \delta_{n,m},

which provides the desired Fourier inversion relation

\label{eqn:discreteFourier:140}
\boxed{
X_m = \inv{2 N + 1} \sum_{k = -N}^N x(t_k) e^{-j m \omega_0 t_k}.
}

# References

[1] Peeter Joot. Condensed matter physics., appendix: Discrete Fourier transform. 2013. URL http://peeterjoot.com/archives/math2013/phy487.pdf. [Online; accessed 02-December-2014].

[2] Giannini and Leuzzi Nonlinear Microwave Circuit Design. Wiley Online Library, 2004.