## 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

## 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.

## Laplace transform refresher

November 27, 2014 ece1254 No comments , , , ,

Laplace transforms were used to solve the MNA equations for time dependent systems, and to find the moments used to in MOR.

For the record, the Laplace transform is defined as:

\label{eqn:laplaceTransformVec:20}
\boxed{
\LL( f(t) ) =
\int_0^\infty e^{-s t} f(t) dt.
}

The only Laplace transform pair used in the lectures is that of the first derivative

\label{eqn:laplaceTransformVec:40}
\begin{aligned}
\LL(f'(t)) &=
\int_0^\infty e^{-s t} \ddt{f(t)} dt \\
&=
\evalrange{e^{-s t} f(t)}{0}{\infty} – (-s) \int_0^\infty e^{-s t} f(t) dt \\
&=
-f(0) + s \LL(f(t)).
\end{aligned}

Here it is loosely assumed that the real part of $$s$$ is positive, and that $$f(t)$$ is “well defined” enough that $$e^{-s \infty } f(\infty) \rightarrow 0$$.

Where used in the lectures, the laplace transforms were of vectors such as the matrix vector product $$\LL(\BG \Bx(t))$$. Because such a product is linear, observe that it can be expressed as the original matrix times a vector of Laplace transforms

\label{eqn:laplaceTransformVec:60}
\begin{aligned}
\LL( \BG \Bx(t) )
&=
\LL {\begin{bmatrix}
G_{i k} x_k(t)
\end{bmatrix}}_i \\
&=
{\begin{bmatrix}
G_{i k} \LL x_k(t)
\end{bmatrix}}_i \\
&=
\BG
{\begin{bmatrix}
\LL x_i(t)
\end{bmatrix}}_i.
\end{aligned}

The following notation was used in the lectures for such a vector of Laplace transforms

\label{eqn:laplaceTransformVec:80}
\BX(s) = \LL \Bx(t) =
{\begin{bmatrix}
\LL x_i(t)
\end{bmatrix}}_i.

## ECE1254H Modeling of Multiphysics Systems. Lecture 19: Model order reduction (cont).. Taught by Prof. Piero Triverio

November 26, 2014 ece1254 No comments , , ,

## Disclaimer

Peeter’s lecture notes from class. These may be incoherent and rough.

## Model order reduction (cont).

An approximation of the following system is sought

\label{eqn:multiphysicsL19:20}
\BG \Bx(t) + C \dot{\Bx}(t) = B \Bu(t)

\label{eqn:multiphysicsL19:40}
\By(t) = \BL^\T \Bx(t).

The strategy is to attempt to find a $$N \times q$$ projector $$\BV$$ of the form

\label{eqn:multiphysicsL19:60}
\BV =
\begin{bmatrix}
\Bv_1 & \Bv_2 & \cdots & \Bv_q
\end{bmatrix}

so that the solution of the constrained q-variable state vector $$\Bx_q$$ is sought after letting

\label{eqn:multiphysicsL19:80}
\Bx(t) = \BV \Bx_q(t).

## Moment matching

\label{eqn:multiphysicsL19:100}
\begin{aligned}
\BF(s)
&= \lr{ \BG + s \BC }^{-1} \BB \\
&= \BM_0 + \BM_1 s + \BM_2 s^2 + \cdots + \BM_{q-1} s^{q-1} + M_q s^q + \cdots
\end{aligned}

The reduced model is created such that

\label{eqn:multiphysicsProblemSet3b:120}
\BF_q(s)
=
\BM_0 + \BM_1 s + \BM_2 s^2 + \cdots + \BM_{q-1} s^{q-1} + \tilde{\BM}_q s^q.

Form an $$N \times q$$ projection matrix

\label{eqn:multiphysicsL19:140}
\BV_q \equiv
\begin{bmatrix}
\BM_0 & \BM_1 & \cdots & \BM_{q-1}
\end{bmatrix}

With the substitution of fig. 1, becomes

This is a system of $$N$$ equations, in $$q$$ unknowns. A set of moments from the frequency domain have been used to project the time domain system. This relies on the following unproved theorem (references to come)

### Theorem

If $$\text{span}\{ \Bv_q \} = \text{span} \{ \BM_0, \BM_1, \cdots, \BM_{q-1} \}$$, then the reduced model will match the first $$q$$ moments.

Left multiplication by $$\BV_q^\T$$ yields fig. 2.

This is now a system of $$q$$ equations in $$q$$ unknowns.

With

\label{eqn:multiphysicsL19:160}
\BG_q = \BV_q^\T \BG \BV_q

\label{eqn:multiphysicsL19:180}
\BC_q = \BV_q^\T \BC \BV_q

\label{eqn:multiphysicsL19:200}
\BB_q = \BV_q^\T \BB

\label{eqn:multiphysicsL19:220}
\BL_q^\T = \BL^\T \BV_q

the system is reduced to

\label{eqn:multiphysicsL19:240}
\BG_q \Bx_q(t) + \BC_q \dot{\Bx}_q(t) = \BB_q \Bu(t).

\label{eqn:multiphysicsL19:260}
\By(t) \approx \BL_q^\T \Bx_q(t)

### Moments calculation

Using
\label{eqn:multiphysicsL19:300}
\lr{ \BG + s \BC }^{-1} \BB = \BM_0 + \BM_1 s + \BM_2 s^2 + \cdots

thus
\label{eqn:multiphysicsL19:320}
\begin{aligned}
\BB &=
\lr{ \BG + s \BC }
\BM_0 +
\lr{ \BG + s \BC }
\BM_1 s +
\lr{ \BG + s \BC }
\BM_2 s^2 + \cdots \\
&=
\BG \BM_0
+ s \lr{ C \BM_0 + \BG \BM_1 }
+ s^2 \lr{ C \BM_1 + \BG \BM_2 }
+ \cdots
\end{aligned}

Since $$\BB$$ is a zeroth order matrix, setting all the coefficients of $$s$$ equal to zero provides a method to solve for the moments

\label{eqn:multiphysicsL19:340}
\begin{aligned}
\BB &= \BG \BM_0 \\
-\BC \BM_0 &= \BG \BM_1 \\
-\BC \BM_1 &= \BG \BM_2 \\
\end{aligned}

The moment $$\BM_0$$ can be found with LU of $$\BG$$, plus the forward and backward substitutions. Proceeding recursively, using the already computed LU factorization, each subsequent moment calculation requires only one pair of forward and backward substitutions.

Numerically, each moment has the exact value

\label{eqn:multiphysicsL19:360}
\BM_q = \lr{- \BG^{-1} \BC }^q \BM_0.

As $$q \rightarrow \infty$$, this goes to some limit, say $$\Bw$$. The value $$\Bw$$ is related to the largest eigenvalue of $$-\BG^{-1} \BC$$. Incidentally, this can be used to find the largest eigenvalue of $$-\BG^{-1} \BC$$.

The largest eigenvalue of this matrix will dominate these factors, and can cause some numerical trouble. For this reason it is desirable to avoid such explicit moment determination, instead using implicit methods.

The key is to utilize the theorem above, and look instead for an alternate basis $$\{ \Bv_q \}$$ that also spans $$\{ \BM_0, \BM_1, \cdots, \BM_q \}$$.

### Space generate by the moments

Write

\label{eqn:multiphysicsL19:380}
\BM_q = \BA^q \BR,

where in this case

\label{eqn:multiphysicsL19:400}
\begin{aligned}
\BA &= – \BG^{-1} \BC \\
\BR &= \BM_0 = -\BG \BB
\end{aligned}

The span of interest is

\label{eqn:multiphysicsL19:420}
\text{span} \{ \BR, \BA \BR, \BA^2 \BR, \cdots, \BA^{q-1} \BR \}.

Such a sequence is called a Krylov subspace. One method to compute such a basis, the Arnoldi process, relies on Gram-Schmidt orthonormalization methods:

Some numerical examples and plots on the class slides.

## Stability of discretized linear differential equations

In class today was a highlight of stability methods for linear multistep methods. To motivate the methods used, it is helpful to take a step back and review stability concepts for LDE systems.

By way of example, consider a second order LDE homogeneous system defined by

\label{eqn:stabilityLDEandDiscreteTime:20}
\frac{d^2 x}{dt^2} + 3 \frac{dx}{dt} + 2 = 0.

Such a system can be solved by assuming an exponential solution, say

\label{eqn:stabilityLDEandDiscreteTime:40}
x(t) = e^{s t}.

Substitution gives

\label{eqn:stabilityLDEandDiscreteTime:60}
e^{st} \lr{ s^2 + 3 s + 2 } = 0,

The polynomial part of this equation, the characteristic equation has roots $$s = -2, -1$$.

The general solution of \ref{eqn:stabilityLDEandDiscreteTime:20} is formed by a superposition of solutions for each value of $$s$$

\label{eqn:stabilityLDEandDiscreteTime:80}
x(t) = a e^{-2 t} + b e^{-t}.

Independent of any selection of the superposition constants $$a, b$$, this function will not blow up as $$t \rightarrow \infty$$.

This “stability” is due to the fact that both of the characteristic equation roots lie in the left hand Argand plane.

Now consider a discretized form of this LDE. This will have the form

\label{eqn:stabilityLDEandDiscreteTime:100}
\begin{aligned}
0 &=
\inv{\lr{\Delta t}^2}
\lr{ x_{n+2} – 2 x_{n-1} + x_n } + \frac{3}{\Delta t} \lr{ x_{n+1} – x_n } + 2
x_n \\
&=
x_{n+2} \lr{
\inv{\lr{\Delta t}^2}
}
+
x_{n+1} \lr{
\frac{3}{\Delta t}
-\frac{2}{\lr{\Delta t}^2}
}
+
x_{n} \lr{
\frac{1}{\lr{\Delta t}^2}
-\frac{3}{\Delta t}
+ 2
},
\end{aligned}

or

\label{eqn:stabilityLDEandDiscreteTime:220}
0
=
x_{n+2}
+
x_{n+1} \lr{
3 \Delta t – 2
}
+
x_{n} \lr{
1 – 3 \Delta t + 2 \lr{ \Delta t}^2
}.

Note that after discretization, each subsequent index corresponds to a time shift. Also observe that the coefficients of this discretized equation are dependent on the discretization interval size $$\Delta t$$. If the specifics of those coefficients are ignored, a general form with the following structure can be observed

\label{eqn:stabilityLDEandDiscreteTime:120}
0 =
x_{n+2} \gamma_0
+
x_{n+1} \gamma_1
+
x_{n} \gamma_2.

It turns out that, much like the LDE solution by characteristic polynomial, it is possible to attack this problem by assuming a solution of the form

\label{eqn:stabilityLDEandDiscreteTime:140}
x_n = C z^n.

A time shift index change $$x_n \rightarrow x_{n+1}$$ results in a power adjustment in this assumed solution. This substitution applied to \ref{eqn:stabilityLDEandDiscreteTime:120} yields

\label{eqn:stabilityLDEandDiscreteTime:160}
0 =
C z^n
\lr{
z^{2} \gamma_0
+
z \gamma_1
+
1 \gamma_2
},

Suppose that this polynomial has roots $$z \in \{z_1, z_2\}$$. A superposition, such as

\label{eqn:stabilityLDEandDiscreteTime:180}
x_n = a z_1^n + b z_2^n,

will also be a solution since insertion of this into the RHS of \ref{eqn:stabilityLDEandDiscreteTime:120} yields

\label{eqn:stabilityLDEandDiscreteTime:200}
a z_1^n
\lr{
z_1^{2} \gamma_0
+
z_1 \gamma_1
+
\gamma_2
}
+
b
z_2^n
\lr{
z_2^{2} \gamma_0
+
z_2 \gamma_1
+
\gamma_2
}
=
a z_1^n \times 0
+b z_2^n \times 0.

The zero equality follows since $$z_1, z_2$$ are both roots of the characteristic equation for this discretized LDE.
In the discrete $$z$$ domain stability requires that the roots satisfy the bound $$\Abs{z} < 1$$, a different stability criteria than in the continuous domain. In fact, there is no a-priori guarantee that stability in the continuous domain will imply stability in the discretized domain. Let's plot those z-domain roots for this example LDE, using $$\Delta t \in \{ 1/2, 1, 2 \}$$. The respective characteristic polynomials are $$\label{eqn:stabilityLDEandDiscreteTime:260} 0 = z^2 - \inv{2} z = z \lr{ z - \inv{2} }$$ $$\label{eqn:stabilityLDEandDiscreteTime:240} 0 = z^2 + z = z\lr{ z + 1 }$$ $$\label{eqn:stabilityLDEandDiscreteTime:280} 0 = z^2 + 4 z + 3 = (z + 3)(z + 1).$$ These have respective roots $$\label{eqn:stabilityLDEandDiscreteTime:300} z = 0, \inv{2}$$ $$\label{eqn:stabilityLDEandDiscreteTime:320} z = 0, -1$$ $$\label{eqn:stabilityLDEandDiscreteTime:340} z = -1, -3$$ Only the first discretization of these three yields stable solutions in the z domain, although it appears that $$\Delta t = 1$$ is right on the boundary.

## ECE1254H Modeling of Multiphysics Systems. Lecture 16: LMS systems and stability. Taught by Prof. Piero Triverio

November 17, 2014 ece1254 No comments , , ,

## Disclaimer

Peeter’s lecture notes from class. These may be incoherent and rough.

## Residual for LMS methods

### Mostly on slides:

12_ODS.pdf

Residual is illustrated in fig. 1, assuming that the iterative method was accurate until $$t_{n}$$

fig. 1. Residual illustrated

### Summary

• [FE]: $$R_{n+1} \sim \lr{ \Delta t}^2$$. This is of order $$p = 1$$.
• [BE]: $$R_{n+1} \sim \lr{ \Delta t}^2$$. This is of order $$p = 1$$.
• [TR]: $$R_{n+1} \sim \lr{ \Delta t}^3$$. This is of order $$p = 2$$.
• [BESTE]: $$R_{n+1} \sim \lr{ \Delta t}^4$$. This is of order $$p = 3$$.

## Global error estimate

Suppose $$t \in [0, 1] s$$, with $$N = 1/{\Delta t}$$ intervals. For a method with local error of order $$R_{n+1} \sim \lr{ \Delta t}^2$$ the global error is approximately $$N R_{n+1} \sim \Delta t$$.

## Stability

Recall that a linear multistep method (LMS) was a system of the form

\label{eqn:multiphysicsL16:20}
\sum_{j=-1}^{k-1} \alpha_j x_{n-j} = \Delta t \sum_{j=-1}^{k-1} \beta_j f( x_{n-j}, t_{n-j} )

Consider a one dimensional test problem

\label{eqn:multiphysicsL16:40}
\dot{x}(t) = \lambda x(t)

where as in fig. 2, $$\Re(\lambda) < 0$$ is assumed to ensure stability.

fig. 2. Stable system

Linear stability theory can be thought of as asking the question: “Is the solution of \ref{eqn:multiphysicsL16:40} computed by my LMS method also stable?”

Application of \ref{eqn:multiphysicsL16:20} to \ref{eqn:multiphysicsL16:40} gives

\label{eqn:multiphysicsL16:60}
\sum_{j=-1}^{k-1} \alpha_j x_{n-j} = \Delta t \sum_{j=-1}^{k-1} \beta_j \lambda x_{n-j},

or
\label{eqn:multiphysicsL16:80}
\sum_{j=-1}^{k-1} \lr{ \alpha_j – \Delta \beta_j \lambda }
x_{n-j} = 0.

With

\label{eqn:multiphysicsL16:100}
\gamma_j = \alpha_j – \Delta \beta_j \lambda,

this expands to
\label{eqn:multiphysicsL16:120}
\gamma_{-1} x_{n+1}
+
\gamma_{0} x_{n}
+
\gamma_{1} x_{n-1}
+
\cdots
+
\gamma_{k-1} x_{n-k} .

This can be seen as a

• discrete time system
• FIR filter

The numerical solution $$x_n$$ will be stable if \ref{eqn:multiphysicsL16:120} is stable.

A characteristic equation associated with \ref{eqn:multiphysicsL16:120} can be defined as

\label{eqn:multiphysicsL16:140}
\gamma_{-1} z^k
+
\gamma_{0} z^{k-1}
+
\gamma_{1} z^{k-2}
+
\cdots
+
\gamma_{k-1} = 0.

This is a polynomial with roots $$z_n$$ (poles). This is stable if the poles satisfy $$\Abs{z_n} < 1$$, as illustrated in fig. 3

Stability

Observe that the $$\gamma’s$$ are dependent on $$\Delta t$$.

FIXME: There’s a lot of handwaving here that could use more strict justification. Check if the text covers this in more detail.

### Example: Forward Euler stability

For $$k = 1$$ step.

\label{eqn:multiphysicsL16:180}
x_{n+1} – x_n = \Delta t f( x_n, t_n ),

the coefficients are $$\alpha_{-1} = 1, \alpha_0 = -1, \beta_{-1} = 0, \beta_0 =1$$. For the simple function above

\label{eqn:multiphysicsL16:200}
\gamma_{-1} = \alpha_{-1} – \Delta t \lambda \beta_{-1} = 1

\label{eqn:multiphysicsL16:220}
\gamma_{0} = \alpha_{0} – \Delta t \lambda \beta_{0} = -1 – \Delta t \lambda.

The stability polynomial is

\label{eqn:multiphysicsL16:240}
1 z + \lr{ -1 – \Delta t \lambda} = 0,

or

\label{eqn:multiphysicsL16:260}
\boxed{
z = 1 + \delta t \lambda.
}

This is the root, or pole.

For stability we must have

\label{eqn:multiphysicsL16:280}
\Abs{ 1 + \Delta t \lambda } < 1,

or
\label{eqn:multiphysicsL16:300}
\Abs{ \lambda – \lr{ -\inv{\Delta t} } } < \inv{\Delta t},

This inequality is illustrated roughly in fig. 4.

fig. 4. Stability region of FE

All poles of my system must be inside the stability region in order to get stable $$\gamma$$.