Numerical Solution of Stieltjes Differential Equations

This work is devoted to the obtaining of a new numerical scheme based on quadrature formulae for the Lebesgue–Stieltjes integral for the approximation of Stieltjes ordinary differential equations. This novel method allows us to numerically approximate models based on Stieltjes ordinary differential equations for which no explicit solution is known. We prove several theoretical results related to the consistency, convergence, and stability of the numerical method. We also obtain the explicit solution of the Stieltjes linear ordinary differential equation and use it to validate the numerical method. Finally, we present some numerical results that we have obtained for a realistic population model based on a Stieltjes differential equation and a system of Stieltjes differential equations with several derivators.


Introduction
In this work, we present a numerical method in order to approximate the solution of a Stieltjes differential equation of the type:    x g (t) = f (t, x(t)), for g almost every t ∈ [0, T), where x g is the Stieltjes derivative with respect to a left-continuous non-decreasing function g. That is, given x : [0, T] → R, we define, for each t ∈ [0, T] \ C g , x g (t) as the following limit in the case it exists: , if t / ∈ D g , where D g denotes the set of discontinuities of g. In this particular case, and C g = {s ∈ R : g is constant on(s − ε, s + ε) for some ε ∈ R + }.
Defining Equation (1) for "g almost every t ∈ [0, T)", we are implicitly considering the Lebesgue-Stieltjes measure space ([0, T], M g , µ g ), where M g is the σ-algebra and µ g the measure constructed in an analogous fashion to the classical Lebesgue measure, where the length of [a, b) is given by µ g ([a, b)) = g(b) − g(a). The interested reader may refer to [1] for details concerning this measure space. The theoretical study of these kinds of derivatives and their applications appears, for instance, in [1][2][3][4][5], but can be traced back in the more general setting of fractional derivatives to [6]. More recent works in this direction include [7], where the authors studied Ψ-fractional derivatives. As stated in [2]  Hypothesis 1. f (·, x) is g-measurable for every x ∈ R; Hypothesis 2. f (·, x 0 ) ∈ L 1 g ([0, T)), where L 1 g ([0, T)) is the space of the equivalence classes of absolutely integrable functions with respect to the g Lebesgue-Stieltjes integral; Hypothesis 3. there exists L ∈ L 1 g ([0, T), [0, ∞)) such that for g almost every t ∈ [0, T) and every x, y ∈ R, we have that: Then, Problem (1) has a unique solution in the space BC g ([0, T]) of bounded g-continuous functions u : [0, T] → R, that is the solution u satisfies, for every t 0 ∈ [0, T], BC g ([0, T]) is a Banach space with the supremum norm-see [2] (Theorem 3.4). Furthermore, the solution of Problem (1) is the unique fixed point of the operator where given t ∈ [0, T], f (s, x(s)) d µ g ; (8) that is, the solution of Problem (1) is such that Furthermore, from [2] (Lemma 7.2) and [1] (Theorem 5.4), the solution will belong to the space AC g ([0, T]) of g-absolutely continuous functions, that is, of those functions u : [0, T] → R such that, for every > 0, there exists δ > 0 satisfying that, if {(a n , b n )} n∈N is a collection of pairwise-disjoint open intervals such that N ∑ n=1 |g(b n ) − g(a n )| < δ, (10) then, It is precisely Expression (9) that motivates the approximation based on the quadrature formulae for the Lebesgue-Stieltjes integral, which we introduce in Section 2. We will see that, in order to obtain error bounds, it will be necessary to impose additional conditions on the regularity of the function f and the solution of Problem (1).
In order to conveniently organize this work, in Section 2, we obtain some numerical quadrature formulae for approximating the Lebesgue-Stieltjes integral. In Section 3, we present a predictor-corrector method based on the quadrature formulae obtained in Section 2. In Section 4, we analyze mathematically the consistency, convergence, and stability of the numerical method derived in Section 3. In order to validate the numerical method, in Section 5, we obtain the explicit solution of the general linear equation of the Stieltjes type. Finally, in Section 6, we present some numerical results that we have obtained for the general linear equation and for realistic population models based on a Stieltjes differential equation.

Quadrature Formulae for the Lebesgue-Stieltjes Integral
We now introduce some convenient notation. Given an increasing left-continuous function g : [a, b] → R, we define ∆ + g : [a, b) → R as ∆ + g(t) = g(t + ) − g(t). In the same way, we define ∆ − h(t) = h(t − ) − h(t) whenever the left limit of h exists at t. Clearly, g is continuous at t 0 ∈ [a, b) if and only if ∆ + g(t 0 ) = 0. It is important to point out the fact that the interval [a, b] contains its right extremity, implying that the function g is bounded. We have that so g has a countable number of discontinuities, say those in D g = {d k } k∈Λ , where Λ ⊂ N. If we define the bounded increasing function g B : [a, b) → [0, +∞) as it is clear that g C : [a, b) → R, given by g C (t) := g(t) − g B (t), is bounded, increasing, and continuous.
We say g C is the continuous part of g and g B is the jump part of g.
As we mentioned in the previous section, the numerical method we propose to approximate the solution of the differential problem (1) in its integral form (9) will be based on the approximation of the Lebesgue-Stieltjes integral. We start this section by proving a result that will allow us to interpret the integral in (9) in terms of a Kurzweil-Stieltjes integral for which it will be possible to establish quadrature formulae. The Kurzweil-Stieltjes integral generalizes the Lebesgue-Stieltjes integral. On the one hand, it uses the gauge of the Henstock-Kurzweil integral in order to obtain an integral that is more powerful than that of Riemann (in fact, a function is Lebesgue integrable if and only if the function and its absolute value are Henstock-Kurzweil integrable); on the other hand, it uses the Stieltjes distribution function g in order to put a weight in the integral. The reader will find a comprehensive study of the Kurzweil-Stieltjes integral in [8].
where on the right-hand side, we consider a Kurzweil-Stieltjes integral. Furthermore, if {d k } k∈Λ is the set of discontinuities of g in (a, b), we have that Proof. This lemma is an immediate consequence of [8] (Theorems 6.12.3 and 6.3.13). That is, since f is Lebesgue-Stieltjes integrable, by [8] (Theorem 6.12.3), it is Kurzweil-Stieltjes integrable as well, and furthermore, Now, due to the fact that g is left-continuous, ∆ − g(b) = 0, and since µ g ({a}) = ∆ + g(a), we obtain: Finally, by [8] Since g is left-continuous, we have, in particular, that ∆ − g(b) = ∆ − g(d k ) = 0 for every k ∈ N, and the desired result follows.
In Lemma 2, we will see that, under certain regularity hypotheses on f and g, we can obtain error estimates for the quadrature formula for a point and the trapeze formula.
where H > 0 and p ∈ (0, 1]. Then, and and Indeed, we can adapt the techniques in [9] for the Riemann-Stieltjes integral to the case of the Kurzweil-Stieltjes'. On the one hand, by [8] (Theorem 6.3.6), given that h : On the other, thanks to [8] (Theorem 6.4.2) (integration by parts), from which, given that g C is continuous, In particular, We have that for every t ∈ [a, b]. Then, thanks to the bound (24), we obtain the bound (23). In order to prove (22), we can proceed in an analogous fashion by integrating by parts where we have already canceled out the terms concerning the sum. From the previous expression, we obtain As we will see later on, it will be of special interest to consider the case when D f ⊂ D g and f C behaves in a similar way to g C . In such a case, we can sharpen the previous quadrature formulae to obtain Lemma 5.  Proof. It is clear that f is g-continuous. Since g is bounded and | f (t) − f (s)| ≤ H|g(t) − g(s)| for every t, s ∈ [a, b], f is bounded as well.
The g-integrability is a straightforward consequence of the definition of the Riemann-Stieltjes integral and the fact that f is g-continuous.
), so f is of bounded variation.
The error estimates obtained in the previous formulae are not enough for our proposes, that is proving the convergence of the numerical approximation to the solution of Problem (1). In order to improve the previous estimations, we must add some extra requirements to the continuous part of functions g and f .
In the next lemma and corollary, we will prove that if f is a g-Lipschitz continuous function, some properties of g C and g B are transferred to f C and f B , respectively. In particular, we will see that f is a g-Lipschitz continuous function and g C is Lipschitz continuous, then f C is also Lipschitz continuous. This property will be fundamental in order to improve the previous quadrature formula.
For the next lemma, we denote by C(X) the set of connected components (in the usual topology) of X ⊂ R. Furthermore, if f is increasing, then f C is g C -Lipschitz continuous with Lipschitz constant H.
Thus, taking the limit when s tends to t from the right, Therefore, f B is g B -Lipschitz and g-Lipschitz with constant H.
We know that ∑ t∈[a,b) ∆ + g(t) < ∞. This implies, on the one hand, that D g = {t n } n∈Λ with Λ ⊂ N is countable and, on the other, that ∆ + g(t n ) → 0. Observe that g is continuous at b, and either g is continuous at a or a = t k for some k ∈ Λ. In this last case, we will assume, without loss of generality, that a = t 1 .
Thus, consider, for n ∈ N, the functions Given A ∈ C([a, b]\{t k } n k=1 ) (that is, an interval that is a connected component of [a, b]\{t k } n k=1 )) and t, s ∈ A, t < s, since there are no jumps of Since F n (t, s) ≥ 0 for every t, s ∈ A and f n and g n are continuous at the points of ∂A, it also holds for t, s ∈ A. Furthermore, F n (t, x) We conclude that F n ≥ 0.
Observe now that f n converges uniformly to f C and g n converges uniformly to g C , so F n converges uniformly to Since F n ≥ 0 for every n ∈ N, F ≥ 0, and thus, f C is g C -Lipschitz continuous with Lipschitz constant H. Proof. Since f is g-Lipschitz continuous, it is g-absolutely continuous and, by [1] (Theorem 5.4), there exists f g µ g -a.e. and Thus, by the definition of the Stieltjes g-derivative, Clearly, both f 1 and f 2 are g-Lipschitz continuous with Lipschitz constant H and increasing, so f C 1 and f C 2 are g C -Lipschitz with Lipschitz constant H. Thus, Proof. This is a direct consequence of the previous corollary since for all In order to simplify the notation, from now on, we will assume, when necessary, that both the continuous part of f and the continuous part of g are Lipschitz continuous with the same Lipschitz constant H (if necessary, we redefine H to be max{H, H 2 }). We also assume that f C and g C are Lipschitz continuous with Lipschitz constant H. Then, and where {d k } k∈Λ is the set of discontinuities of g in (a, b).
Proof. First observe that, since f is g-continuous, D f ⊂ D g -see [2] (Proposition 3.2). Separating the jump part from the continuous part in both f and g, we have that: where the three first integrals correspond to series and the fourth can be approximated using the previous quadrature formulae. Indeed, using analogous reasoning as that in [8] (Theorems 6.3.12 and 6.3 Now, using the same argumentation as before, since . The proof of identity (44) is analogous.

Description of the Numerical Method
In this section, we present a predictor-corrector method based on the previous quadrature formulae. We will assume g : [0, T] → R is continuous at t 0 = 0 and that D g is finite. Let x ∈ BC g ([0, T]) be a solution of problem (1), and in what follows, let f * (x)(t) := f (t, x(t)), Consider now a set {t k } N+1 k=0 ⊂ [0, T] satisfying: Hypothesis 5. f * (x) is g-continuous, and in particular, D f * (x) ⊂ D g and f C * (x) and g C are Lipschitz-continuous with the same Lipschitz constant H.
Let us define x k := x(t k ) and x + k := x(t + k ), for k = 0, . . . , N + 1. By the definition of the Stieltjes derivative, and in the particular case t k ∈ D g , ∆ + g(t k ) = 0, so we have that x + k = x k . Then, for every k = 0, . . . , N + 1, where the integral is of the Kurzweil-Stieltjes type. Using (44) on each interval, Thus, in the case that we use (45), we have: Observe that the condition D g ⊂ {t k } N+1 k=0 implies that, on each interval, the quadrature formulae lose the terms related to the interior jumps. Restricting to [t k , t k+1 ], we have: Hence, and : Taking into account the previous formulae, it is clear that, if we want to use (54) to approximate (52), the a-priori ignorance of the value x k+1 forces us to estimate it. In order to do this, we just have to use (53), for which no further estimation is needed. That is, we will use (53) as the predictor and (54) as the corrector. Thus, the method will be as follows.

Error Analysis
In this section, we analyze the numerical method introduced in the previous one. As happens with those numerical methods based on quadrature formulae-cf. [10,11]-it will be crucial at this point to study the error of approximating the integral in this way.
As before, we will need certain regularity hypotheses on the derivator g, as well as on the function f when composed with the solution of the problem. Thus, we will assume Hypotheses 1-3 since they are necessary to guarantee the existence of Problem (1)-see page 2; Hypotheses 4-6, established in the previous sections in order to formulate the numerical method and the following additional hypotheses for proving the convergence of the method: Hypothesis 8. f (t + , ·) ∈ C 1 (R) for every t ∈ R, and there exits K 3 ∈ R such that We must emphasize that the above hypotheses are not independent. For example, Hypothesis 6 implies Hypotheses 1 and 2, and Hypothesis 7 implies Hypothesis 3. Therefore, for our purposes, it is sufficient that Hypotheses 4-8 are fulfilled. We will now establish the basic notions related to the truncating error associated with the quadrature formula of the predictor-corrector method.
Definition 2 (Local truncating error associated with the quadrature formula). Given a partition P = {t k } N+1 k=0 satisfying (H4), we define: • The local error associated with Formula (44) as: with k = 0, . . . , N. • The local error associated with Formula (45) as: , we define the local truncating error of the predictor-corrector method associated with (58), in terms of the exact solution, as: with k = 0, . . . , N.

Remark 2.
It is usual to find in the literature local truncating errors defined relative to the discretization step h, that is σ * [10,11]. In those cases, for k = 0, . . . , N, We opted for the first set of errors in order to simplify the notation. In any case, the relation between both definitions is clear. Now, we present some bounds of the errors. Lemma 6. For every k = 0, . . . , N, we have the following bounds: Proof. The two first assertions are a direct consequence of Lemma 5. In order to obtain the third one, we manipulate the definition of τ k+1 , leaving: from which we obtain, using the definition of σ k+1 , By the mean value theorem of differential calculus, there exists c k+1 in the open interval of extremities x k+1 , x * k+1 such that Then, taking the absolute value, Using the bounds obtained for |σ k+1 | and |σ * k+1 |, Corollary 4 (Consistency of the numerical method). In the functional framework in which Lemma 6 is valid, the method is consistent.

Corollary 5.
Indeed, thanks to the bounds provided by Lemma 6, we obtain: from which we deduce the consistency of the method in the classical sense.

Remark 3.
In view of the bounds in Lemma 6, we must observe that the introduction of a predictor in the quadrature formula does not penalize its convergence order. This is due to the fact that σ * k+1 (which is the predictor term in the formula) appears multiplied by h in (66).
In our case, due to the regularity of the terms involved, we are not capable of improving the order of convergence of the two-point formula with respect to the one-point one. This is not usually the case, as in the literature, we can see examples-for instance [10,11]-where the two-point quadrature formula has a better convergence order-without the predictor penalizing the global order of the method-than the one-point one.
Definition 3 (Local error of the algorithm). We define the following errors associated with the numerical algorithm.
• e + k := u + k − x + k , with k = 0, . . . , N, is the local error of the corrector regarding the limit from the right at t k . It is clear that it does not make sense to consider this error for k = N + 1.
. . , N + 1, is the local error of the predictor at the point t k . • e k = u k − x k , with k = 0, . . . , N + 1 the local error of the predictor at the point t k , and e 0 is the error associated with the initial condition.
In Lemma 7, we obtain bounds for the previous lemmata based on recurrence formulae, which, afterwards, we will analyze in order to obtain the bounds of the error at each of the points of the temporal discretization.

Lemma 7.
Under the hypotheses of Lemma 6, we derive the following formulae for e * k+1 , e k+1 and e + k , with k = 0, . . . , N, Proof. We compute each of the error bounds separately.
• Local error of the corrector regarding the limit from the right: We have that: where c k belongs to the open interval of extremities u k and x k . Taking the absolute value on both sides, • Local error of the predictor at the point: We have: where c + k belongs to the open interval of extremities u + k and x + k . Taking the absolute value on both sides, • Local error of the predictor at the point: We have that where c + k belongs to the open interval of extremities u + k , x + k and c * k+1 belongs to the open interval of extremities u * k+1 and x * k+1 . Taking the absolute value on both sides, where: Remark 4. Observe that previous error formulae can be simplified in the case t k / ∈ D g . In this situation, those errors concerning the limit from the right coincide with the ones of the corrector at the point, and we recover the classical error formulae.
From the formulae in Lemma 7, we can prove the following result concerning the error of the numerical method.
where τ = max{|τ k | : k = 1, . . . , N + 1}. Proof. Using the notation of Lemma 7, we have that Thus, applying the previous bound recursively, Accounting for the number discontinuities of the derivator (which we denote by #D g ), we obtain: Now, taking into account that, for a given number G ≥ 0, and that 1 + G ≤ exp(G), we have: Now, we will prove the main theorem of this section. In it, we will see that, in the framework of the previous results, we can guarantee the convergence of the method introduced in the previous section.
Theorem 1 (Convergence of the predictor-corrector method). Under the hypotheses of Lemma 5, if we assume e 0 = 0, we have, for a given t j ∈ [0, T], that Furthermore, we get the following error bounds: for every j = 1, . . . , N + 1, where Proof. We analyze each case separately.
• Errors associated with the corrector. From the previous lemma: where when h → 0. Hence, given t j ∈ [0, T], we get Then, given that |e 0 | = 0, we have the convergence of the corrector to the solution of the problem: • Errors associated with the predictor: Using the bounds in Lemmas 7 and 5, we have that, given t j ∈ [0, T], from which we obtain the convergence. • Errors associated with the right limit: Using the bounds in Lemmas 7 and 5, we have that, given t j ∈ [0, T], We obtain the same kind of convergence as in the previous case. It is worth noting that, in the case t j / ∈ D g , the error associated with the right limit coincides with the error of the predictor.

Remark 5.
Observe that the order of convergence of the method equals the order of τ minus one, that is of the order of τ. In the case that we deal with functions with extra regularity, we may be able to improve the order of τ, which would better the order of convergence of the numerical method. Last, we would like to mention that the method we presented generalizes the classical order two Runge-Kutta. This assertion is motivated by the fact that the usual derivative is a particular instance of the Stieltjes derivative in the case g(x) = x.
Last, we analyze the stability of the method with the intention of evaluating its sensitivity towards the perturbations generated by the rounding errors produced while evaluating the different elements of Scheme (58). We omit the proof of the following result, as it is essentially a modification of that of Theorem 1.

The General Linear Equation
In order to validate the numerical approximation of the solution of Problem (1), we will consider the following general linear equation as a test problem: Under (97) and (98), we know there is a unique solution of (96), which can be computed explicitly-see [2]-as the unique solution of the problem where are L 1 g ([0, T)) functions thanks to [2] (Proposition 6.8). Therefore, by [2] (Proposition 6.7), the solution of Problem (99) is given by where given an element c ∈ L 1 g ([0, T)), {s 1 , . . . , s N } = {t ∈ [0, T) ∩ D g : 1 + c(t)∆ + g(t) < 0} being the set of points such that 1 + c(t)∆ + g(t) < 0 and s N+1 = T. This set has finite cardinality-see [2] (Lemma 6.4). In our case, c = d; thus, if and only if 1 < d(t)∆ + g(t). We will still denote by and s N+1 = T, so As we can see above, the general expression of the exponential e d (t), and therefore, of the solution of the general linear Equation (102), has a convoluted statement. This expression can be simplified if we consider the particular case where d is constant and d∆ + g(t) < 1, ∀t ∈ [a, b) ∩ D g -the case that we will consider in the numerical experiments. In that case, It is also remarkable that, in the case of g(t) = t, we recover the classical exponential. Furthermore, we have the following direct result for the problem with constant coefficients. Theorem 3. Let g : [0, T] → R be increasing, left-continuous, and such that g(0) = 0; x 0 ∈ R, h ∈ R, and d ∈ R such that Then, the solution of the problem is given by the following expression: If we assume that the set of discontinuities of function g is finite and we consider a time discretization {t k } N+1 k=0 ⊂ [0, T] satisfying hypothesis (H4) with d∆ + g(t) < 1, ∀t ∈ [a, b) ∩ D g , then the condition ∑ t∈[a,b)∩D g |ln (1 − d∆ + g(t))| < ∞ is trivially satisfied. Therefore, we have the following corollary for the homogenous case. The proof is straightforward. Corollary 6. Let g : [0, T] → R be increasing and left-continuous, such that g(0) = 0 with a set of discontinuity points that we can assume equal to the discretization points, that is D g = {t 1 , . . . , t N } ⊂ (0, T), where t k < t k+1 for k = 1, . . . , N − 1. Then, the solution of the problem where x 0 , d ∈ R and d∆ + g(t) < 1, ∀t ∈ [a, b) ∩ D g , is given by: (1 − d∆ + g(t k )) exp(d∆ + g(t k )), t ∈ (t n , t n+1 ]; n = 1, . . . , N − 1, Finally, it is remarkable that, in the previous corollary, we can change the hypothesis d∆ + g(t) < 1, ∀t ∈ [a, b) ∩ D g by d∆ + g(t) = 1, ∀t ∈ [a, b) ∩ D g , and obtain a similar expression for the solution taking into account the general Formula (102). The last hypothesis is more general than the previous one, but in order to present the results in a clear way, we will assume that the first hypothesis is fulfilled.

Numerical Simulations
In this section we will present some numerical results that we reached using the scheme (58) for approximating the solution of the homogeneous linear Equation (113) with constant coefficients. We will also compare the numerical solution with the explicit solution (113) that we obtained in the previous section. To test the robustness of the method, we will use the numerical scheme to approximate the solution of a silkworm population model based on the example presented in [3]. Finally, we approximate the solution of a system of Stieltjes differential equations with several derivators related to a fish population model subject to fishing with open and closed seasons-cf. [12]. All numerical simulations were performed using MATLAB 2019b software on a 2.4 GHz Intel Core i5 CPU with 8 GB of RAM.

Approximation of the General Linear Equation
In order to validate the scheme (58) for different numbers of discontinuities in the derivator g (the main difficulty of the problem), we will consider an increasing regular continuous part g C , and we will obtain several g test functions adding to the previous one the jump part g B associated with several choices of jumps. We consider the following function: where α > 0. We have that ϕ is an increasing C ∞ (R) function, and we can use it to construct a more sophisticated increasing function g C that will be constant in some intervals. For instance, for α = 4, we can consider the following function g C in the time interval [0, 10] and, from it, build the derivator g by adding the jump function g B .
In Figure 1a, we observe that we have concatenated three times the function ϕ, and in order to obtain the derivator function g, we added four jumps at the timest k = 2k, k = 1, 2, 3, with ∆ + g(t k ) = 1. In Figure 2a, we plot the solution for d = −0.5 and in Figure 2b the solution for d = 0.5. As we can see in both figures, we have inactivity periods where the function g is constant and impulses at the times where the function g presents discontinuities.
(a) Continuous part g C .
(b) Derivator g.  We summarize the results obtained for different values of time step h taking x 0 = 1, d = −0.5, and for different values of #D g , with ∆ + g(s) = 1, ∀s ∈ D g . CPU times were calculated by averaging the CPU times of ten simulations in order to minimize the influence that other processes may have on the measurement.
From the Table 1, we can observe that numerical errors grow as the number of discontinuities in the derivator increases. This behavior is consistent with the error bounds obtained in Theorem 1 in which the term [1 + G 2 ] #D g appears by multiplying the error expressions. In Figure 3a, we can observe the error evolution for the predictor and, in Figure 3b, the error evolution for the corrector. We realize that the global behavior in terms of h for the predictor is O(h) and O(h 2 ) for the corrector. This improvement in the order of convergence with respect to the one predicted in theory is a consequence of the fact that, thanks to the regularity of the solution, the trapezoidal formula is more accurate.
(b) Corrector error in log 10 − log 10 scale. Finally, in Figure 4a, we can see the exact solution and the predictor for #D g = 2 using as the time step h = 1.e − 01 and, in Figure 4b, h = 1.e − 05. In Figure 4c,d, we have the analogous graphical representations for #D g = 4.

Approximation of a Silkworm Population Model
We present in this section the numerical approximation of a realistic case that corresponds to a silkworm population model based on the example presented in [3], which we will briefly summarize for the convenience of the reader. In this example, the authors considered that the life cycle of silkworms has three stages: worm, cocoon, and moth. Moths lay eggs and die soon after, then eggs hatch and produce a completely new colony of silkworms.

Eggs
(5k + 4, 5k + 5], k = 0, 1, 2, . . . In order to take into account the previous behavior, they considered the following derivator g : [0, ∞) → R: and they solved the following Stieltjes equation: with c > 0, λ > 0. In [3] (Proposition 5.1), the authors obtained the explicit solution of the previous model: Now, we approximate the solution of the Stieltjes differential Equation (118) using the scheme (58) in the time interval [0, 15], for λ = 1.1, c = 1.2, and x 0 = 8. We realize that in order to evaluate function (119), we have to approximate the integral value using a classical quadrature formulae, so the convergence order of the full scheme will be penalized by this approximation. In our case, we considered a composite trapezoidal rule. In the following Table 2, we summarize the numerical results that we obtained in this case (we omit the errors for the predictor and the limits from the right): Finally, in Figure 5a, we can see the exact solution and the predictor using as the time step h = 1.e − 01 and, in Figure 6, h = 1.e − 05.
In Table 3, we can see the successive approximations in some selected times. We can see that the largest errors correspond to the points where the impulses appear.   In this section, we present the application of the algorithm previously deduced to a system of Stieltjes differential equations with several derivators: , x(t)), for g i almost every t ∈ [0, T), where x(t) = (x 1 (t), . . . , x n (t)). As stated in [12] (Theorem 4.3), in the case g : [0, T] → R n + , g = (g 1 , . . . , g n ), such that for each i ∈ {1, . . . , n}, g i is nondecreasing and left-continuous, and f : [0, T] × R n → R n satisfying the following hypotheses: Hypothesis 9. f i (·, x) is g i -measurable for every x ∈ R and i = 1, . . . , n; Hypothesis 10. f i (·, x 0 ) ∈ L 1 g i ([0, T)) for every i = 1, . . . , n; Hypothesis 11. there exists L : [0, T] → [0, ∞) such that for every i = 1, . . . , n, L ∈ L 1 g i ([0, T)), and Then, there exists an unique solution x ∈ BC g ([0, T]), where: such that x is the unique fixed point of the operator: where given t ∈ [0, T], In order to adapt the scheme (58) to this problem (we leave out the details as they are completely analogous to those of the case of a single equation), we have to assume the following additional hypothesis: T] is such that t 0 = 0, t N+1 = T, t k+1 − t k = h > 0, for every k = 0, . . . , N and D g i ⊂ {t k } N+1 k=0 , for every i = 1, . . . , n.
Let us see now how we can use the previous scheme to approximate the solution of the model for a fish population subject to fishing with open and closed seasons that appears in [12] (Section 5.1). In this work, the authors identified each year with an interval of the form (4k − 4, 4k], k ∈ N; they supposed that at times t = 4k − 1, k ∈ N, new fishes are born, and they assumed that intervals of the form (4k − 2, 4k], k ∈ N, correspond to the closed season. Finally, they considered the following system for the population size of a fish species p(t) at time t and for the number of captured individuals h(t) at time t (cf. [12] (Section 5.1) for more details): where: and g 1 (t) = 5 + g 1 (t − 4), g 2 (t) = 4 + g 2 (t − 4), for t > 4. The parameter α > 0 represents the death rate of the fish population caused by the environment they live in; β > 0 is the fishing rate; and γ > 0 is the decay rate for h. Parameters 0 < c p , c h < 1 represent the number of newborns of species p and the number of fished individuals at the beginning of the open season. As far as we know, there is no known explicit solution for the system; therefore, the use of numerical techniques is necessary to study the behavior of the system (127). Now, we present some numerical results that we have obtained for T = 10, p 0 = 4, h 0 = 3, α = 0.05, β = 0.1, γ = 0.05, c p = 0.9, and c h = 0.9. In Figure 6, we can see the evolution of the number of a fish species and captured individuals in each time step. In view of the results obtained, we can observe that: • Between two consecutive new generations of fish, p decays. This decrease is due to mortality from natural causes and the effect of fishing.
• As the open season progresses, the number of fished individuals (h) also decreases as a consequence of less resources being available. • Once closed season starts, the number of fished individuals collapses and then remains constant during the closed season.
As there is no known explicit solution, we cannot compare the numerical results with the exact solution.
However, the behavior that we observe in the numerical solution agrees with the interpretation of the system that appears in [12] (Section 5.1).

Conclusions
In this work, we obtained a new numerical scheme based on quadrature formulae for the Lebesgue-Stieltjes integral for the approximation of Stieltjes ordinary differential equations. We proved several theoretical results related to the consistency, convergence, and stability of the numerical method. The techniques that we developed can be generalized to obtain higher order Runge-Kutta-like methods.
In order to validate the numerical method, we compared the numerical and the exact solution for the linear equation. We realize that the order of convergence that we obtained in the numerical experiments is higher than we proved in the theory. This fact may mean that it is possible to improve the estimates obtained under more stringent regularity assumptions regarding the functions involved.
We also compared the numerical and the exact solution of a realistic silkworm population model (118) based on a Stieltjes differential equation. We obtained a good behavior of the numerical method even though we approximated the right-hand side by a quadrature formula. The behavior observed in the solution of Equation (118) (0D model) can also be appreciated in the 2D model based on a parabolic partial differential equation of the Stieltjes type (cf. [13] (Section 5)): where Ω ⊂ R 2 , η > 0, u 0 ∈ L 2 (Ω), and f is given by Formula (119). Finally, we proposed a numerical method to approximate the solution of systems of Stieltjes differential equations with several derivators, and we applied it to a system related to a fish population model for which no explicit solution is known. We observed that the behavior of the numerical solution is compatible with the interpretation of the model given in [12].
The ideas we developed in this work can be used to design an algorithm that improves the numerical approximation of the linear Stieltjes parabolic problems given in [13]. In this work, the authors obtained a numerical approximation of (130) by truncating the Fourier series of the solution and using the exact solution (120) to obtain each of the coefficients of the aforementioned series. Instead of truncating the Fourier series, we can perform a temporal discretization similar to the one that we developed in this work for Stieltjes differential equations. Indeed, if we assume Hypotheses 1-6 in the scalar case, we can consider the following modification of the approximation (53): In this case, we obtain the following implicit numerical method: given u 0 = x 0 , we compute {(u + k−1 , u k )} N+1 k=1 as    u + k = u k + f (t k , u k )∆ + g(t k ), For example, if we follow the same scheme as in this work and consider the following linear Stieltjes parabolic problem: in Ω,