phasor

Polarization review

January 24, 2015 ece1229 , , , , , , ,

[Click here for a PDF of this post with nicer formatting]

It seems worthwhile to review how a generally polarized field phasor leads to linear, circular, and elliptic geometries.

The most general field polarized in the \( x, y \) plane has the form

\begin{equation}\label{eqn:polarizationReview:20}
\BE
= \lr{ \xcap a e^{j \alpha} + \ycap b e^{j \beta} } e^{j \lr{ \omega t -k z }}
= \lr{ \xcap a e^{j \lr{\alpha – \beta}/2} + \ycap b e^{j \lr{ \beta – \alpha}/2} } e^{j \lr{ \omega t -k z + \lr{\alpha + \beta}/2 }}.
\end{equation}

Knowing to factor out the average phase angle above is only because I tried initially without that and things got ugly and messy. I guessed this would help (it does).

Let \( \boldsymbol{\mathcal{E}} = \text{Re} \BE = \xcap x + \ycap y \), \( \theta = \omega t + (\alpha + \beta)/2 \), and \( \phi = (\alpha – \beta)/2 \), so that

\begin{equation}\label{eqn:polarizationReview:40}
\BE
= \lr{ \xcap a e^{j \phi} + \ycap b e^{-j \phi} } e^{j \theta }.
\end{equation}

The coordinates can now be read off

\begin{equation}\label{eqn:polarizationReview:60}
\frac{x}{a} = \cos\phi \cos\theta – \sin\phi \sin\theta
\end{equation}
\begin{equation}\label{eqn:polarizationReview:80}
\frac{y}{b} = \cos\phi \cos\theta + \sin\phi \sin\theta,
\end{equation}

or in matrix form

\begin{equation}\label{eqn:polarizationReview:100}
\begin{bmatrix}
x/a \\
y/b \\
\end{bmatrix}
=
\begin{bmatrix}
\cos\phi & – \sin\phi \\
\cos\phi & \sin\phi
\end{bmatrix}
\begin{bmatrix}
\cos\theta \\
\sin\theta
\end{bmatrix}.
\end{equation}

The goal is to eliminate all the \( \theta \) (i.e. time dependence), converting the parametric relationship into a conic form.
Assuming that neither \( \cos\theta \), nor \( \sin\theta \) are zero for now (those are special cases and lead to linear polarization), inverting the matrix will allow the \( \theta \) dependence to be eliminated

\begin{equation}\label{eqn:polarizationReview:120}
\inv{\sin\lr{ 2\phi }}
\begin{bmatrix}
\sin\phi & \sin\phi \\
– \cos\phi & \cos\phi
\end{bmatrix}
\begin{bmatrix}
x/a \\
y/b \\
\end{bmatrix}
=
\begin{bmatrix}
\cos\theta \\
\sin\theta
\end{bmatrix}.
\end{equation}

Squaring and summing both rows of these equation gives
\begin{equation}\label{eqn:polarizationReview:140}
\begin{aligned}
1
&=
\inv{\sin^2 \lr{ 2\phi}}
\lr{
\sin^2\phi
\lr{
\frac{x}{a}
+\frac{y}{b}
}^2
+
\cos^2\phi
\lr{
-\frac{x}{a}
+\frac{y}{b}
}^2
} \\
&=
\inv{\sin^2 \lr{ 2\phi}}
\lr{
\frac{x^2}{a^2}
+\frac{y^2}{b^2}
+2 \frac{x y}{a b} \lr{ \sin^2\phi – \cos^2\phi }
} \\
&=
\inv{\sin^2 \lr{ 2\phi}}
\lr{
\frac{x^2}{a^2}
+\frac{y^2}{b^2}
-2 \frac{x y}{a b} \cos \lr{2\phi}
}
\end{aligned}
\end{equation}

Time to summarize and handle the special cases.

  1. To have \( \cos\phi = 0 \), the phase angles must satisfy \( \alpha – \beta = \lr{ 1 + 2 k } \pi, \, k \in \mathbb{Z} \).

    For this case \ref{eqn:polarizationReview:50} reduces to

    \begin{equation}\label{eqn:polarizationReview:160}
    -\frac{x}{a} = \frac{y}{b},
    \end{equation}

    which is just a line.

    Example.

    Let \( \alpha = 0, \beta = -\pi \), so that the phasor has the value

    \begin{equation}\label{eqn:polarizationReview:260}
    \BE = \lr{ \xcap a – \ycap b } e^{j \omega t}
    \end{equation}

  2. For have \( \sin\phi = 0 \), the phase angles must satisfy \( \alpha – \beta = 2 \pi k, \, k \in \mathbb{Z} \).

    For this case \ref{eqn:polarizationReview:60} and \ref{eqn:polarizationReview:80} reduce to

    \begin{equation}\label{eqn:polarizationReview:180}
    \frac{x}{a} = \frac{y}{b},
    \end{equation}

    also just a line.

    Example.

    Let \( \alpha = \beta = 0 \), so that the phasor has the value

    \begin{equation}\label{eqn:polarizationReview:280}
    \BE = \lr{ \xcap a + \ycap b } e^{j \omega t}
    \end{equation}

  3. Last is the circular and elliptically polarized case. The system is clearly elliptically polarized if \( \cos(2 \phi) = 0\), or \( \alpha – \beta = (\pi/2)( 1 + 2 k ), k \in \mathbb{Z}\). When that is the case and \( a = b \) also holds, the ellipse is a circle.

    When the \( \cos( 2 \phi) = 0 \) condition does not hold, a rotation of coordinates

    \begin{equation}\label{eqn:polarizationReview:200}
    \begin{bmatrix}
    x \\
    y
    \end{bmatrix}
    =
    \begin{bmatrix}
    \cos\mu & \sin\mu \\
    -\sin\mu & \cos\mu
    \end{bmatrix}
    \begin{bmatrix}
    u \\
    v
    \end{bmatrix}
    \end{equation}

    where

    \begin{equation}\label{eqn:polarizationReview:220}
    \mu = \inv{2} \tan^{-1} \lr{ \frac{ 2 \cos (\alpha – \beta)}{b – a}}
    \end{equation}

    puts the trajectory into a standard (but messy) conic form

    \begin{equation}\label{eqn:polarizationReview:240}
    1 = \frac{u^2}{ab} \lr{
    \frac{b}{a} \cos^2 \mu
    + \frac{a}{b} \sin^2 \mu
    + \inv{2} \sin\lr{2 \mu + \alpha – \beta}
    }
    +
    \frac{v^2}{ab} \lr{
    \frac{b}{a} \sin^2 \mu
    + \frac{a}{b} \cos^2 \mu
    – \inv{2} \sin\lr{2 \mu + \alpha – \beta}
    }
    \end{equation}

    It isn’t obvious to me that the factors of the \( u^2, v^2 \) terms are necessarily positive, which is required for the conic to be an ellipse and not a hyperbola.

    Circular polarization example.

    With \( a = b = E_0 \), \( \alpha = 0 \), \( \beta = \pm \pi/2 \), all the circular polarization conditions are met, leaving the phasor with values

    \begin{equation}\label{eqn:polarizationReview:300}
    \BE = E_0 \lr{ \xcap \pm j \ycap } e^{j \omega t}
    \end{equation}

Fundamental parameters of antennas

January 22, 2015 ece1229 , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,

[Click here for a PDF of this post with nicer formatting]

This is my first set of notes for the UofT course ECE1229, Advanced Antenna Theory, taught by Prof. Eleftheriades, covering ch. 2 [1] content.

Unlike most of the other classes I have taken, I am not attempting to take comprehensive notes for this class. The class is taught on slides that match the textbook so closely, there is little value to me taking notes that just replicate the text. Instead, I am annotating my copy of textbook with little details instead. My usual notes collection for the class will contain musings of details that were unclear, or in some cases, details that were provided in class, but are not in the text (and too long to pencil into my book.)

Poynting vector

The Poynting vector was written in an unfamiliar form

\begin{equation}\label{eqn:chapter2Notes:560}
\boldsymbol{\mathcal{W}} = \boldsymbol{\mathcal{E}} \cross \boldsymbol{\mathcal{H}}.
\end{equation}

I can roll with the use of a different symbol (i.e. not \(\BS\)) for the Poynting vector, but I’m used to seeing a \( \frac{c}{4\pi} \) factor ([6] and [5]). I remembered something like that in SI units too, so was slightly confused not to see it here.

Per [3] that something is a \( \mu_0 \), as in

\begin{equation}\label{eqn:chapter2Notes:580}
\boldsymbol{\mathcal{W}} = \inv{\mu_0} \boldsymbol{\mathcal{E}} \cross \boldsymbol{\mathcal{B}}.
\end{equation}

Note that the use of \( \boldsymbol{\mathcal{H}} \) instead of \( \boldsymbol{\mathcal{B}} \) is what wipes out the requirement for the \( \frac{1}{\mu_0} \) term since \( \boldsymbol{\mathcal{H}} = \boldsymbol{\mathcal{B}}/\mu_0 \), assuming linear media, and no magnetization.

Typical far-field radiation intensity

It was mentioned that

\begin{equation}\label{eqn:advancedantennaL1:20}
U(\theta, \phi)
=
\frac{r^2}{2 \eta_0} \Abs{ \BE( r, \theta, \phi) }^2
=
\frac{1}{2 \eta_0} \lr{ \Abs{ E_\theta(\theta, \phi) }^2 + \Abs{ E_\phi(\theta, \phi) }^2},
\end{equation}

where the intrinsic impedance of free space is

\begin{equation}\label{eqn:advancedantennaL1:480}
\eta_0 = \sqrt{\frac{\mu_0}{\epsilon_0}} = 377 \Omega.
\end{equation}

(this is also eq. 2-19 in the text.)

To get an understanding where this comes from, consider the far field radial solutions to the electric and magnetic dipole problems, which have the respective forms (from [3]) of

\begin{equation}\label{eqn:chapter2Notes:740}
\begin{aligned}
\boldsymbol{\mathcal{E}} &= -\frac{\mu_0 p_0 \omega^2 }{4 \pi } \frac{\sin\theta}{r} \cos\lr{w t – k r} \thetacap \\
\boldsymbol{\mathcal{B}} &= -\frac{\mu_0 p_0 \omega^2 }{4 \pi c} \frac{\sin\theta}{r} \cos\lr{w t – k r} \phicap \\
\end{aligned}
\end{equation}
\begin{equation}\label{eqn:chapter2Notes:760}
\begin{aligned}
\boldsymbol{\mathcal{E}} &= \frac{\mu_0 m_0 \omega^2 }{4 \pi c} \frac{\sin\theta}{r} \cos\lr{w t – k r} \phicap \\
\boldsymbol{\mathcal{B}} &= -\frac{\mu_0 m_0 \omega^2 }{4 \pi c^2} \frac{\sin\theta}{r} \cos\lr{w t – k r} \thetacap \\
\end{aligned}
\end{equation}

In neither case is there a component in the direction of propagation, and in both cases (using \( \mu_0 \epsilon_0 = 1/c^2\))

\begin{equation}\label{eqn:chapter2Notes:780}
\Abs{\boldsymbol{\mathcal{H}}}
= \frac{\Abs{\boldsymbol{\mathcal{E}}}}{\mu_0 c}
= \Abs{\boldsymbol{\mathcal{E}}} \sqrt{\frac{\epsilon_0}{\mu_0}}
= \inv{\eta_0}\Abs{\boldsymbol{\mathcal{E}}} .
\end{equation}

A superposition of the phasors for such dipole fields, in the far field, will have the form

\begin{equation}\label{eqn:chapter2Notes:800}
\begin{aligned}
\BE &= \inv{r} \lr{ E_\theta(\theta, \phi) \thetacap + E_\phi(\theta, \phi) \phicap } \\
\BB &= \inv{r c} \lr{ E_\theta(\theta, \phi) \thetacap – E_\phi(\theta, \phi) \phicap },
\end{aligned}
\end{equation}

with a corresponding time averaged Poynting vector

\begin{equation}\label{eqn:chapter2Notes:820}
\begin{aligned}
\BW_{\textrm{av}}
&= \inv{2 \mu_0} \BE \cross \BB^\conj \\
&=
\inv{2 \mu_0 c r^2}
\lr{ E_\theta \thetacap + E_\phi \phicap } \cross
\lr{ E_\theta^\conj \thetacap – E_\phi^\conj \phicap } \\
&=
\frac{\thetacap \cross \phicap}{2 \mu_0 c r^2}
\lr{ \Abs{E_\theta}^2 + \Abs{E_\phi}^2 } \\
&=
\frac{\rcap}{2 \eta_0 r^2}
\lr{ \Abs{E_\theta}^2 + \Abs{E_\phi}^2 },
\end{aligned}
\end{equation}

verifying \ref{eqn:advancedantennaL1:20} for a superposition of electric and magnetic dipole fields. This can likely be shown for more general fields too.

Field plots

We can plot the fields, or intensity (or log plots in dB of these).
It is pointed out in [3] that when there is \( r \) dependence these plots are done by considering the values of at fixed \( r \).

The field plots are conceptually the simplest, since that vector parameterizes
a surface. Any such radial field with magnitude \( f(r, \theta, \phi) \) can
be plotted in Mathematica in the \( \phi = 0 \) plane at \( r = r_0 \), or in
3D (respectively, but also at \( r = r_0\)) with code like that of the
following listing

ParametricPlotListing

Intensity plots can use the same code, with the only difference being the interpretation. The surface doesn’t represent the value of a vector valued radial function, but is the magnitude of a scalar valued function evaluated at \( f( r_0, \theta, \phi) \).

The surfaces for \( U = \sin\theta, \sin^2\theta \) in the plane are parametrically plotted in fig. 2, and for cosines in fig. 1 to compare with textbook figures.

CoSineAndCoSineSqFig1pn

fig 1. Cosinusoidal radiation intensities

SineAndSinSqFig3pn

fig 2. Sinusoidal radiation intensities

 

Visualizations of \( U = \sin^2 \theta\) and \( U = \cos^2 \theta\) can be found in fig. 3 and fig. 4 respectively. Even for such simple functions these look pretty cool.

SineSq3DFig4pn

fig 3. Square sinusoidal radiation intensity

 

CoSineSq3DFig2pn

fig 4. Square cosinusoidal radiation intensity

 

dB vs dBi

Note that dBi is used to indicate that the gain is with respect to an “isotropic” radiator.
This is detailed more in [2].

Trig integrals

Tables 1.1 and 1.2 produced with tableOfTrigIntegrals.nb have some of the sine and cosine integrals that are pervasive in this chapter.

trigIntegralsUpToPiBy2

trigIntegralsUpToPi

Polarization vectors

The text introduces polarization vectors \( \rhocap \) , but doesn’t spell out their form. Consider a plane wave field of the form

\begin{equation}\label{eqn:chapter2Notes:840}
\BE
=
E_x e^{j \phi_x} e^{j \lr{ \omega t – k z }} \xcap
+
E_y e^{j \phi_y} e^{j \lr{ \omega t – k z }} \ycap.
\end{equation}

The \( x, y \) plane directionality of this phasor can be written

\begin{equation}\label{eqn:chapter2Notes:860}
\Brho =
E_x e^{j \phi_x} \xcap
+
E_y e^{j \phi_y} \ycap,
\end{equation}

so that

\begin{equation}\label{eqn:chapter2Notes:880}
\BE = \Brho e^{j \lr{ \omega t – k z }}.
\end{equation}

Separating this direction and magnitude into factors

\begin{equation}\label{eqn:chapter2Notes:900}
\Brho = \Abs{\BE} \rhocap,
\end{equation}

allows the phasor to be expressed as

\begin{equation}\label{eqn:chapter2Notes:920}
\BE = \rhocap \Abs{\BE} e^{j \lr{ \omega t – k z }}.
\end{equation}

As an example, suppose that \( E_x = E_y \), and set \( \phi_x = 0 \). Then

\begin{equation}\label{eqn:chapter2Notes:940}
\rhocap = \xcap + \ycap e^{j \phi_y}.
\end{equation}

Phasor power

In section 2.13 the phasor power is written as

\begin{equation}\label{eqn:chapter2Notes:620}
I^2 R/2,
\end{equation}

where \( I, R \) are the magnitudes of phasors in the circuit.

I vaguely recall this relation, but had to refer back to [4] for the details.
This relation expresses average power over a period associated with the frequency of the phasor

\begin{equation}\label{eqn:chapter2Notes:640}
\begin{aligned}
P
&= \inv{T} \int_{t_0}^{t_0 + T} p(t) dt \\
&= \inv{T} \int_{t_0}^{t_0 + T} \Abs{\BV} \cos\lr{ \omega t + \phi_V }
\Abs{\BI} \cos\lr{ \omega t + \phi_I} dt \\
&= \inv{T} \int_{t_0}^{t_0 + T} \Abs{\BV} \Abs{\BI}
\lr{
\cos\lr{ \phi_V – \phi_I } + \cos\lr{ 2 \omega t + \phi_V + \phi_I}
}
dt \\
&= \inv{2} \Abs{\BV} \Abs{\BI} \cos\lr{ \phi_V – \phi_I }.
\end{aligned}
\end{equation}

Introducing the impedance for this circuit element

\begin{equation}\label{eqn:chapter2Notes:660}
\BZ = \frac{ \Abs{\BV} e^{j\phi_V} }{ \Abs{\BI} e^{j\phi_I} } = \frac{\Abs{\BV}}{\Abs{\BI}} e^{j\lr{\phi_V – \phi_I}},
\end{equation}

this average power can be written in phasor form

\begin{equation}\label{eqn:chapter2Notes:680}
\BP = \inv{2} \Abs{\BI}^2 \BZ,
\end{equation}

with
\begin{equation}\label{eqn:chapter2Notes:700}
P = \textrm{Re} \BP.
\end{equation}

Observe that we have to be careful to use the absolute value of the current phasor \( \BI \), since \( \BI^2 \) differs in phase from \( \Abs{\BI}^2 \). This explains the conjugation in the [4] definition of complex power, which had the form

\begin{equation}\label{eqn:chapter2Notes:720}
\BS = \BV_{\textrm{rms}} \BI^\conj_{\textrm{rms}}.
\end{equation}

Radar cross section examples

Flat plate.

\begin{equation}\label{eqn:chapter2Notes:960}
\sigma_{\textrm{max}} = \frac{4 \pi \lr{L W}^2}{\lambda^2}
\end{equation}

RCSsquareGeometryFig1

fig. 6. Square geometry for RCS example.

 

Sphere.

In the optical limit the radar cross section for a sphere

RCSsphereGeometryFig3

fig. 7. Sphere geometry for RCS example.

 

\begin{equation}\label{eqn:chapter2Notes:980}
\sigma_{\textrm{max}} = \pi r^2
\end{equation}

Note that this is smaller than the physical area \( 4 \pi r^2 \).

Cylinder.

RCScylinderGeometryFig1

fig. 8. Cylinder geometry for RCS example.

 

\begin{equation}\label{eqn:chapter2Notes:1000}
\sigma_{\textrm{max}} = \frac{ 2 \pi r h^2}{\lambda}
\end{equation}

Tridedral corner reflector

trihedralCornerReflectorFig6

fig. 9. Trihedral corner reflector geometry for RCS example.

 

\begin{equation}\label{eqn:chapter2Notes:1020}
\sigma_{\textrm{max}} = \frac{ 4 \pi L^4}{3 \lambda^2}
\end{equation}

Scattering from a sphere vs frequency

Frequency dependence of spherical scattering is sketched in fig. 10.

  • Low frequency (or small particles): Rayleigh\begin{equation}\label{eqn:chapter2Notes:1040}
    \sigma = \lr{\pi r^2} 7.11 \lr{\kappa r}^4, \qquad \kappa = 2 \pi/\lambda.
    \end{equation}
  • Mie scattering (resonance),\begin{equation}\label{eqn:chapter2Notes:1060}
    \sigma_{\textrm{max}}(A) = 4 \pi r^2
    \end{equation}
    \begin{equation}\label{eqn:chapter2Notes:1080}
    \sigma_{\textrm{max}}(B) = 0.26 \pi r^2.
    \end{equation}
  • optical limit ( \(r \gg \lambda\) )\begin{equation}\label{eqn:chapter2Notes:1100}
    \sigma = \pi r^2.
    \end{equation}
sphericalScatteringFig5

fig 10. Scattering from a sphere vs frequency (from Prof. Eleftheriades’ class notes).

FIXME: Do I have a derivation of this in my optics notes?

Notation

  • Time average.
    Both Prof. Eleftheriades
    and the text [1] use square brackets \( [\cdots] \) for time averages, not \( <\cdots> \). Was that an engineering convention?
  • Prof. Eleftheriades
    writes \(\Omega\) as a circle floating above a face up square bracket, as in fig. 1, and \( \sigma \) like a number 6, as in fig. 1.
  • Bold vectors are usually phasors, with (bold) calligraphic script used for the time domain fields. Example: \( \BE(x,y,z,t) = \ecap E(x,y) e^{j \lr{\omega t – k z}}, \boldsymbol{\mathcal{E}}(x, y, z, t) = \textrm{Re} \BE \).
greekStyleOmegaFig1

fig. 11. Prof. handwriting decoder ring: Omega

sigmaFig1

fig 12. Prof. handwriting decoder ring: sigma

 

References

[1] Constantine A Balanis. Antenna theory: analysis and design. John Wiley \& Sons, 3rd edition, 2005.

[3] David Jeffrey Griffiths and Reed College. Introduction to electrodynamics. Prentice hall Upper Saddle River, NJ, 3rd edition, 1999.

[4] J.D. Irwin. Basic Engineering Circuit Analysis. MacMillian, 1993.

[5] JD Jackson. Classical Electrodynamics. John Wiley and Sons, 2nd edition, 1975.

[6] L.D. Landau and E.M. Lifshitz. The classical theory of fields. Butterworth-Heinemann, 1980. ISBN 0750627689.

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

December 9, 2014 ece1254 , , , , , ,

[Click here for a PDF of this post with nicer formatting]

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

\begin{equation}\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}) \\
\end{bmatrix}, \qquad
\Bu_a =
\begin{bmatrix}
u_b(t_{-N}) \\
u_b(t_{1-N}) \\
\vdots \\
u_b(t_{N-1}) \\
u_b(t_{N}) \\
\end{bmatrix}, \qquad
\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},
\end{equation}

and Fourier component vectors

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

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

\begin{equation}\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},
\end{equation}

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

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

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

\begin{equation}\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}
\end{equation}

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

\begin{equation}\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}
\end{equation}

With a vector of fourier coeffient vectors

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

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

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

the complete block diagonalization is

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

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

\begin{equation}\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}
\end{equation}

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

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

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

\begin{equation}\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}.
\end{equation}

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

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

If \( N = 2 \) was used, then we would have instead

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

Representing non-linear sources

The time domain MNA fig. 1.

diodeWithResistorFig1

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

\begin{equation}\label{eqn:diodeRLCSample:580}
\BG =
\begin{bmatrix}
Z_s & -Z_s \\
-Z_s & Z_s + Z_g \\
\end{bmatrix}, \qquad
\Bx =
\begin{bmatrix}
v_1(t) \\
v_2(t) \\
\end{bmatrix}
\end{equation}
\begin{equation}\label{eqn:diodeRLCSample:600}
\BB =
\begin{bmatrix}
I^{(s)}/2 & -I^{(0)} & I^{(s)}/2 \\
-I^{(s)}/2 & I^{(0)} & -I^{(s)}/2
\end{bmatrix}, \qquad
\Bu(t) =
\begin{bmatrix}
e^{-j \omega_0 t} \\
1 \\
e^{j \omega_0 t}
\end{bmatrix}
\end{equation}
\begin{equation}\label{eqn:diodeRLCSample:620}
\BD
=
\begin{bmatrix}
I^{(0)} \\
-I^{(0)}
\end{bmatrix}, \qquad
\Bw(t) = e^{(v_1(t) – v_2(t))/V_T}.
\end{equation}

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

\begin{equation}\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}
\end{equation}

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

\begin{equation}\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}.
\end{equation}

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

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

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

\begin{equation}\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},
\end{equation}

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

\begin{equation}\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}.
\end{equation}

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

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

The DFT inverse is

\begin{equation}\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 }.
\end{equation}

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

\begin{equation}\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}.
\end{equation}

With the introduction a term by term exponentiation operator

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

this one Fourier coefficient vector is

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

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

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

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

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

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

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

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

\begin{equation}\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 }.
\end{equation}

The partials now follow immediately

\begin{equation}\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}
}.
\end{equation}

Generalization

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

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

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

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

so

\begin{equation}\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 }.
\end{equation}

The Jacobian partials are

\begin{equation}\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}
}.
\end{equation}

Note that this Jacobian

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

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

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

\begin{equation}\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}.
\end{equation}

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

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

References

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

A sample diode RLC circuit

December 5, 2014 ece1254 , , , , , ,

[Click here for a PDF of this post with nicer formatting]

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.

diodeRLCFig1

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

\begin{equation}\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}
\end{equation}

Let’s write this as

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

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

\begin{equation}\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}
\end{equation}

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

\begin{equation}\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}
}
\end{equation}

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.