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

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