A Procedure for Constructing the Solution of a Nonlinear Fredholm Integro-Differential Equation of Second Order

Abstract

In this work, a large class of integro-differential equations, arising from the description of heat transfer problems, is considered, particularly the nonlinear equations. We propose a procedure for constructing their solution in a very simple and reliable way in which the only needed tool is the same one employed to solve a linear second-order ordinary differential equation, subject to Robin boundary conditions. Proofs of the convergence, existence, and uniqueness are presented. Some special cases are simulated to illustrate the proposed tools.

1. Introduction

In this work, the solution of a Fredholm integro-differential equation is represented as the limit of a sequence whose elements can be obtained from the minimization of quadratic functionals.

The integro-differential problem to be considered here has, as a particular case, the mathematical description of the heat transfer process in symmetrical sets of fins subjected to thermal radiation heat exchange.

In fact, any heat transfer problem involving nonconvex fins (or symmetrical sets of fins) at high temperature levels needs to consider thermal radiant heat transfer from/to the fin. The amount of reflected or emitted thermal radiation from a fin, directly reaching this same fin (or the same set of fins), is represented by an integral operator in the governing equation, giving rise to a second order integro-differential equation. The problems involving thermal radiation heat transfer are inherently nonlinear.

The procedure to be employed for constructing the exact solution may be used for obtaining approximations, for instance, by means of a finite difference scheme or by means of a finite element approximation (taking advantage of the quadratic functional).

Due to its applications in several areas of physics, mathematics, and engineering, Fredholm integro-differential equations continue to be an area of interest.

In the last two decades, their numerical simulations and mathematical analysis have been found with great frequency in scientific articles. Several procedures for solving integro-differential equations have been used, for instance, the Taylor polynomial approach [1,2], block-pulse functions [2,3], the CAS wavelet operational matrix [4,5,6], the Tau method [7], the Spectral Homotopy Analysis method [8,9], the Legendre collocation method [10], the Chebyshev finite difference method [11], the Decomposition Method [12], the Pade approximant [13], and other procedures [14,15].

The main contribution of this work lies in the extreme simplicity of the proposed procedure. It is presented a simple and reliable way to construct the exact solution for a given class of nonlinear Fredholm integro-differential equations, subject to Robin boundary conditions. The mentioned procedure can be also used to carry out numerical simulations for the considered equations. The required tools are available to any undergraduate engineering student.

2. The Considered Problem

The main subject of this work is the problem represented in Equation (1), which generalizes the mathematical description of nonconvex sets of cylindrical fins [16,17,18,19,20,21], in which there is direct thermal radiant interchange among points (far positioned) of these fins. The main objective is to find the function $u$, the solution of

$$\begin{array}{l}\frac{d^{2}u}{d{x}^2}-A\left\{\widehat{f}\left(u\right)\right\}-B\left\{\widehat{h}\left(u\right)-{\displaystyle {\int }_0^1\widehat{h}\left(u\right)K\left(x,\xi \right)d\xi }\right\}+C=00<x<1\\ \frac{du}{dx}={\gamma }_L\left(u-u_L\right),x=0\\ -\frac{du}{dx}={\gamma }_R\left(u-u_R\right),x=1\end{array},$$

(1)

in which ${\gamma }_L$ and ${\gamma }_R$ are positive constants, while $u_L$ and $u_R$ are nonnegative constants. In addition,

$$\begin{array}{l}f=\widehat{f}\left(u\right)\to \mathrm{nondecreasing}\mathrm{function}\mathrm{of}u,\mathrm{with}\widehat{f}\left(0\right)\le 0\\ h=\widehat{h}\left(u\right)\to \mathrm{nondecreasing}\mathrm{function}\mathrm{of}u,\mathrm{with}\widehat{h}\left(0\right)=0\\ A=\tilde{A}\left(x\right)\ge 0\to \mathrm{known}\mathrm{function}\\ B=\tilde{B}\left(x\right)>0\to \mathrm{known}\mathrm{function}\\ C=\tilde{C}\left(x\right)\ge 0\to \mathrm{known}\mathrm{function}\\ K\left(x,\xi \right)\ge 0,0\le x\le 1,0\le \xi \le 1\\ 0\le {\displaystyle {\int }_0^{1}K\left(x,\xi \right)d\xi }\le \mu <1,0\le x\le 1,\mu \to \mathrm{constant}\end{array},$$

(2)

by means of a sequence whose elements can be easily obtained.

For instance, when problem (1) represents the heat transfer process in a set of two parallels fins, the kernel is given by [22].

$$K\left(x,\xi \right)=\frac{{\left(d/H\right)}^2}{4{\left({\left(x-\xi \right)}^2+{\left(d/H\right)}^2\right)}^{3/2}},$$

(3)

in which $d$ is the distance between the fins, and $H$ is the length of each fin.

When problem (1) represents the heat transfer process in a set of two fins (with an angle $2\theta$), the kernel is given by [22].

$$K\left(x,\xi \right)=\frac{x\xi {\sin }^2\left(2\theta \right)}{4{\left({\left(4x\xi \sin \theta \right)}^2+{\left(x-\xi \right)}^2\right)}^{3/2}}.$$

(4)

When the surfaces are assumed to be black, the function $h=\widehat{h}\left(u\right)$ is given by [23,24,25,26].

$$h=\widehat{h}\left(u\right)=B{\left|T\right|}^{3}T,B=\mathrm{constant}.$$

(5)

For a porous fin, the function $f=\widehat{f}\left(u\right)$ is usually given by [20].

$$f=\widehat{f}\left(u\right)=A\left|T-T_{\infty }\right|\left(T-T_{\infty }\right),A=\mathrm{constant},T_{\infty }=\mathrm{constant}.$$

(6)

For a solid cylindrical fin, we usually have [25,26].

$$f=\widehat{f}\left(u\right)=A\left(T-T_{\infty }\right),A=\mathrm{constant},T_{\infty }=\mathrm{constant}.$$

(7)

The function $C=\tilde{C}\left(x\right)$ plays the role of an external source. Many times, it is assumed to be zero everywhere [27].

3. The Solution $\mathit{u}$ Is Nonnegative

Aiming to prove that the solution of (1), denoted here by $u$, is nonnegative everywhere, let us begin by assuming that $u$ has a minimum within the interval $\left(0,1\right)$. Denoting by $\overline{x}$ the point at which the minimum is reached, we have in a small vicinity [28,29] of this point that

$$-A\left\{\widehat{f}\left(u\right)\right\}-B\left\{\widehat{h}\left(u\right)-{\displaystyle {\int }_0^1\widehat{h}\left(u\right)K\left(x,\xi \right)d\xi }\right\}+C\le 0,\overline{x}-{\delta }_1<x<\overline{x}+{\delta }_2.$$

(8)

Taking into account (2), we may write

$$A\left\{\widehat{f}\left(u\right)\right\}+B\left\{\widehat{h}\left(u\right)-{\displaystyle {\int }_0^1\widehat{h}\left(u\right)K\left(x,\xi \right)d\xi }\right\}\ge 0,\overline{x}-{\delta }_1<x<\overline{x}+{\delta }_2.$$

(9)

Therefore, denoting by $u_{MIN}$ the minimum value of $u$, we have (within the considered vicinity)

$$A\left\{\widehat{f}\left(u_{MIN}\right)\right\}+B\left\{\widehat{h}\left(u_{MIN}\right)\right\}\ge 0\Rightarrow u_{MIN}\ge 0.$$

(10)

On the other hand, let us assume that $u$ does not assume a minimum for $x\in \left(0,1\right)$. In this case, we must have the minimum at $x=0$ or at $x=1$. If the minimum is reached at $x=0$, then the derivative of $u$ is nonnegative at $x=0$; hence,

$${\gamma }_L\left(u_{MIN}-u_L\right)\ge 0\mathrm{at}x=0\Rightarrow u_{MIN}\ge 0,$$

(11)

and if the minimum is reached at $x=1$, then the derivative of $u$ is nonpositive at $x=1$; thus,

$${\gamma }_R\left(u_{MIN}-u_R\right)\ge 0\mathrm{at}x=1\Rightarrow u_{MIN}\ge 0.$$

(12)

Hence, it is ensured that the solution of (1), denoted by $u$, is nonnegative everywhere.

4. An Upper Bound for the Solution $u$

Let us assume that $u$ assumes its maximum at the (interior) point $x=\overline{\overline{x}}$. In a sufficiently small neighborhood of this point, we have

$$-A\left\{\widehat{f}\left(u\right)\right\}-B\left\{\widehat{h}\left(u\right)-{\displaystyle {\int }_0^L\widehat{h}\left(u\right)K\left(x,\xi \right)d\xi }\right\}+C\ge 0,\overline{\overline{x}}-{\delta }_1<x<\overline{\overline{x}}+{\delta }_2.$$

(13)

Hence,

$$\left\{\underset{0\le x\le 1}{\max }C\right\}\ge \left\{\underset{0\le x\le 1}{\min }A\right\}\left\{\widehat{f}\left(u_{MAX}\right)\right\}+\left\{\underset{0\le x\le 1}{\min }B\right\}\left(1-\mu \right)\left\{\widehat{h}\left(u_{MAX}\right)\right\}.$$

(14)

The above inequality consists of an upper bound for $u$, provided it assumes a maximum for $x\in \left(0,1\right)$.

On the other hand, let us assume that $u$ does not reach a maximum for $x\in \left(0,1\right)$. In this case, we must have the maximum at $x=0$ or at $x=1$. If the maximum is reached at $x=0$, then the derivative of $u$ is nonpositive at $x=0$; so,

$${\gamma }_L\left(u_{MAX}-u_L\right)\le 0\mathrm{at}x=0\Rightarrow u_{MAX}\le u_L,$$

(15)

and if the maximum is reached at $x=1$, then the derivative of $u$ is nonnegative at $x=1$; hence,

$${\gamma }_R\left(u_{MAX}-u_R\right)\le 0\mathrm{at}x=1\Rightarrow u_{MAX}\le u_R.$$

(16)

Therefore, we are able to evaluate an upper bound for the solution of (1), denoted by $u$. This value will be the larger one among $u_L$, $u_R$, and $u_{\ast }$, in which $u_{\ast }$ is the (unique) root of the equation below.

$$\left\{\underset{0\le x\le 1}{\max }C\right\}=\left\{\underset{0\le x\le 1}{\min }A\right\}\left\{\widehat{f}\left(u_{\ast }\right)\right\}+\left\{\underset{0\le x\le 1}{\min }B\right\}\left(1-\mu \right)\left\{\widehat{h}\left(u_{\ast }\right)\right\}.$$

(17)

5. Solution Construction

The solution of problem (1), denoted here by $u$, is given by the limit of the sequence $\left[{\mathsf{\Psi }}^0,{\mathsf{\Psi }}^1,{\mathsf{\Psi }}^2,{\mathsf{\Psi }}^3,\dots \right]$, whose elements are given by

$$\begin{array}{l}\frac{d^2{\mathsf{\Psi }}^{i+1}}{d{x}^2}-\alpha {\mathsf{\Psi }}^{i+1}+{\beta }^i=0,0<x<1\\ {\beta }^i=\alpha {\mathsf{\Psi }}^i-A\widehat{f}\left({\mathsf{\Psi }}^i\right)-B\left\{\widehat{h}\left({\mathsf{\Psi }}^i\right)-{\displaystyle {\int }_0^1\widehat{h}\left({\mathsf{\Psi }}^i\right)K\left(x,\xi \right)d\xi }\right\}+C\\ \frac{d{\mathsf{\Psi }}^{i+1}}{dx}={\gamma }_L\left({\mathsf{\Psi }}^{i+1}-u_L\right),x=0\\ -\frac{d{\mathsf{\Psi }}^{i+1}}{dx}={\gamma }_R\left({\mathsf{\Psi }}^{i+1}-u_R\right),x=1\end{array},$$

(18)

in which ${\mathsf{\Psi }}^0\equiv 0$ and $\alpha$ is a (large) positive constant. In other words,

$$u\equiv \underset{i\to \infty }{\lim }{\mathsf{\Psi }}^i.$$

(19)

6. The Behavior of the Sequence and the Constant $\alpha$

In order to show that $\left[{\mathsf{\Psi }}^0,{\mathsf{\Psi }}^1,{\mathsf{\Psi }}^2,{\mathsf{\Psi }}^3,\dots \right]$ is a non-decreasing sequence, the first step is to show that ${\mathsf{\Psi }}^1$ is nonnegative everywhere. For this, let us consider $i=0$ and write

$$\begin{array}{l}\frac{d^2{\mathsf{\Psi }}^1}{d{x}^2}-\alpha {\mathsf{\Psi }}^1+{\beta }^0=0,0<x<1\\ {\beta }^0=\alpha {\mathsf{\Psi }}^0-A\widehat{f}\left({\mathsf{\Psi }}^0\right)-B\left\{\widehat{h}\left({\mathsf{\Psi }}^0\right)-{\displaystyle {\int }_0^1\widehat{h}\left({\mathsf{\Psi }}^0\right)K\left(x,\xi \right)d\xi }\right\}+C=-A\widehat{f}\left(0\right)+\tilde{C}\left(x\right)\ge 0\\ \frac{d{\mathsf{\Psi }}^1}{dx}={\gamma }_L\left({\mathsf{\Psi }}^1-u_L\right),x=0\\ -\frac{d{\mathsf{\Psi }}^1}{dx}={\gamma }_R\left({\mathsf{\Psi }}^1-u_R\right),x=1\end{array}.$$

(20)

Suppose that ${\mathsf{\Psi }}^1$ assumes a minimum for $x\in \left(0,1\right)$. In this case, within a sufficiently small vicinity of the point $x=\overline{x}$, in which ${\mathsf{\Psi }}^1$ reaches its minimum, we must have (denoting the minimum by ${\mathsf{\Psi }}_{MIN}^1$).

$$-\alpha {\mathsf{\Psi }}_{MIN}^1+{\beta }^0\le 0\Rightarrow {\mathsf{\Psi }}_{MIN}^1\ge \frac{{\beta }^0}{\alpha }\ge 0.$$

(21)

On the other hand, if ${\mathsf{\Psi }}^1$ assumes its minimum at $x=0$, we must have

$${\gamma }_L\left({\mathsf{\Psi }}_{MIN}^1-u_L\right)\ge 0\Rightarrow {\mathsf{\Psi }}_{MIN}^1\ge u_L\ge 0,$$

(22)

while, if the minimum is reached at $x=1$, we have

$${\gamma }_R\left({\mathsf{\Psi }}_{MIN}^1-u_R\right)\ge 0\Rightarrow {\mathsf{\Psi }}_{MIN}^1\ge u_R\ge 0.$$

(23)

Therefore, the minimum of ${\mathsf{\Psi }}^1$ is nonnegative. In other words, the function ${\widehat{\mathsf{\Psi }}}^1\left(x\right)$ is nonnegative everywhere. In addition, we have proven that ${\mathsf{\Psi }}^1\ge {\mathsf{\Psi }}^0\equiv 0$.

Now, let us consider (18) for two consecutive values of $i$. The difference yields

$$\begin{array}{l}\frac{d^2\left({\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\right)}{d{x}^2}-\alpha \left({\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\right)+\left({\beta }^i-{\beta }^{i-1}\right)=0,0<x<1\\ {\beta }^i-{\beta }^{i-1}=\alpha \left({\mathsf{\Psi }}^i-{\mathsf{\Psi }}^{i-1}\right)-A\left\{\widehat{f}\left({\mathsf{\Psi }}^i\right)-\widehat{f}\left({\mathsf{\Psi }}^{i-1}\right)\right\}-\\ -B\left\{\widehat{h}\left({\mathsf{\Psi }}^i\right)-\widehat{h}\left({\mathsf{\Psi }}^{i-1}\right)-{\displaystyle {\int }_0^1\left(\widehat{h}\left({\mathsf{\Psi }}^i\right)-\widehat{h}\left({\mathsf{\Psi }}^{i-1}\right)\right)K\left(x,\xi \right)d\xi }\right\}\\ \frac{d\left({\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\right)}{dx}={\gamma }_L\left({\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\right),x=0\\ -\frac{d\left({\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\right)}{dx}={\gamma }_R\left({\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\right),x=1\end{array}.$$

(24)

If the difference ${\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i$ assumes its minimum at the point $\overline{x}\in \left(0,1\right)$, then, in a sufficiently small neighborhood of this point, we must have

$$\alpha \left({\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\right)\ge \left({\beta }^i-{\beta }^{i-1}\right).$$

(25)

If the minimum is reached at $x=0$, we must have

$${\gamma }_L\left({\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\right)\ge 0\Rightarrow {\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\ge 0,$$

(26)

while, if the minimum is reached at $x=1$, we have

$${\gamma }_R\left({\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\right)\ge 0\Rightarrow {\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\ge 0.$$

(27)

In this way, in order to ensure that ${\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\ge 0$ everywhere, we must ensure that ${\beta }^i-{\beta }^{i-1}\ge 0$ for all $x\in \left(0,1\right)$.

This condition is always satisfied when the constant $\alpha$ is chosen in such a way that

$$\begin{array}{l}\alpha \left({\mathsf{\Psi }}^i-{\mathsf{\Psi }}^{i-1}\right)\ge \\ \ge A\left\{\widehat{f}\left({\mathsf{\Psi }}^i\right)-\widehat{f}\left({\mathsf{\Psi }}^{i-1}\right)\right\}+B\left\{\widehat{h}\left({\mathsf{\Psi }}^i\right)-\widehat{h}\left({\mathsf{\Psi }}^{i-1}\right)-{\displaystyle {\int }_0^1\left(\widehat{h}\left({\mathsf{\Psi }}^i\right)-\widehat{h}\left({\mathsf{\Psi }}^{i-1}\right)\right)K\left(x,\xi \right)d\xi }\right\}\end{array}$$

(28)

for all $x\in \left(0,1\right)$.

Hence, if ${\mathsf{\Psi }}^i\ge {\mathsf{\Psi }}^{i-1}$, a sufficient (not necessary) condition for ensuring that ${\mathsf{\Psi }}^{i+1}\ge {\mathsf{\Psi }}^i$ is the following

$$\alpha \ge \frac{A\left\{\widehat{f}\left({\mathsf{\Psi }}^i\right)-\widehat{f}\left({\mathsf{\Psi }}^{i-1}\right)\right\}+B\left\{\widehat{h}\left({\mathsf{\Psi }}^i\right)-\widehat{h}\left({\mathsf{\Psi }}^{i-1}\right)\right\}}{{\mathsf{\Psi }}^i-{\mathsf{\Psi }}^{i-1}},x\in \left(0,1\right),i=1,2,3,\dots$$

(29)

Therefore, since ${\mathsf{\Psi }}^1\ge {\mathsf{\Psi }}^0$, we ensure that ${\mathsf{\Psi }}^2\ge {\mathsf{\Psi }}^1$ provided

$$\alpha \ge \frac{A\left\{\widehat{f}\left({\mathsf{\Psi }}^1\right)-\widehat{f}\left({\mathsf{\Psi }}^0\right)\right\}-B\left\{\widehat{h}\left({\mathsf{\Psi }}^1\right)-\widehat{h}\left({\mathsf{\Psi }}^0\right)\right\}}{{\mathsf{\Psi }}^1-{\mathsf{\Psi }}^0}.$$

(30)

Repeating this procedure we have, for sufficiently large $\alpha$, that ${\mathsf{\Psi }}^{i+1}\ge {\mathsf{\Psi }}^i$.

7. The Solution $\mathit{u}$ as an Upper Bound for ${\mathbf{\Psi }}^{\mathit{i}}$

From (1) and (18), we may write

$$\begin{array}{l}\frac{d^2\left(u-{\mathsf{\Psi }}^{i+1}\right)}{d{x}^2}+\alpha \left({\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\right)-\\ -A\left\{\widehat{f}\left(u\right)-\widehat{f}\left({\mathsf{\Psi }}^i\right)\right\}-B\left\{\widehat{h}\left(u\right)-\widehat{h}\left({\mathsf{\Psi }}^i\right)-{\displaystyle {\int }_0^1\left(\widehat{h}\left(u\right)-\widehat{h}\left({\mathsf{\Psi }}^i\right)\right)K\left(x,\xi \right)d\xi }\right\}=0\\ \frac{d\left(u-{\mathsf{\Psi }}^{i+1}\right)}{dx}={\gamma }_L\left(u-{\mathsf{\Psi }}^{i+1}\right),x=0\\ -\frac{d\left(u-{\mathsf{\Psi }}^{i+1}\right)}{dx}={\gamma }_R\left(u-{\mathsf{\Psi }}^{i+1}\right),x=1\end{array}.$$

(31)

Let us consider that $u-{\mathsf{\Psi }}^{i+1}$ assumes its minimum at the interior point $x=\overline{x}$. So, in a sufficiently small neighborhood of this point, we must have

$$\alpha \left({\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\right)-A\left\{\widehat{f}\left(u\right)-\widehat{f}\left({\mathsf{\Psi }}^i\right)\right\}-B\left\{\widehat{h}\left(u\right)-\widehat{h}\left({\mathsf{\Psi }}^i\right)-{\displaystyle {\int }_0^1\left(\widehat{h}\left(u\right)-\widehat{h}\left({\mathsf{\Psi }}^i\right)\right)K\left(x,\xi \right)d\xi }\right\}\le 0$$

(32)

or, in a more convenient form, assuming that (29) holds,

$$\begin{array}{l}\alpha \left({\mathsf{\Psi }}^{i+1}-{\mathsf{\Psi }}^i\right)-A\left\{\widehat{f}\left({\mathsf{\Psi }}^{i+1}\right)-\widehat{f}\left({\mathsf{\Psi }}^i\right)\right\}-\\ -B\left\{\widehat{h}\left({\mathsf{\Psi }}^{i+1}\right)-\widehat{h}\left({\mathsf{\Psi }}^i\right)-{\displaystyle {\int }_0^1\left(\widehat{h}\left(u\right)-\widehat{h}\left({\mathsf{\Psi }}^i\right)\right)K\left(x,\xi \right)d\xi }\right\}\le \\ \le A\left\{\widehat{f}\left(u\right)-\widehat{f}\left({\mathsf{\Psi }}^{i+1}\right)\right\}+B\left\{\widehat{h}\left(u\right)-\widehat{h}\left({\mathsf{\Psi }}^{i+1}\right)\right\}\end{array}.$$

(33)

Since $u$ is nonnegative everywhere, and ${\mathsf{\Psi }}^0\equiv 0$ (and taking into account (29)—the definition of $\alpha$), we conclude that

$$A\left\{\widehat{f}\left(u\right)-\widehat{f}\left({\mathsf{\Psi }}^1\right)\right\}+B\left\{\widehat{h}\left(u\right)-\widehat{h}\left({\mathsf{\Psi }}^1\right)\right\}\ge 0$$

(34)

in the neighborhood of $x=\overline{x}$. Therefore, in this case, $u\ge {\mathsf{\Psi }}^1$ everywhere.

On the other hand, if the difference $u-{\mathsf{\Psi }}^{i+1}$ assumes its minimum at $x=0$, we must have

$${\gamma }_L\left(u-{\mathsf{\Psi }}^{i+1}\right)\ge 0\Rightarrow u-{\mathsf{\Psi }}^{i+1}\ge 0,\mathrm{at}x=0,$$

(35)

while, if the minimum is reached at $x=1$, we have

$${\gamma }_R\left(u-{\mathsf{\Psi }}^{i+1}\right)\ge 0\Rightarrow u-{\mathsf{\Psi }}^{i+1}\ge 0,\mathrm{at}x=1.$$

(36)

Repeating this procedure, we can conclude that the minimum of $u-{\mathsf{\Psi }}^{i+1}$ is nonnegative. Therefore,

$$u\ge \dots \ge {\mathsf{\Psi }}^{i+1}\ge {\mathsf{\Psi }}^i\ge \dots \ge {\mathsf{\Psi }}^3\ge {\mathsf{\Psi }}^2\ge {\mathsf{\Psi }}^1\ge {\mathsf{\Psi }}^0\equiv 0,\mathrm{for}x\in \left[0,1\right].$$

(37)

In other words, the solution of the original problem represents an upper bound for the sequence $\left[{\mathsf{\Psi }}^0,{\mathsf{\Psi }}^1,{\mathsf{\Psi }}^2,{\mathsf{\Psi }}^3,\dots \right]$. This fact ensures the convergence.

Since the solution $u$ is nonnegative and has a known upper bound, the constant $\alpha$ may be chosen from the following formula (this is not a necessary choice)

$$\alpha =\left\{\underset{0\le x\le 1}{\max }A\right\}\left\{\underset{0\le \theta \le u_{MAX}}{\max }\frac{df}{d\theta }\right\}+\left\{\underset{0\le x\le 1}{\max }B\right\}\left\{\underset{0\le \theta \le u_{MAX}}{\max }\frac{dh}{d\theta }\right\},$$

(38)

provided the derivatives of $f$ and of $h$ are bounded.

It is remarkable that the convergence may be reached for lower values of the constant $\alpha$.

8. Solution Uniqueness

In order to demonstrate that the limit of the sequence $\left[{\mathsf{\Psi }}^0,{\mathsf{\Psi }}^1,{\mathsf{\Psi }}^2,{\mathsf{\Psi }}^3,\dots \right]$ is the unique solution of problem (1), let us assume that $u$ is different from ${\mathsf{\Psi }}^{\infty }$ (the limit).

Since (37) holds, we must have

$$u\ge {\mathsf{\Psi }}^{\infty },0\le x\le 1,$$

(39)

and we only need to show that the maximum of $u-{\mathsf{\Psi }}^{\infty }$ is not positive.

For this, let us assume that $u-{\mathsf{\Psi }}^{\infty }$ assumes its maximum at the point $\overline{\overline{x}}\in \left(0,1\right)$. So, in a sufficiently small vicinity of this point, we must have

$$-A\left\{\widehat{f}\left(u\right)-\widehat{f}\left({\mathsf{\Psi }}^{\infty }\right)\right\}-B\left\{\widehat{h}\left(u\right)-\widehat{h}\left({\mathsf{\Psi }}^{\infty }\right)-{\displaystyle {\int }_0^1\left(\widehat{h}\left(u\right)-\widehat{h}\left({\mathsf{\Psi }}^{\infty }\right)\right)K\left(x,\xi \right)d\xi }\right\}\ge 0.$$

(40)

So, from the definition of the constant $\alpha$, we have, in the considered vicinity,

$$\alpha \left(u-{\mathsf{\Psi }}^{\infty }\right)\le 0\Rightarrow u\le {\mathsf{\Psi }}^{\infty }.$$

(41)

Now, suppose that $u-{\mathsf{\Psi }}^{\infty }$ assumes its maximum at the point $x=0$. In this case, we have

$${\gamma }_L\left(u-{\mathsf{\Psi }}^{\infty }\right)\le 0\Rightarrow u\le {\mathsf{\Psi }}^{\infty },\mathrm{at}x=0.$$

(42)

On the other hand, if $u-{\mathsf{\Psi }}^{\infty }$ assumes its maximum at the point $x=1$, we have

$${\gamma }_R\left(u-{\mathsf{\Psi }}^{\infty }\right)\le 0\Rightarrow u\le {\mathsf{\Psi }}^{\infty },\mathrm{at}x=1.$$

(43)

Consequently, it is demonstrated that the solution is unique and is represented by the limit of the sequence $\left[{\mathsf{\Psi }}^0,{\mathsf{\Psi }}^1,{\mathsf{\Psi }}^2,{\mathsf{\Psi }}^3,\dots \right]$.

9. Variational Formulation

The solution of problem (18), denoted by ${\mathsf{\Psi }}^{i+1}$, is the function that minimizes the (quadratic) functional below

$$I\left[v\right]=\frac{1}{2}{\displaystyle \underset{0}{\overset{1}{\int }}\left\{{\left(\frac{dv}{dx}\right)}^2+\alpha v^2\right\}dx}-{\displaystyle \underset{0}{\overset{1}{\int }}{\beta }^{i}vdx}+\frac{1}{2}{\gamma }_L{\left[{\left(v-u_L\right)}^2\right]}_{x=0}+\frac{1}{2}{\gamma }_R{\left[{\left(v-u_R\right)}^2\right]}_{x=1}.$$

(44)

The first variation of $I\left[v\right]$ is given by [30]

$$\delta I\left[v\right]={\displaystyle \underset{0}{\overset{1}{\int }}\left\{\left(\frac{dv}{dx}\right)\frac{d\delta v}{dx}+\alpha v\delta v\right\}dx}-{\displaystyle \underset{0}{\overset{1}{\int }}{\beta }^i\delta vdx}+{\gamma }_L{\left[\left(v-u_L\right)\delta v\right]}_{x=0}+{\gamma }_R{\left[\left(v-u_R\right)\delta v\right]}_{x=1},$$

(45)

or in a more convenient form by

$$\delta I\left[v\right]={\displaystyle \underset{0}{\overset{1}{\int }}\left\{\frac{d}{dx}\left(\left(\frac{dv}{dx}\right)\delta v\right)-\frac{d^{2}v}{d{x}^2}\delta v+\alpha v\delta v\right\}dx}-{\displaystyle \underset{0}{\overset{1}{\int }}{\beta }^i\delta vdx}+{\gamma }_L{\left[\left(v-u_L\right)\delta v\right]}_{x=0}+{\gamma }_R{\left[\left(v-u_R\right)\delta v\right]}_{x=1}$$

(46)

in which $\delta v$ is any admissible variation.

Considering that

$${\displaystyle \underset{0}{\overset{1}{\int }}\frac{d}{dx}\left(\left(\frac{dv}{dx}\right)\delta v\right)dx}={\left[\left(\frac{dv}{dx}\right)\delta v\right]}_{x=1}-{\left[\left(\frac{dv}{dx}\right)\delta v\right]}_{x=0},$$

(47)

we may write

$$\begin{array}{l}\delta I\left[v\right]={\displaystyle \underset{0}{\overset{1}{\int }}\left\{-\frac{d^{2}v}{d{x}^2}\delta v+\alpha v\delta v\right\}dx}-{\displaystyle \underset{0}{\overset{1}{\int }}{\beta }^i\delta vdx}+{\gamma }_L{\left[\left(v-u_L\right)\delta v\right]}_{x=0}+\\ +{\gamma }_R{\left[\left(v-u_R\right)\delta v\right]}_{x=1}+{\left[\left(\frac{dv}{dx}\right)\delta v\right]}_{x=1}-{\left[\left(\frac{dv}{dx}\right)\delta v\right]}_{x=0}\end{array}.$$

(48)

The extremum of the functional is obtained making $\delta I\left[v\right]=0$. Considering that $\delta v$ is arbitrary, we obtain the Euler–Lagrange equation and the natural boundary conditions as follows

$$\begin{array}{l}-\frac{d^{2}v}{d{x}^2}+\alpha v-{\beta }^i=0,0<x<1\\ {\gamma }_L\left(v-u_L\right)-\left(\frac{dv}{dx}\right)=0,x=0\\ {\gamma }_R\left(v-u_R\right)+\left(\frac{dv}{dx}\right)=0,x=1\end{array}.$$

(49)

Problems (44) and (18) are the same problem, since ${\beta }^i$ is known.

Since the functional $I\left[v\right]$ is strictly convex (see Appendix A), its extremum is unique and corresponds to a minimum [29,30]. Thus, each element ${\mathsf{\Psi }}^{i+1}$ may be obtained from the minimization of $I\left[v\right]$.

The existence of this minimum is ensured, since the functional is coercive [31].

10. A Numerical Approximation

The procedure proposed for constructing the solution of (1) may be used for obtaining numerical approximations too.

To illustrate this fact, let us consider the following (piecewise linear) approximation for the element ${\mathsf{\Psi }}^i$, given by

$${\mathsf{\Psi }}^i=\left({\mathsf{\Psi }}_{j+1}^i-{\mathsf{\Psi }}_j^i\right)\left(\frac{x-x_j}{\Delta x}\right)+{\mathsf{\Psi }}_j^i,x_j\le x\le x_{j+1},x_j=\left(j-1\right)\Delta x,j=1,2,\dots ,N,\Delta x=\frac{1}{N}$$

(50)

in which ${\mathsf{\Psi }}_j^i$ represents the approximation for ${\mathsf{\Psi }}^i$ at the spatial point $x_j$, and $\Delta x=1/N$.

In this case, the functional $I\left[v\right]$ becomes, for each $i$, the following function

$$\begin{array}{l}{F}^{i+1}\left(v_1,v_2,\dots ,v_N,v_{N+1}\right)=\frac{1}{2}{\gamma }_L{\left(v_1-u_L\right)}^2+\frac{1}{2}{\gamma }_R{\left(v_N-u_R\right)}^2+\\ +\frac{1}{2}{\displaystyle \sum _{j=1}^N{\displaystyle \underset{x_j}{\overset{x_{j+1}}{\int }}\left\{{\left(\frac{v_{j+1}-v_j}{\Delta x}\right)}^2+\alpha {\left(\left(v_{j+1}-v_j\right)\left(\frac{x-x_j}{\Delta x}\right)+v_j\right)}^2-2{\beta }^i\left(\left(v_{j+1}-v_j\right)\left(\frac{x-x_j}{\Delta x}\right)+v_j\right)\right\}dx}}\end{array}$$

(51)

in which ${\beta }^i$ is given by (considering the piecewise approximation (50))

$$\begin{array}{l}{\beta }^i=\alpha \left(\left({\mathsf{\Psi }}_{j+1}^i-{\mathsf{\Psi }}_j^i\right)\left(\frac{x-x_j}{\Delta x}\right)+{\mathsf{\Psi }}_j^i\right)-A\widehat{f}\left(\left({\mathsf{\Psi }}_{j+1}^i-{\mathsf{\Psi }}_j^i\right)\left(\frac{x-x_j}{\Delta x}\right)+{\mathsf{\Psi }}_j^i\right)-\\ -B\left\{\widehat{h}\left(\left({\mathsf{\Psi }}_{j+1}^i-{\mathsf{\Psi }}_j^i\right)\left(\frac{x-x_j}{\Delta x}\right)+{\mathsf{\Psi }}_j^i\right)-{\displaystyle {\int }_0^1\widehat{h}\left(\left({\mathsf{\Psi }}_{j+1}^i-{\mathsf{\Psi }}_j^i\right)\left(\frac{\xi -x_j}{\Delta x}\right)+{\mathsf{\Psi }}_j^i\right)K\left(x,\xi \right)d\xi }\right\}+C\end{array}.$$

(52)

The values of $v_j$ that minimize the function $F^{i+1}\left(v_1,v_2,\dots ,v_N,v_{N+1}\right)$ defined in (46) are exactly the values of ${\mathsf{\Psi }}_j^{i+1}$, obtained from the following system (linear),

$$\frac{\partial }{\partial v_j}\left\{F^{i+1}\left(v_1,v_2,\dots ,v_N,v_{N+1}\right)\right\}=0,j=1,2,\dots ,N+1.$$

(53)

11. Two Examples (With Known Exact Solutions)

Let us consider the following problem

$$\begin{array}{l}\frac{d^{2}u}{d{x}^2}-u-\left\{{\left|u\right|}^{3}u-{\displaystyle {\int }_0^1{\left|u\right|}^{3}u\frac{5}{31}d\xi }\right\}+x+{\left(1+x\right)}^4=0,0<x<1\\ \frac{du}{dx}=u,x=0\\ -\frac{du}{dx}=u-3,x=1\end{array},$$

(54)

whose exact solution is given by

$$u=1+x,0\le x\le 1.$$

(55)

Comparing (53) with (1) we have that

$$\begin{array}{l}A=B=1,C=x+{\left(1+x\right)}^4,{\gamma }_L={\gamma }_R=1,u_L=0,u_R=3\\ K(x,\xi )=5/31,\widehat{f}\left(u\right)=u,\widehat{h}\left(u\right)={\left|u\right|}^{3}u\end{array}.$$

(56)

Some elements of the sequence $\left[{\mathsf{\Psi }}^0,{\mathsf{\Psi }}^1,{\mathsf{\Psi }}^2,{\mathsf{\Psi }}^3,\dots \right]$, approximated by (50), are shown in Figure 1, as well as the exact solution.

It is to be noticed that the speed of convergence was strongly affected by the value of $\alpha$. As $\alpha$ increased, the speed of convergence decreased. Nevertheless, $\alpha$ cannot be very small, since this may give rise to nonincreasing sequences, and the convergence may not be achieved.

To illustrate the role of the constant $\alpha$ in an explicit way, let us consider the following problem:

$$\begin{array}{l}\frac{d^{2}u}{d{x}^2}-\frac{1}{2}u-\left\{{\left|u\right|}^{3}u-{\displaystyle {\int }_0^1{\left|u\right|}^{3}u\frac{1}{2}d\xi }\right\}+1=0,0<x<1\\ \frac{du}{dx}=u-1,x=0\\ -\frac{du}{dx}=u-1,x=1\end{array},$$

(57)

whose exact solution is given by

$$u=1,0\le x\le 1.$$

(58)

The elements of the sequence are obtained from

$$\begin{array}{l}\frac{d^2{\mathsf{\Psi }}^{i+1}}{d{x}^2}-\alpha {\mathsf{\Psi }}^{i+1}+{\beta }^i=0,0<x<1\\ {\beta }^i=\alpha {\mathsf{\Psi }}^i-\frac{1}{2}{\mathsf{\Psi }}^i-\left\{{\left|{\mathsf{\Psi }}^i\right|}^3{\mathsf{\Psi }}^i-{\displaystyle {\int }_0^1{\left|{\mathsf{\Psi }}^i\right|}^3{\mathsf{\Psi }}^i\frac{1}{2}d\xi }\right\}+1\\ \frac{d{\mathsf{\Psi }}^{i+1}}{dx}=\left({\mathsf{\Psi }}^{i+1}-1\right),x=0\\ -\frac{d{\mathsf{\Psi }}^{i+1}}{dx}=\left({\mathsf{\Psi }}^{i+1}-1\right),x=1\end{array}.$$

(59)

Let us assume that for each $i$, ${\mathsf{\Psi }}^i$ is a constant for $0\le x\le 1$. This is equivalent to imposing ${\mathsf{\Psi }}_2^i={\mathsf{\Psi }}_1^i$, assuming $N=1$. In this case, the functional becomes the function

$$F^{i+1}\left(v\right)=\frac{1}{2}{\left(v-1\right)}^2+\frac{1}{2}{\left(v-1\right)}^2+\frac{1}{2}{\displaystyle \underset{0}{\overset{1}{\int }}\left\{\alpha {\left(v\right)}^2-2{\beta }^i\left(v\right)\right\}dx}.$$

(60)

So, the constant value for ${\mathsf{\Psi }}^{i+1}$ will be obtained from the minimization of the function $F^{i+1}\left(v\right)$. In other words, it will be the root of

$$\frac{d}{dv}{F}^{i+1}\left(v\right)=\left(v-1\right)+\left(v-1\right)+\alpha v-{\displaystyle \underset{0}{\overset{1}{\int }}{\beta }^{i}dx}.$$

(61)

The root of the above expression will be exactly the value of ${\mathsf{\Psi }}^{i+1}$ (assumed a constant for each $i$), given by

$$\left({\mathsf{\Psi }}^{i+1}-1\right)+\left({\mathsf{\Psi }}^{i+1}-1\right)+\alpha {\mathsf{\Psi }}^{i+1}-{\displaystyle \underset{0}{\overset{1}{\int }}{\beta }^{i}dx}=0\Rightarrow {\mathsf{\Psi }}^{i+1}=\frac{{\displaystyle \underset{0}{\overset{1}{\int }}{\beta }^{i}dx}+2}{2+\alpha },$$

(62)

or considering the expression for ${\beta }^i$, we have

$${\mathsf{\Psi }}^{i+1}=\frac{2+{\beta }^i}{2+\alpha }=\frac{3+\left(\alpha -\frac{1}{2}\right){\mathsf{\Psi }}^i-\frac{1}{2}{\left({\mathsf{\Psi }}^i\right)}^4}{2+\alpha }.$$

(63)

Table 1. The constant ${\mathsf{\Psi }}^i$ obtained for nine different values of $\alpha$.

$\mathit{\alpha }\to$	100	50	20	10	5	2	0.5	0.1	0.01
i = 1	0.0294	0.0577	0.1364	0.2500	0.4286	0.7500	1.200	1.4286	1.4925
i = 2	0.0581	0.1126	0.2572	0.4478	0.7017	0.9917	0.7853	0.1648	−0.1058
i = 3	0.0861	0.1649	0.3643	0.6028	0.8623	1.0010	1.1239	1.3970	1.5183
i = 4	0.1134	0.2146	0.4588	0.7217	0.9434	0.9999	0.8808	0.2556	−0.1995
i = 5	0.1400	0.2620	0.5420	0.8101	0.9785	1.0000	1.0796	1.3789	1.5408
i = 6	0.1660	0.3070	0.6149	0.8734	0.9921	1.0000	0.9283	0.3053	−0.285
i = 7	0.1913	0.3499	0.6781	0.9172	0.9972	1.0000	1.0515	1.3684	1.5604
i = 8	0.2161	0.3906	0.7326	0.9466	0.9990	1.0000	0.9555	0.3332	−0.3625
i = 9	0.2402	0.4293	0.7792	0.9659	0.9996	1.0000	1.0333	1.3622	1.5766
i = 10	0.2637	0.4660	0.8186	0.9784	0.9999	1.0000	0.9720	0.3494	−0.4288
i = 11	0.2866	0.5009	0.8517	0.9864	1.0000	1.0000	1.0215	1.3585	1.5887
i = 12	0.3090	0.5339	0.8794	0.9915	1.0000	1.0000	0.9823	0.3589	−0.4793
i = 13	0.3307	0.5651	0.9022	0.9946	1.0000	1.0000	1.0138	1.3563	1.5963
i = 14	0.3520	0.5947	0.9210	0.9966	1.0000	1.0000	0.9887	0.3646	−0.5116
i = 15	0.3727	0.6226	0.9363	0.9979	1.0000	1.0000	1.0089	1.3549	1.6002
i = 16	0.3929	0.6489	0.9488	0.9987	1.0000	1.0000	0.9928	0.3681	−0.5287
i = 17	0.4125	0.6737	0.9590	0.9992	1.0000	1.0000	1.0057	1.3541	1.6020
i = 18	0.4317	0.6970	0.9671	0.9995	1.0000	1.0000	0.9954	0.3702	−0.5364
i = 19	0.4504	0.7189	0.9737	0.9997	1.0000	1.0000	1.0036	1.3536	1.6027
i = 20	0.4685	0.7395	0.9790	0.9998	1.0000	1.0000	0.9971	0.3715	−0.5395
i = 21	0.4862	0.7587	0.9832	0.9999	1.0000	1.0000	1.0023	1.3533	1.6030
i = 22	0.5035	0.7768	0.9866	0.9999	1.0000	1.0000	0.9981	0.3723	−0.5407
i = 23	0.5202	0.7936	0.9893	1.0000	1.0000	1.0000	1.0015	1.3531	1.6031
i = 24	0.5365	0.8093	0.9915	1.0000	1.0000	1.0000	0.9988	0.3727	−0.5411
i = 25	0.5524	0.8240	0.9932	1.0000	1.0000	1.0000	1.0010	1.3530	1.6031
i = 26	0.5678	0.8376	0.9946	1.0000	1.0000	1.0000	0.9992	0.3730	−0.5413
i = 27	0.5828	0.8503	0.9957	1.0000	1.0000	1.0000	1.0006	1.3529	1.6031
i = 28	0.5973	0.8621	0.9966	1.0000	1.0000	1.0000	0.9995	0.3732	−0.5414
i = 29	0.6115	0.8730	0.9973	1.0000	1.0000	1.0000	1.0004	1.3529	1.6031
i = 30	0.6252	0.8832	0.9978	1.0000	1.0000	1.0000	0.9997	0.3733	−0.5414
i = 31	0.6386	0.8926	0.9983	1.0000	1.0000	1.0000	1.0003	1.3528	1.6031
i = 32	0.6515	0.9012	0.9986	1.0000	1.0000	1.0000	0.9998	0.3734	−0.5414
i = 33	0.6641	0.9093	0.9989	1.0000	1.0000	1.0000	1.0002	1.3528	1.6031
i = 34	0.6762	0.9167	0.9991	1.0000	1.0000	1.0000	0.9999	0.3734	−0.5414
i = 35	0.6881	0.9235	0.9993	1.0000	1.0000	1.0000	1.0001	1.3528	1.6031
i = 36	0.6995	0.9298	0.9995	1.0000	1.0000	1.0000	0.9999	0.3734	−0.5414
i = 37	0.7106	0.9356	0.9996	1.0000	1.0000	1.0000	1.0001	1.3528	1.6031
i = 38	0.7213	0.9409	0.9997	1.0000	1.0000	1.0000	0.9999	0.3735	−0.5414
i = 39	0.7318	0.9459	0.9997	1.0000	1.0000	1.0000	1.0000	1.3528	1.6031
i = 40	0.7418	0.9504	0.9998	1.0000	1.0000	1.0000	1.0000	0.3735	−0.5414
i = 41	0.7516	0.9545	0.9998	1.0000	1.0000	1.0000	1.0000	1.3528	1.6031
i = 42	0.7610	0.9584	0.9999	1.0000	1.0000	1.0000	1.0000	0.3735	−0.5414
i = 43	0.7701	0.9619	0.9999	1.0000	1.0000	1.0000	1.0000	1.3528	1.6031
i = 44	0.7789	0.9651	0.9999	1.0000	1.0000	1.0000	1.0000	0.3735	−0.5414
i = 45	0.7874	0.9680	0.9999	1.0000	1.0000	1.0000	1.0000	1.3528	1.6031
i = 46	0.7957	0.9707	0.9999	1.0000	1.0000	1.0000	1.0000	0.3735	−0.5414
i = 47	0.8036	0.9732	1.0000	1.0000	1.0000	1.0000	1.0000	1.3528	1.6031
i = 48	0.8113	0.9755	1.0000	1.0000	1.0000	1.0000	1.0000	0.3735	−0.5414
i = 49	0.8187	0.9776	1.0000	1.0000	1.0000	1.0000	1.0000	1.3528	1.6031
i = 50	0.8258	0.9795	1.0000	1.0000	1.0000	1.0000	1.0000	0.3735	−0.5414

illustrates the convergence for several values of the constant $\alpha$. Notice that for $\alpha =0.1$ and for $\alpha =0.01$, the convergence was not reached. For $\alpha =2$, we reached convergence, but the sequence was not a nondecreasing one. The speed of convergence, for $\alpha \ge 5$, decreased as $\alpha$ increased.

In other words, each column of Table 1 represents the elements: ${\mathsf{\Psi }}^1$, ..., ${\mathsf{\Psi }}^{50}$, obtained with nine different values of the constant $\alpha$.

12. Conclusions

The procedure proposed here differs from others due to its simplicity. Only basic tools are needed, which are available for undergraduate students. This is the main novelty and contribution of this work.

In addition, since the exact solution is represented by the limit of a sequence, there is no limit of accuracy when employing a numerical approximation.

The proposed procedure may be used for problems involving Dirichlet and some Neumann boundary conditions. This may be performed by employing very large or very small values for ${\gamma }_L$ and ${\gamma }_R$. When ${\gamma }_L$ (and/or ${\gamma }_R$) is very large, we approximate a Dirichlet boundary condition. When ${\gamma }_L$ (and/or ${\gamma }_R$) is very small, we approximate a Neumann boundary condition (in this case representing an insulated edge).

When modeling a heat transfer problem involving high temperature levels, thermal radiant heat transfer plays a meaningful role, since the thermal interaction among far-positioned points becomes more significant as the temperature levels become larger. This fact gives rise to an integral operator, which represents the amount of thermal radiant energy arriving at each point of the body.

Besides the heat transfer phenomena, in which the thermal radiation plays a nonnegligible role, integro-differential equations are found in several other branches, such as optimal control problems [32]. Integro-differential equations involving fractional derivatives also consist of a potential issue to be explored due to their increasing interest and applicability [33,34,35,36,37,38,39].

The heat transfer process in a nonsymmetrical system of fins, as well as in multiphase bodies [40], gives rise to systems of second-order integro-differential equations. It seems logical that, with some adjustments, the procedure proposed here may be extended to such a class of problems.