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

## Disclaimer

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

## Nonlinear differential equations

Assume that the relationships between the zeroth and first order derivatives has the form

\label{eqn:multiphysicsL15:20}
F\lr{ x(t), \dot{x}(t) } = 0

\label{eqn:multiphysicsL15:40}
x(0) = x_0

The backward Euler method where the derivative approximation is

\label{eqn:multiphysicsL15:60}
\dot{x}(t_n) \approx \frac{x_n – x_{n-1}}{\Delta t},

can be used to solve this numerically, reducing the problem to

\label{eqn:multiphysicsL15:80}
F\lr{ x_n, \frac{x_n – x_{n-1}}{\Delta t} } = 0.

This can be solved with Newton’s method. How do we find the initial guess for Newton’s? Consider a possible system in fig. 1.

fig. 1. Possible solution points

One strategy for starting each iteration of Newton’s method is to base the initial guess for $$x_1$$ on the value $$x_0$$, and do so iteratively for each subsequent point. One can imagine that this may work up to some sample point $$x_n$$, but then break down (i.e. Newton’s diverges when the previous value $$x_{n-1}$$ is used to attempt to solve for $$x_n$$). At that point other possible strategies may work. One such strategy is to use an approximation of the derivative from the previous steps to attempt to get a better estimate of the next value. Another possibility is to reduce the time step, so the difference between successive points is reduced.

## Analysis, accuracy and stability ($$\Delta t \rightarrow 0$$)

Consider a differential equation

\label{eqn:multiphysicsL15:100}
\dot{x}(t) = f(x(t), t)

\label{eqn:multiphysicsL15:120}
x(t_0) = x_0

A few methods of solution have been considered

• (FE) $$x_{n+1} – x_n = \Delta t f(x_n, t_n)$$
• (BE) $$x_{n+1} – x_n = \Delta t f(x_{n+1}, t_{n+1})$$
• (TR) $$x_{n+1} – x_n = \frac{\Delta t}{2} f(x_{n+1}, t_{n+1}) + \frac{\Delta t}{2} f(x_{n}, t_{n})$$

A common pattern can be observed, the generalization of which are called
\textit{linear multistep methods}
(LMS), which have the form

\label{eqn:multiphysicsL15:140}
\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} )

The FE (explicit), BE (implicit), and TR methods are now special cases with

• (FE) $$\alpha_{-1} = 1, \alpha_0 = -1, \beta_{-1} = 0, \beta_0 = 1$$
• (BE) $$\alpha_{-1} = 1, \alpha_0 = -1, \beta_{-1} = 1, \beta_0 = 0$$
• (TR) $$\alpha_{-1} = 1, \alpha_0 = -1, \beta_{-1} = 1/2, \beta_0 = 1/2$$

Here $$k$$ is the number of timesteps used. The method is explicit if $$\beta_{-1} = 0$$.

### Definition: Convergence

With

• $$x(t)$$ : exact solution
• $$x_n$$ : computed solution
• $$e_n$$ : where $$e_n = x_n – x(t_n)$$, is the global error

The LMS method is convergent if

\begin{equation*}%\label{eqn:multiphysicsL15:180}
\max_{n, \Delta t \rightarrow 0} \Abs{ x_n – t(t_n) } \rightarrow 0 %\xrightarrow[t \rightarrow 0 ]{} 0
\end{equation*}

Convergence: zero-stability and consistency (small local errors made at each iteration),

where zero-stability is “small sensitivity to changes in initial condition”.

### Definition: Consistency

A local error $$R_{n+1}$$ can be defined as

\begin{equation*}%\label{eqn:multiphysicsL15:220}
R_{n+1} = \sum_{j = -1}^{k-1} \alpha_j x(t_{n-j}) – \Delta t \sum_{j=-1}^{k-1} \beta_j f(x(t_{n-j}), t_{n-j}).
\end{equation*}

The method is consistent if

\begin{equation*}%\label{eqn:multiphysicsL15:240}
\lim_{\Delta t} \lr{
\max_n \Abs{ \inv{\Delta t} R_{n+1} } = 0
}
\end{equation*}

or $$R_{n+1} \sim O({\Delta t}^2)$$

## Disclaimer

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

## Struts and Joints, Node branch formulation

Let’s consider the simple strut system of fig. 1 again.

fig. 1. Simple strut system

Our unknowns are

1. Forces at each of the points we have a force with two components
\label{eqn:multiphysicsL12:20}
\Bf_A = \lr{ f_{A,x}, f_{A,y} }
We construct a total force vector\label{eqn:multiphysicsL12:40}
\Bf =
\begin{bmatrix}
f_{A,x} \\
f_{A,y} \\
f_{B,x} \\
f_{B,y} \\
\vdots
\end{bmatrix}
2. Positions of the joints\label{eqn:multiphysicsL12:60}
\Br =
\begin{bmatrix}
x_1 \\
y_1 \\
y_1 \\
y_2 \\
\vdots
\end{bmatrix}

Our given variables are

1. The load force $$\Bf_L$$.
2. The joint positions at rest.
3. parameter of struts.

### Conservation laws

The conservation laws are

\label{eqn:multiphysicsL12:80}
\Bf_A + \Bf_B + \Bf_C = 0

\label{eqn:multiphysicsL12:100}
-\Bf_C + \Bf_D + \Bf_L = 0

which we can put in matrix form

\label{eqn:multiphysicsL12:120}
\begin{bmatrix}
1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\
0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & -1 & 0 & 1 & 0 \\
0 & 1 & 0 & 1 & 0 & -1 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix}
f_{A,x} \\
f_{A,y} \\
f_{B,x} \\
f_{B,y} \\
f_{C,x} \\
f_{C,y} \\
f_{D,x} \\
f_{D,y}
\end{bmatrix}
=
\begin{bmatrix}
0 \\
0 \\
-f_{L,x} \\
-f_{L,y}
\end{bmatrix}

Here the block matrix is called the incidence matrix $$\BA$$, and we write

\label{eqn:multiphysicsL12:160}
A \Bf = \Bf_L.

### Constitutive laws

Given a pair of nodes as in fig. 2.

fig. 2. Strut spanning nodes

each component has an equation relating the reaction forces of that strut based on the positions

\label{eqn:multiphysicsL12:180}
f_{c,x} = S_x \lr{ x_1 – x_2, y_1 – y_2, p_c }

\label{eqn:multiphysicsL12:200}
f_{c,y} = S_y \lr{ x_1 – x_2, y_1 – y_2, p_c },

where $$p_c$$ represent the parameters of the system. We write

\label{eqn:multiphysicsL12:220}
\Bf =
\begin{bmatrix}
f_{A,x} \\
f_{A,y} \\
f_{B,x} \\
f_{B,y} \\
\vdots
\end{bmatrix}
=
\begin{bmatrix}
S_x \lr{ x_1 – x_2, y_1 – y_2, p_c } \\
S_y \lr{ x_1 – x_2, y_1 – y_2, p_c } \\
\vdots
\end{bmatrix},

or

\label{eqn:multiphysicsL12:260}
\Bf = S(\Br)

### Putting the pieces together

The node branch formulation is

\label{eqn:multiphysicsL12:280}
\begin{aligned}
A \Bf – \Bf_L &= 0 \\
\Bf – S(\Br) &= 0
\end{aligned}

We’ll want to approximate this system using the Jacobian methods discussed, and can expect the cost of that Jacobian calculation to potentially be expensive. To move to the nodal formulation we eliminate forces (the equivalent of currents in this system)

\label{eqn:multiphysicsL12:320}
A S(\Br) – \Bf_L = 0

We cannot use this nodal formulation when we have struts that are so stiff that the positions of some of the nodes are fixed, but can work around that as before by introducing an additional unknown for each component of such a strut.

## Damped Newton’s method

We want to be able to deal with the oscillation that we can have in examples like that of fig. 3.

fig. 3. Oscillatory Newton’s iteration

Large steps can be dangerous. We want to modify Newton’s method as follows

Our algorithm is

• Guess $$\Bx^0, k = 0$$.
• REPEAT
• Compute $$F$$ and $$J_F$$ at $$\Bx^k$$
• Solve linear system $$J_F(\Bx^k) \Delta \Bx^k = – F(\Bx^k)$$
• $$\Bx^{k+1} = \Bx^k + \alpha^k \Delta \Bx^k$$
• $$k = k + 1$$
• UNTIL converged

with $$\alpha^k = 1$$ we have standard Newton’s method. We want to pick $$\alpha^k$$ so that we minimize

## Continuation parameters

Newton’s method converges given a close initial guess. We can generate a sequence of problems where the previous problem generates a good initial guess for the next problem.

An example is a heat conducting bar, with a final heat distribution. We can start the numeric iteration with $$T = 0$$, and gradually increase the temperatures until we achieve the final desired heat distribution.

Suppose that we want to solve

\label{eqn:multiphysicsL12:340}
F(\Bx) = 0.

We modify this problem by introducing

\label{eqn:multiphysicsL12:360}
\tilde{F}(\Bx(\lambda), \lambda) = 0,

where

• $$\tilde{F}(\Bx(0), 0) = 0$$ is easy to solve
• $$\tilde{F}(\Bx(1), 1) = 0$$ is equivalent to $$F(\Bx) = 0$$.
• (more on slides)

The source load stepping algorithm is

• Solve $$\tilde{F}(\Bx(0), 0) = 0$$, with $$\Bx(\lambda_{\text{prev}} = \Bx(0)$$
• (more on slides)

This can still have problems, for example, when the parameterization is multivalued as in fig. 4.

fig. 4. Multivalued parameterization

We can attempt to adjust $$\lambda$$ so that we move along the parameterization curve.

## ECE1254H Modeling of Multiphysics Systems. Lecture 10: Nonlinear systems. Taught by Prof. Piero Triverio

October 21, 2014 ece1254 ,

## Disclaimer

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

## Nonlinear systems

On slides, some examples to motivate:

• struts
• fluids
• diode (exponential)Example in fig. 1.

fig. 1. Diode circuit

\label{eqn:multiphysicsL10:20}
I_d = I_s \lr{ e^{V_d/V_t} – 1 } = \frac{10 – V_d}{10}.

## Richardson and Linear Convergence

Seeking the exact solution $$x^\conj$$ for

\label{eqn:multiphysicsL10:40}
f(x^\conj) = 0,

Suppose that

\label{eqn:multiphysicsL10:60}
x^{k + 1} = x^k + f(x^k)

If $$f(x^k) = 0$$ then we have convergence, so $$x^k = x^\conj$$.

### Convergence analysis

Write the iteration equations at a sample point and the solution as

\label{eqn:multiphysicsL10:80}
x^{k + 1} = x^k + f(x^k)

\label{eqn:multiphysicsL10:100}
x^{\conj} = x^\conj +
\underbrace{
f(x^\conj)
}_{=0}

Taking the difference we have

\label{eqn:multiphysicsL10:120}
x^{k+1} – x^\conj = x^k – x^\conj + \lr{ f(x^k) – f(x^\conj) }.

The last term can be quantified using the mean value theorem \ref{thm:multiphysicsL10:140}, giving

\label{eqn:multiphysicsL10:140}
x^{k+1} – x^\conj
= x^k – x^\conj +
\evalbar{\PD{x}{f}}{\tilde{x}} \lr{ x^k – x^\conj }
=
\lr{ x^k – x^\conj }
\lr{
1 + \evalbar{\PD{x}{f}}{\tilde{x}} }.

The absolute value is thus

\label{eqn:multiphysicsL10:160}
\Abs{x^{k+1} – x^\conj } =
\Abs{ x^k – x^\conj }
\Abs{
1 + \evalbar{\PD{x}{f}}{\tilde{x}} }.

We have convergence provided $$\Abs{ 1 + \evalbar{\PD{x}{f}}{\tilde{x}} } < 1$$ in the region where we happen to iterate over. This could easily be highly dependent on the initial guess.

Stated more accurately we have convergence provided

\label{eqn:multiphysicsL10:180}
\Abs{
1 + \PD{x}{f} }
\le \gamma < 1
\qquad \forall \tilde{x} \in [ x^\conj – \delta, x^\conj + \delta ],

and $$\Abs{x^0 – x^\conj } < \delta$$. This is illustrated in fig. 3.

fig. 3. Convergence region

It could very easily be difficult to determine the convergence regions.

We have some problems

• Convergence is only linear
• $$x, f(x)$$ are not in the same units (and potentially of different orders). For example, $$x$$ could be a voltage and $$f(x)$$ could be a circuit current.
• (more on slides)

Examples where we may want to use this:

• Spice Gummal Poon transistor model. Lots of diodes, …
• Mosfet model (30 page spec, lots of parameters).

## Newton’s method

The core idea of this method is sketched in fig. 4. To find the intersection with the x-axis, we follow the slope closer to the intersection.

fig. 4. Newton’s method

To do this, we expand $$f(x)$$ in Taylor series to first order around $$x^k$$, and then solve for $$f(x) = 0$$ in that approximation

\label{eqn:multiphysicsL10:200}
f( x^{k+1} ) \approx f( x^k ) + \evalbar{ \PD{x}{f} }{x^k} \lr{ x^{k+1} – x^k } = 0.

This gives

\label{eqn:multiphysicsL10:220}
\boxed{x^{k+1} = x^k – \frac{f( x^k )}{\evalbar{ \PD{x}{f} }{x^k}}}

### Example: Newton’s method

For the solution of

\label{eqn:multiphysicsL10:260}
f(x) = x^3 – 2,

it was found (table 1.1)

The error tails off fast as illustrated roughly in fig. 6.

fig. 6. Error by iteration

### Convergence analysis

The convergence condition is

\label{eqn:multiphysicsL10:280}
0 = f(x^k) + \evalbar{ \PD{x}{f} }{x^k} \lr{ x^{k+1} – x^k }.

The Taylor series for $$f$$ around $$x^k$$, using a mean value formulation is

\label{eqn:multiphysicsL10:300}
f(x)
= f(x^k)
+ \evalbar{ \PD{x}{f} }{x^k} \lr{ x – x^k }.
+ \inv{2} \evalbar{ \PDSq{x}{f} }{\tilde{x} \in [x^\conj, x^k]} \lr{ x – x^k }^2.

Evaluating at $$x^\conj$$ we have

\label{eqn:multiphysicsL10:320}
0 = f(x^k)
+ \evalbar{ \PD{x}{f} }{x^k} \lr{ x^\conj – x^k }.
+ \inv{2} \evalbar{ \PDSq{x}{f} }{\tilde{x} \in [x^\conj, x^k]} \lr{ x^\conj – x^k }^2.

and subtracting this from \ref{eqn:multiphysicsL10:280} we are left with

\label{eqn:multiphysicsL10:340}
0 = \evalbar{\PD{x}{f}}{x^k} \lr{ x^{k+1} – x^k – x^\conj + x^k }
– \inv{2} \evalbar{\PDSq{x}{f}}{\tilde{x}} \lr{ x^\conj – x^k }^2.

Solving for the difference from the solution, the error is

\label{eqn:multiphysicsL10:360}
x^{k+1} – x^\conj
= \inv{2} \lr{ \PD{x}{f} }^{-1} \evalbar{\PDSq{x}{f}}{\tilde{x}} \lr{ x^k – x^\conj }^2,

or in absolute value

\label{eqn:multiphysicsL10:380}
\Abs{ x^{k+1} – x^\conj }
= \inv{2} \Abs{ \PD{x}{f} }^{-1} \Abs{ \PDSq{x}{f} } \Abs{ x^k – x^\conj }^2.

We see that convergence is quadratic in the error from the previous iteration. We will have trouble if the derivative goes small at any point in the iteration region, for example in fig. 7, we could easily end up in the zero derivative region.

fig. 7. Newton’s method with small derivative region

### When to stop iteration

One way to check is to look to see if the difference

\label{eqn:multiphysicsL10:420}
\Norm{ x^{k+1} – x^k } < \epsilon_{\Delta x},

however, when the function has a very step slope

this may not be sufficient unless we also substitute our trial solution and see if we have the match desired.

Alternatively, if the slope is shallow as in fig. 7, then checking for just $$\Abs{ f(x^{k+1} } < \epsilon_f$$ may also mean we are off target.

Finally, we may also need a relative error check to avoid false convergence. In fig. 10, we may have both

fig. 10. Possible relative error difference required

\label{eqn:multiphysicsL10:440}
\Abs{x^{k+1} – x^k} < \epsilon_{\Delta x}

\label{eqn:multiphysicsL10:460}
\Abs{f(x^{k+1}) } < \epsilon_{f},

however, we may also want a small relative difference

\label{eqn:multiphysicsL10:480}
\frac{\Abs{x^{k+1} – x^k}}{\Abs{x^k}}
< \epsilon_{f,r}.

This can become problematic in real world engineering examples such as to diode of fig. 11, where we have shallow regions and fast growing or dropping regions.

fig. 11. Diode current curve

## Theorems

### Mean value theorem

For a continuous and differentiable function $$f(x)$$, the difference can be expressed in terms of the derivative at an intermediate point

\begin{equation*}
f(x_2) – f(x_1)
= \evalbar{ \PD{x}{f} }{\tilde{x}} \lr{ x_2 – x_1 }
\end{equation*}

where $$\tilde{x} \in [x_1, x_2]$$.

This is illustrated (roughly) in fig. 2.

fig. 2. Mean value theorem illustrated