multiphysicsProblemSet2bProblem2.tex

%
% Copyright © 2014 Peeter Joot.  All Rights Reserved.
% Licenced as described in the file LICENSE under the root directory of this GIT repository.
%
\index{Finite element method}
\makeproblem{Finite element methods.}{multiphysics:problemSet2b:2}{
When modeling distributions of charged particles governed by drift-diffusion equations, equilibrium analysis typically leads to the following nonlinear Poisson equation
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:20}
- \PDSq{x}{\psi} = - \lr{ e^{\psi(x)} - e^{-\psi(x)} },
\end{equation}
where we will consider the interval \( x \in [0,1] \) with the boundary conditions \( \psi(0) = -V \) and \( \psi(1) = V\) .
If a simple finite-difference scheme is used to solve the nonlinear Poisson equation on an N-node grid, the discrete equations are
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:40}
2 \psi_i - \psi_{i+1} - \psi_{i-1} + \Delta x^2 \lr{ e^{\psi_i} - e^{-\psi_i} } = 0,
\end{equation}
for \( i \in [2, \cdots, N - 1] \),
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:60}
2 \psi_1 - \psi_{2} - \lr{ -V } + \Delta x^2 \lr{ e^{\psi_1} - e^{-\psi_1} } = 0,
\end{equation}
and
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:80}
2 \psi_N - \psi_{N-1} - V + \Delta x^2 \lr{ e^{\psi_N} - e^{-\psi_N} } = 0.
\end{equation}

Note, \( \Delta x = 1/(N +1) \) and not \( 1/(N -1) \). The nodes at \( x = 0 \) and \( x = 1 \) are not included in the discretization, but rather enter through the boundary conditions.
\makesubproblem{}{multiphysics:problemSet2b:2a}
Prove that the Jacobian associated with the above discretized equations is nonsingular regardless of the values for the \( \psi_i \)'s.
\makesubproblem{}{multiphysics:problemSet2b:2b}
Based on this observation, how do you expect damped Newton methods will perform, when applied to solving this problem?
\makesubproblem{}{multiphysics:problemSet2b:2c}
Solve the above equations with a multidimensional Newton's method using a zero initial guess for the \( \psi_i\)'s, and \( N = 100 \).  Demonstrate that your program achieves quadratic convergence.
\makesubproblem{}{multiphysics:problemSet2b:2d}
and determine the number of Newton iterations required to insure one part in \( 10^6 \) accuracy in the solution for two cases: when \( V = 1 \) and when \( V = 20 \).
\makesubproblem{}{multiphysics:problemSet2b:2e}
What happens when \( V = 100 \)?
} % makeproblem
%
\makeanswer{multiphysics:problemSet2b:2}{
\withproblemsetsParagraph{
\makeSubAnswer{}{multiphysics:problemSet2b:2a}
%
Before starting consider \cref{fig:ps2bFiniteElement:ps2bFiniteElementFig3}.  that illustrates the partitioning of a possible solution.  This shows why \( \Delta x = 1/(N+1) \).
%
\imageFigure{../figures/ece1254-multiphysics/ps2bFiniteElementFig3}{Partitioning intervals.}{fig:ps2bFiniteElement:ps2bFiniteElementFig3}{0.3}
%
In matrix form the nodal equations are
%
%2 \psi_i - \psi_{i+1} - \psi_{i-1} + \Delta x^2 \lr{ e^{\psi_i} - e^{-\psi_i} } = 0,
%2 \psi_2 - \psi_{3} - \psi_{1} + \Delta x^2 \lr{ e^{\psi_2} - e^{-\psi_2} } = 0,
%2 \psi_N-1 - \psi_{N} - \psi_{N-2} + \Delta x^2 \lr{ e^{\psi_N-1} - e^{-\psi_N-1} } = 0,
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:100}
\begin{bmatrix}
 2 & -1 &         &    & \\
-1 &  2 & -1      &    & \\
   &    & \ddots  &    & \\
   &    &    -1   & 2  &  -1  \\
   &    &         & -1 &  2  \\
\end{bmatrix}
\begin{bmatrix}
\psi_1 \\
\psi_2 \\
\vdots \\
\psi_{N-1} \\
\psi_N \\
\end{bmatrix}
+
\lr{\Delta x}^2
\begin{bmatrix}
e^{\psi_1} - e^{-\psi_1} \\
e^{\psi_2} - e^{-\psi_2} \\
\vdots \\
e^{\psi_{N-1}} - e^{-\psi_{N-1}} \\
e^{\psi_N} - e^{-\psi_N} \\
\end{bmatrix}
=
\begin{bmatrix}
-V \\
0 \\
\vdots \\
0 \\
V
\end{bmatrix}.
\end{equation}
%
Let
%
\begin{subequations}
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:160}
G =
\begin{bmatrix}
 2 & -1 &         &    & \\
-1 &  2 & -1      &    & \\
   &    & \ddots  &    & \\
   &    &    -1   & 2  &  -1  \\
   &    &         & -1 &  2  \\
\end{bmatrix}
\end{equation}
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:180}
\Bx =
\begin{bmatrix}
\psi_1 \\
\psi_2 \\
\vdots \\
\psi_{N-1} \\
\psi_N \\
\end{bmatrix}
\end{equation}
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:200}
H = 2 \lr{ \Delta x }^2
{\begin{bmatrix}
\sinh \psi_i
\end{bmatrix}}_i
\end{equation}
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:240}
\Bb =
\begin{bmatrix}
-V \\
0 \\
\vdots \\
0 \\
V
\end{bmatrix},
\end{equation}
\end{subequations}
%
and
%
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:220}
F(\Bx) = G \Bx + H(\Bx) - \Bb.
\end{equation}
%
The zeros of \( F(\Bx) \) are sought.  Expansion of \( F(\Bx) \) in Taylor series produces the Jacobian term that is to be shown non-singular.  That Taylor expansion, around \( \Bx^0 \) is
%
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:260}
\begin{aligned}
F(\Bx^1)
&= G \Bx^0 + H(\Bx^0) - \Bb +
{\begin{bmatrix}
\evalbar{\PD{\psi_j}{} \lr{G_{ik} \psi_k + H_i} }{\psi_j = \psi_j^0} \lr{ \psi_j^1 - \psi_j^0}
\end{bmatrix}}_i
\\ &= G \Bx^0 + H(\Bx^0) - \Bb +
J(\Bx^0) \lr{ \Bx^1 - \Bx^0 },
\end{aligned}
\end{equation}
%
where the Jacobian evaluated at \( \Bx^0 \) is
%
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:280}
J(\Bx^0) =
G +
\evalbar{
{\begin{bmatrix}
\PD{\psi_j}{H_i}
\end{bmatrix}}_{ij}
}{\psi_j = \psi_j^0}.
\end{equation}
%
In this case the Jacobian is the sum of a constant matrix and one with hyperbolic cosines along the diagonal
%
\begin{subequations}
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:520}
J_H(\Bx^0) = 2 \lr{ \Delta x }^2
{\begin{bmatrix}
\cosh\psi_i \delta_{ij}
\end{bmatrix}}_{ij}.
\end{equation}
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:320}
J(\Bx^0) = G + J_H(\Bx^0)
\end{equation}
\end{subequations}
%
This Jacobian is symmetric, so the eigenvalues are real.  Gershgorin's theorem for the interior rows states that
%
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:440}
\Abs{\lambda_i - \lr{ 2 + 2 \lr{\Delta x}^2 \cosh(\psi_i) }} < 2,
\end{equation}
%
or
%
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:460}
2 \lr{\Delta x}^2 \cosh(\psi_i) < \lambda_i
< 4 + 2 \lr{\Delta x}^2 \cosh(\psi_i).
\end{equation}
%
For the first and last rows Gershgorin's theorem requires
%
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:480}
\Abs{\lambda_i - \lr{ 2 + 2 \lr{\Delta x}^2 \cosh(\psi_i) }} < 1,
\end{equation}
%-\lambda_i + \lr{ 1 + 2 \lr{\Delta x}^2 \cosh(\psi_i) } < 0
%
or
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:500}
1 + 2 \lr{\Delta x}^2 \cosh(\psi_i) < \lambda_i
< 3 + 2 \lr{\Delta x}^2 \cosh(\psi_i).
\end{equation}
%
\Cref{eqn:multiphysicsProblemSet2bProblem2:460} shows that the lower bound on any of the eigenvalues is \( 2 \lr{ \Delta x}^2 \), so while not strictly singular anywhere, the Jacobian can potentially be nearly singular.
%
\makeSubAnswer{}{multiphysics:problemSet2b:2b}
%
Despite the Jacobian being non-singular, there is the possibility that it can have near zero eigenvalues.  This could trigger oscillatory behavior if the iteration changes too much at once (motivating the use of damping).  To apply the damping method outlined in the class slides, the bounds \( beta \) and \( \rho \) of
%
\begin{subequations}
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:540}
\Norm{ J^{-1} } \le \beta
\end{equation}
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:560}
\Norm{ J(\Bx') - J(\Bx) } \le \rho \Norm{ \Bx' - \Bx},
\end{equation}
\end{subequations}
%
are required.  The Jacobian \( J = G + J_H \) is the sum of a tridiagonal matrix \( G \) and a positive definite diagonal matrix \( J_H \).  Since \( J_H \) only increases the magnitude of all the diagonals, it seems reasonable to make the approximation
%
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:580}
\Norm{ J^{-1} } \le \Norm{ G },
\end{equation}
%
a value that is constant and can be precomputed before any iteration.  For the Lipshitz bound
%
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:600}
\begin{aligned}
\Norm{ J(\Bx') - J(\Bx) }
&=
\Norm{ G + J_H(\Bx') - \lr{ G + J_H(\Bx)} } \\
&=
2 \lr{ \Delta x }^2
\Norm{
{\begin{bmatrix}
\lr{\cosh \psi_i' - \cosh \psi_i} \delta_{ij}
\end{bmatrix}}_{ij}
} \\
&=
2 \lr{ \Delta x }^2
\Norm{
{\begin{bmatrix}
\lr{ \inv{2} \lr{\psi_i'}^2 - \inv{2} \psi_i^2 }\delta_{ij}
\end{bmatrix}}_{ij}
} \\
&\lesssim
\lr{ \Delta x }^2 \lr{ 2 V }
\Norm{
{\begin{bmatrix}
\lr{ \psi_i' - \psi_i }\delta_{ij}
\end{bmatrix}}_{ij}
} \\
&\approx
\lr{ \Delta x }^2 V \Norm{ \Bx' - \Bx }.
\end{aligned}
\end{equation}
%
This yields a pair of crude approximations that can be used to determine the damping constant
%
\begin{subequations}
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:620}
\beta = \Norm{ G^{-1} }
\end{equation}
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:640}
\rho = 2 \lr{ \Delta x }^2 V.
\end{equation}
\end{subequations}
%
The damping algorithm of the class slides was
%
\index{damping}
\makealgorithm{Damping algorithm.}{alg:multiphysicsProblemSet2bProblem2:1}{
\IF{\( \Norm{F} \le \inv{\beta^2 \rho} \)}
\STATE \( \alpha = 1 \)
\ELSE
\STATE \( \alpha = \inv{\beta^2 \rho \Norm{F}} \)
\ENDIF
}
%
%%% Newton's method is make the approximation \( F(\Bx^1) \approx 0 \), and to solve for \( \Bx^1 \) in
%%%
%%% \begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:300}
%%% \Bx^1 = \Bx^0 - \lr{ J(\Bx^0) }^{-1} F(\Bx^0),
%%% \end{equation}
%%%
%%% however, rather than calculate the inverse of the Jacobian, LU factorization can be used to solve for \( \Bx^1 - \Bx^0 \).  With a constant Jacobian such a pre-factorization should be very effective in this case.
%%%
%%% %%%The difference of exponentials can be expanded in Taylor series
%%% %%%
%%% %%%\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:120}
%%%\begin{aligned}
%%% %%%e^{\psi_i} - e^{-\psi_i}
%%% %%%&=
%%% %%%\lr{ 1 + \psi_i + \psi_i^3/2 + \psi_i^3/6 + \cdots }
%%% %%%-
%%% %%%\lr{ 1 - \psi_i + \psi_i^3/2 - \psi_i^3/6 + \cdots }
%%% %%%\\ &=
%%% %%%2 \psi_i + \inv{3} \psi_i^3 + \cdots,
%%% %%%
%%%\end{aligned}
%%%\end{equation}
%%% %%%
%%% %%%so the Jacobian is just \( 2 I \), which is never singular.  The linear approximation of the nodal equations is
%%% %%%
%%% %%%\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:140}
%%% %%%\begin{bmatrix}
%%% %%% 2 \lr{1 + \lr{\Delta x}^2} & -1 &         &    & \\
%%% %%%-1 &  2 \lr{1 + \lr{\Delta x}^2} & -1      &    & \\
%%% %%%   &    & \ddots  &    & \\
%%% %%%   &    &    -1   & 2 \lr{ 1 + \lr{\Delta x}^2 } &  -1  \\
%%% %%%   &    &         & -1 &  2 \lr{ 1 + \lr{\Delta x}^2 } \\
%%% %%%\end{bmatrix}
%%% %%%\begin{bmatrix}
%%% %%%\psi_1 \\
%%% %%%\psi_2 \\
%%% %%%\vdots \\
%%% %%%\psi_{N-1} \\
%%% %%%\psi_N \\
%%% %%%\end{bmatrix}
%%% %%%=
%%% %%%\begin{bmatrix}
%%% %%%-V \\
%%% %%%0 \\
%%% %%%\vdots \\
%%% %%%0 \\
%%% %%%V
%%% %%%\end{bmatrix}.
%%% %%%\end{equation}
%%% %%%
%%% %%%FIXME: does he want us to show that this complete Jacobian is non-singular (which is what matters since that's what has to be inverted to apply newton's method).
%%%
%%% In order to apply damped Newton's method, we are looking to minimize \( \Norm{F \lr{ \Bx^0 + \alpha \Delta} }^2 \) with respect to \( \alpha \).  Here \( \Delta = -J^{-1} F \), the direction of the iteration adjustment.  That comes from solving
%%%
%%% \begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:380}
%%%\begin{aligned}
%%% 0
%%% &= \PD{\alpha}{} \Norm{ F( \Bx^0 + \alpha \Delta ) }^2
%%% \\ &= 2 F \cdot \PD{ \alpha }{F}
%%% \\ &= 2 \lr{ G \lr{ \Bx^0 + \alpha \Delta } - \Bb + H } \cdot \lr{
%%% G \Delta +
%%% 2 \lr{ \Delta x }^2
%%% {\begin{bmatrix}
%%% \Delta_i \cosh\lr{ \psi_i + \alpha \Delta_i }
%%% \end{bmatrix}}_i
%%% }
%%%
%%%\end{aligned}
%%%\end{equation}
%%%
%%% Writing
%%%
%%% \begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:400}
%%%\begin{aligned}
%%% H' &=
%%% 2 \lr{ \Delta x }^2
%%% {\begin{bmatrix}
%%% \Delta_i \cosh\lr{ \psi_i + \alpha \Delta_i }
%%% \end{bmatrix}}_i
%%% \approx
%%% 2 \lr{ \Delta x }^2
%%% {\begin{bmatrix}
%%% \Delta_i \cosh \psi_i
%%% \end{bmatrix}}_i
%%% \\ &= J_H(\Bx^0) \Delta
%%% ,
%%%
%%%\end{aligned}
%%%\end{equation}
%%%
%%% then to zeroth order in \( \alpha \) this has solution
%%%
%%% \begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:420}
%%%\begin{aligned}
%%% \alpha
%%% &= -\frac{F(\Bx^0) \cdot \lr{ G \Delta + H' }}{
%%% (G \Delta) \cdot \lr{ G \Delta + H' }
%%% }
%%% \\ &=
%%% - \frac{F^\T J \Delta}{\Delta^\T G J \Delta}
%%% .
%%%
%%%\end{aligned}
%%%\end{equation}
%%%
%%% % BOGUS thought, the 2's along the diagonal will dominate.
%%% %With \( 2 ( 1 + \lr{ \Delta x}^2 ) \approx \) along the diagonal of the The inverse Jacobian being \( O(N^2) \)
%%% %
%%% %\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:360}
%%% %\lr{ J(\Bx^0) }^{-1} = \frac{(N+1)^2}{2} I,
%%% %\end{equation}
%%% %
%%% %the step sizes could be large unless damping is used.
%
\makeSubAnswer{}{multiphysics:problemSet2b:2c}
%
To demonstrate quadratic convergence, let
%
\begin{subequations}
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:660}
u(k) = \Norm{ \Bx^k - \Bx^\conj }
\end{equation}
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:680}
y(k) = u^{k+1},
\end{equation}
\end{subequations}
%
A log-log plot can be used to test for a monomial fit of the form
%
\begin{equation}\label{eqn:multiphysicsProblemSet2bProblem2:700}
y(k) = L u^m(k).
\end{equation}
%
For \( V = 10 \) the plot of \cref{fig:demonstrateQuadraticV10:demonstrateQuadraticV10Fig1}
generated with \matlabFuncPath{plotConvergence()}{ps2b:plotConvergence.m}
%
%V=10 ;
%ee=1e-6 ;
%\matlabFuncPath{[x, r] = newtonsDiffusion( 100, V, ee, ee, ee, 50, 1, 1 )}{ps2b:newtonsDiffusion.m}
%\matlabFuncPath{[L, m] = quadraticDifferences( x, r.allX(:,2:end-1) )}{ps2b:quadraticDifferences.m}
%
shows an almost linear relation, with a monomial least squares fitting of \( L = 0.394, m = 2.004 \).
%
\mathImageFigure{../figures/ece1254-multiphysics/demonstrateQuadraticV10Fig1.pdf}{Quadratic convergence \(V = 10\).}{fig:demonstrateQuadraticV10:demonstrateQuadraticV10Fig1}{0.3}{ps2b:quadraticDifferences.m}
%
For \( V = 20 \), the polynomial fitting divergences slightly \( L = 0.386, m = 1.889 \) from quadratic, with most of that divergence occurring at the end of the iteration as seen in
the plot of \cref{fig:demonstrateQuadraticV20:demonstrateQuadraticV20Fig2}.
%
These iterations used no damping (standard Newton's method), nor used any stepping.
%
\mathImageFigure{../figures/ece1254-multiphysics/demonstrateQuadraticV20Fig2.pdf}{Quadratic convergence \(V = 20\).}{fig:demonstrateQuadraticV20:demonstrateQuadraticV20Fig2}{0.3}{ps2b:quadraticDifferences.m}
%
\makeSubAnswer{}{multiphysics:problemSet2b:2d}
%
Plots with \( V = 1, 10, 20, 100 \) are given in \cref{fig:ps2bP2PlotOfPsiV1StepSize1NoDamping:ps2bP2PlotOfPsiV1StepSize1NoDampingFig1}, \cref{fig:ps2bP2PlotOfPsiV10StepSize1NoDamping:ps2bP2PlotOfPsiV10StepSize1NoDampingFig5}, \cref{fig:ps2bP2PlotOfPsiV20StepSize1NoDamping:ps2bP2PlotOfPsiV20StepSize1NoDampingFig9}, \cref{fig:ps2bP2PlotOfPsiV100StepSize1WithDamping:ps2bP2PlotOfPsiV100StepSize1WithDampingFig15} respectively.
%
%\cref{fig:ps2bP2PlotOfPsiV100StepSize1NoDamping:ps2bP2PlotOfPsiV100StepSize1NoDampingFig13}.
%\cref{fig:ps2bP2PlotOfPsiV100StepSize5NoDamping:ps2bP2PlotOfPsiV100StepSize5NoDampingFig14}.
%\cref{fig:ps2bP2PlotOfPsiV100StepSize5WithDamping:ps2bP2PlotOfPsiV100StepSize5WithDampingFig16}.
%\cref{fig:ps2bP2PlotOfPsiV10StepSize1WithDamping:ps2bP2PlotOfPsiV10StepSize1WithDampingFig7}.
%\cref{fig:ps2bP2PlotOfPsiV10StepSize5NoDamping:ps2bP2PlotOfPsiV10StepSize5NoDampingFig6}.
%\cref{fig:ps2bP2PlotOfPsiV10StepSize5WithDamping:ps2bP2PlotOfPsiV10StepSize5WithDampingFig8}.
%\cref{fig:ps2bP2PlotOfPsiV1StepSize1WithDamping:ps2bP2PlotOfPsiV1StepSize1WithDampingFig3}.
%\cref{fig:ps2bP2PlotOfPsiV1StepSize5NoDamping:ps2bP2PlotOfPsiV1StepSize5NoDampingFig2}.
%\cref{fig:ps2bP2PlotOfPsiV1StepSize5WithDamping:ps2bP2PlotOfPsiV1StepSize5WithDampingFig4}.
%\cref{fig:ps2bP2PlotOfPsiV20StepSize1WithDamping:ps2bP2PlotOfPsiV20StepSize1WithDampingFig11}.
%\cref{fig:ps2bP2PlotOfPsiV20StepSize5NoDamping:ps2bP2PlotOfPsiV20StepSize5NoDampingFig10}.
%\cref{fig:ps2bP2PlotOfPsiV20StepSize5WithDamping:ps2bP2PlotOfPsiV20StepSize5WithDampingFig12}.
\mathImageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV1StepSize1NoDampingFig1.pdf}{\( \psi(x), V = 1 \).}{fig:ps2bP2PlotOfPsiV1StepSize1NoDamping:ps2bP2PlotOfPsiV1StepSize1NoDampingFig1}{0.3}{ps2b:countItersAndPlot.m}
%\imageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV1StepSize1WithDampingFig3}{CAPTION}{fig:ps2bP2PlotOfPsiV1StepSize1WithDamping:ps2bP2PlotOfPsiV1StepSize1WithDampingFig3}{0.3}
%\imageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV1StepSize5NoDampingFig2}{CAPTION}{fig:ps2bP2PlotOfPsiV1StepSize5NoDamping:ps2bP2PlotOfPsiV1StepSize5NoDampingFig2}{0.3}
%\imageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV1StepSize5WithDampingFig4}{CAPTION}{fig:ps2bP2PlotOfPsiV1StepSize5WithDamping:ps2bP2PlotOfPsiV1StepSize5WithDampingFig4}{0.3}
%
\mathImageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV10StepSize1NoDampingFig5.pdf}{\( \psi(x), V = 10 \).}{fig:ps2bP2PlotOfPsiV10StepSize1NoDamping:ps2bP2PlotOfPsiV10StepSize1NoDampingFig5}{0.3}{ps2b:countItersAndPlot.m}
%\imageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV10StepSize1WithDampingFig7}{CAPTION}{fig:ps2bP2PlotOfPsiV10StepSize1WithDamping:ps2bP2PlotOfPsiV10StepSize1WithDampingFig7}{0.3}
%\imageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV10StepSize5NoDampingFig6}{CAPTION}{fig:ps2bP2PlotOfPsiV10StepSize5NoDamping:ps2bP2PlotOfPsiV10StepSize5NoDampingFig6}{0.3}
%\imageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV10StepSize5WithDampingFig8}{CAPTION}{fig:ps2bP2PlotOfPsiV10StepSize5WithDamping:ps2bP2PlotOfPsiV10StepSize5WithDampingFig8}{0.3}
%
\mathImageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV20StepSize1NoDampingFig9.pdf}{\( \psi(x), V = 20 \).}{fig:ps2bP2PlotOfPsiV20StepSize1NoDamping:ps2bP2PlotOfPsiV20StepSize1NoDampingFig9}{0.3}{ps2b:countItersAndPlot.m}
%\imageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV20StepSize1WithDampingFig11}{CAPTION}{fig:ps2bP2PlotOfPsiV20StepSize1WithDamping:ps2bP2PlotOfPsiV20StepSize1WithDampingFig11}{0.3}
%\imageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV20StepSize5NoDampingFig10}{CAPTION}{fig:ps2bP2PlotOfPsiV20StepSize5NoDamping:ps2bP2PlotOfPsiV20StepSize5NoDampingFig10}{0.3}
%\imageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV20StepSize5WithDampingFig12}{CAPTION}{fig:ps2bP2PlotOfPsiV20StepSize5WithDamping:ps2bP2PlotOfPsiV20StepSize5WithDampingFig12}{0.3}
%
%\imageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV100StepSize1NoDampingFig13}{\( \psi(x), V = 100 \) }{fig:ps2bP2PlotOfPsiV100StepSize1NoDamping:ps2bP2PlotOfPsiV100StepSize1NoDampingFig13}{0.3}
\mathImageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV100StepSize1WithDampingFig15.pdf}{\( \psi(x), V = 100 \).}{fig:ps2bP2PlotOfPsiV100StepSize1WithDamping:ps2bP2PlotOfPsiV100StepSize1WithDampingFig15}{0.3}{ps2b:countItersAndPlot.m}
%\imageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV100StepSize5NoDampingFig14}{CAPTION}{fig:ps2bP2PlotOfPsiV100StepSize5NoDamping:ps2bP2PlotOfPsiV100StepSize5NoDampingFig14}{0.3}
%\imageFigure{../figures/ece1254-multiphysics/ps2bP2PlotOfPsiV100StepSize5WithDampingFig16}{CAPTION}{fig:ps2bP2PlotOfPsiV100StepSize5WithDamping:ps2bP2PlotOfPsiV100StepSize5WithDampingFig16}{0.3}
%
The number of iterations required are tabulated in \cref{tab:multiphysicsProblemSet2bProblem2:720} (first 3 columns damped).  Introducing stepping on \( V \) slows convergence for small \( V \), but helps markedly for larger \( V \), more so than the damping algorithm.
%
\captionedTable{Number of iterations to converge on solution.}{tab:multiphysicsProblemSet2bProblem2:720}{
\begin{tabular}{|l|l|l|l|l|}
\hline
%V & \multicolumn{2}{ |c| }{number of steps (No damping)} & \multicolumn{2}{ |c| }{number of steps (With damping)} \\
%\hline
V   & 0 steps& 5 steps& 0 steps& 5 steps(damped)   \\ \hline
1   & 3  & 15 & 3  & 15  \\ \hline
10  & 8  & 23 & 8  & 23  \\ \hline
20  & 16 & 28 & 8  & 28  \\ \hline
100 & 91 & 37 & 89 & 100 \\ \hline
\end{tabular}
}
%
\makeSubAnswer{}{multiphysics:problemSet2b:2e}
%
At \( V = 100 \) without damping or stepping to find the solution, Matlab produces a number of ``Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.'', but does eventually converge.  It is likely that there is a great deal of oscillation reaching the solution.
%
Use of either stepping (5 step intervals) or damping was sufficient to avoid the ill conditioning issues with the Jacobian.  Stepping was more effective than damping in this case to converge on the solution with only 37 iterations required instead of 89.  Using both stepping and damping slowed convergence more than using either separately.
%
Note that the damping that used the \( 1/\beta^2 \rho \) factor was too extreme, perhaps due to the simplifications made computing these \( \beta, \rho \) factors.  A heuristic adjustment of the damping factor was used, first using a damping threshold of \( 1/{\beta \rho} \) and then \( 1/{\sqrt{\beta} \rho} \).  That damping factor was enough to handle the ill conditioning problem, but also allow convergence in a not unreasonable number of iterations.
}
}