On Some Iterative Numerical Methods for Mixed Volterra–Fredholm Integral Equations

Abstract

In this paper, we propose a class of simple numerical methods for approximating solutions of one-dimensional mixed Volterra–Fredholm integral equations of the second kind. These methods are based on fixed point results for the existence and uniqueness of the solution (results which also provide successive iterations of the solution) and suitable cubature formulas for the numerical approximations. We discuss in detail a method using Picard iteration and the two-dimensional composite trapezoidal rule, giving convergence conditions and error estimates. The paper concludes with numerical experiments and a discussion of the methods proposed.

1. Introduction

Integral equations play an important role in the area of applied mathematics, as they arise from a variety of physical, engineering and biological problems.

In this work, we consider mixed Volterra–Fredholm integral equations (MVFIEs) of the form

$$\begin{array}{ccc}\hfill u\left(x\right)& =& \underset{a}{\overset{x}{{\displaystyle \int }}}\underset{a}{\overset{b}{{\displaystyle \int }}}K\left(x,t,s,u\left(s\right)\right)dsdt+f\left(x\right),x\in [a,b],\hfill \end{array}$$

(1)

where $K\in C\left([a,b]\times [a,b]\times [a,b]\times \mathbb{R}\right)$ and $f\in C[a,b]$. More assumptions are made on K and f below. Such equations arise in many applications in areas of physics, fluid dynamics, electrodynamics, and biology. Various formulations of boundary value problems, with Neumann, Dirichlet or mixed boundary conditions are reduced to such integral equations. They also provide mathematical models for the development of an epidemic and numerous other physical and biological problems.

Over time, especially recently, there have been many papers devoted to studying these equations and their properties. Procedures for approximating their solutions numerically have been developed via collocation methods ([1,2,3]), CAS wavelets ([4]), Taylor expansion methods ([5]), block-pulse functions ([6]), linear programming ([7]), spectral methods ([8]), etc. For more considerations on mixed integral equations, see, e.g., [9,10].

The aim of this paper is to present a class of simple, yet reliable numerical methods for approximating the solution of MVFIE’s, using a combination of fixed point results for the existence and uniqueness of the solution and a suitable cubature formula for the numerical approximation of the successive iterates.

The rest of the paper is organized as follows. In Section 2, we analyze the solvability of Equation (1), using results from fixed point theory. In Section 3, we develop a numerical procedure for approximating the solution of Equation (1), when an appropriate cubature formula is used, discussing conditions for convergence and giving error estimates. In particular, we analyze in detail the convergence and give an error bound for the case when the two-dimensional composite trapezoidal rule is used. Section 4 contains numerical experiments that show the applicability of the described method. In Section 5, we draw conclusions and discuss future research ideas.

2. Solvability of the MVFIE in Equation (1)

We study the solvability of Equation (1) using fixed point theory. We recall the main results for fixed points on a Banach space.

Definition 1.

Let $(X,||·|\left|\right)$ be a Banach space. A mapping $F:X\to X$ is called a$\mathit{q}—$contractionif there exists $0<q<1$ such that

$$\left|\right|Fu-Fv\left|\right|\le q\left|\right|u-v\left|\right|,$$

for all $u,v\in X$.

We have the classical result, the contraction principle on a Banach space.

Theorem 1.

Let $(X,||·|\left|\right)$ be a Banach space and $F:X\to X$ a $q—$contraction. Then,

(a) 

equation $u=Fu$ has exactly one solution $u^{*}\in X$;

(b) 

the sequence of successive approximations $u_{n+1}=F{u}_n,n\in \mathbb{N},$ converges to the solution $u^{*}$, for any arbitrary choice of initial point $u_0\in X$; and

(c) 

the error estimate

$$\begin{array}{ccc}\hfill \left|\right|u_n-u^{*}\left|\right|& \le & {\displaystyle \frac{q^n}{1-q}}\left|\right|u_1-u_0\left|\right|\hfill \end{array}$$

holds for every $n\in \mathbb{N}$.

Remark 1.

Theorem 1 still holds if X is replaced by any closed subset $Y\subseteq X$ satisfying $F\left(Y\right)\subseteq Y$. This is important below (see Remark 2), since, to approximate the solution of the MVFIE, we want to apply this fixed point result locally (on such a subset Y), not globally, on the entire space X.

We apply Theorem 1 to the MVFIE in Equation (1). To this end, we define the integral operator $F:C[a,b]\to C[a,b]$ associated with Equation (1), by

$$\begin{array}{ccc}\hfill Fu\left(x\right)& :=& \underset{a}{\overset{x}{{\displaystyle \int }}}\underset{a}{\overset{b}{{\displaystyle \int }}}K\left(x,t,s,u\left(s\right)\right)dsdt+f\left(x\right).\hfill \end{array}$$

(2)

Then, finding a solution of the integral in Equation (1) is equivalent to finding a fixed point for the associated operator F:

$$u=Fu.$$

(3)

To use Theorem 1 (actually, Remark 1), we consider the space $X=C[a,b]$ equipped with the Chebyshev norm

$$\begin{array}{ccc}\hfill \left|\right|u\left|\right|& :=& \underset{x\in [a,b]}{max}\left|u\left(x\right)\right|,u\in X\hfill \end{array}$$

and the ball $B_R:=\{u\in C[a,b]|\left|\right|u-f\left|\right|\le R\}$, for some $R>0$. Then, $(X,||·|\left|\right)$ is a Banach space and $B_R\subseteq X$ is a closed convex subset. Thus, we have:

Theorem 2.

Let $K\in C\left([a,b]\times [a,b]\times [a,b]\times \mathbb{R}\right)$ and $f\in C[a,b]$. Assume that:

(i)

there exists a constant $L>0$ such that

$$\begin{array}{ccc}\hfill \left|K\right(x,t,s,u)-K(x,t,s,v\left)\right|& \le & L|u-v|,\hfill \end{array}$$

(4)

for all $x,t,s\in [a,b]$ and all $u,v\in [R_1-R,R_2+R]$, where ${\displaystyle R_1:=\underset{x\in [a,b]}{min}f\left(x\right),R_2:=\underset{x\in [a,b]}{max}f\left(x\right);}$

(ii)

$$\begin{array}{ccc}\hfill M_K{(b-a)}^2& \le & R,\hfill \end{array}$$

(5)

where ${\displaystyle M_K:=max\left|K(x,t,s,u)\right|}$ over all $x,t,s\in [a,b]$ and all $u,v\in [R_1-R,R_2+R]$

(iii)

$$q:=L{(b-a)}^2<1;$$

(6)

Then,

(a) 

Equation (3) has exactly one solution $u^{*}\in B_R$;

(b) 

the sequence of successive approximations

$$u_{n+1}=F{u}_n,n=0,1,\dots$$

(7)

converges to the solution $u^{*},$ for any arbitrary initial point $u_0\in B_R$; and

(c) 

the error estimate

$$\begin{array}{ccc}\hfill \left|\right|u_n-u^{*}\left|\right|& \le & {\displaystyle \frac{q^n}{1-q}}\left|\right|u_1-u_0\left|\right|\hfill \end{array}$$

(8)

holds for every $n\in \mathbb{N}$.

Proof. 

Let $u\in B_R$, arbitrary. Since $\left|\right|u-f\left|\right|\le R$, it follows that $u\left(x\right)\in [R_1-R,R_2+R]$, for all $x\in [a,b]$. Then, for every fixed $x\in [a,b]$, we have

$$\begin{array}{ccc}\hfill \left|Fu\right(x)-f(x\left)\right|& \le & \underset{a}{\overset{x}{{\displaystyle \int }}}\underset{a}{\overset{b}{{\displaystyle \int }}}|K\left(x,t,s,u\left(s\right)\right)|dsdt\le M_K{(b-a)}^2.\hfill \end{array}$$

Thus, by Equation (5), $Fu\in B_R$.

Next, again for every fixed $x\in [a,b]$, by Equation (4), we have:

$$\begin{array}{ccc}\hfill \left|Fu\right(x)-Fv(x\left)\right|& \le & \underset{a}{\overset{x}{{\displaystyle \int }}}\underset{a}{\overset{b}{{\displaystyle \int }}}|K\left(x,t,s,u\left(s\right)\right)-K\left(x,t,s,v\left(s\right)\right)|dsdt\hfill \\ & \le & L\left|\right|u-v\left|\right|\underset{a}{\overset{x}{{\displaystyle \int }}}\underset{a}{\overset{b}{{\displaystyle \int }}}dsdt\hfill \\ & \le & q\left|\right|x-y\left|\right|.\hfill \end{array}$$

Then,

$$\begin{array}{ccc}\hfill \left|\right|Fu-Fv\left|\right|& \le & q\left|\right|u-v\left|\right|\hfill \end{array}$$

and by Equation (6), all the conclusions follow from Theorem 1 and Remark 1. □

Remark 2.

Let us discuss the Lipschitz condition in Equation (4). One case where it is relatively easy to verify it is if $\left|{\displaystyle \frac{\partial K}{\partial u}}\right|$ is bounded. Then, any bound of the derivative is a Lipschitz constant L for the function. Still, this condition can be quite restrictive, as well, if imposed on the entire space. However, it is much more relaxed when the partial derivative has to be bounded only locally, which is the reason we use a local fixed point result, instead of a global one. We use this in our numerical examples. Of course, there exist functions with unbounded derivatives that still satisfy Equation (4). Even some non-differentiable kernels (for instance, some kernels with weak singularities) can verify a Lipschitz condition (see, e.g., [1]).

For more considerations on fixed point results, see, e.g., [11].

3. Numerical Approximation of the Solution

Now, we have an iterative procedure, given by Equation (7), for approximating the true solution $u^{*}$. However, to use it, we have to approximate the integrals numerically and give approximations of the iterative solution at some mesh points. We consider a cubature formula

$$\underset{a}{\overset{b}{{\displaystyle \int }}}\underset{c}{\overset{d}{{\displaystyle \int }}}\phi (t,s)dsdt={\displaystyle \sum }_{i=0}^{m_1}{\displaystyle \sum }_{j=0}^{m_2}{a}_{ij}\phi (t_i,s_j)+R_{\phi },$$

(9)

with nodes $a=t_0<t_1<\cdots <t_{m_1}=b,c=s_0<s_1<\cdots <s_{m_2}=d$, coefficients $a_{ij}\in \mathbb{R},i=0,1,\dots ,m_1,j=0,1,\dots ,m_2$ and for which the remainder satisfies

$$|R_{\phi }|\le M,$$

(10)

for some $M>0$, with $M\to 0$ as $m_1,m_2\to \infty$.

Since we need to integrate on $[a,x]\times [a,b],x\in [a,b]$, we use Equation (9) the following way: Let $a=x_0<x_1<\cdots <x_m=b$ be a partition of $[a,b]$ and let $u_0={\tilde{u}}_0\equiv f$ be the initial approximation. Then, we use the iteration in Equation (7) and the numerical integration scheme in Equation (9) to approximate $u_n\left(x_k\right)$ by ${\tilde{u}}_n\left(x_k\right)$, for $k=\overline{0,m}$ and $n=0,1,\dots$. We have the following approximations:

$$\begin{array}{ccc}\hfill u_1\left(x_k\right)& =& \underset{a}{\overset{x_k}{{\displaystyle \int }}}\underset{a}{\overset{b}{{\displaystyle \int }}}K\left(x_k,t,s,f\left(s\right)\right)dsdt+f\left(x_k\right)\hfill \\ & =& {\displaystyle \sum }_{i=0}^k{\displaystyle \sum }_{j=0}^m{a}_{ij}K\left(x_k,x_i,x_j,f\left(x_j\right)\right)+R_K+f\left(x_k\right)\hfill \\ & =& {\tilde{u}}_1\left(x_k\right)+{\tilde{R}}_1,\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill {\tilde{u}}_1\left(x_k\right)& =& {\displaystyle \sum }_{i=0}^k{\displaystyle \sum }_{j=0}^m{a}_{ij}K\left(x_k,x_i,x_j,f\left(x_j\right)\right)+f\left(x_k\right).\hfill \end{array}$$

Then, denoting the maximum error at the nodes by

$$\begin{array}{ccc}\hfill e(u_1,{\tilde{u}}_1)& :=& \underset{x_k\in [a,b]}{max}|u_1\left(x_k\right)-{\tilde{u}}_1\left(x_k\right)|,\hfill \end{array}$$

we have, by Equation (10),

$$e(u_1,{\tilde{u}}_1)\le \left|{\tilde{R}}_1\right|\le M.$$

(11)

We proceed further in a similar way to get:

$$\begin{array}{ccc}\hfill u_2\left(x_k\right)& =& \underset{a}{\overset{x_k}{{\displaystyle \int }}}\underset{a}{\overset{b}{{\displaystyle \int }}}K\left(x_k,t,s,u_1\left(s\right)\right)dsdt+f\left(x_k\right)\hfill \\ \hfill & =& {\displaystyle \sum }_{i=0}^k{\displaystyle \sum }_{j=0}^m{a}_{ij}K\left(x_k,x_i,x_j,u_1\left(x_j\right)\right)+R_K+f\left(x_k\right)\hfill \\ & =& {\displaystyle \sum }_{i=0}^k{\displaystyle \sum }_{j=0}^m{a}_{ij}K\left(x_k,x_i,x_j,{\tilde{u}}_1\left(x_j\right)+{\tilde{R}}_1\right)+R_K+f\left(x_k\right)\hfill \\ \hfill & =& {\displaystyle \sum }_{i=0}^k{\displaystyle \sum }_{j=0}^m{a}_{ij}K\left(x_k,x_i,x_j,{\tilde{u}}_1\left(x_j\right)\right)+f\left(x_k\right)+{\tilde{R}}_2\hfill \\ \hfill & =& {\tilde{u}}_2\left(x_k\right)+{\tilde{R}}_2,\hfill \end{array}$$

(12)

where

$$\begin{array}{ccc}\hfill {\tilde{u}}_2\left(x_k\right)& =& {\displaystyle \sum }_{i=0}^k{\displaystyle \sum }_{j=0}^m{a}_{ij}K\left(x_k,x_i,x_j,{\tilde{u}}_1\left(x_j\right)\right)+f\left(x_k\right).\hfill \end{array}$$

Thus, the values ${\tilde{u}}_2\left(x_k\right)$ can be computed from the previous iteration.

Now, let us estimate the error. To this end, denote by ${\displaystyle \gamma :=L{\displaystyle \sum }_{i=0}^m{\displaystyle \sum }_{j=0}^m|a_{ij}|}$. Then, by Equation (11), we have

$$\begin{array}{ccc}\hfill e(u_2,{\tilde{u}}_2)& \le & |{\tilde{R}}_2|\hfill \\ \hfill & \le & {\displaystyle {\displaystyle \sum }_{i=0}^k{\displaystyle \sum }_{j=0}^m|a_{ij}\left|L\right|{\tilde{R}}_1|+|R_K|}\hfill \\ & \le & LM{\displaystyle \sum }_{i=0}^m{\displaystyle \sum }_{j=0}^m|a_{ij}|+M\hfill \\ \hfill & =& M(1+\gamma ).\hfill \end{array}$$

(13)

Similarly, denoting by

$$\begin{array}{ccc}\hfill {\tilde{u}}_n\left(x_k\right)& =& {\displaystyle \sum }_{i=0}^k{\displaystyle \sum }_{j=0}^m{a}_{ij}K\left(x_k,x_i,x_j,{\tilde{u}}_{n-1}\left(x_j\right)\right)+f\left(x_k\right),k=0,1,\dots ,m,\hfill \end{array}$$

(14)

inductively, we get

$$\begin{array}{ccc}\hfill e(u_n,{\tilde{u}}_n)& \le & |{\tilde{R}}_n|\hfill \\ & \le & \gamma |{\tilde{R}}_{n-1}|+M\hfill \\ \hfill & \le & M\left(1+\gamma +\cdots +{\gamma }^{n-1}\right).\hfill \end{array}$$

(15)

Now, we can give an error estimate for our approximations.

Theorem 3.

Assume the conditions of Theorem 2 hold. Further, assume that the coefficients in the cubature in Equation (9) satisfy the condition

$${\displaystyle \gamma =L{\displaystyle \sum }_{i=0}^m{\displaystyle \sum }_{j=0}^m|a_{ij}|<1.}$$

(16)

Then, for the true solution $u^{*}$ of Equation (3) and the approximations ${\tilde{u}}_n$ given by Equation (14), the error estimate

$$\begin{array}{ccc}\hfill e({\tilde{u}}_n,u^{*})& \le & {\displaystyle \frac{q^n}{1-q}}\left|\right|u_1-u_0\left|\right|+{\displaystyle \frac{M}{1-\gamma }}\hfill \end{array}$$

(17)

holds for every $n\in \mathbb{N}$.

Proof. 

For every $k=0,\dots ,m$ and $n=0,1,\dots$, we have

$$\begin{array}{ccc}\hfill |u^{*}\left(x_k\right)-{\tilde{u}}_n\left(x_k\right)|& \le & |u^{*}\left(x_k\right)-u_n\left(x_k\right)|+|u_n\left(x_k\right)-{\tilde{u}}_n\left(x_k\right)|.\hfill \end{array}$$

The assertion then follows from Equations (15) and (16) and Theorem 2. □

Using the Trapezoidal Rule

Thus, to approximate the iterates $u_n\left(x_k\right)$, we can choose any numerical integration formula that satisfies condition in Equation (16). Next, we propose a simple such formula, the two-dimensional composite trapezoidal rule:

$$\begin{array}{ccc}\hfill \underset{a}{\overset{b}{{\displaystyle \int }}}\underset{c}{\overset{d}{{\displaystyle \int }}}\phi (t,s)dsdt& =& {\displaystyle \frac{(b-a)(d-c)}{4m_1{m}_2}}[\phi (a,c)+\phi (b,c)+\phi (a,d)+\phi (b,d)\hfill \\ & +& 2{\displaystyle \sum }_{i=1}^{m_1-1}\left(\phi (t_i,c)+\phi (t_i,d)\right)\hfill \\ \hfill & +& 2{\displaystyle \sum }_{j=1}^{m_2-1}\left(\phi (a,s_j)+\phi (b,s_j)\right)+4{\displaystyle \sum }_{i=1}^{m_1-1}{\displaystyle \sum }_{j=1}^{m_2-1}\phi \left(t_i,s_j)\right)]+R_{\phi },\hfill \end{array}$$

(18)

with equidistant nodes $t_i=a+{\displaystyle \frac{b-a}{m_1}}i,s_j=c+{\displaystyle \frac{d-c}{m_2}}j,i=\overline{0,m_1},j=\overline{0,m_2}$. The remainder is given by

$$\begin{array}{c}\hfill R_{\phi }=-[{\displaystyle \frac{{(b-a)}^3(d-c)}{12m_1^2{m}_2}}{\phi }^{(2,0)}(\xi ,{\eta }_1)+{\displaystyle \frac{(b-a){(d-c)}^3}{12m_1{m}_2^2}}{\phi }^{(0,2)}({\xi }_1,\eta )\\ \hfill +{\displaystyle \frac{{(b-a)}^3{(d-c)}^3}{144m_1^2{m}_2^2}}{\phi }^{(2,2)}(\xi ,\eta )],\xi ,{\xi }_1\in (a,b),\eta ,{\eta }_1\in (c,d),\end{array}$$

(19)

where the notation ${\phi }^{(\alpha ,\beta )}(t,s)={\displaystyle \frac{{\partial }^{\alpha +\beta }\phi }{\partial t^{\alpha }\partial s^{\beta }}}(t,s)$ is used.

For a fixed m, we consider the nodes $x_k=a+{\displaystyle \frac{b-a}{m}}k,k=\overline{0,m}$. Then, we have

$$\begin{array}{ccc}\hfill \underset{a}{\overset{x_k}{{\displaystyle \int }}}\underset{a}{\overset{b}{{\displaystyle \int }}}K\left(x_k,t,s,u_n\left(s\right)\right)dsdt& =& {\displaystyle \frac{(x_k-a)(b-a)}{4km}}[K\left(x_k,a,a,u_n\left(a\right)\right)\hfill \\ & \hfill & +K\left(x_k,x_k,a,u_n\left(a\right)\right)+K\left(x_k,a,b,u_n\left(b\right)\right)+K\left(x_k,x_k,b,u_n\left(b\right)\right)\hfill \\ & & +2{\displaystyle \sum }_{i=1}^{k-1}\left(K\left(x_k,x_i,a,u_n\left(a\right)\right)+K\left(x_k,x_i,b,u_n\left(b\right)\right)\right)\hfill \\ & \hfill & +2{\displaystyle \sum }_{j=1}^{m-1}\left(K\left(x_k,a,x_j,u_n\left(x_j\right)\right)+K\left(x_k,x_k,x_j,u_n\left(x_j\right)\right)\right)\hfill \\ \hfill & & +4{\displaystyle \sum }_{i=1}^{k-1}{\displaystyle \sum }_{j=1}^{m-1}K\left(x_k,x_i,x_j,u_n\left(x_j\right)\right)]+R_K,\hfill \end{array}$$

(20)

for each $k=0,1,\dots ,m$.

Now, ${\displaystyle \frac{x_k-a}{k}}={\displaystyle \frac{b-a}{m}}$, thus, in this case, $\gamma \le L{(b-a)}^2=q$, which is already assumed to be strictly less than 1, by Equation (6).

To find M from Equation (10), we have to find bounds for $K^{(2,0)}\left(x,t,s,u_n\left(s\right)\right)$, $K^{(0,2)}\left(x,t,s,u_n\left(s\right)\right)$ and $K^{(2,2)}\left(x,t,s,u_n\left(s\right)\right)$, as functions of t and s. We have

$$\begin{array}{ccc}\hfill K^{(2,0)}\left(x,t,s,u_n\left(s\right)\right)& =& {\displaystyle \frac{{\partial }^{2}K}{\partial t^2}}\left(x,t,s,u_n\left(s\right)\right),\hfill \\ \hfill K^{(0,2)}\left(x,t,s,u_n\left(s\right)\right)& =& {\displaystyle \frac{{\partial }^{2}K}{\partial s^2}}\left(x,t,s,u_n\left(s\right)\right)+2{\displaystyle \frac{{\partial }^{2}K}{\partial s\partial u}}\left(x,t,s,u_n\left(s\right)\right)u_n^{\prime }\left(s\right)\hfill \\ & +& {\displaystyle \frac{{\partial }^{2}K}{\partial u^2}}\left(x,t,s,u_n\left(s\right)\right){\left(u_n^{\prime }\left(s\right)\right)}^2+{\displaystyle \frac{\partial K}{\partial u}}\left(x,t,s,u_n\left(s\right)\right)u_n^{″}\left(s\right),\hfill \\ \hfill K^{(2,2)}\left(x,t,s,u_n\left(s\right)\right)& =& {\displaystyle \frac{{\partial }^{4}K}{\partial t^2\partial s^2}}\left(x,t,s,u_n\left(s\right)\right)+2{\displaystyle \frac{{\partial }^{4}K}{\partial t^2\partial s\partial u}}\left(x,t,s,u_n\left(s\right)\right)u_n^{\prime }\left(s\right)\hfill \\ & +& {\displaystyle \frac{{\partial }^{4}K}{\partial t^2\partial u^2}}\left(x,t,s,u_n\left(s\right)\right){\left(u_n^{\prime }\left(s\right)\right)}^2+{\displaystyle \frac{{\partial }^{3}K}{\partial t^2\partial u}}\left(x,t,s,u_n\left(s\right)\right)u_n^{″}\left(s\right)\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill u_n\left(x\right)& =& \underset{a}{\overset{x}{{\displaystyle \int }}}\underset{a}{\overset{b}{{\displaystyle \int }}}K\left(x,t,s,u_{n-1}\left(s\right)\right)dsdt+f\left(x\right),\hfill \\ \hfill u_n^{\prime }\left(x\right)& =& \underset{a}{\overset{b}{{\displaystyle \int }}}K\left(x,x,s,u_{n-1}\left(s\right)\right)ds+f^{\prime }\left(x\right),\hfill \\ \hfill u_n^{″}\left(x\right)& =& \underset{a}{\overset{b}{{\displaystyle \int }}}\left({\displaystyle \frac{\partial K}{\partial x}}\left(x,x,s,u_{n-1}\left(s\right)\right)+{\displaystyle \frac{\partial K}{\partial t}}\left(x,x,s,u_{n-1}\left(s\right)\right)\right)ds+f^{″}\left(x\right).\hfill \end{array}$$

From the calculations above, it is clear that, if K is a $C^4$ function with bounded fourth-order partial derivatives and f is a $C^2$ function with bounded second-order derivatives, then there exists a constant $\mathcal{C}>0$, independent of n or m, such that

$$\begin{array}{ccc}\hfill |R_K|& \le & {\displaystyle \frac{{(b-a)}^4}{12m^2}}\mathcal{C}=:M.\hfill \end{array}$$

(21)

Thus, when the trapezoidal rule is used, we have the following approximation result:

Theorem 4.

Assume the conditions of Theorem 2 hold with $K\in C^4\left([a,b]\times [a,b]\times [a,b]\times \mathbb{R}\right)$ and $f\in C^2[a,b]$. Then, for the true solution $u^{*}$ of Equation (3) and the approximations ${\tilde{u}}_n$ given by Equation (14), for all $n=1,2,\dots$ and any $m\in {\mathbb{N}}^{*}$, we have the error estimate

$$\begin{array}{ccc}\hfill e({\tilde{u}}_n,u^{*})& \le & {\displaystyle \frac{q^n}{1-q}}\left|\right|u_1-u_0\left|\right|+{\displaystyle \frac{{(b-a)}^4}{12m^2}}{\displaystyle \frac{\mathcal{C}}{1-\gamma }},\hfill \end{array}$$

(22)

with $\mathcal{C}$ given in Equation (21).

4. Numerical Experiments

Example 1.

Let us start with the nonlinear MVFIE

$$\begin{array}{ccc}\hfill u\left(x\right)& =& {\displaystyle \frac{1}{2}}\underset{0}{\overset{x}{{\displaystyle \int }}}\underset{0}{\overset{\pi /4}{{\displaystyle \int }}}sint\left(1-u^2\left(s\right)\right)sdsdt+{\displaystyle \frac{\pi -2}{16}}\left[1+\left({\displaystyle \frac{16\sqrt{2}}{\pi -2}}-1\right)cosx\right],x\in \left[0,\frac{\pi }{4}\right],\hfill \end{array}$$

(23)

with exact solution $u^{*}\left(x\right)=\sqrt{2}cosx$.

We have $K(x,t,s,u)={\displaystyle \frac{s}{2}}sint\left(1-u^2\right)$ and $f\left(x\right)={\displaystyle \frac{\pi -2}{16}}\left[1+\left({\displaystyle \frac{16\sqrt{2}}{\pi -2}}-1\right)cosx\right]$. Then, on $\left[0,\frac{\pi }{4}\right]$, f is a decreasing function, thus $R_1=f(\pi /4)\approx 1.0209$ and $R_2=f\left(0\right)=\sqrt{2}$. In addition, on $\left[0,\frac{\pi }{4}\right]$, $0\le sint\le \frac{1}{\sqrt{2}}$. If we choose $R=1$, then

$${\displaystyle \underset{[R_1-R,R_2+R]}{max}|1-u^2|=|1-{(R_2+R)}^2|=2(1+\sqrt{2}).}$$

Thus, for $x,t,s\in \left[0,\frac{\pi }{4}\right]$ and $u\in [R_1-R,R_2+R]$, we have

$$\begin{array}{ccc}\hfill M_K{(b-a)}^2& \le & {\displaystyle \frac{\pi }{8}}{\displaystyle \frac{1}{\sqrt{2}}}2(1+\sqrt{2}){\left({\displaystyle \frac{\pi }{4}}\right)}^2={\left(\frac{\pi }{4}\right)}^3\left(1+{\displaystyle \frac{1}{\sqrt{2}}}\right)\approx 0.827<R.\hfill \end{array}$$

Now, to find the Lipschitz constant L, we compute ${\displaystyle \frac{\partial K}{\partial u}}=-susint$. Then, for $x,t,s\in \left[0,\frac{\pi }{4}\right]$ and $u\in [R_1-R,R_2+R]$, we get

$$\begin{array}{ccc}\hfill L{(b-a)}^2& \le & {\displaystyle \frac{\pi }{4}}{\displaystyle \frac{1}{\sqrt{2}}}(R_2+R){\left({\displaystyle \frac{\pi }{4}}\right)}^2={\left(\frac{\pi }{4}\right)}^3\left(1+{\displaystyle \frac{1}{\sqrt{2}}}\right)\approx 0.827<1.\hfill \end{array}$$

Thus, the conditions of Theorem 4 are satisfied.

We use the trapezoidal rule with $m=12$ and $m=24$, with the corresponding nodes $x_k={\displaystyle \frac{\pi }{4m}}k,k=\overline{0,m}$.

Table 1. Errors for Example 1.

	m	12	24
n
1		$1.061540\times {10}^{-1}$	$1.923056\times {10}^{-2}$
5		$9.798965\times {10}^{-4}$	$2.386057\times {10}^{-4}$
10		$1.318075\times {10}^{-4}$	$3.071171\times {10}^{-5}$

contains the errors $e(u^{*},{\tilde{u}}_n)$, with initial approximation $u_0\left(x\right)=f\left(x\right)={\displaystyle \frac{\pi -2}{16}}\left[1+\left({\displaystyle \frac{16\sqrt{2}}{\pi -2}}-1\right)cosx\right]$.

Example 2.

Next, consider the nonlinear MVFIE

$$\begin{array}{ccc}\hfill u\left(x\right)& =& {\displaystyle \frac{1}{4}}\underset{0}{\overset{x}{{\displaystyle \int }}}\underset{0}{\overset{1}{{\displaystyle \int }}}t{e}^{x-2s}{u}^2\left(s\right)dsdt+{\displaystyle \frac{1}{600}}x(120-x)e^x,x\in [0,1],\hfill \end{array}$$

(24)

whose true solution is $u^{*}\left(x\right)={\displaystyle \frac{1}{5}}x{e}^x.$

In this case, $K(x,t,s,u)={\displaystyle \frac{1}{4}}{e}^{x}t{e}^{-2s}{u}^2$, $f\left(x\right)={\displaystyle \frac{1}{600}}x(120-x)e^x$ and ${\displaystyle \frac{\partial K}{\partial u}}={\displaystyle \frac{1}{2}}{e}^{x}t{e}^{-2s}u$. Now $R_1=0$ and $R_2\approx 0.5391$.

Choose again $R=1$. Then, for $x,t,s\in [0,1]$ and $u\in [R_1-R,R_2+R]$, we have

$$\begin{array}{ccc}\hfill M_K{(b-a)}^2& \le & {\displaystyle \frac{1}{4}}{(R_2+R)}^2\approx 0.5922\le R\mathrm{and}\hfill \\ \hfill L{(b-a)}^2& \le & {\displaystyle \frac{1}{2}}(R_2+R)\approx 0.7695<1,\hfill \end{array}$$

thus Theorem 4 applies.

For $m=12$ and $m=24$, we use the trapezoidal rule with nodes $x_k={\displaystyle \frac{1}{m}}k,k=\overline{0,m}$. In

Table 2. Errors for Example 2.

	m	12	24
n
1		$4.110239\times {10}^{-2}$	$1.260717\times {10}^{-2}$
5		$2.253881\times {10}^{-4}$	$5.973695\times {10}^{-5}$
10		$1.521347\times {10}^{-5}$	$3.273994\times {10}^{-6}$

, we give the errors $e(u^{*},{\tilde{u}}_n)$, with initial approximation $u_0\left(x\right)=f\left(x\right)={\displaystyle \frac{1}{600}}x(120-x)e^x$.

5. Conclusions

In this paper, we develop a class of numerical iterative methods for the solutions of one-dimensional mixed Volterra–Fredholm integral equations of the second kind. These methods use successive iterations of the analytical solution (given by fixed point theory results) and numerical approximations of the integrals involved in the iterative solutions. Many results in fixed point theory can be used, as long as they provide an iterative procedure that converges to the true analytic solution. In addition, we describe the conditions that a quadrature rule has to satisfy in order to guarantee the convergence of the approximations of the solution at some mesh points to the true values at those points. Again, many such numerical schemes can be employed. We have chosen the contraction principle for the first part and the two-dimensional composite trapezoidal rule for the numerical integration of the iterates. These choices were made in order to develop a fairly simple numerical method, simple in the proof of the convergence, the conditions to be met and, especially, in implementation. Most mathematical software have the trapezoidal rule built in, thus they are very handy, very little further coding being required. Another advantage of these types of methods, compared to, e.g., collocation, Galerkin, Nyström methods—or other numerical methods that search the solution in a certain form, substitute it into the equation and then force the equation to be true at some mesh points—is that they do not lead to a system of algebraic equations, which oftentimes is ill-conditioned and, thus, difficult to solve. Although the method presented is fairly easy to use and implement, as the numerical examples show, it gives good approximations even with a relatively small number of iterations and of quadrature nodes.

Similar ideas can be used for other types of mixed integral equations, with more complicated kernel functions, such as kernels with some type of singularity or kernels with modified arguments. In addition, other fixed point successive approximation results can be considered (such as Mann iteration or Krasnoselskii iteration), which, under certain conditions, may converge faster than Picard iteration. To increase the speed of convergence of the method, other, higher-order numerical integration schemes can also be used (such as spline quadratures, see, e.g., [12,13,14]), as long as they satisfy the convergence conditions of Theorem 3.