From 69896e27f0b410582a18b0b3a6424d2028655c5d Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 31 May 2012 15:17:02 +0000 Subject: [PATCH] Slightly more. git-svn-id: https://svn.dealii.org/trunk@25587 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-15/doc/intro.dox | 104 ++++++++++++++++++++--- deal.II/examples/step-15/doc/results.dox | 9 ++ deal.II/examples/step-15/step-15.cc | 5 +- 3 files changed, 101 insertions(+), 17 deletions(-) diff --git a/deal.II/examples/step-15/doc/intro.dox b/deal.II/examples/step-15/doc/intro.dox index 90103af62d..b9cabbb931 100644 --- a/deal.II/examples/step-15/doc/intro.dox +++ b/deal.II/examples/step-15/doc/intro.dox @@ -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. -

Summary

- -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}$. +

Choice of step length and globalization

+ +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 results section. + + +

Summary of the algorithm and testcase

+ +Overall, the algorithm we use in this program works as follows: +
    +
  1. + 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 + MinimalSurfaceProblem::set_boundary_values). Set + $n=0$. +
    2. + +
  2. + 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$. +
    3. + +
  3. + 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 + MinimalSurfaceProblem::determine_step_length. +
    4. + +
  4. + The new approximation of the solution is given by + $u^{n+1}=u^{n}+\alpha^n \delta u^{n}$. +
    5. + +
  5. + 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 + MinimalSurfaceProblem::set_boundary_values). 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. +
    6. + +
  6. + Set $n\leftarrow n+1$ and go to step 2. +
    7. +
diff --git a/deal.II/examples/step-15/doc/results.dox b/deal.II/examples/step-15/doc/results.dox index 730f329641..b85402b9ae 100644 --- a/deal.II/examples/step-15/doc/results.dox +++ b/deal.II/examples/step-15/doc/results.dox @@ -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 + + + + +

Possibilities for extensions

+ +

Step length control

+ +

Integrating mesh refinement and nonlinear solvers

diff --git a/deal.II/examples/step-15/step-15.cc b/deal.II/examples/step-15/step-15.cc index 6ab31f48b4..41d36fa8b8 100644 --- a/deal.II/examples/step-15/step-15.cc +++ b/deal.II/examples/step-15/step-15.cc @@ -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 triangulation; @@ -150,8 +150,7 @@ namespace Step15 double BoundaryValues::value (const Point &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 MinimalSurfaceProblem class implementation} -- 2.39.5