multiphysicsL5.tex

%
% Copyright © 2014 Peeter Joot.  All Rights Reserved.
% Licenced as described in the file LICENSE under the root directory of this GIT repository.
%
%\input{../blogpost.tex}
%\renewcommand{\basename}{multiphysicsL5}
%\renewcommand{\dirname}{notes/ece1254/}
%\newcommand{\keywords}{ECE1254H}
%\input{../peeter_prologue_print2.tex}
%
%\usepackage{kbordermatrix}
%
%\beginArtNoToc
%
%\generatetitle{ECE1254H Modeling of Multiphysics Systems.  Lecture 5: Numerical error and conditioning.  Taught by Prof.\ Piero Triverio}
%\chapter{Numerical errors and conditioning}
\label{chap:multiphysicsL5}
%
%\section{Disclaimer}
%
%Peeter's lecture notes from class.  These may be incoherent and rough.
%
\section{Strict diagonal dominance.}
\index{strict diagonal dominance}
%
Related to a theorem on one of the slides:
%
\index{strict diagonal dominance}
\makedefinition{Strictly diagonally dominant.}{dfn:multiphysicsL5:20}{
%
A matrix \( [ M_{ij} ] \) is strictly diagonally dominant if
%
\begin{equation}\label{eqn:multiphysicsL5:40}
\Abs{M_{ii}} > \sum_{j \ne i} \Abs{ M_{ij} } \quad \forall i
\end{equation}
}
%
For example, the stamp matrix
%
\begin{equation}\label{eqn:multiphysicsL5:60}
\kbordermatrix{
    & i      & j \\
i & \inv{R}  & -\inv{R} \\
j & -\inv{R} & \inv{R}
}
\end{equation}
%
is not strictly diagonally dominant.  For row \(i\) this strict dominance can be achieved by adding a reference resistor
%
\begin{equation}\label{eqn:multiphysicsL5:80}
\kbordermatrix{
    & i      & j \\
i & \inv{R_0} + \inv{R}  & -\inv{R} \\
j & -\inv{R} & \inv{R}
}
\end{equation}
%
However, even with strict dominance, there will be trouble with ill posed (perturbative) systems.
%
Round off error examples with \textAndIndex{double precision}
%
\begin{equation}\label{eqn:multiphysicsL5:100}
\lr{ 1 - 1 } + \pi 10^{-17} = \pi 10^{-17},
\end{equation}
%
vs.
%
\begin{equation}\label{eqn:multiphysicsL5:120}
\lr{ 1 + \pi 10^{-17} } -1 = 0.
\end{equation}
%
This is demonstrated by
%
\begin{verbatim}
#include <stdio.h>
#include <math.h>
%
// produces:
// 0 3.14159e-17
int main()
{
   double d1 = (1 + M_PI * 1e-17) - 1 ;
   double d2 = M_PI * 1e-17 ;
%
   printf( "%g %g\n", d1, d2 ) ;
%
   return 0 ;
}
\end{verbatim}
%
Note that a union and bitfield \citep{peeterjRoundingCpp} can be useful for exploring double precision representation.
%
\section{Exploring uniqueness and existence.}
\index{uniqueness}
\index{existence}
%
For a matrix system \( \BM x = \Bb \) in column format, with
%
\begin{equation}\label{eqn:multiphysicsL5:140}
\begin{bmatrix}
\BM_1 & \BM_2 & \cdots & \BM_N
\end{bmatrix}
\begin{bmatrix}
   x_1 \\
   x_2 \\
   \vdots \\
   x_N
\end{bmatrix}
= \Bb.
\end{equation}
%
This can be written as
%
\begin{equation}\label{eqn:multiphysicsL5:160}
\begin{aligned}
\mathLabelBox
{
   x_1
}
{weight}
   \BM_1 +
\mathLabelBox
[
   labelstyle={below of=m\themathLableNode, below of=m\themathLableNode}
]
{
   x_2
}
{weight}
   \BM_2 + \cdots x_N \BM_N \\ &= \Bb.
\end{aligned}
\end{equation}
%
Linear dependence means
%
\begin{equation}\label{eqn:multiphysicsL5:180}
y_1 \BM_1 + y_2 \BM_2 + \cdots y_N \BM_N = 0,
\end{equation}
%
or \( M \By = 0 \).
%
\captionedTable{Solution space.}{tab:lecture5:1}{
\begin{tabular}{|l|l|l|}
\hline
& \( \Bb \in \Span \{M\_i\} \)
& \( \Bb \notin \Span \{M\_i\} \) \\ \hline
col(M) lin. indep. & \(\Bx\) exists and is unique & No solution \\ \hline
col(M) lin. dep. & \(\Bx\) exists.  \( \infty \) solutions & No solution \\ \hline
\end{tabular}
}
%
With a linear dependency an additional solution, given solution \( \Bx \) is \( \Bx^1 = \Bx + \alpha y \).  This becomes relevant for numerical processing since for a system \( M \Bx^1 = \Bb \) a \( \alpha M \By \) can often be found, for which
%
\begin{equation}\label{eqn:multiphysicsL5:200}
M \Bx + \alpha M \By = \Bb,
\end{equation}
%
where \( \alpha M \By \) is of order \( 10^{-20} \).
%
\section{Perturbation and norms.}
\index{perturbation}
\index{vector norm}
%
Consider a perturbation to the system \( M \Bx = \Bb \)
%
\begin{equation}\label{eqn:multiphysicsL5:240}
\lr{ M + \delta M } \lr{ \Bx + \delta \Bx } = \Bb.
\end{equation}
%
Some vector norms
%
\begin{itemize}
\item \(L_1\) norm
%
\begin{equation}\label{eqn:multiphysicsL5:260}
\Norm{ \Bx }_1 = \sum_i \Abs{ x_i }
\end{equation}
%
\item \(L_2\) norm
%
\begin{equation}\label{eqn:multiphysicsL5:280}
\Norm{ \Bx }_2 = \sqrt{ \sum_i \Abs{ x_i }^2 }
\end{equation}
%
\item \(L_\infty\) norm
%
\begin{equation}\label{eqn:multiphysicsL5:300}
\Norm{ \Bx }_\infty = \max_i \Abs{ x_i }.
\end{equation}
\end{itemize}
%
These are illustrated  for \( \Bx = (x_1, x_2) \) in \cref{fig:lecture5:lecture5Fig2}.
%
\imageFigure{../figures/ece1254-multiphysics/lecture5Fig2}{Some vector norms.}{fig:lecture5:lecture5Fig2}{0.15}
%
\section{Matrix norm.}
\index{matrix norm}
%
A matrix operation \( \By = M \Bx \) can be thought of as a transformation as in \cref{fig:lecture5:lecture5Fig3}.
%
\imageFigure{../figures/ece1254-multiphysics/lecture5Fig3}{Matrix as a transformation.}{fig:lecture5:lecture5Fig3}{0.1}
%
The 1-norm for a Matrix is defined as
%
\begin{equation}\label{eqn:multiphysicsL5:340}
\Norm{ M } = \max_{ \Norm{\Bx}_1 = 1 } \Norm{ M \Bx },
\end{equation}
%
and the matrix 2-norm is defined as
%
\begin{equation}\label{eqn:multiphysicsL5:360}
\Norm{ M }_2 = \max_{ \Norm{\Bx}_2 = 1 } \Norm{ M \Bx }_2.
\end{equation}
%
%\EndArticle
%\EndNoBibArticle