Finite difference method for a Volterra equation with a power-type nonlinearity

TL;DR

This paper proves the convergence of explicit finite-difference methods for nonlinear Volterra equations with power-type nonlinearities, overcoming challenges posed by non-Lipschitz kernels using discrete Gronwall's lemmas.

Contribution

It introduces a convergent explicit finite-difference scheme for nonlinear Volterra equations with non-Lipschitz kernels, including convergence analysis and order bounds.

Findings

Convergent explicit finite-difference methods for the problem.

Upper bounds on the convergence order.

Application to example problems.

Abstract

In this work we prove that a family of explicit numerical finite-difference methods is convergent when applied to a nonlinear Volterra equation with a power-type nonlinearity. In that case the kernel is not of Lipschitz type, therefore the classical analysis cannot be applied. We indicate several difficulties that arise in the proofs and show how they can be remedied. The tools that we use consist of variations on discreet Gronwall's lemmas and comparison theorems. Additionally, we give an upper bound on the convergence order. We conclude the paper with a construction of a convergent method and apply it for solving some examples.

Tables

Table 1. Table 1: Simulated order of convergence of the midpoint rule ( 55 ) applied to ( 56 ). Calculation has been done at x = 0.001 𝑥 0.001 x=0.001 .

m	$1$	$1.5$	$2$	$10$	${10}^2$	${10}^3$
Order	$\infty$	1.000	0.999	0.98	0.927	0.918
Est. order	0	$0.333$	$0.500$	$0.900$	$0.990$	$0.999$

Table 2. Table 2: Simulated order of convergence of the Midpoint rule ( 55 ) applied to ( 58 ). Calculation has been done at x = 0.001 𝑥 0.001 x=0.001 .

m	$1$	$1.5$	$2$	$10$	${10}^2$	${10}^3$
Order	0.999	1.000	1.000	0.980	0.927	0.918
Est. order	$--$	0.332	0.499	0.899	0.990	0.999

Equations

$$u(x)={\int }_0^{x}K(x,t)u(t{)}^{\frac{1}{m+1}}dt,m>0,x\in [0,X].$$

$$y(x{)}^{m+1}={\int }_0^{x}K(x,t)y(t)dt,m>0,x\in [0,X].$$

$$0<C\le K(x,t)\le D,$$

$$x_n=nh,h=\frac{X}{N},n=0,1,2,...,N.$$

$${\int }_0^{x_n}K(x_n,t)y(t)dt=h\sum _{i=1}^{n-1}{w}_{n,i}K(x_n,x_i)y(x_i)+{\delta }_n(h),$$

$$\delta (h):=\underset{1\le n\le N}{\max }|{\delta }_n(h)|,$$

$$0<w_{n,i}\le W,n,i=0,1,2,...,N,$$

$$y_n^{m+1}=h\sum _{i=1}^{n-1}{w}_{n,i}{K}_{n,i}{y}_i,n=2,3,...,N,$$

$$\begin{array}{rl}y((n+1)h{)}^{m+1}& ={\int }_0^{(n+1)h}K((n+1)h,t)y(t)dt=h\sum _{i=1}^n{w}_{n+1,i}{K}_{n+1,i}y(ih)+{\delta }_{n+1}(h)\\ & \ge h\sum _{i=1}^n{w}_{n+1,i}{K}_{n+1,i}{y}_i+{\delta }_{n+1}(h)\ge y_{n+1}^{m+1}.\end{array}$$

$$y(x{)}^{m+1}\lessgtr C{x}^{\mu }{\int }_0^x{t}^{\nu }y(t)dt,$$

$$y(x)\lessgtr {\left(C\frac{m}{m+1}\frac{1}{1+\nu +\frac{\mu }{m+1}}\right)}^{\frac{1}{m}}{x}^{\frac{1+\nu +\mu }{m}},$$

$$z(x):={\int }_0^x{t}^{\nu }y(t)dt,$$

$$y(x)\ge C^{\frac{1}{m+1}}{x}^{\frac{\mu }{m+1}}z(x{)}^{\frac{1}{m+1}},$$

$$z^{\prime }(x)=x^{\nu }y(x)\ge C^{\frac{1}{m+1}}{x}^{\nu +\frac{\mu }{m+1}}z(x{)}^{\frac{1}{m+1}}.$$

$${\int }_0^{x}z(t{)}^{-\frac{1}{m+1}}{z}^{\prime }(t)dt\ge C^{\frac{1}{m+1}}{\int }_0^x{t}^{\nu +\frac{\mu }{m+1}}dt,$$

$$\frac{m+1}{m}z(x{)}^{\frac{m}{m+1}}\ge C^{\frac{1}{m+1}}\frac{1}{1+\nu +\frac{\mu }{m+1}}{x}^{1+\nu +\frac{\mu }{m+1}}.$$

$$y(x)\ge {\left(C\frac{m}{m+1}\frac{1}{1+\nu +\frac{\mu }{m+1}}\right)}^{\frac{1}{m}}{x}^{\frac{1+\nu +\mu }{m}}.$$

$$K(x,t)\sim C{x}^{\mu }\text{as}x\to 0^{+}\text{and}t\le x.$$

$$y(x)\sim {\left(C\frac{m}{m+1}\frac{1}{1+\frac{\mu }{m+1}}\right)}^{\frac{1}{m}}{x}^{\frac{1+\mu }{m}},x\to 0^{+}.$$

$$y(x{)}^{m+1}=K(x,\tau (x)){\int }_0^{x}y(t)dt,$$

$$z(x):={\int }_0^{x}y(t)dt,$$

$$z^{\prime }(x)=y(x)={\left(K(x,\tau (x))z(x)\right)}^{\frac{1}{m+1}}.$$

$$z(x)={\left(\frac{m}{m+1}{\int }_0^{x}K(t,\tau (t){)}^{\frac{1}{m+1}}dt\right)}^{\frac{m+1}{m}}.$$

$$y(x)=K(x,\tau (x){)}^{\frac{1}{m+1}}{\left(\frac{m}{m+1}{\int }_0^{x}K(t,\tau (t){)}^{\frac{1}{m+1}}dt\right)}^{\frac{1}{m}}.$$

$$y(x)\sim {\left(C{x}^{\mu }\right)}^{\frac{1}{m+1}}{\left(C^{\frac{1}{m+1}}\frac{m}{m+1}\frac{1}{1+\frac{\mu }{m+1}}{x}^{1+\frac{\mu }{m+1}}\right)}^{\frac{1}{m}}.$$

$$y_1\ge \left(1-\frac{m+1}{m}\frac{ϵ(h)}{h^{1+\frac{1}{m}}}\right){\left(C\frac{m}{m+1}h\right)}^{\frac{1}{m}},$$

$$y_n\ge \left(1-\frac{m+1}{m}\frac{ϵ(h)}{h^{1+\frac{1}{m}}}\right){\left(C\frac{m}{m+1}nh\right)}^{\frac{1}{m}},n\ge 2,$$

$$y_{n+1}^{m+1}\ge C\left(1-\frac{m+1}{m}\frac{ϵ(h)}{h^{1+\frac{1}{m}}}\right){\left(C\frac{m}{m+1}\right)}^{\frac{1}{m}}h\sum _{i=1}^n{w}_{n,i}(ih{)}^{\frac{1}{m}}.$$

$$\begin{array}{rl}{y}_{n+1}^{m+1}& \ge C\left(1-\frac{m+1}{m}\frac{ϵ(h)}{h^{1+\frac{1}{m}}}\right){\left(C\frac{m}{m+1}\right)}^{\frac{1}{m}}\left({\int }_0^{(n+1)h}{t}^{\frac{1}{m}}dt-{ϵ}_{n+1}(h)\right)\\ & =\left(1-\frac{m+1}{m}\frac{ϵ(h)}{h^{1+\frac{1}{m}}}\right){\left(C\frac{m}{m+1}(n+1)h\right)}^{\frac{m+1}{m}}\left(1-\frac{m}{m+1}\frac{{ϵ}_{n+1}(h)}{{\left((n+1)h\right)}^{\frac{m+1}{m}}}\right),\end{array}$$

$$1-\frac{m}{m+1}\frac{{ϵ}_{n+1}(h)}{{\left((n+1)h\right)}^{\frac{m+1}{m}}}\ge {\left(1-\frac{m+1}{m}\frac{ϵ(h)}{h^{1+\frac{1}{m}}}\right)}^m.$$

Full text

Keywords: Volterra equation, nonlinearity, power-type, finite difference method

1 Introduction

Integral equations are one of most useful tools used in mathematical analysis and modelling. Many essential problems formulated as differential equations can be transformed into the equivalent integral equations which, many times, are more useful for showing existence, uniqueness and deriving estimates. This field has reached its maturity and there is a wealth of literature reviewing many subjects such as the theory for Volterra [5, 18, 13] and Fredholm equations [40, 37, 21], numerical methods [2, 22, 5, 19, 4, 3] and topics in integro-differential equations [12, 31] and fractional calculus [26, 16, 23, 36].

In this work we consider a family of numerical finite-difference methods for solving a certain class of nonlinear Volterra equations. These, in particular, arise as models of dynamics in porous media [41, 29, 30], heat transfer [24], propagation of shock-waves in gas filled tubes [20] and anomalous diffusion [32, 34]. The latter application is the main motivation behind our investigations (more details can be found in [33]). The general form of the family of equations that we consider can be written as follows

[TABLE]

Notice that the nonlinearity is not Lipschitz so that we cannot use the most of the classical theory. However, there are many results concerning existence and uniqueness. In the seminal papers [9, 17] the problem was stated and several conditions for existence and uniqueness were given. Later, in a number of work those results were strengthened and generalized (see [29, 10, 11, 27, 28]). Abel integral equations with power-type nonlinearities have also been recently investigated [8, 42]. For further details see other works of cited authors.

According to the best of authors’ knowledge, the literature consists of only very few papers concerning numerical methods for this kind of equations. In [6] a iterative fixed-point way of solving (1) was proposed. Also, in [7] a short review article concerning theoretical and numerical concepts of nonlinear Volterra equation has been published.

In what follows we consider a family of explicit methods for solving (1) and prove that under certain boundedness assumptions on the kernel they are convergent. Moreover, we find a bound on the convergence rate and illustrate the theory with the numerical simulations.

2 Finite difference method

Consider the following nonlinear Volterra equation which arises from (1) via the transformation $y(x)^{m+1}=u(x)$

[TABLE]

We will need the boundedness assumption on the kernel $K\in C([0,X]^{2})$

[TABLE]

from which only the lower one is a nontrivial requirement. It can be shown that the positive solution of (2) with (3) exists and is unique (see [9] and other papers mentioned in the Introduction).

In order to construct the numerical method fix a natural number $N$ and introduce the grid

[TABLE]

The next step is to discretize the integral in (2), say

[TABLE]

where $\delta_{n}(h)$ is the local consistency error. Next, we define

[TABLE]

and further assume that

[TABLE]

for some $W>0$. Denoting $y_{n}$ as a numerical approximation to $y(nh)$ we may then propose the following explicit finite difference scheme for solving (2)

[TABLE]

where $K_{n,i}:=K(x_{n},x_{i})$. Since the equation (2) has a trivial solution, it is necessary to start the above iteration with a value which will force the convergence to the nontrivial one. In the next section we will show one way of making that choice.

Below we will prove that (8) is convergent to the unique positive solution of (2). Before that, however, we need some auxiliary results. The first is a simple observation concerning the iteration scheme (8).

Proposition 1**.**

Let $y$ be a solution of (2) and construct $y_{n}$ via the iteration (8). If $\delta_{n}(h)$ in (5) is positive (negative) for $n=1,2,...,N$, then $y(nh)\geq y_{n}$ ($y(nh)\leq y_{n}$) provided that $y(h)\geq y_{1}$ ($y(h)\leq y_{1}$).

Proof.

It suffices to consider only the case $\delta_{n}\geq 0$ for all $n=1,2,...,N$. Since by assumption we have $y(h)\geq y_{1}$ then, by induction, we can assume that the assertion holds for the $n$ first terms. The inductive step can be conducted as follows

[TABLE]

The first inequality is a consequence of the inductive assumption, while in the second we use the fact that $\delta_{n+1}\geq 0$ and (8). It follows that $y((n+1)h)\geq y_{n+1}$ what ends the proof. ∎

Next, we prove a comparison theorem that can be thought as a generalization of the Gronwall-Bellman’s Lemma.

Lemma 1**.**

Let $y\in C{[0,1]}$ be a positive function satisfying

[TABLE]

with $\nu,\mu,m>0$. Then

[TABLE]

where the sign $\lessgtr$ denotes either less or greater than.

Proof.

The proof will proceed for the $\geq$-case. The reasoning for the other follows exactly the same route. Let $y$ be as in the assumption, and denote

[TABLE]

with $z(0)=0$. We have then

[TABLE]

and therefore

[TABLE]

After division by $z^{\frac{1}{m+1}}$ we can integrate to obtain

[TABLE]

whence

[TABLE]

Using the above estimate in (13) implies

[TABLE]

This ends the proof. ∎

In a similar vein we can prove the following asymptotic result.

Lemma 2**.**

Let $y\in C{[0,1]}$ be a positive solution of (2) with

[TABLE]

Then

[TABLE]

Proof.

We start by deriving an representation of the solution to (2). By the Mean Value Theorem we have

[TABLE]

because $y$ is positive while $0\leq\tau(x)\leq x$. Similarly as before we substitute

[TABLE]

and therefore

[TABLE]

Solving this differential equation yields

[TABLE]

Finally,

[TABLE]

If now $K(x,t)\sim Cx^{\mu}$ as $x\rightarrow 0^{+}$ with $t\leq x$ we have

[TABLE]

The proof is complete. ∎

The second lemma is a discrete analogue of the above.

Lemma 3**.**

Let $\left\{y_{n}\right\}$, $n=1,2,...,N$ be a sequence of positive numbers satisfying the recurrence (8) with (3) and (5). Assume that the initial value is chosen to satisfy

[TABLE]

then we have

[TABLE]

where $\epsilon(h)$ is the maximal local consistency error (6) of the quadrature (5) applied to the function $t^{\frac{1}{m}}$.

Proof.

Due to the assumption we can proceed to the inductive step. To this end, let $y_{i}$ satisfy (27) for every $i=1,2,...n$. Observe that from (8) with (3) we have

[TABLE]

Because the sum above is a discretization of the integral (5) we can further write

[TABLE]

where $\epsilon_{n+1}$ is the $n+1$th local consistency error for the quadrature. Now, in order to complete the induction we have to show that

[TABLE]

Since by definition, $\epsilon_{n+1}\leq\epsilon$ and obviously $(n+1)h\geq h$ we further have

[TABLE]

Now, due to convexity we have the inequality $(1-x)^{m}\leq 1-x$ for $m\geq 1$ which grants the validity of the above estimate. This ends the proof. ∎

Finally, we prove a variant of the discreet generalization of the Gronwall-Bellman’s lemma.

Lemma 4**.**

Let $\left\{e_{n}\right\}$, $n=1,2,...$ be a sequence of positive numbers satisfying

[TABLE]

for $A>0$. If we define

[TABLE]

then

[TABLE]

Proof.

First, notice that $-1<A\zeta(1-A)<0$ for $0<A<1$ [39]. Therefore, the above bounds are well-defined. To prove them we proceed by induction. Assume (33) and notice that obviously we have $e_{1}\leq M$.

Next, we make the inductive step and assume that (34) holds for all terms up to a fixed $n$. We will show that this assertion is also true for the $(n+1)$th term. First, assume that $A\geq 1$ and use the inductive assumption

[TABLE]

Now, we need an estimate for the above sum. Since the function $x^{A-1}$ is increasing, we can bound the series by an integral

[TABLE]

therefore

[TABLE]

since $B\leq M$. On the other hand, for the case $0<A<1$ we have

[TABLE]

The appropriate bound for the sum comes from the Riemann-Siegel formula [39], which gives

[TABLE]

whence

[TABLE]

where the last inequality follows from the fact that $B\leq-MA\zeta(1-A)$. This concludes the reasoning. ∎

As we noted, the above two lemmas are generalizations of the Gronwall-Bellman’s results. Similar theorems arise in many situations when considering integral equations and numerical methods used in solving them. There is a large collection of papers about various generalizations of the aforementioned inequalities from which the monograph [1] gives a thorough exposition of the continuous case. Many other results concerning the discreet version can be found for instance in [43, 38, 15, 25].

We can now proceed to the main theorem. We claim that the numerical method (8) is locally convergent to the unique positive solution of (2). We difficulty of the proof lies in the non-Lipschitz behaviour of the nonlinearity in (2). The usual method do not work since we cannot make estimate which will invoke the standard Gronwall-Bellman’s lemma. To overcome this difficulty we use a comparison result which helps us to bound the solution from below with a known function for which the numerical method is convergent.

Theorem 1**.**

Let $y=y(x)$ be the nontrivial solution of (2) with (3). Moreover, by $e_{n}:=y(nh)-y_{n}$ denote the error of the numerical scheme (5)-(8) and assume that $\epsilon(h)h^{-\frac{m+1}{m}}<1$ as $h\rightarrow 0^{+}$, where $\epsilon(h)$ is defined in Lemma 3. Then, if $y_{1}$ satisfies (26) with $\left|e_{1}\right|=O(h^{p})$ for some $p>0$, we have

[TABLE]

where $E(h):=1-\frac{m+1}{m}\epsilon(h)h^{-\frac{m+1}{m}}=O(1)$ as $h\rightarrow 0^{+}$.

Proof.

First, we can use the Langrange’s mean-value theorem

[TABLE]

where $\xi$ is between $y(nh)$ and $y_{n}$. Now, we use the above estimate with (2) and (5) to obtain

[TABLE]

Now, we have to consider two cases depending on the relative size of $y(nh)$ and $y_{n}$. If $y_{n}\leq y(nh)$, then $\xi\geq y_{n}$ and from Lemma 3 it follows that

[TABLE]

Eventually, defining

[TABLE]

we can write

[TABLE]

On the other hand, if $y(nh)\leq y_{n}$, then $\xi\geq y(nh)$ and

[TABLE]

where in the first inequality we have used Lemma 1. Hence

[TABLE]

since $E(h)<1$ due to the assumption. It follows that regardless of the relative size of $y(nh)$ and $y_{n}$, the error term satisfies the above recurrence inequality. Now, we can invoke the Lemma 4 to conclude that since $E(h)=O(1)$ as $h\rightarrow 0^{+}$ we have

[TABLE]

where $A:=\frac{1}{m}\frac{WD}{CE(h)}$. Now, when $h\rightarrow 0^{+}$ with $nh\rightarrow\text{const.}$ we obtain the assertion. The proof is complete. ∎

The analysis of the above proof can lead to a stronger convergence result for a specific family of quadratures.

Corollary 1**.**

Let $y=y(x)$ be a solution of (2) with (3). Moreover, choose a quadrature with $\delta_{n}\leq 0$ for each $n=1,2,...,N$. If the starting value $y_{1}$ has $|e_{1}|=O(h^{p})$ for some $p>0$, then

[TABLE]

Proof.

If $\delta_{n}\leq 0$, then by Proposition 1 we have $y(nh)\leq y_{n}$. The proof follows the same steps as in Theorem 1 with the difference that we always have $\xi\geq y(nh)$ in (42). Hence the conclusion can be drawn from (47) without any subsequent estimates. ∎

3 Construction of a quadrature and numerical examples

In this section we will construct a convergent numerical method for solving (8). By Lemma 19, we can identify several problems in numerically solving this integral equation. For example, when $m>1$ the non-smooth behaviour of the solution at $x=0$ can severely slow down the convergence. Recall that for a smooth function, its Newton-Cotes quadrature has an error bounded by the value of a sufficiently high derivative. If the considered function is not smooth enough, the order of the method can be reduced [35, 14]. On the other hand, when $0<m<1$ the Corollary 50 can give very weak estimates on the convergence rate. Having these remarks in mind, in what follows, we will illustrate the above results by constructing a numerical method for the case $m\geq 1$ leaving the other case for the future work.

To maximize the theoretical convergence rate we will use a quadrature with $\delta_{n}(h)\leq 0$ and use Corollary 50. Notice that $y(x)\propto x^{\frac{1}{m}}$ so that for $m\geq 1$ the solution is locally concave. As a quadrature we choose the midpoint method which is of second order for sufficiently smooth functions. However, near $x=0$ the method will be only $O(h^{1+\frac{1}{m}})$. In order to see this fact we use the result from [14] stating that for any function with bounded variation the local consistency error is the following

[TABLE]

Now, due to the monotonicity we have $V_{0}^{h}(t^{\frac{1}{m}})=h^{\frac{1}{m}}$, which states that the quadrature has order $1+\frac{1}{m}$. Note that this is in contrast with usual estimate for smooth functions.

Moreover, we can restrict ourselves only to an arbitrary and fixed interval $[0,\xi]$, so that the derivative singularity at the origin would be integrated. For the rest of the interval any other quadrature can be used since the nonlinearity of the equation (2) is of the Lipschitz class.

In order to convince oneself that the midpoint quadrature will yield a greater area under a curve than the exact integral consider an arbitrary interval $[x_{i},x_{i+1}]$. Construct a line tangent to $x_{i+\frac{1}{2}}$ and notice that due to concavity, the area of a trapezium generated by it is greater than the are area under $y$. On the other hand, the trapezium theorem states that its area is equal to the area of a rectangle spanned by its height and the midline. This precisely is the midpoint quadrature. These considerations can be confirmed by the standard error analysis but this reasoning requires existence of a second derivative.

Before we construct the quadrature we have to choose the starting value for the iteration (8). We have several ways of doing this. Keeping in mind that we have to choose $y_{1}\geq y(h)$, we can try to bound the defining integral equation by some simple, not necessarily open, quadrature. The simplest choice is to use the rectangle rule

[TABLE]

because $y$ is increasing and $K$ positive. We can thus define the initial step for the iteration as

[TABLE]

This choice has a disadvantage of producing a large starting error $O(h^{\frac{1}{m}})$, since by Lemma 19 we have

[TABLE]

Therefore we can put $p=\frac{1}{m}$ in the Corollary 50. In order to increase the accuracy we can always use Richardson’s extrapolation and we implicitly assume that it is utilized in our quadrature.

Having made all the preparations we can propose the following iteration scheme to solve (2)

[TABLE]

The first term $y_{1}$ is calculated by the formula (53) and $y_{0}=0$. This quadrature uses intervals of length $2h$ and approximates the integral by the midpoint rule. Therefore, we have to distinguish between the even and odd steps ($k=0,1$ respectively) with regard to whether or not the terminal point is a midpoint of the interval. The first term on the right-hand side of the above formula is nonzero only in the odd term and is a result of a trapezoid rule applied to the interval $[0,h]$. We will illustrate the above results with some numerical examples.

Example 1. As a first example we consider the following equation

[TABLE]

which can be readily solved with

[TABLE]

Here, $K(x,t)\equiv C=D=1$ and hence, thanks to (55) we have $W=2$. From (41) we immediately deduce that the scheme is convergent with order at least equal to $1+\frac{1}{m}-\frac{2}{m}=1=\frac{1}{m}$. Recall, that since the solution is proportional to $x^{\frac{1}{m}}$, the midpoint method will converge with an order of $1+\frac{1}{m}$.

The results for numerical simulations are given in Tab. 1. We can see that the real order of convergence depends very weakly on $m$ and is close to $1$. Note that we are comparing the exact with numerical solution at a fixed point $x=0.001$. In order to compute the order we use the standard technique of linear regression on a log-log plot of the error versus the iteration step which, in our simulation, changes from $10^{-1}$ to $10^{-5}$. Notice that the integration is exact for $m=1$ which is not surprising since midpoint rule is exact for the linear functions.

Example 2. The next equation we consider is the following equation with a convolution kernel

[TABLE]

which is solved by

[TABLE]

Here, $C=1$, $D=e^{X}$ and $W=2$, hence the order of convergence can be estimated as equal to at least $1+\frac{1}{m}-\frac{2}{m}e^{X}>1-\frac{1}{m}$. Choosing a sufficiently small $X$ assures that the method is convergent (specifically $X<\ln\frac{m+1}{2}$). Results of the simulations are given in Tab. 2. Once again, the order of convergence is very close to $1$.

4 Conclusion

We have shown that an explicit method for solving the Volterra integral equation with a power-type nonlinearity is convergent. The subjects of the future work consist of generalizing our results to the implicit schemes for integral equations with kernels bounded from below by an arbitrary power functions.

Acknowledgements

This research was supported by the National Science Centre, Poland under the project with a signature NCN $2015/17/D/ST1/00625$.

Ł. P. would like to thank Prof. W. Okrasiński for a lot of fruitful talks that shaped the final form and quality of the paper.