]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Slightly more.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 31 May 2012 15:17:02 +0000 (15:17 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 31 May 2012 15:17:02 +0000 (15:17 +0000)
git-svn-id: https://svn.dealii.org/trunk@25587 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-15/doc/intro.dox
deal.II/examples/step-15/doc/results.dox
deal.II/examples/step-15/step-15.cc

index 90103af62db27e7ee331fa440bf86881aa43c473..b9cabbb93115f2f7b6206a88738b39e308622d09 100644 (file)
@@ -48,11 +48,11 @@ space. In this example, we choose $\Omega$ as the unit disk.
 
 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[
@@ -181,28 +181,104 @@ $B$ is symmetric, and so $A$ is symmetric as well.
 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>
index 730f329641caad04b39d498f0c27b7402e035b29..b85402b9aec1f30aa97e8f926a4d76104214eefc 100644 (file)
@@ -97,3 +97,12 @@ because there isn't much change in the solution. The final solution
 and mesh are shown here:
 
 @image html step-15.grid.png
+
+
+
+<a name="extensions"></a>
+<h3>Possibilities for extensions</h3>
+
+<h4> Step length control </h4>
+
+<h4> Integrating mesh refinement and nonlinear solvers </h4>
index 6ab31f48b49844a93fa0086881cdb5d64bda7fc2..41d36fa8b86667ed177a011a41d8e1b3b31c8799 100644 (file)
@@ -97,7 +97,7 @@ namespace Step15
       void refine_grid ();
       void set_boundary_values ();
       double compute_residual (const double alpha) const;
-      double determine_step_length() const;
+      double determine_step_length () const;
 
       Triangulation<dim>   triangulation;
 
@@ -150,8 +150,7 @@ namespace Step15
   double BoundaryValues<dim>::value (const Point<dim> &p,
                                     const unsigned int /*component*/) const
   {
-    double return_value=sin(2*M_PI*(p[0]+p[1]));
-    return return_value;
+    return std::sin(2 * numbers::PI * (p[0]+p[1]));
   }
 
                                   // @sect3{The <code>MinimalSurfaceProblem</code> class implementation}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.