Programming with Chebfun. Case study: Richards equation

Abstract

The Chebfun software system is a Matlab extension based on representations of (piece-wise) smooth one-variable functions by expansions in Chebyshev polynomials. One of Chebfun's attractive features is the ability to provide solutions to nonlinear boundary value problems (BVP) with accuracy close to the machine precision. This is done by the chebop class which provides automatic solutions by performing linearizations of the nonlinear BVP with a Newton method in function spaces and automatic differentiation. A drawback of chebop automatic approach is the possible lack of convergence of the Newton method if the initial guess is not close enough to the exact solution. An explicit functional Newton-Frechet linearization done for each particular shape of the differential operator (i.e. without automatic differentiation) proves to be more robust than the chebop class and allows an enlargement of the range of convergence. Another alternative is the implicit L-scheme (quasi-Newton approach with derivatives replaced by suitable positive constants L), with a much simpler implementation and globally convergent. Mixed methods L-scheme/chebop and L-scheme/Frechet may be used to speed up the computation while preserving the robustness and global convergence of the L-scheme. While chebop is the easiest way to solve BVPs, provided that it converges, the alternative approaches overcome the convergence issue, yielding accurate solutions to a wide class of steady-state one-dimensional problems governed by Richards equation. Chebfun2 and Chebfun3, which currently cannot solve BVPs, provide efficient tools for accuracy assessments of the unsteady solutions in one or two spatial dimensions obtained by classical discretization schemes.