multiphysicsL10.tex

%
% Copyright © 2014 Peeter Joot.  All Rights Reserved.
% Licenced as described in the file LICENSE under the root directory of this GIT repository.
%
% for template copy, run:
%
% ~/bin/ct multiphysicsL1  multiphysicsLectureN tl1
%
%\input{../blogpost.tex}
%\renewcommand{\basename}{multiphysicsL10}
%\renewcommand{\dirname}{notes/ece1254/}
%\newcommand{\keywords}{ECE1254H}
%\input{../peeter_prologue_print2.tex}
%
%%\usepackage{kbordermatrix}
%
%\beginArtNoToc
%\generatetitle{ECE1254H Modeling of Multiphysics Systems.  Lecture 10: Nonlinear systems.  Taught by Prof.\ Piero Triverio}
%\chapter{Nonlinear systems}
\label{chap:multiphysicsL10}
%\section{Disclaimer}
%
%Peeter's lecture notes from class.  These may be incoherent and rough.
%
\section{Nonlinear systems.}
\index{nonlinear systems}
On slides, some examples to motivate:
\begin{itemize}
\item struts
\item fluids
\item diode (exponential)

Example in \cref{fig:lecture10:lecture10Fig1}.
%\imageFigure{../figures/ece1254-multiphysics/lecture10Fig1}{Diode circuit}{fig:lecture10:lecture10Fig1}{0.2}
\imageFigure{../figures/ece1254-multiphysics/diode-circuit-voltage-source.pdf}{Diode circuit.}{fig:lecture10:lecture10Fig1}{0.2}

\begin{equation}\label{eqn:multiphysicsL10:20}
I_{\textrm{d}} = I_{\textrm{s}} \lr{ e^{V_{\textrm{d}}/V_{\textrm{T}}} - 1 } = \frac{10 - V_{\textrm{d}}}{10}.
\end{equation}
\end{itemize}
\section{Richardson and Linear Convergence.}
Seeking the exact solution \( x^\conj \) for
\begin{equation}\label{eqn:multiphysicsL10:40}
f(x^\conj) = 0,
\end{equation}

Suppose that
\begin{equation}\label{eqn:multiphysicsL10:60}
x^{k + 1} = x^k + f(x^k).
\end{equation}

If \( f(x^k) = 0 \) then the iteration has converged at \( x^k = x^\conj \).
\paragraph{Convergence analysis}
\index{Richardson!convergence}
Write the iteration equations at a sample point and the solution as
\begin{subequations}
\begin{equation}\label{eqn:multiphysicsL10:80}
x^{k + 1} = x^k + f(x^k)
\end{equation}
\begin{equation}\label{eqn:multiphysicsL10:100}
x^{\conj} = x^\conj +
\mathLabelBox
{
f(x^\conj)
}
{\(=0\)}.
\end{equation}
\end{subequations}

The difference is
\begin{equation}\label{eqn:multiphysicsL10:120}
x^{k+1} - x^\conj = x^k - x^\conj + \lr{ f(x^k) - f(x^\conj) }.
\end{equation}

The last term can be quantified using the mean value theorem \ref{thm:multiphysicsL10:140}, giving
\begin{equation}\label{eqn:multiphysicsL10:140}
\begin{aligned}
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}}  }.
\end{aligned}
\end{equation}

The absolute value is thus
\begin{equation}\label{eqn:multiphysicsL10:160}
\Abs{x^{k+1} - x^\conj } =
\Abs{ x^k - x^\conj }
\Abs{
1 + \evalbar{\PD{x}{f}}{\tilde{x}}  }.
\end{equation}

Convergence is obtained when \( \Abs{ 1 + \evalbar{\PD{x}{f}}{\tilde{x}}  } < 1 \) in the iteration region.  This could easily be highly dependent on the initial guess.

A more accurate statement is that convergence is obtained provided
\begin{equation}\label{eqn:multiphysicsL10:180}
\Abs{
1 + \PD{x}{f} }
\le \gamma < 1
\qquad \forall \tilde{x} \in [ x^\conj - \delta, x^\conj + \delta ],
\end{equation}
and \( \Abs{x^0 - x^\conj } < \delta \).  This is illustrated in \cref{fig:lecture10:lecture10Fig3}.
\imageFigure{../figures/ece1254-multiphysics/lecture10Fig3}{Convergence region.}{fig:lecture10:lecture10Fig3}{0.2}

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

There are some problems
\begin{itemize}
\item Convergence is only linear
\item \( 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.
\item (more on slides)
\end{itemize}

Examples where this may be desirable include
\begin{itemize}
\item Spice Gummal Poon transistor model.  Lots of diodes, ...
\item Mosfet model (30 page spec, lots of parameters).
\end{itemize}
\section{Newton's method.}
\index{Newton's method!single variable}
The core idea of this method is sketched in \cref{fig:lecture10:lecture10Fig4.}.  To find the intersection with the x-axis, follow the slope closer to the intersection.
\imageFigure{../figures/ece1254-multiphysics/lecture10Fig4}{Newton's method.}{fig:lecture10:lecture10Fig4}{0.2}

To do this, expand \( f(x) \) in Taylor series to first order around \( x^k \), and then solve for \( f(x) = 0 \) in that approximation
\begin{equation}\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.
\end{equation}

This gives
\boxedEquation{eqn:multiphysicsL10:220}{
x^{k+1} = x^k - \frac{f( x^k )}{\evalbar{ \PD{x}{f} }{x^k}}.
}
\index{Newton's method}
\makeexample{Newton's method.}{example:multiphysicsL10:240}{
For the solution of
\begin{equation}\label{eqn:multiphysicsL10:260}
f(x) = x^3 - 2,
\end{equation}
it was found

\captionedTable{Numerical Newton's method example.}{tab:multiphysicsL10:400}{
\begin{tabular}{|l|l|l|}
\hline
k      & \( x^k \) & \( \Abs{ x^k - x^\conj} \) \\
\hline
0      & 10                        & 8.74                                                      \\
1      & \( \cdot \)                    & 5.4                                                       \\
\( \vdots \) & \( \cdot \)                  & \( \vdots \)                                                 \\
8      & \( \cdot \)                  & \( 1.7 \times 10^{-3}  \)                  \\
9      & \( \cdot  \)                 & \( 2.4 \times 10^{-6}   \)                 \\
10     & \( \cdot   \)                & \( 6.6 \times 10^{-12}   \)                \\
\hline
\end{tabular}
}

The error tails off fast as illustrated roughly in \cref{fig:lecture10:lecture10Fig6}.
\imageFigure{../figures/ece1254-multiphysics/lecture10Fig6}{Error by iteration.}{fig:lecture10:lecture10Fig6}{0.1}
}
%
\paragraph{Convergence analysis}
\index{Newton's method!convergence}
%
The convergence condition is
%
\begin{equation}\label{eqn:multiphysicsL10:280}
0 = f(x^k) + \evalbar{ \PD{x}{f} }{x^k} \lr{ x^{k+1} - x^k }.
\end{equation}
%
The Taylor series for \( f \) around \( x^k \), using a mean value formulation is
%
\begin{equation}\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.
\end{equation}
%
Evaluating at \( x^\conj \) gives
%
\begin{equation}\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,
\end{equation}
%
and subtracting this from \cref{eqn:multiphysicsL10:280} leaves
%
\begin{equation}\label{eqn:multiphysicsL10:340}
0 = \evalbar{\PD{x}{f}}{x^k} \lr{ x^{k+1} - \cancel{x^k} - x^\conj + \cancel{x^k} }
- \inv{2} \evalbar{\PDSq{x}{f}}{\tilde{x}} \lr{ x^\conj - x^k }^2.
\end{equation}
%
Solving for the difference from the solution, the error is
%
\begin{equation}\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,
\end{equation}
%
or in absolute value
%
\begin{equation}\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.
\end{equation}
%
Convergence is quadratic in the error from the previous iteration.  There will be trouble if the derivative goes small at any point in the iteration region. For example in \cref{fig:lecture10:lecture10Fig7}, the iteration could easily end up in the zero derivative region.
%
\imageFigure{../figures/ece1254-multiphysics/lecture10Fig7}{Newton's method with small derivative region.}{fig:lecture10:lecture10Fig7}{0.3}
%
\paragraph{When to stop iteration}
%
One way to check is to look to see if the difference
%
\begin{equation}\label{eqn:multiphysicsL10:420}
\Norm{ x^{k+1} - x^k } < \epsilon_{\Delta x},
\end{equation}
%
however, when the function has a very step slope,
%
%F8
%
this may not be sufficient unless the trial solution is also substituted to see if the desired match has been achieved.

Alternatively, if the slope is shallow as in \cref{fig:lecture10:lecture10Fig7}, then checking for just \( \Abs{ f(x^{k+1} } < \epsilon_f \) may also put the iteration off target.

Finally, a relative error check to avoid false convergence may also be required.  \Cref{fig:lecture10:lecture10Fig10} shows both absolute convergence criteria
%
\imageFigure{../figures/ece1254-multiphysics/lecture10Fig10}{Possible relative error difference required.}{fig:lecture10:lecture10Fig10}{0.3}
%
\begin{subequations}
\begin{equation}\label{eqn:multiphysicsL10:440}
\Abs{x^{k+1} - x^k} < \epsilon_{\Delta x}
\end{equation}
\begin{equation}\label{eqn:multiphysicsL10:460}
\Abs{f(x^{k+1}) } < \epsilon_{f}.
\end{equation}
\end{subequations}
%
A small relative difference may also be required
%
\begin{equation}\label{eqn:multiphysicsL10:480}
\frac{\Abs{x^{k+1} - x^k}}{\Abs{x^k}}
< \epsilon_{f,r}.
\end{equation}

This can become problematic in real world engineering examples such as to diode of \cref{fig:lecture10:lecture10Fig11}, containing both shallow regions and fast growing or dropping regions.
%
\imageFigure{../figures/ece1254-multiphysics/lecture10Fig11}{Diode current curve.}{fig:lecture10:lecture10Fig11}{0.3}
%
%\section{Theorems}
%\EndArticle
%\EndNoBibArticle