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

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