Damping for Picard fixed point iterations

Part of the !5714 (merged)

Test case

The test case is a heat conduction problem:

Governing PDE (transient heat conduction, 1D, x \in [0,1]) reduces to:

\frac{\partial T}{\partial t} - \frac{\partial^2 T}{\partial x^2} = 0

Initial condition T(x,0) = 0. Boundary conditions:

  • Left, x=0 — Dirichlet: T(0,t) = 0
  • Right, x=1 — nonlinear Robin flux: q(T) = \alpha(T),(T - T_\infty), \qquad \alpha(T) = \frac{h_0}{1 + T^2} with T_\infty = 20, h_0 = 100.

Damping

Picard treats the nonlinear coefficient with a lag: at iterate x_k it freezes \alpha at the old temperature, solves one linear system, gets raw output g(x_k). Fixpoint iteration:

x_{k+1}=g(x_k)

converges only if the map g is a contraction near the solution, i.e. spectral radius \rho(g'(x^*))<1.

Here the coefficient is violently nonlinear:

\alpha(0)=h_0=100,\qquad \alpha(T_\infty)=\frac{100}{1+400}\approx 0.25

Two orders of magnitude swing over the operating range, plus large driving gap (T_\infty=20 vs IC 0). The frozen-coefficient lag over-corrects each step: guess low T gives huge \alpha, huge flux, overshoots high; high T gives tiny \alpha, undershoots back. So g'(x^*) is negative with magnitude >1 - iteration oscillates and diverges. Plain Picard never converges (that is why undamped test caps at 200 iters yet still fails).

Under-relaxation blends old and raw:

x_{k+1}=(1-\beta)\,x_k+\beta\,g(x_k),\qquad \beta\in(0,1]

Effective map derivative becomes \;1-\beta+\beta\,g'\;. Picking \beta small enough pulls the magnitude below 1:

|1-\beta(1-g')|<1

For oscillatory divergence (g'<-1) this needs \beta<\tfrac{2}{1-g'}. \beta=0.39 sits inside that band; \beta=1 does not. Damping trades speed for a guaranteed contraction.


  1. Tests covering your feature were added?
  2. Any new feature or behaviour change was documented?
Edited by Dmitri Naumov

Merge request reports

Loading
Loading