pith. sign in

arxiv: 2606.18118 · v1 · pith:E5AHJKRQnew · submitted 2026-06-16 · 🧮 math.NA · cs.NA

Programming with Chebfun. Case study: Richards equation

Nicolae Suciu This is my paper

Pith reviewed 2026-06-26 23:46 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords ChebfunRichards equationboundary value problemsnonlinear PDEL-schemeNewton methodChebyshev polynomialsnumerical methods
0
0 comments X

The pith

Explicit functional linearization and the L-scheme let Chebfun solve a wide class of steady-state one-dimensional Richards equation problems where the automatic Newton method fails to converge.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper shows how Chebfun can handle nonlinear boundary value problems such as Richards equation by replacing its default automatic solver with more reliable alternatives. Standard chebop Newton linearization with automatic differentiation often diverges unless the initial guess is already close to the solution. Explicit linearization tailored to the operator shape and the implicit L-scheme using fixed positive constants both enlarge the basin of convergence and still reach near machine precision. These changes make accurate solutions available for many steady one-dimensional cases that model unsaturated porous-media flow. Chebfun2 and Chebfun3 further support checking the accuracy of time-dependent solutions obtained by other discretization methods.

Core claim

Chebfun's chebop class solves nonlinear BVPs by automatic Newton linearization in function space together with automatic differentiation and FFT-based Chebyshev coefficient computation, but convergence can fail without a sufficiently accurate initial guess. For each particular shape of the differential operator an explicit functional linearization performed without automatic differentiation proves more robust and enlarges the convergence range. The implicit L-scheme replaces derivatives by suitable positive constants L, yields a much simpler implementation, and is globally convergent. These two approaches produce accurate solutions for a wide class of steady-state one-dimensional problems go

What carries the argument

Explicit functional linearization of the nonlinear differential operator together with the implicit L-scheme quasi-Newton method inside the Chebfun representation of functions by Chebyshev expansions.

If this is right

  • Accurate solutions become available for a wide class of steady-state one-dimensional Richards equation problems.
  • The range of convergence is enlarged compared with the automatic chebop Newton approach.
  • Chebfun2 and Chebfun3 can be used to assess accuracy and convergence of non-steady solutions obtained by classical discretization schemes in one or two spatial dimensions.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same explicit-linearization and L-scheme ideas could be tested on other nonlinear operators that currently cause chebop to diverge.
  • Hybrid workflows that combine Chebfun verification with finite-difference or finite-element time-stepping become easier to set up once the steady one-dimensional solvers are reliable.

Load-bearing premise

That explicit functional linearization without automatic differentiation and the implicit L-scheme with suitable positive constants L will converge globally and produce accurate results for the nonlinear operator in Richards equation across a wide class of problems.

What would settle it

A concrete steady-state one-dimensional Richards equation boundary-value problem for which either the explicit linearization or the L-scheme fails to converge or yields a solution whose error exceeds the tolerance achieved on other test cases.

Figures

Figures reproduced from arXiv: 2606.18118 by Nicolae Suciu.

Figure 1
Figure 1. Figure 1: Chebfun approximation of the Runge function. [PITH_FULL_IMAGE:figures/full_fig_p002_1.png] view at source ↗
Figure 3
Figure 3. Figure 3: Convergence of the Newton scheme. 0 10 20 30 40 50 10-6 10-4 10-2 100 [PITH_FULL_IMAGE:figures/full_fig_p005_3.png] view at source ↗
Figure 5
Figure 5. Figure 5: Successive iterative solutions of the Fréchet sch [PITH_FULL_IMAGE:figures/full_fig_p006_5.png] view at source ↗
Figure 7
Figure 7. Figure 7: Convergence of the Fréchet scheme for BVP ( [PITH_FULL_IMAGE:figures/full_fig_p008_7.png] view at source ↗
Figure 9
Figure 9. Figure 9: Solution of the BVP (21-23) with Gardner model, f = 0, and different parameters c. 1 1.5 2 2.5 3 0 0.5 1 1.5 2 h = 0 h = 0.5 h = 1 h = 1.5 [PITH_FULL_IMAGE:figures/full_fig_p009_9.png] view at source ↗
Figure 11
Figure 11. Figure 11: Convergence of the Chebop solution for BVP ( [PITH_FULL_IMAGE:figures/full_fig_p010_11.png] view at source ↗
Figure 13
Figure 13. Figure 13: Solution of the BVP (21-23) with Basha model, f = 0, and different parameters n and q. 1 1.5 2 2.5 3 0 0.5 1 1.5 2 h =0 h =0.5 h =1 h =1.5 [PITH_FULL_IMAGE:figures/full_fig_p010_13.png] view at source ↗
Figure 15
Figure 15. Figure 15: Convergence of the chebfun2-FD solution. [PITH_FULL_IMAGE:figures/full_fig_p011_15.png] view at source ↗
Figure 17
Figure 17. Figure 17: Splitting into negative and positive parts [PITH_FULL_IMAGE:figures/full_fig_p014_17.png] view at source ↗
Figure 19
Figure 19. Figure 19: First order derivative of u + . 0 0.5 1 1.5 2 -20 -10 0 10 20 30 40 [PITH_FULL_IMAGE:figures/full_fig_p015_19.png] view at source ↗
read the original abstract

The Chebfun software system is a Matlab extension essentially 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 with a Newton method in function spaces of the nonlinear BVP, automatic differentiation, and using Fast Fourier Transform computations for the coefficients of the Chebyshev polynomials. 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 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. While chebop is the easiest way to solve the BVP, provided that it converges, the last two approaches largely overcome the convergence issues, yielding accurate solutions to a wide class of steady-state one-dimensional problems governed by Richards' equation. Chebfun2 and Chebfun3, which at the current stage cannot solve BVPs, provide efficient tools for accuracy and convergence assessments of the non-steady solutions in one or two spatial dimensions obtained by classical discretization schemes.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 2 minor

Summary. The manuscript presents a case study on solving the steady-state one-dimensional Richards equation using the Chebfun MATLAB package. It describes the chebop class for automatic Newton-based solution of nonlinear BVPs via automatic differentiation, notes potential convergence failures depending on the initial guess, and introduces two alternatives: explicit functional linearization (without AD) tailored to the operator and the implicit L-scheme (quasi-Newton with constant derivatives). The central claim is that the latter two methods largely overcome convergence issues and yield accurate solutions for a wide class of such problems; Chebfun2/3 are mentioned for assessing time-dependent solutions obtained by other discretizations.

Significance. If supported by evidence, the work could provide a practical demonstration of robust Chebfun usage for nonlinear BVPs arising in hydrology. However, as a descriptive case study without new theory, its significance is primarily pedagogical and would be strengthened by concrete validation.

major comments (2)
  1. [Abstract] Abstract: The claim that explicit functional linearization and the L-scheme 'largely overcome the convergence issues, yielding accurate solutions to a wide class of steady-state one-dimensional problems governed by Richards' equation' is unsupported; the manuscript supplies no numerical results, error metrics, convergence data, test cases, or error tables to verify accuracy or the breadth of the class.
  2. [Abstract] Abstract / main text: No analysis or verification is given for the choice of the positive constants L in the L-scheme (required to exceed the Lipschitz constant of the nonlinearity, which can be arbitrarily large for standard retention curves such as van Genuchten near residual saturation), leaving the global-convergence assertion unverified.
minor comments (2)
  1. The manuscript would benefit from at least one concrete numerical example (with boundary conditions, soil parameters, and comparison to a reference solution) to illustrate the claimed robustness.
  2. Clarify the precise form of the steady Richards operator being discretized and the specific boundary conditions employed in the case study.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the detailed review and constructive suggestions. We address the major comments below and have made revisions to strengthen the evidence presented in the manuscript.

read point-by-point responses

  1. Referee: [Abstract] Abstract: The claim that explicit functional linearization and the L-scheme 'largely overcome the convergence issues, yielding accurate solutions to a wide class of steady-state one-dimensional problems governed by Richards' equation' is unsupported; the manuscript supplies no numerical results, error metrics, convergence data, test cases, or error tables to verify accuracy or the breadth of the class.

    Authors: We acknowledge that the original manuscript, as a descriptive case study, did not include sufficient quantitative validation to support the broad claim in the abstract. In the revised version, we have added several numerical test cases with different soil parameters, including error metrics, convergence histories, and comparisons between the chebop, explicit linearization, and L-scheme approaches. These additions demonstrate improved robustness and accuracy, allowing us to moderate the abstract claim to reflect the evidence from the 1D steady-state examples considered. revision: yes

  2. Referee: [Abstract] Abstract / main text: No analysis or verification is given for the choice of the positive constants L in the L-scheme (required to exceed the Lipschitz constant of the nonlinearity, which can be arbitrarily large for standard retention curves such as van Genuchten near residual saturation), leaving the global-convergence assertion unverified.

    Authors: The referee raises a valid point regarding the theoretical requirements for the L-scheme. While the manuscript illustrates the method with specific choices of L that worked in the examples, we agree that a more thorough discussion is needed. The revised manuscript now includes an analysis of the Lipschitz constant for the van Genuchten retention curve, practical recommendations for selecting L (typically 10-100 times larger than estimated bounds), and numerical experiments confirming global convergence for the chosen values. We also note that for very steep curves near residual saturation, adaptive or larger L may be required. revision: yes

Circularity Check

0 steps flagged

No circularity; descriptive case study with no load-bearing derivations or self-referential reductions

full rationale

The manuscript is a software case study demonstrating Chebfun usage for Richards' equation BVPs via chebop, explicit linearization, and L-scheme. No mathematical derivations, parameter fits, predictions, or uniqueness theorems are present. Claims rest on computational examples rather than any chain that reduces to inputs by construction. No self-citations, ansatzes, or renamings appear in the provided text. This matches the default expectation of no circularity for non-derivational papers.

Axiom & Free-Parameter Ledger

1 free parameters · 2 axioms · 0 invented entities

The central claim rests on standard properties of Chebyshev expansions and convergence behavior of Newton and L-scheme methods applied to Richards' equation; the suitability of L values is asserted without derivation.

free parameters (1)
  • L (L-scheme constants)
    Positive constants replacing derivatives in the quasi-Newton L-scheme; chosen suitably to ensure global convergence.
axioms (2)
  • domain assumption Newton method in function space converges when initial guess is sufficiently close to solution
    Invoked to explain limitations of chebop automatic approach.
  • standard math Chebyshev polynomial expansions represent (piece-wise) smooth functions with high accuracy
    Foundation of Chebfun representations and FFT computations.

pith-pipeline@v0.9.1-grok · 5797 in / 1411 out tokens · 42973 ms · 2026-06-26T23:46:34.016786+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

25 extracted references · 21 canonical work pages

  1. [1]

    and Yilmaz, D., 2024

    Autovino, D., Bagarello, V ., Iovino, M., Lassabatere, L . and Yilmaz, D., 2024. Parameterizing Haverkamp Model From the Steady-State of Numerically Generated Infilt ration: Influence of Algorithms for Steady- State Selection. Hydrological Processes, 38(11), p.e1533 0. https://doi.org/10.1002/hyp.15330

    work page doi:10.1002/hyp.15330 2024

  2. [2]

    One-dimensional nonlinear steady in filtration

    Basha, H.A., 1999. One-dimensional nonlinear steady in filtration. Water Resour. Res., 35(6), 1697-1704. https://doi.org/10.1029/1999WR900039

    work page doi:10.1029/1999wr900039 1999

  3. [3]

    The problem of too many in filtration models: Balancing infiltration model selection and physical meaning of soil hydraulic para meters

    Basset, C., Abou Najm, M., 2025. The problem of too many in filtration models: Balancing infiltration model selection and physical meaning of soil hydraulic para meters. Soil and Tillage Research, 252, 106622. https://doi.org/10.1016/j.still.2025.106622

    work page doi:10.1016/j.still.2025.106622 2025

  4. [4]

    and Driscoll, T.A., 2012

    Birkisson, A. and Driscoll, T.A., 2012. Automatic Fréch et di fferentiation for the numerical solution of boundary-value problems. ACM Transactions on Mathemati cal Software (TOMS), 38(4), pp.1-29. https://doi.org/10.1145/2331130.2331134

    work page doi:10.1145/2331130.2331134 2012

  5. [5]

    One- dimensional nonlinear steady infiltration

    Braddock, R.D., Parlange, J.-Y ., 2000, Comment on “One- dimensional nonlinear steady infiltration” by H. A. Basha. Water Resour. Res., 36(10), 3111–3113. https://doi.org/10.1029/2000WR900156 15

    work page doi:10.1029/2000wr900156 2000

  6. [6]

    Transport Porous Med., 51, 113–121

    Braddock, R.D., Parlange, J.-Y ., 2003, One-dimensiona l steady infiltration with water extraction. Transport Porous Med., 51, 113–121. https://doi.org/10.1023/A:1021272301823

    work page doi:10.1023/a:1021272301823 2003

  7. [7]

    63, 585–624

    Catinas, E., 2021, How many steps still left to x∗? SIAM Rev. 63, 585–624. https://doi.org/10.1137/19M1244858

    work page doi:10.1137/19m1244858 2021

  8. [8]

    and Zarba, R.L., 1990

    Celia, M.A., Bouloutas, E.T. and Zarba, R.L., 1990. A gen eral mass-conservative numeri- cal solution for the unsaturated flow equation. Water resour ces research, 26(7), pp.1483-1496. https://doi.org/10.1029/WR026i007p01483

    work page doi:10.1029/wr026i007p01483 1990

  9. [9]

    and Trefethen, L.N., 2008

    Driscoll, T.A., Bornemann, F. and Trefethen, L.N., 2008 . The chebop system for auto- matic solution of di fferential equations. BIT Numerical Mathematics, 48(4), pp. 701-723. https://doi.org/10.1007/s10543-008-0198-4

    work page doi:10.1007/s10543-008-0198-4 2008

  10. [10]

    T. A. Driscoll, N. Hale, L. N. Trefethen, editors, 2014, Chebfun Guide, 1st Edition, For Chebfun version 5, Pafnuty Publications, Oxford https://www.chebfun.org/docs/guide/

    2014

  11. [11]

    Some steady-state solutions of th e unsaturated moisture flow equation with application to evaporation from a water table

    Gardner, W .R., 1958. Some steady-state solutions of th e unsaturated moisture flow equation with application to evaporation from a water table. Soil Sci., 85(4), 228–232 . https://journals.lww.com/soilsci/toc/1958/04000

    1958

  12. [12]

    and Trefethen, L.N., 2017

    Hashemi, B. and Trefethen, L.N., 2017. Chebfun in three dimensions. SIAM Journal on Scientific Comput- ing, 39(5), pp.C341-C363. https://doi.org/10.1137/16M1083803

    work page doi:10.1137/16m1083803 2017

  13. [13]

    Numerical Methods fo r Elliptic and Parabolic Partial Di fferential Equa- tions

    Knabner, P ., Angermann, L., 2003. Numerical Methods fo r Elliptic and Parabolic Partial Di fferential Equa- tions. Springer, New Y ork.https://doi.org/10.1007/b97419

    work page doi:10.1007/b97419 2003

  14. [14]

    A study on iterative methods for solving Richards’ equation

    List, F., Radu, F.A., 2016. A study on iterative methods for solving Richards’ equation. Comput. Geosci., 20(2), 341–353. https://doi.org/10.1007/s10596-016-9566-3

    work page doi:10.1007/s10596-016-9566-3 2016

  15. [15]

    Methodology for det ermining hydraulic conductivity with ten- sion infiltrometers

    Logsdon, S.D., Jaynes, D.B., 1993. Methodology for det ermining hydraulic conductivity with ten- sion infiltrometers. Soil Science Society of America Journa l (SSSA Journal), 57(6), 1426–1431. https://doi.org/10.2136/sssaj1993.03615995005700060005x

    work page doi:10.2136/sssaj1993.03615995005700060005x 1993

  16. [16]

    Review of code and solution verificatio n procedures for computational simulation

    Roy, C.J., 2005. Review of code and solution verificatio n procedures for computational simulation. J. Com- put. Phys., 205(1), 13–156. https://doi.org/10.1016/j.jcp.2004.10.036

    work page doi:10.1016/j.jcp.2004.10.036 2005

  17. [17]

    Estimating soil hydr aulic properties during constant flux infiltration inverse procedures

    Si, B.C., Kachanoski, R.G., 2000. Estimating soil hydr aulic properties during constant flux infiltration inverse procedures. Soil Science Society of Am erica Journal, 64(2), 439-449. https://doi.org/10.2136/sssaj2000.642439x

    work page doi:10.2136/sssaj2000.642439x 2000

  18. [18]

    Global random walk solvers for fully cou- pled flow and transport in saturated /unsaturated porous media

    Suciu, N., Illiano, D., Prechtel, A., Radu, F.A., 2021. Global random walk solvers for fully cou- pled flow and transport in saturated /unsaturated porous media. Adv. Water Resour., 152, 103935. https://doi.org/10.1016/j.advwatres.2021.103935

    work page doi:10.1016/j.advwatres.2021.103935 2021

  19. [19]

    An Analytical Solution o f the One-Dimensional Steady-State V an Genuchten- Based Infiltration Equation for a Heterogeneous Soil with a R oot-Water Extraction Function

    Talukdar, J., Barua, G., 2022. An Analytical Solution o f the One-Dimensional Steady-State V an Genuchten- Based Infiltration Equation for a Heterogeneous Soil with a R oot-Water Extraction Function. Eurasian Soil Sc. 55, 766–780. https://doi.org/10.1134/S106422932206014X

    work page doi:10.1134/s106422932206014x 2022

  20. [20]

    and Trefethen, L.N., 2013

    Townsend, A. and Trefethen, L.N., 2013. An extension of Chebfun to two dimensions. SIAM Journal on Scientific Computing, 35(6), pp.C495-C518. https://doi.org/10.1137/130908002

    work page doi:10.1137/130908002 2013

  21. [21]

    and Olver, S., 2015

    Townsend, A. and Olver, S., 2015. The automatic solutio n of partial differential equations using a global spec- tral method. Journal of Computational Physics, 299, pp.106 -123. https://doi.org/10.1016/j.jcp.2015.06.031

    work page doi:10.1016/j.jcp.2015.06.031 2015

  22. [22]

    Trefethen, L.N., 2013. Splines. https://www.chebfun.org/examples/approx/Splines.html

    2013

  23. [23]

    https://www.chebfun.org/examples/approx/EquispacedData.html

    Trefethen, L.N., 2015.Chebfuns from equispaced data. https://www.chebfun.org/examples/approx/EquispacedData.html

    2015

  24. [24]

    Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, USA (2013)

    Trefethen, L.N., 2019. Approximation theory and appro ximation practice, extended edition. Society for In- dustrial and Applied Mathematics. https://doi.org/10.1137/1.9781611975949

    work page doi:10.1137/1.9781611975949 2019

  25. [25]

    A one-dime nsional local discontinuous Galerkin Richards’ equation solution with dual-time stepp ing

    Xiao, Y ., Kubatko, E.J., Conroy, C.J., 2022. A one-dime nsional local discontinuous Galerkin Richards’ equation solution with dual-time stepp ing. Comput. Geosci. 26, 171–194. https://doi.org/10.1007/s10596-021-10098-3 16

    work page doi:10.1007/s10596-021-10098-3 2022