Analysis of a hybrid numerical method – decomposing leaf hydraulic conductance

Abstract We describe a hybrid numerical method to solve a boundary value problem where an unknown parameter of the model is chosen to satisfy an additional boundary condition. After the solution of the differential equation is approximated using a one-step method, a secant method is used to update the value of the unknown parameter. The model is a generalization of a model first used to describe water flow through roots, which was later used to describe water flow through the tank bromeliad Guzmania lingulata. In both cases, identification of the unknown parameter represents the decomposition of overall plant conductance into components in the radial and axial directions. We describe convergence of the one-step and secant portions of the method in a base case corresponding to previous applications of the model and in an intermediate case corresponding to a first approximation of the geometry of the leaf. We demonstrate that in the more general case, which better represents the geometry of G. lingulata, the one-step method also converges as expected. Finally, we discuss the implications of including a better description of the geometry of the leaf in context of radial conductance and show that our modeling of the leaf geometry increases the component of the overall leaf conductance in the radial direction by as much as 25%.


Introduction
We present a hybrid numerical method to approximate the solution of the second order, linear differential equation, together with boundary data d 2 u dx 2 = a 2 g(x)u, u(0) = 1, on the interval 0 ≤ x ≤ 1, where the non-negative function g satisfies 0 ≤ g(x) ≤ 1, and the parameter a is chosen to satisfy the boundary condition u (1) = 0.
A variant of this model was first considered by Landsberg and Fowkes (1978) in the context of water transport through a root and later used by North, Lynch, Maharaj, Phillips, and Woodside (2013) to describe water transport through leaves of the tank bromeliad Guzmania lingulata. Both of these modeling efforts seek to better understand how plants transport water by decomposing overall hydraulic conductance into components into the axial (parallel to growth of the plant) and radial (perpendicular to the surface of the plant organ, such as a leaf) directions. The ability to compare conductance in the axial and radial directions allows plant physiologists to categorize a plant's proclivity to conserve water (North et al., 2013), to compare profiles of water potential within a leaf (Frensch & Steudle, 1989), to understand rates of gas exchange within individual leaves (Ocheltree, Nippert, & Prasad, 2012), and to examine the impact of radial water movement on photosynthesis and growth (McKown, Cochard, & Sack, 2010).
In the context of G. lingulata, water is absorbed by the tank region of the leaf and is then transported throughout the blade region, where photosynthesis and water loss occur. The dimensionless spatial variable x represents length from the tank-blade boundary (x = 0) toward the tip of the leaf (x = 1), as shown in Figure 1, and the unknown function u represents scaled water potential within the leaf. Therefore, conditions on the derivative u model the flux of water, so that the condition u (0) = −v 0 indicates an inward flow at the tank-blade boundary and the additional condition u (1) = 0 specifies that there is no flux through the end of the leaf.
A novel component of this modeling approach is the addition of the function g, which describes the width of the plant as a function of the distance along the blade of the plant. We seek to build upon the model described by North et al. (2013) which assumes that the width of G. lingulata is constant along the blade of the plant (equivalently, in the system (1), g(x) ≡ 1) through incorporation of the nonconstant function g. Under the constant width assumption (i.e. g(x) ≡ 1), the exact solution of the system (1) can be found using elementary methods, and the parameter a can be shown to satisfy a transcendental equation (Landsberg & Fowkes, 1978). In general, the exact solution may not be available, which motivates development of the hybrid numerical method.
In Section 2, we describe how a system of the form (1) together with the additional condition u (1) = 0 model properties of water flow in G. lingulata. Here, we also give an estimate for v 0 and describe how an increase in the potential difference which drives the flow corresponds to a decrease in v 0 .
In Section 3, we describe the hybrid numerical method used to approximate the solution of the model system. We begin with the observation that the system (1), although dependent on a spatial variable, is mathematically equivalent to an initial value problem. The hybrid method employs a one-step method to approximate the solution of the system (1) with a guess for the unknown parameter a and a root-finding method to determine the value of the unknown parameter a satisfying the additional condition u (1) = 0. Alternating between an improved guess to the unknown parameter and an approximation for the solution to the system (1) using that improved guess, we expect to identify a reference value of the unknown parameter a and an approximate solution to the system (1) simultaneously.
Our method is analogous to a shooting method for a boundary value problem, which instead incorporates a guess for an artificial initial value to approximate the solution of an  (North et al., 2013). initial value problem, then refines that initial value with a root-finding method in order to satisfy the desired boundary value. Such a shooting method may be found in many introductory textbooks on numerical analysis, see for instance Bradie (2006). However, the literature lacks convergence analysis of a two-state method where the root-finding method is used to determine an unknown parameter of the differential equation, which is a second novel component of this approach.
In Section 4, we first make use of the constant width assumption to analyze the hybrid numerical method. In this base case, we compare the approximate solution obtained using the hybrid method to the exact solution determined using elementary methods and show that the one-step method and the root-finding method both converge as expected, demonstrating proof of concept. As an intermediate case, we discuss a specific curvilinear form for g for which an exact solution to the system (1) is available in terms of special functions. We similarly show that the one-step method and the root-finding method both converge as expected. Finally, we apply the hybrid numerical method to the system (1) with a form for g obtained by fitting the shape of the leaf with a cubic spline. We show that the one-step method again converges as expected for a broad range of the flux condition v 0 specified by the boundary condition u (0) = −v 0 which corresponds to a potential difference up to a theoretical value of O 10 4 MPa.
We close by returning to the context of the stated goal of decomposing total leaf conductance of G. lingulata into components in the axial and radial directions, specifically using a cubic spline representation for g. Our results suggest that improved accounting for the geometry of the leaf with a cubic spline representation gives a value for the decomposed radial conductance that is significantly higher in comparison to the base case obtained previously by North et al. (2013), perhaps by as much as 25%.

Modeling equations
Here, we derive the equations that are used to decompose the leaf conductance into components in the axial and radial directions. Following Ogée, Cuntz, Peylin, and Bariac (2007), we model the xylem as parallel veins of length L. Similar to previous modeling efforts by Landsberg and Fowkes (1978) and by North et al. (2013), we divide water flow within the leaf into components in the axial and radial directions and consider flow from the tank-blade boundary (z = 0) to the tip of the leaf (z = L). In the axial direction we assume that the flux J A is driven by the the gradient of the water potential, as where R A is axial resistance. Our assumptions on the xylem allow the axial conductance (the inverse of axial resistance) to be approximated via a measurement of the radius of an individual vein (North et al., 2013). In the radial direction, we assume that the flux J R is driven by a difference between the water potential in the vein and that in the mesophyll, given by mes , as where R R is the unknown radial resistance. In the blade, the potential difference − mes > 0, which gives rise to a positive outward flux, J R > 0. In a cross-section of the leaf on the interval [z, z + z], conservation of flow gives that where wg(z) gives the width of the leaf atz ∈ [z, z + z]. As z → 0,z → z and Therefore, If the water potential at the tank-blade boundary is known and the flux through the tip of the leaf is zero, the boundary conditions are (0) = 0 and (L) = 0. At the base of the leaf, the total flux J leaf is equal to the flux in the axial direction, so that R A J A (0) = R A J leaf , or equivalently (0) = −R A J leaf . The nondimensional scaling where x = z/L gives rise to the differential equation, together with boundary data where the flux condition is described by the dimensionless value the parameter a satisfies and a is chosen to satisfy u (1) = 0. Since the only portion of a that is unknown is the axial resistance, determination of a is equivalent to decomposing leaf conductance into components in the axial and radial directions.
Since u (0) = −v 0 , the physical interpretation of v 0 is scaled inward flux at x = 0. In many numerical experiments that follow, we use the value v 0 = .8. That value is derived from order of magnitude estimates described by North et al. (2013), summarized here: In other numerical experiments, a range of values is used for v 0 . Since v 0 is inversely proportional to the potential difference 0 − mes , the range of values O 10 −4 < v 0 < O 1 corresponds to an increase in the potential difference on the order of 10 4 MPa.
In comparison to similar derivations by Landsberg and Fowkes (1978) and North et al. (2013), a benefit of this nondimensional approach is that it reduces the number of parameters in the model to just two: the parameter a, representing the ratio of radial to axial resistance, and the flux condition u (0) = −v 0 .

Numerical method
We seek to identify the value of the parameter a such that the approximate solution to the system (1) satisfies u (1) = 0. The value of a is found using two iterative schemes: a onestep method to approximate the solution to the corresponding initial value problem and a root-finding method to ensure that the solution satisfies u (1) = 0. Although description 'initial value problem' is typically reserved for a differential equation with time as the independent variable, the system (1) is mathematically equivalent.
Given any guess as to the unknown parameter a, we can approximate the solution to the system (1) with a one-step method by treating it as an initial value problem. We divide the interval 0 ≤ x ≤ 1 into N subintervals with x j = hj, where the width of each subinterval is h = 1 N . We identify the approximate solution to the system (1) with guess a k using N subintervals as Given an approximate solution v k, N corresponding to the unknown parameter a k , we identify the error in the first derivative at x = 1 as e k, N : is obtained directly from the approximate solution of the differential equation. Since we seek a such that u (1) = 0, the representation of the error quickly reduces to We then update the parameter a k with the usual secant method We expect that e k, N → 0 for fixed N as k → ∞. Table 1. The algorithm requires two initial guesses of the unknown parameter and combines a one-step method with a root-finding method until error is below a stated tolerance.
Step Input Scheme Output In practice, we begin with a fixed discretization of the spatial interval and two initial guesses of the unknown parameter represented by a 0 and a 1 . As depicted in Table 1, we initialize the algorithm by finding an approximation to the solution corresponding to the guesses of the unknown parameter given by a 0 and a 1 , and compute the corresponding error terms e 0, N and e 1, N . The algorithm continues in a two-state process (i.e. obtain a k+1 via the secant method, then perform the one-step method with a k+1 to obtain the corresponding error term e k+1, N ) until the error falls below a stated tolerance, at which point the method is terminated.
Examined alone, convergence properties of a given one-step approach to an initial value problem are well understood. Similarly, the rate of convergence of the secant method is well known. Below, we demonstrate numerically that these methods employed in series, where the secant method changes a parameter of the differential equation in question, exhibit the expected convergence properties.

Numerical results
We begin with the observation that the assumption g(x) ≡ 1 reduces the system (1) to one equivalent to those considered by Landsberg and Fowkes (1978) and by North et al. (2013). Under the assumption g(x) ≡ 1, an elementary solution to the system (1) in terms of the unknown parameter a is readily available.
In this so-called base case, we compare the results of the hybrid numerical method to the elementary solution and show that both the root-finding method and the one-step method converge as expected. As an intermediate case, we use a linearly decreasing function to represent the geometry of the leaf for which an elementary solution to the system (1) can be expressed in terms of special functions. In this case, we also show that the root-finding method and the one-step method converge as expected.
To better incorporate the nonconstant width of the leaf into the model, we approximate the width as a function of the distance along the blade with a cubic spline. In this case, an elementary solution is no longer available. Nevertheless, we show that the one-step method converges as expected.
In many of the numerical experiments below, we use v 0 = .8, which was derived in Section 2 as a representative value for the dimensionless inward flux. In other numerical experiments, we use a range of values of v 0 that correspond to a broad range of values of the potential difference that drives the flow.  Figure 1 is oriented such that x = 0 corresponds to the tank-blade boundary and x = 1 corresponds to the end of the leaf.
Note: In the base case, the width w of the leaf is assumed to be constant.

Base case
We assume that the width of the leaf is constant along the length of the blade so that g(x) ≡ 1 as illustrated in Figure 2.
In this case, the solution to the system (1) is given by The value of a for which u (1) = 0 satisfies the transcendental equation Since the flux at the boundary x = 0 is inward, the boundary condition u (0) = −v 0 has v 0 > 0, so there is a unique solution to Equation (15) satisfying a > 0 for all v 0 . Given v 0 , the reference value of a for which u (1) = 0 is quickly determined in practice by Newton's method.
To examine the convergence properties of the root-finding method within the hybrid numerical scheme, we fix the grid size N and define the error between the desired reference value of a and the kth approximation determined by the root-finding method as e k = a−a k . By definition, a pth order root-finding method satisfies where the constant C describes the rate of convergence. Therefore, for a pth order method, the ratio |e k |/|e k−1 | p should approach a constant as k gets large, which reveals the order of the method.
In particular, the root-finding method we use within the hybrid numerical method is the secant method, which is well-known to be an order φ method, where φ = 1+ √ 5 2 is the golden ratio. The data in Table 2 show that the secant method employed within the hybrid numerical scheme indeed converges with order φ where C ≈ .9 is the rate of convergence. We note that rounding error affects the computation of |e k |/|e k−1 | φ when k = 9, leading to the seemingly inconsistent value shown in the final ratio. These data are obtained using a fourth-order Runge-Kutta method on a grid with N = 640 and the flux condition u (0) = −.8. Since the convergence of the root-finding method depends on the random initial guesses, this behavior is often difficult to illustrate. In order to examine the convergence properties of the one-step method, we make use of the elementary solution (14). If the one-step method is an nth order method (i.e. the local truncation error is O h n as h → 0), then we expect the approximate solution v k, N corresponding to guess a k to satisfy v k, N (x j ) = u(x j , a k ) + O h n . We note that the structure of the system (1) ensures that oscillations are impossible. Since u ≥ 0 and u (0) < 0, the additional constraint on the parameter a (namely, identify a such that u (1) = 0) guarantees that u is a monotonically decreasing function of x. Since the approximate solution is determined by a one-step method, the error in the approximate solution at x = 1 is the largest error on the interval 0 ≤ x ≤ 1.
Suppose that the root-finding method (as described in Section 3) converges to a stated tolerance in K steps. Let U N := v K, N represent the approximate solution to the system (1) using the parameter a K and N grid points. If U N is obtained through an nth order one-step method, then we expect that E N , the difference between the approximate solution and the elementary solution evaluated at x = 1, will have the property that The data in Table 3 show the results of the two-level refinement study designed to illustrate nth order convergence of three separate one-step methods. The refinement study compares the error of the numerical approximation using N grid points to the error of the numerical approximation using 2N grid points. When using an nth order one-step method, the ratio of absolute values of successive errors given by R N satisfies Therefore, the base-two logarithm of the ratio of absolute values of successive errors, applied to an nth order method should reveal the order of the method, namely log 2 R N = n.
To examine convergence properties of the one-step portion of the hybrid numerical scheme, we use three different one-step methods. Euler's method is a first order method, while the two-point Runge-Kutta (RK2) and four-point Runge-Kutta (RK4) methods are second order and fourth order, respectively. The data in Table 3 under the column 1.645e-05 1.398e-10 a The RK4 method is near machine precision at N = 2560, which makes further refinements in the two-level approach inaccurate.
heading log 2 R N illustrate first, second, and fourth order methods, respectively. Although the behavior shown in Table 3 is for the flux condition u (0) = −.8, this behavior is representative of the hybrid numerical method for a wide range of v 0 .

Intermediate case
An immediate constitutive form for the function modeling the geometry of the leaf in the model system (1) is g(x) = 1 − x. This geometry, depicted in Figure 3, models the width of the leaf as a linearly decreasing function of the length along the blade, with zero width at the tip of the blade. In this case, the elementary solution to the model system (1) can be written as the linear combination where Ai and Bi are Airy functions of the first and second kind, and the parameter α satisfies α 3 = a 2 . The Wronskian of the Airy functions is equal to 1/π, so the constants that satisfy the boundary conditions at x = 0 are To satisfy the no-flux condition at x = 1, we choose the parameter α to satisfy In practice, a reference value for α (equivalently, for a) is determined quickly by Newton's method.  Convergence results similar to that found in Section 4.1 are shown in Tables 4 and 5. Specifically, Table 4 demonstrates that the secant method, operating within the two-state hybrid method, converges with order φ where C ≈ .7 is the rate of convergence. Again, we see that some rounding error exists in the computation of |e k |/|e k−1 | φ when k = 9, which affects the final ratio in the table. Table 5 shows the results of a two-level refinement study applied to the algorithm for u (0) = −.8. The data shown in the columns of Table 5 under the heading log 2 R N illustrate first, second, and fourth order convergence of the Euler, RK2, and RK4 methods, respectively.

General case
To better mimic the shape of the leaf with a nonconstant function g, we obtain numerical values for the function g by reading ten data points into WebPlotDigitizer (Rohatgi, 2017), as shown in Figure 4, and interpolating these points with a cubic spline, which gives a numerical representation for the function g. In this case, the elementary solution to the system (1) is no longer available, so a reference value of the parameter a that satisfies the condition u (1) = 0 cannot be determined. Therefore examination of the convergence properties of the root-finding method is not possible.
Nevertheless, examination of the convergence properties of the one-step method can be achieved by using approximate solutions of the (1) in a three-level refinement study. When Table 5. The results of the two-level refinement test show that the various one-step methods, operating within the hybrid numerical scheme, converge to the exact solution of the differential equation as expected in the intermediate case. 2.298e-05 6.896e-10 a The RK4 method is near machine precision at N = 1280, which makes further refinements in the two-level approach inaccurate.
x 1 0 Figure 4. In order to determine a function g that represents the scaled width of the leaf as a function of length along the leaf, data points representing the boundary are captured (left) by WebPlotDigitizer.
Note: These points are fit with a cubic spline (right) and used in the numerical approximation of the solution of the model system.
the value of the parameter a (again, determined by the root-finding method described in Section 3) has converged to a given tolerance in K iterations, we again define U N := v K, N as the approximate solution using that parameter value and N grid points. For an nth order method, we expect the difference D N between the approximate solution using N and 2N 1.251e-05 5.983e-10 a The RK4 method is near machine precision at N = 1280, which makes further refinements in the three-level approach inaccurate.
grid points to satisfy In this case, when R N represents the ratio of absolute values of successive differences, we expect it to satisfy Again, the base-two logarithm of R N reveals the order of the method, log 2 R N = n. The data in Table 6 show the results of the three-level refinement test designed to illustrate convergence of the various one-step methods in the general case incorporating the geometry of the leaf. Similar to the previous refinement test, the data in Table 6 under the column heading log 2 R N illustrate first, second, and fourth order methods, respectively.
Once again, these data are obtained for the flux condition given by u (0) = −.8. Nevertheless, the method exhibits the same behavior for a robust range of values of the parameter v 0 . This robust behavior is illustrated in in Figure 5, which shows the base-two logarithm of R N once the method has converged against a variety of flux conditions, v 0 .

Discussion
The numerical method discussed in Sections 3 and 4 was motivated by the observation that the width of the leaf of G. lingulata is a nonconstant function of the length along the blade, suggesting a nonconstant function g in the system (1). When g(x) ≡ 1 (the so-called base case), the solution of the system (1) can be found using elementary methods so that a reference value of the parameter a that satisfies the additional condition u (1) = 0 can be determined by a basic root-finding method. As such, with a nonconstant g in mind,  Figure 6. The graph shows numerical approximations to the scaled potential u as a function of the scaled spatial variable x where x = 0 represents the tank-blade boundary and x = 1 represents the tip of the leaf.
Notes: For the same value of the flux condition v 0 , the scaled potential in the model which considers the geometry of the leaf is higher than the so-called base case, with the case where g is linearly decreasing along the length of the blade showing the largest difference. Note: For a robust range of the flux condition v 0 , that ratio is greater than one. the base case gave proof of concept for the numerical method used in the more general case. In an intermediate case where g is a linearly decreasing function of the length along the blade, the solution of the system (1) can be found in terms of Airy functions. Once again, a reference value of the parameter a that satisfies the additional condition u (1) = 0 can be determined by a basic root-finding method. In a more general case, the function g representing the geometry of the leaf is determined numerically through a cubic spline, so those elementary solutions are no longer available. Nevertheless, the one-step portion of the hybrid numerical method converged as expected, yielding approximations to the desired solution u and reference value a, which achieves the stated goal of decomposing overall plant conductance into components in the axial and radial directions.
In the context of the model plant, G. lingulata, using the cubic spline geometry yields a value of the scaled water potential u that is between the value found by the base geometry and the value found by the intermediate geometry. Figure 6 illustrates this behavior for the parameter v 0 = .8, which is representative of a robust range of that parameter. Since v 0 is a proxy for axial resistance (see Equation (9)), an increase in the leaf water potential for the same flux condition v 0 in the cubic spline case gives rise to a greater potential difference between the xylem and the mesophyll, suggesting that the radial conductance may be higher in the cubic spline case than the base case.
Indeed, in comparing the value of the parameter a found in the base case to the value in the general case, we see agreement with this analysis. Since a 2 is inversely proportional to radial resistance (see Equation (10)), it serves as a proxy for radial conductance. Figure 7 illustrates radial conductance, as determined by the general case and the base case, showing that the ratio of the value found in the general case to the base case is greater than one for a robust range of flux condition v 0 . This ratio, which is valid for a large range of v 0 , illustrates that accounting for the leaf geometry increases the component of the total leaf conductance in the radial direction by as much as 25%.

Disclosure statement
No potential conflict of interest was reported by the authors.