As described above, we solve this equation using Newton's method in which we
compute the $n$th approximate solution from the $n-1$st one, and use
-a damping parameter $\lambda^n$ to get better global convergence behavior:
+a damping parameter $\alpha^n$ to get better global convergence behavior:
@f{align*}
F'(u^{n},\delta u^{n})&=- F(u^{n})
\\
- u^{n+1}&=u^{n}+\lambda^n \delta u^{n}
+ u^{n+1}&=u^{n}+\alpha^n \delta u^{n}
@f}
with
@f[
On the other hand, $B$ is also positive definite, which confers the same
property onto $A$. This can be seen by noting that the vector $v_1 =
\frac{\nabla u^n}{|\nabla u^n|}$ is an eigenvector of $B$ with eigenvalue
-$\lambda_1=1-\frac{1}{1+|\nabla u^n|^2} > 0$ while all vectors $v_2\ldots v_d$
+$\lambda_1=a_n \left(1-\frac{|\nabla u^n|^2}{1+|\nabla u^n|^2}\right) > 0$ while all vectors $v_2\ldots v_d$
that are perpendicular to $v_1$ and each other are eigenvectors with
-eigenvalue $1$. Since all eigenvalues are positive, $B$ is positive definite
+eigenvalue $a_n$. Since all eigenvalues are positive, $B$ is positive definite
and so is $A$. We can thus use the CG method for solving the Newton steps.
It is worth noting, however, that the positive definiteness degenerates for
problems where $\nabla u$ becomes large. In other words, if we simply multiply
all boundary values by 2, then to first order $u$ and $\nabla u$ will also be
multiplied by two, but as a consequence the smallest eigenvalue of $B$ will
-become smaller and the matrix will become more ill-conditioned. It is simple
+become smaller and the matrix will become more ill-conditioned. (More
+specifically, for $|\nabla u^n|\rightarrow\infty$ we have that
+$\lambda_1 \propto a_n \frac{1}{|\nabla u^n|^2}$ whereas
+$\lambda_2\ldots \lambda_d=a_n$; thus, the condition number of $B$,
+which is a multiplicative factor in the condition number of $A$ grows
+like ${\cal O}(|\nabla u^n|^2)$.) It is simple
to verify with the current program that indeed multiplying the boundary values
used in the current program by larger and larger values results in a problem
that will ultimately no longer be solvable using the simple preconditioned CG
method we use here.
-<h3>Summary</h3>
-
-Starting with the function $u^{0}\equiv 0$, the first Newton update is computed by
-solving the system $A^{0}\;\delta U^{0}=b^{0}$ with boundary condition $\delta u^{0}=g$ on
- $\partial \Omega$. The new approximation of the solution is given by
- $u^{1}=u^{0}+\lambda^0 \delta u^{0}$. The next updates are given as solution of
- the linear system $A^{n}\;\delta U^{n}=b^{n}$ with boundary condition $\delta u^{n}=0$ on
- $\partial \Omega$ and the new approximation given by $u^{n+1}=u^{n}+\lambda^n
- \delta u^{n}$.
+<h3> Choice of step length and globalization </h3>
+
+As stated above, Newton's method works by computing a direction
+$\delta u^n$ and then performing the update $u^{n+1} = u^{n}+\alpha^n
+\delta u^{n}$ with a step length $0 < \alpha^n \le 1$. It is a common
+observation that for strongly nonlinear models, Newton's method does
+not converge if we always choose $\alpha^n=1$ unless one starts with
+an initial guess $u^0$ that is sufficiently close to the solution $u$
+of the nonlinear problem. In practice, we don't always have such an
+initial guess, and consequently taking full Newton steps (i.e., using
+$\alpha=1$) does frequently not work.
+
+A common strategy therefore is to use a smaller step length for the
+first few steps while the iterate $u^n$ is still far away from the
+solution $u$ and as we get closer use larger values for $\alpha^n$
+until we can finally start to use full steps $\alpha^n=1$ as we are
+close enough to the solution. The question is of course how to choose
+$\alpha^n$. There are basically two widely used approaches: line
+search and trust region methods.
+
+In this program, we simply always choose the step length equal to
+0.1. This makes sure that for the testcase at hand we do get
+convergence although it is clear that by not eventually reverting to
+full step lengths we forego the rapid, quadratic convergence that
+makes Newton's method so appealing. Obviously, this is a point one
+eventually has to address if the program was made into one that is
+meant to solve more realistic problems. We will comment on this issue
+some more in the <a href="#Results">results section</a>.
+
+
+<h3> Summary of the algorithm and testcase </h3>
+
+Overall, the algorithm we use in this program works as follows:
+<ol>
+<li>
+ Start with the function $u^{0}\equiv 0$ and modify it in such a way
+ that the values of $u^0$ along the boundary equal the correct
+ boundary values $g$ (this happens in
+ <code>MinimalSurfaceProblem::set_boundary_values</code>). Set
+ $n=0$.
+</li>
+
+<li>
+ Compute the Newton update by solving the system $A^{n}\;\delta
+ U^{n}=b^{n}$
+ with boundary condition $\delta u^{n}=0$ on $\partial \Omega$.
+</li>
+
+<li>
+ Compute a step length $\alpha^n$. In this program, we always set
+ $\alpha^n=0.1$. To make things easier to extend later on, this
+ happens in a function of its own, namely in
+ <code>MinimalSurfaceProblem::determine_step_length</code>.
+</li>
+
+<li>
+ The new approximation of the solution is given by
+ $u^{n+1}=u^{n}+\alpha^n \delta u^{n}$.
+</li>
+
+<li>
+ If $n$ is a multiple of 5 then refine the mesh, transfer the
+ solution $u^{n+1}$ to the new mesh and set the values of $u^{n+1}$
+ in such a way that along the boundary we have
+ $u^{n+1}|_{\partial\Gamma}=g$ (again in
+ <code>MinimalSurfaceProblem::set_boundary_values</code>). Note that
+ this isn't automatically
+ guaranteed even though by construction we had that before mesh
+ refinement $u^{n+1}|_{\partial\Gamma}=g$ because mesh refinement
+ adds new nodes to the mesh where we have to interpolate the old
+ solution to the new nodes upon bringing the solution from the old to
+ the new mesh. The values we choose by interpolation may be close to
+ the exact boundary conditions but are, in general, nonetheless not
+ the correct values.
+</li>
+
+<li>
+ Set $n\leftarrow n+1$ and go to step 2.
+</li>
+</ol>
projects / dealii-svn.git / commitdiff