Galerkin-Bernstein Approximations for the System of Third-Order Nonlinear Boundary Value Problems

Snigdha Dhar, Md. Shafiqul Islam^*

Department of Applied Mathematics, University of Dhaka, Dhaka, Bangladesh.

DOI: 10.4236/jamp.2024.126127                

Abstract

This paper is devoted to find the numerical solutions of one dimensional general nonlinear system of third-order boundary value problems (BVPs) for the pair of functions using Galerkin weighted residual method. We derive mathematical formulations in matrix form, in detail, by exploiting Bernstein polynomials as basis functions. A reasonable accuracy is found when the proposed method is used on few examples. At the end of the study, a comparison is made between the approximate and exact solutions, and also with the solutions of the existing methods. Our results converge monotonically to the exact solutions. In addition, we show that the derived formulations may be applicable by reducing higher order complicated BVP into a lower order system of BVPs, and the performance of the numerical solutions is satisfactory.

Keywords

System of Third-Order BVP, Galerkin Method, Bernstein Polynomials, Nonlinear BVP, Higher-Order BVP

1. Introduction

Ordinary differential systems have received a lot of interest in studies as a result of their frequent occurrence in numerous applications. Second-order ordinary differential systems can simulate a number of natural phenomena. For instance, while studying chemically reacting systems, fluid mechanics, relativistic mechanics, gas dynamics, and nuclear physics, the known Emden-Fowler equations appear. Numerous authors have made an effort to increase accuracy in the literature on numerical analysis. Linear second-order boundary value problems can be effectively solved using the finite difference approach. However, when it comes to deal with nonlinear equations, solving the corresponding boundary value problems can be quite challenging. The Galerkin method [1] has been used to solve the two point BVP by the authors Bhatti and Bracken [2]. Bernoulli polynomials [3] have been used to solve second-order both linear and nonlinear BVPs with Dirichlet, Neumann and Robin boundary conditions. Spline functions [4] [5] have been extensively researched, because piece-wise polynomials [6] can be differentiated, integrated, and approximate to any function with any desired accuracy.

Second order ordinary differential systems are used to describe a variety of problems in biology, engineering, and physics. A lot of works have been done for obtaining numerical solutions for the linear system of second-order boundary value problems, for example, Geng and Cui [7] represented the approximate solutions of system of second-order linear and nonlinear BVP in the form of series in the reproducing kernel space. For solving a nonlinear system of second-order BVP, the variational iteration technique [8] was introduced. A nonlinear system of second order BVPs was solved by Dehghan and Saadatmandi using the sinc-collocation method [9]. The exact and approximate solutions were expressed in the reproducing kernel space by Du and Cui [10]. A novel strategy to solve the nonlinear systems of second order BVPs is provided, depending on the homotopy perturbation method and the reproducing kernel method. The advantages of both of these techniques are combined in the homotopy perturbation-reproducing kernel method [11], which can be utilised to effectively handle systems of nonlinear boundary value problems. A few examples of linear and nonlinear systems were solved using the Galerkin weighted residual technique [6]. In order to solve systems of second-order BVPs, a family of boundary value methods was used in [12], as a block unification approach. In order to solve systems of singular boundary value problems, an optimization approach [13] was presented and then solved via continuous genetic algorithm [14].

Recently, Galerkin finite element method [15] has been used to determine the approximate solutions of a coupled second-order BVPs. In recent years, coupled Lane-Emden [16] [17] boundary value problems have received a lot of attention. The singularity is the primary barrier to solutions, and many authors are working to solve it. In order to efficiently solve the system of Lane-Emden type equations that arise in physics, star structure, and astrophysics, Ala’yed, Saadeh and Qazza [18] proposed a collocation approach based on cubic B-spline functions.

While second-order BVPs are well-studied, third-order BVPs present unique challenges due to their complexity. The study of third-order BVPs is crucial to understanding more complex physical systems and phenomena. For instance, third-order BVPs arise in Nano boundary layer fluid flows over a stretching surface [19]. Third-order BVPs are used at resonance on the half-line [20].

Solving a nonlinear system of third-order boundary value problems is quite difficult. Consequently, a few researchers have attempted to solve numerically. Ezzati and Aqhamohamadi [21] used He’s Homotopy Perturbation method to solve the nonlinear system of third-order BVPs. Block method [22] was used to solve the nonlinear system of third-order boundary value problems. We observe that a little attention has been given for solving third-order system of boundary value problems. Thus, we are motivated to solve the system of third-order boundary value problems by Modified Galerkin weighted residual method with Bernstein polynomials as trial functions. However, we organize this research work as follows.

First of all, we give basic idea on Bernstein polynomials in Section 2. Mathematical formulations are described elaborately, in matrix form, in Section 3. Section 4 is devoted to numerical experiments and the discussion of the obtained results in tabular form and graphically. Then in Section 5, we provide an application that sixth order BVP may be modeled into a system of third order BVPs which can be solved by the proposed technique. Finally, conclusions and references are amended.

2. Bernstein Polynomials

The general form of the Bernstein polynomials of degree n over the interval $\left[a\mathrm{,}b\right]$ is defined by [3] [6] [23]

${\varphi }_{i\mathrm{,}n}\left(x\right)=\left(\begin{array}{c}n\\ i\end{array}\right)\frac{{\left(x-a\right)}^i{\left(b-x\right)}^{n-i}}{{\left(b-a\right)}^n}\mathrm{,<mi>\; </mi>}a\le x\le b\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}i=\mathrm{0,1,2,}\cdots \mathrm{,}n$

It is important to observe that each of these $n+1$ polynomials, with a degree of n, fulfills the following properties:

(i) ${\varphi }_{i,n}\left(x\right)=0$ if $i<0$ or $i>n$

(ii) ${\displaystyle {\sum }_{i\mathrm{<mo>=</mo>0}}^n}{\varphi }_{i\mathrm{,}n}\left(x\right)=1$

(iii) ${\varphi }_{i,n}\left(a\right)={\varphi }_{i,n}\left(b\right)=0$ , $1\le i<n$

For simplification, we denote ${\varphi }_{i\mathrm{,}n}\left(x\right)$ by ${\varphi }_i\left(x\right)$ . The Bernstein polynomials of degree 3, 4 and 5 are given below, respectively

Degree 3: ${\varphi }_0\left(x\right)={\left(1-x\right)}^3$ ${\varphi }_1\left(x\right)=3x{\left(1-x\right)}^2$ ${\varphi }_2\left(x\right)=3x^2\left(1-x\right)$

${\varphi }_3\left(x\right)=x^3$

Degree 4: ${\varphi }_0\left(x\right)={\left(1-x\right)}^4$ ${\varphi }_1\left(x\right)=4x{\left(1-x\right)}^3$ ${\varphi }_2\left(x\right)=6x^2{\left(1-x\right)}^2$

${\varphi }_3\left(x\right)=4x^3\left(1-x\right)$ ${\varphi }_4\left(x\right)=x^4$

Degree 5: ${\varphi }_0\left(x\right)={\left(1-x\right)}^5$ ${\varphi }_1\left(x\right)=5x{\left(1-x\right)}^4$ ${\varphi }_2\left(x\right)=10x^2{\left(1-x\right)}^3$

${\varphi }_3\left(x\right)=10x^3{\left(1-x\right)}^2$ ${\varphi }_4\left(x\right)=5x^4\left(1-x\right)$ ${\varphi }_5\left(x\right)=x^5$

To solve a BVP using the Galerkin method, it is required that each of these polynomials satisfies the homogeneous representation of the essential boundary conditions. Here, ${\varphi }_0\left(x\right)$ and ${\varphi }_n\left(x\right)$ do not satisfy the homogeneous boundary conditions. Therefore, in order to satisfy the homogeneous boundary conditions, we use only ${\varphi }_i\left(x\right)$ for $1\le i\le \left(n-1\right)$ . We use Bernstein polynomials of degree 3, 4 and 5 throughout this paper.

3. Mathematical Formulation of System of Third-Order BVPs

In recent times, the interest in boundary value problems of a system of ordinary differential equations has been sparked among researchers in mathematics, physics, engineering, biology, and other fields.

The general linear system of two third order differential equations in two unknown functions $p\left(x\right)$ and $q\left(x\right)$ of the form below is taking into consideration [6]

$\{\begin{array}{l}{a}_1\left(x\right)p^{‴}+a_2\left(x\right)p^{″}+a_3\left(x\right)p^{\prime }+a_4\left(x\right)p+a_5\left(x\right)q^{‴}+a_6\left(x\right)q^{″}\\ +a_7\left(x\right)q^{\prime }+a_8\left(x\right)q=g_1\left(x\right)\\ b_1\left(x\right)p^{‴}+b_2\left(x\right)p^{″}+b_3\left(x\right)p^{\prime }+b_4\left(x\right)p+b_5\left(x\right)q^{‴}+b_6\left(x\right)q^{″}\\ +b_7\left(x\right)q^{\prime }+b_8\left(x\right)q=g_2\left(x\right)\end{array}$ (1)

where $a_j\left(x\right)$ , $b_j\left(x\right)$ , $g_1\left(x\right)$ , $g_2\left(x\right)$ are given functions, and $a_j\left(x\right)$ , $b_j\left(x\right)$ are continuous for $j=\mathrm{1,2,}\cdots \mathrm{,8}$ .

Let’s have a look into system of third-order linear ordinary boundary value problems in one dimension for the pair of functions $p\left(x\right)$ and $q\left(x\right)$ in the following form [15]

$\{\begin{array}{l}{p}^{‴}+a_1\left(x\right)p^{″}+a_2\left(x\right)p^{\prime }+a_3\left(x\right)p+a_4\left(x\right)q^{″}+a_5\left(x\right)q^{\prime }\\ +a_6\left(x\right)q=f\left(x\right),a\le x\le b\\ q^{‴}+b_1\left(x\right)q^{″}+b_2\left(x\right)q^{\prime }+b_3\left(x\right)q+b_4\left(x\right)p^{″}+b_5\left(x\right)p^{\prime }\\ +b_6\left(x\right)p=g\left(x\right),a\le x\le b.\end{array}$ (2)

Here, both of the equations are of third-order, therefore three boundary conditions are needed for both $p\left(x\right)$ and $q\left(x\right)$ . Let us assume the boundary conditions at the ends, i.e.

$p\left(a\right)={\alpha }_1,\mathrm{<mi>\; </mi>}{p}^{\prime }\left(a\right)={\gamma }_1,\mathrm{<mi>\; </mi>}{p}^{\prime }\left(b\right)={\beta }_1,\mathrm{<mi>\; </mi>}q\left(a\right)={\alpha }_2,\mathrm{<mi>\; </mi>}{q}^{\prime }\left(a\right)={\gamma }_2,\mathrm{<mi>\; </mi>}{q}^{\prime }\left(b\right)={\beta }_2.$ (3)

The terms $a_i\left(x\right)$ and $b_i\left(x\right)$ are continuous functions for $i=\mathrm{1,2,}\cdots \mathrm{,}n$ and, $f\left(x\right)$ and $g\left(x\right)$ are non-homogeneous.

The trial solutions for the two functions $p\left(x\right)$ and $q\left(x\right)$ of system (2) can be written as

$\{\begin{array}{l}\tilde{p}\left(x\right)={\displaystyle \sum _{i=1}^{n-1}}{a}_i{\varphi }_i\left(x\right),n\ge 1\hfill \\ \tilde{q}\left(x\right)={\displaystyle \sum _{i=1}^{n-1}}{b}_i{\varphi }_i\left(x\right),n\ge 1\hfill \end{array}$(4)

where $a_i$ , $b_i$ are unknown parameters and ${\varphi }_i\left(x\right)$ are basis functions which must satisfy the boundary conditions (3).

Using these approximations, we can define two residual functions ${\epsilon }_p\left(x\right)$ and ${\epsilon }_q\left(x\right)$ :

$\{\begin{array}{l}{\epsilon }_p\left(x\right)={\tilde{p}}^{‴}+a_1\left(x\right){\tilde{p}}^{″}+a_2\left(x\right){\tilde{p}}^{\prime }+a_3\left(x\right)\tilde{p}+a_4\left(x\right){\tilde{q}}^{″}+a_5\left(x\right){\tilde{q}}^{\prime }+a_6\left(x\right)\tilde{q}-f\left(x\right)\\ {\epsilon }_q\left(x\right)={\tilde{q}}^{‴}+b_1\left(x\right){\tilde{q}}^{″}+b_2\left(x\right){\tilde{q}}^{\prime }+b_3\left(x\right)\tilde{q}+b_4\left(x\right){\tilde{p}}^{″}+b_5\left(x\right){\tilde{p}}^{\prime }+b_6\left(x\right)\tilde{p}-g\left(x\right).\end{array}$ (5)

Now applying the Galerkin method, we get weighted residual system of equations

$\{\begin{array}{l}{\displaystyle {\int }_a^b}{\epsilon }_p\left(x\right){\varphi }_i\left(x\right)\text{d}x=0\hfill \\ {\displaystyle {\int }_a^b}{\epsilon }_q\left(x\right){\varphi }_i\left(x\right)\text{d}x=0.\hfill \end{array}$ (6)

Substituting ${\epsilon }_p\left(x\right)$ and ${\epsilon }_q\left(x\right)$ , we get

$\{\begin{array}{l}{\displaystyle {\int }_a^b}({\tilde{p}}^{‴}+a_1\left(x\right){\tilde{p}}^{″}+a_2\left(x\right){\tilde{p}}^{\prime }+a_3\left(x\right)\tilde{p}+a_4\left(x\right){\tilde{q}}^{″}\\ +a_5\left(x\right){\tilde{q}}^{\prime }+a_6\left(x\right)\tilde{q}){\varphi }_i\left(x\right)\text{d}x={\displaystyle {\int }_a^b}f\left(x\right){\varphi }_i\left(x\right)\text{d}x\\ {\displaystyle {\int }_a^b}({\tilde{q}}^{‴}+b_1\left(x\right){\tilde{q}}^{″}+b_2\left(x\right){\tilde{q}}^{\prime }+b_3\left(x\right)\tilde{q}+b_4\left(x\right){\tilde{p}}^{″}\\ +b_5\left(x\right){\tilde{p}}^{\prime }+b_6\left(x\right)\tilde{p}){\varphi }_i\left(x\right)\text{d}x={\displaystyle {\int }_a^b}g\left(x\right){\varphi }_i\left(x\right)\text{d}x.\end{array}$(7)

Applying integration by parts in the first term of (7) and setting ${\varphi }_i\left(x\right)=0$ at the boundary $x=a$ and $x=b$ , we obtain

$\begin{array}{c}{\displaystyle {\int }_a^b}\frac{{\text{d}}^3\tilde{p}}{\text{d}{x}^3}{\varphi }_i\left(x\right)\text{d}x={\left[\frac{\text{d}\tilde{p}}{\text{d}x}{\varphi }_i\left(x\right)\right]}_a^b-{\displaystyle {\int }_a^b}\frac{\text{d}\tilde{p}}{\text{d}x}\frac{\text{d}{\varphi }_i}{\text{d}x}\text{d}x\\ =\frac{\text{d}\tilde{p}}{\text{d}x}\left(x=b\right){\varphi }_i\left(b\right)-\frac{\text{d}\tilde{p}}{\text{d}x}\left(x=a\right){\varphi }_i\left(a\right)-{\displaystyle {\int }_a^b}\frac{\text{d}\tilde{p}}{\text{d}x}\frac{\text{d}{\varphi }_i}{\text{d}x}\text{d}x\\ =-{\left[\frac{\text{d}}{\text{d}x}{\varphi }_i\left(x\right)\frac{\text{d}\tilde{p}}{\text{d}x}\right]}_a^b+{\displaystyle {\int }_a^b}\frac{{\text{d}}^2}{\text{d}{x}^2}{\varphi }_i\left(x\right)\frac{\text{d}\tilde{p}}{\text{d}x}\text{d}x\\ =-{\left[\frac{\text{d}}{\text{d}x}{\varphi }_i(x)\frac{\text{d}\tilde{p}}{\text{d}x}\right]}_{x=b}+{\left[\frac{\text{d}}{\text{d}x}{\varphi }_i\left(x\right)\frac{\text{d}\tilde{p}}{\text{d}x}\right]}_{x=a}+{\displaystyle {\int }_a^b}\frac{{\text{d}}^2}{\text{d}{x}^2}{\varphi }_i\left(x\right)\frac{\text{d}\tilde{p}}{\text{d}x}\text{d}x\end{array}$

$\begin{array}{l}{\displaystyle {\int }_a^b}\frac{{\text{d}}^3\tilde{q}}{\text{d}{x}^3}{\varphi }_i\left(x\right)\text{d}x=-{\left[\frac{\text{d}}{\text{d}x}{\varphi }_i\left(x\right)\frac{\text{d}\tilde{q}}{\text{d}x}\right]}_{x=b}+{\left[\frac{\text{d}}{\text{d}x}{\varphi }_i\left(x\right)\frac{\text{d}\tilde{q}}{\text{d}x}\right]}_{x=a}\\ +{\displaystyle {\int }_a^b}\frac{{\text{d}}^2}{\text{d}{x}^2}{\varphi }_i\left(x\right)\frac{\text{d}\tilde{q}}{\text{d}x}\text{d}x.\end{array}$

Substituting these in (7), the system of weighted simultaneous residual equations is obtained

$\begin{array}{l}-{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}\frac{\text{d}\tilde{p}}{\text{d}x}\right]}_{x=b}+{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}\frac{\text{d}\tilde{p}}{\text{d}x}\right]}_{x=a}+{\displaystyle {\int }_a^b}{{\varphi }^{″}}_i\left(x\right){\tilde{p}}^{\prime }\left(x\right)\text{d}x+{\displaystyle {\int }_a^b}(a_1\left(x\right){\tilde{p}}^{″}+a_2\left(x\right){\tilde{p}}^{\prime }\\ +a_3\left(x\right)\tilde{p}+a_4\left(x\right){\tilde{q}}^{″}+a_5\left(x\right){\tilde{q}}^{\prime }+a_6\left(x\right)\tilde{q}){\varphi }_i\left(x\right)\text{d}x={\displaystyle {\int }_a^b}f\left(x\right){\varphi }_i\left(x\right)\text{d}x\end{array}$(8a)

$\begin{array}{l}-{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}\frac{\text{d}\tilde{q}}{\text{d}x}\right]}_{x=b}+{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}\frac{\text{d}\tilde{q}}{\text{d}x}\right]}_{x=a}+{\displaystyle {\int }_a^b}{{\varphi }^{″}}_i\left(x\right){\tilde{q}}^{\prime }\left(x\right)\text{d}x+{\displaystyle {\int }_a^b}(b_1\left(x\right){\tilde{q}}^{″}+b_2\left(x\right){\tilde{q}}^{\prime }\\ +b_3\left(x\right)\tilde{q}+b_4\left(x\right){\tilde{p}}^{″}+b_5\left(x\right){\tilde{p}}^{\prime }+b_6\left(x\right)\tilde{p}){\varphi }_i\left(x\right)\text{d}x={\displaystyle {\int }_a^b}g\left(x\right){\varphi }_i\left(x\right)\text{d}x\end{array}$(8b)

Applying the boundary conditions (3) in (8a) and (8b), we have respectively

$\begin{array}{l}-{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}{\beta }_1\right]}_{x=b}+{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}{\gamma }_1\right]}_{x=a}+{\displaystyle {\int }_a^b}{{\varphi }^{″}}_i\left(x\right){\tilde{p}}^{\prime }\left(x\right)\text{d}x+{\displaystyle {\int }_a^b}(a_1\left(x\right){\tilde{p}}^{″}+a_2\left(x\right){\tilde{p}}^{\prime }\\ +a_3\left(x\right)\tilde{p}+a_4\left(x\right){\tilde{q}}^{″}+a_5\left(x\right){\tilde{q}}^{\prime }+a_6\left(x\right)\tilde{q}){\varphi }_i\left(x\right)\text{d}x={\displaystyle {\int }_a^b}f\left(x\right){\varphi }_i\left(x\right)\text{d}x\end{array}$(9a)

$\begin{array}{l}-{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}{\beta }_2\right]}_{x=b}+{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}{\gamma }_2\right]}_{x=a}+{\displaystyle {\int }_a^b}{{\varphi }^{″}}_i\left(x\right){\tilde{q}}^{\prime }\left(x\right)\text{d}x+{\displaystyle {\int }_a^b}(b_1\left(x\right){\tilde{q}}^{″}+b_2\left(x\right){\tilde{q}}^{\prime }\\ +b_3\left(x\right)\tilde{q}+b_4\left(x\right){\tilde{p}}^{″}+b_5\left(x\right){\tilde{p}}^{\prime }+b_6\left(x\right)\tilde{p}){\varphi }_i\left(x\right)\text{d}x={\displaystyle {\int }_a^b}g\left(x\right){\varphi }_i\left(x\right)\text{d}x\end{array}$(9b)

Now putting trial solutions (4) into (9a) and (9b), and simplifying these for $i=\mathrm{1,2,}\cdots \mathrm{,}n$ , we obtain

$\begin{array}{l}{\displaystyle \sum _{j=1}^{n-1}}\left({\displaystyle {\int }_a^b}\left[{{\varphi }^{″}}_i\left(x\right){{\varphi }^{\prime }}_j\left(x\right)+\left(a_1\left(x\right){{\varphi }^{″}}_j\left(x\right)+a_2\left(x\right){{\varphi }^{\prime }}_j\left(x\right)+a_3\left(x\right){\varphi }_j\left(x\right)\right){\varphi }_i\left(x\right)\right]\text{d}x\right)a_j\\ +{\displaystyle \sum _{j=1}^{n-1}}\left({\displaystyle {\int }_a^b}\left(a_4\left(x\right){{\varphi }^{″}}_j\left(x\right)+a_5\left(x\right){{\varphi }^{\prime }}_j\left(x\right)+a_6\left(x\right){\varphi }_j\left(x\right)\right){\varphi }_i\left(x\right)\text{d}x\right)b_j\\ ={\displaystyle {\int }_a^b}f\left(x\right){\varphi }_i\left(x\right)\text{d}x+{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}{\beta }_1\right]}_{x=b}-{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}{\gamma }_1\right]}_{x=a},\end{array}$ (10a)

$\begin{array}{l}{\displaystyle \sum _{j=1}^{n-1}}\left({\displaystyle {\int }_a^b}\left[{{\varphi }^{″}}_i\left(x\right){{\varphi }^{\prime }}_j\left(x\right)+\left(b_1\left(x\right){{\varphi }^{″}}_j\left(x\right)+b_2\left(x\right){{\varphi }^{\prime }}_j\left(x\right)+b_3\left(x\right){\varphi }_j\left(x\right)\right){\varphi }_i\left(x\right)\right]\text{d}x\right)b_j\\ +{\displaystyle \sum _{j=1}^{n-1}}\left({\displaystyle {\int }_a^b}\left(b_4\left(x\right){{\varphi }^{″}}_j\left(x\right)+b_5\left(x\right){{\varphi }^{\prime }}_j\left(x\right)+b_6\left(x\right){\varphi }_j\left(x\right)\right){\varphi }_i\left(x\right)\text{d}x\right)a_j\\ ={\displaystyle {\int }_a^b}g\left(x\right){\varphi }_i\left(x\right)\text{d}x+{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}{\beta }_2\right]}_{x=b}-{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}{\gamma }_2\right]}_{x=a}.\end{array}$ (10b)

The matrix form is equivalent to the previously mentioned equations

$\{\begin{array}{l}{\displaystyle \sum _{j=1}^{n-1}}\left(A_{j,i}{a}_j+H_{j,i}{b}_j\right)=F_i\hfill \\ {\displaystyle \sum _{j=1}^{n-1}}\left(C_{j,i}{b}_j+D_{j,i}{a}_j\right)=G_i\hfill \end{array}$ (11)

where,

$A_{j\mathrm{,}i}={\displaystyle {\int }_a^b}\left[{{\varphi }^{″}}_i\left(x\right){{\varphi }^{\prime }}_j\left(x\right)+\left(a_1\left(x\right){{\varphi }^{″}}_j\left(x\right)+a_2\left(x\right){{\varphi }^{\prime }}_j\left(x\right)+a_3\left(x\right){\varphi }_j\left(x\right)\right){\varphi }_i\left(x\right)\right]\text{d}x$

$H_{j\mathrm{,}i}={\displaystyle {\int }_a^b}\left(a_4\left(x\right){{\varphi }^{″}}_j\left(x\right)+a_5\left(x\right){{\varphi }^{\prime }}_j\left(x\right)+a_6\left(x\right){\varphi }_j\left(x\right)\right){\varphi }_i\left(x\right)\text{d}x$

$F_i={\displaystyle {\int }_a^b}f\left(x\right){\varphi }_i\left(x\right)\text{d}x+{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}{\beta }_1\right]}_{x=b}-{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}{\gamma }_1\right]}_{x=a}$

$C_{j\mathrm{,}i}={\displaystyle {\int }_a^b}\left[{{\varphi }^{″}}_i\left(x\right){{\varphi }^{\prime }}_j\left(x\right)+\left(b_1\left(x\right){{\varphi }^{″}}_j\left(x\right)+b_2\left(x\right){{\varphi }^{\prime }}_j\left(x\right)+b_3\left(x\right){\varphi }_j\left(x\right)\right){\varphi }_i\left(x\right)\right]\text{d}x$

$D_{j\mathrm{,}i}={\displaystyle {\int }_a^b}\left(b_4\left(x\right){{\varphi }^{″}}_j\left(x\right)+b_5\left(x\right){{\varphi }^{\prime }}_j\left(x\right)+b_6\left(x\right){\varphi }_j\left(x\right)\right){\varphi }_i\left(x\right)\text{d}x$

$G_i={\displaystyle {\int }_a^b}g\left(x\right){\varphi }_i\left(x\right)\text{d}x+{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}{\beta }_2\right]}_{x=b}-{\left[\frac{\text{d}{\varphi }_i}{\text{d}x}{\gamma }_2\right]}_{x=a}$

$i\mathrm{,}j=\mathrm{1,2,}\cdots \mathrm{,}n$

Now, we are going to solve nonlinear system of two third order differential equations in two unknown functions $p\left(x\right)$ and $q\left(x\right)$ of the form below is taking into consideration [8]

$\{\begin{array}{l}{a}_1\left(x\right)p^{‴}+a_2\left(x\right)p^{″}+a_3\left(x\right)p^{\prime }+a_4\left(x\right)p+a_5\left(x\right)q^{‴}+a_6\left(x\right)q^{″}\\ +a_7\left(x\right)q^{\prime }+a_8\left(x\right)q+M_1\left(p\mathrm{,}q\right)=g_1\left(x\right)\\ b_1\left(x\right)p^{‴}+b_2\left(x\right)p^{″}+b_3\left(x\right)p^{\prime }+b_4\left(x\right)p+b_5\left(x\right)q^{‴}+b_6\left(x\right)q^{″}\\ +b_7\left(x\right)q^{\prime }+b_8\left(x\right)q+M_2\left(p\mathrm{,}q\right)=g_2\left(x\right)\end{array}$ (12)

where $a_j\left(x\right)$ , $b_j\left(x\right)$ , $g_1\left(x\right)$ , $g_2\left(x\right)$ are given functions, $M_1$ , $M_2$ represent nonlinear functions and $a_j\left(x\right)$ , $b_j\left(x\right)$ are continuous for $j=\mathrm{1,2,}\cdots \mathrm{,8}$ . Since the system of third-order boundary value problems consist of nonlinear terms, then we can’t solve the system directly. In this case, mathematical formulation depends on the problem. In order to get the initial values of the parameters, we must neglect the nonlinear terms and solve the linear system. After getting the initial values of the parameters we will be able to solve the system. Then putting the parameters into the trial solutions, we will obtain the approximate solutions for the functions $p\left(x\right)$ and $q\left(x\right)$ . The details are described in the following section.

4. Numerical Results & Discussions

In this study, we use four systems; one linear and three nonlinear, which are available in the existing literature. To verify the effectiveness of the derived formulations, Dirichlet and Neumann boundary conditions are considered. For each case we find the approximate solutions using different number of parameters with Bernstein polynomials, and we compare these solutions with the exact solutions, and graphically which are shown in the same figures. Since the convergence of linear BVP is calculated by

$E=\left|{\tilde{p}}_{n+1}\left(x\right)-{\tilde{p}}_n\left(x\right)\right|<{\delta }_1$

where ${\tilde{p}}_n\left(x\right)$ denotes the approximate solution using n polynomials and ${\delta }_1$ depends on the problem. In this case, ${\delta }_1<{10}^{-8}$ . In addition, the convergence of nonlinear BVP is assumed when the absolute error of two consecutive iterations, ${\delta }_2$ satisfies

$E=\left|{\tilde{p}}_n^{N+1}\left(x\right)-{\tilde{p}}_n^N\left(x\right)\right|<{\delta }_2$

and, in this case, ${\delta }_2<{10}^{-10}$ .

Example 1. Consider the following system of third-order nonlinear boundary value problem [21] [22]

$\{\begin{array}{l}{p}^{‴}\left(x\right)+2p^{\prime }\left(x\right)+xq\left(x\right)=f\left(x\right),0<x<1\\ q^{‴}\left(x\right)+\frac{p^{″}\left(x\right)q^{″}\left(x\right)}{6}=g\left(x\right),0<x<1\\ p\left(0\right)=q\left(0\right)=0,\mathrm{<mi>\; </mi><mi>\; </mi>}p\left(1\right)=q\left(1\right)=0,\mathrm{<mi>\; </mi><mi>\; </mi>}{p}^{\prime }\left(0\right)=q^{\prime }\left(0\right)=0\end{array}$ (13)

where

$f\left(x\right)=x^5-x^3-18x^2+12x-18$

$g\left(x\right)=-36x^3+12x^2+30x-2$

The exact solutions are $p\left(x\right)=3x^2-3x^3$ and $q\left(x\right)=x^4-x^2$ .

Here, we use Bernstein polynomials as trial approximate solution to solve the system (13). Let us consider the trial solution of the form

$\{\begin{array}{l}\tilde{p}\left(x\right)={\theta }_0\left(x\right)+{\displaystyle \sum _{i=1}^{n-1}}{a}_i{\varphi }_i\left(x\right)\hfill \\ \tilde{q}\left(x\right)={\varphi }_0\left(x\right)+{\displaystyle \sum _{i=1}^{n-1}}{b}_i{\varphi }_i\left(x\right)\hfill \end{array}$(14)

where $a_i$ and $b_i$ are parameters and ${\varphi }_i\left(x\right)$ are trial functions (Bernstein polynomials) which satisfy the boundary conditions. Here, we can choose ${\theta }_0\left(x\right)=0$ and ${\varphi }_0\left(x\right)=0$ because the boundary conditions are homogeneous.

Now applying the Galerkin method, we get weighted residual system of equations

$\{\begin{array}{l}{\displaystyle {\int }_0^1}\left({\tilde{p}}^{‴}\left(x\right)+2{\tilde{p}}^{\prime }\left(x\right)+x\tilde{q}\left(x\right)\right){\varphi }_i\left(x\right)\text{d}x={\displaystyle {\int }_0^1}f\left(x\right){\varphi }_i\left(x\right)\text{d}x\hfill \\ {\displaystyle {\int }_0^1}\left({\tilde{q}}^{‴}\left(x\right)+\frac{{\tilde{p}}^{″}\left(x\right){\tilde{q}}^{″}\left(x\right)}{6}\right){\varphi }_i\left(x\right)\text{d}x={\displaystyle {\int }_0^1}g\left(x\right){\varphi }_i\left(x\right)\text{d}x\hfill \end{array}$(15)

Applying integration by parts in the first term of (26), we obtain

$\begin{array}{c}{\displaystyle {\int }_0^1}\frac{{\text{d}}^3\tilde{p}}{\text{d}{x}^3}{\varphi }_i\left(x\right)\text{d}x={\left[\frac{\text{d}\tilde{p}}{\text{d}x}{\varphi }_i(x)\right]}_0^1-{\displaystyle {\int }_0^1}\frac{\text{d}\tilde{p}}{\text{d}x}\frac{\text{d}{\varphi }_i}{\text{d}x}\text{d}x\\ =\frac{\text{d}\tilde{p}}{\text{d}x}\left(x=1\right){\varphi }_i\left(1\right)-\frac{\text{d}\tilde{p}}{\text{d}x}\left(x=0\right){\varphi }_i\left(0\right)-{\displaystyle {\int }_0^1}\frac{\text{d}\tilde{p}}{\text{d}x}\frac{\text{d}{\varphi }_i}{\text{d}x}\text{d}x\\ =-{\left[\frac{\text{d}}{\text{d}x}{\varphi }_i\left(x\right)\frac{\text{d}\tilde{p}}{\text{d}x}\right]}_0^1+{\displaystyle {\int }_0^1}\frac{{\text{d}}^2}{\text{d}{x}^2}{\varphi }_i\left(x\right)\frac{\text{d}\tilde{p}}{\text{d}x}\text{d}x\end{array}$

${\displaystyle {\int }_0^1}\frac{{\text{d}}^3\tilde{q}}{\text{d}{x}^3}{\varphi }_i\left(x\right)\text{d}x=-{\left[\frac{\text{d}}{\text{d}x}{\varphi }_i\left(x\right)\frac{\text{d}\tilde{q}}{\text{d}x}\right]}_0^1+{\displaystyle {\int }_0^1}\frac{{\text{d}}^2}{\text{d}{x}^2}{\varphi }_i\left(x\right)\frac{\text{d}\tilde{q}}{\text{d}x}\text{d}x$

This system can be converted into modified Galerkin form in the usual sense, and using (14) to obtain

$\begin{array}{l}{\displaystyle \sum _{j=1}^{n-1}}\left({\displaystyle {\int }_0^1}\left[{{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{″}}_i\left(x\right)+2{{\varphi }^{\prime }}_j\left(x\right){\varphi }_i\left(x\right)\right]\text{d}x-{\left[{{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{\prime }}_i\left(x\right)\right]}_{x=1}\right)a_j\\ +{\displaystyle \sum _{j=1}^{n-1}}\left({\displaystyle {\int }_0^1}x{\varphi }_j\left(x\right){\varphi }_i\left(x\right)\text{d}x\right)b_j={\displaystyle {\int }_0^1}f\left(x\right){\varphi }_i\left(x\right)\text{d}x,i=1,2,\cdots ,n\end{array}$ (16a)

Figure 1. Exact and approximate solutions of $p\left(x\right)$ for Example 1.

Table 1. Absolute errors $\left|p\left(x\right)-\tilde{p}\left(x\right)\right|$ for Example 1.

x	Present Method (GWRM)			HPM [21]
	Degree 3	Degree 4	Degree 5
0.1	3.010565 × 10−5	8.406867 × 10−10	1.680418 × 10−9	8.74250 × 10−8
0.2	1.976596 × 10−5	3.596601 × 10−8	2.284382 × 10−9	3.41082 × 10−5
0.3	1.836086 × 10−5	9.477548 × 10−8	5.387116 × 10−9	7.35812 × 10−5
0.4	7.161663 × 10−5	1.653393 × 10−7	1.220891 × 10−8	1.23210 × 10−4
0.5	1.273431 × 10−4	2.343983 × 10−7	2.200062 × 10−8	6.24822 × 10−1
0.6	1.728822 × 10−4	2.873641 × 10−7	3.242916 × 10−8	2.31537 × 10−4
0.7	1.955756 × 10−4	3.083190 × 10−7	3.996300 × 10−8	2.75507 × 10−5
0.8	1.827652 × 10−4	2.800157 × 10−7	4.025768 × 10−8	2.92801 × 10−4
0.9	1.217927 × 10−4	1.838780 × 10−7	2.854130 × 10−8	2.39684 × 10−4

$\begin{array}{l}{\displaystyle \sum _{j=1}^{n-1}}\left({\displaystyle {\int }_0^1}{{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{″}}_i\left(x\right)\text{d}x-{\left[{{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{\prime }}_i\left(x\right)\right]}_{x=1}\right)b_j\\ +\frac{1}{6}{\displaystyle {\int }_0^1}\left({\displaystyle \sum _{k=1}^{n-1}}{a}_k{{\varphi }^{″}}_k\left(x\right)\right)\left({\displaystyle \sum _{k=1}^{n-1}}{b}_k{{\varphi }^{″}}_k\left(x\right)\right){\varphi }_i\left(x\right)\text{d}x={\displaystyle {\int }_0^1}g\left(x\right){\varphi }_i\left(x\right)\text{d}x,i=1,2,\cdots ,n\end{array}$ (16b)

The above equations are equivalent to the matrix form

${\displaystyle \sum _{j=1}^{n-1}}\left(A_{j,i}{a}_j+H_{j,i}{b}_j\right)=F_i,i=1,2,\cdots ,n$ (17a)

${\displaystyle \sum _{j=1}^{n-1}}{C}_{j,i}{b}_j=G_i+N_i,i=1,2,\cdots ,n$ (17b)

in which

$A_{j\mathrm{,}i}={\displaystyle {\int }_0^1}\left({{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{″}}_i\left(x\right)+2{{\varphi }^{\prime }}_j\left(x\right){\varphi }_i\left(x\right)\right)\text{d}x-{\left[{{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{\prime }}_i\left(x\right)\right]}_{x=1}$

$H_{j\mathrm{,}i}={\displaystyle {\int }_0^1}x{\varphi }_j\left(x\right){\varphi }_i\left(x\right)\text{d}x$

$F_i={\displaystyle {\int }_0^1}\left(x^5-x^3-18x^2+12x-18\right){\varphi }_i\left(x\right)\text{d}x$

$C_{j\mathrm{,}i}={\displaystyle {\int }_0^1}{{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{″}}_i\left(x\right)\text{d}x-{\left[{{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{\prime }}_i\left(x\right)\right]}_{x=1}$

$G_i={\displaystyle {\int }_0^1}\left(-36x^3+12x^2+30x-2\right){\varphi }_i\left(x\right)\text{d}x$

$N_i=-\frac{1}{6}{\displaystyle {\int }_0^1}\left({\displaystyle \sum _{k=1}^n}{a}_k{{\varphi }^{″}}_k\left(x\right)\right)\left({\displaystyle \sum _{k=1}^n}{b}_k{{\varphi }^{″}}_k\left(x\right)\right){\varphi }_i\left(x\right)\text{d}x$

$i,j=1,2,\cdots ,n$

Neglecting the nonlinear term $N_i$ in (17b), the initial values of $a_j$ and $b_j$ are obtained. Therefore, we obtain initial values solving the system

$\{\begin{array}{l}{\displaystyle \sum _{j=1}^{n-1}}\left(A_{j,i}{a}_j+H_{j,i}{b}_j\right)=F_i,i=1,2,\cdots ,n\hfill \\ {\displaystyle \sum _{j=1}^{n-1}}{C}_{j,i}{b}_j=G_i,i=1,2,\cdots ,n\hfill \end{array}$ (18)

Table 2. Absolute errors $\left|q\left(x\right)-\tilde{q}\left(x\right)\right|$ for Example 1.

x	Present Method (GWRM)			HPM [21]
	Degree 3	Degree 4	Degree 5
0.1	2.305211 × 10−2	1.390676 × 10−8	3.551924 × 10−7	1.34169 × 10−5
0.2	2.226453 × 10−2	7.162991 × 10−6	1.658188 × 10−6	5.22346 × 10−5
0.3	8.856149 × 10−3	1.724459 × 10−5	2.285279 × 10−6	1.12022 × 10−4
0.4	8.354171 × 10−3	2.692795 × 10−5	4.263502 × 10−7	9.61854 × 10−2
0.5	2.294756 × 10−2	3.375417 × 10−5	3.196644 × 10−6	2.00026 × 10−1
0.6	3.090514 × 10−2	3.613628 × 10−5	6.873504 × 10−6	3.30775 × 10−4
0.7	3.060805 × 10−2	3.335918 × 10−5	8.739233 × 10−6	3.78034 × 10−4
0.8	2.283741 × 10−2	2.557966 × 10−5	7.607479 × 10−6	3.86227 × 10−4
0.9	1.077435 × 10−2	1.382641 × 10−5	3.803978 × 10−6	3.11011 × 10−4

Then we substitute these values into (17a) and (17b) and obtain the new values of $a_j$ and $b_j$ . The iterative process keeps going till the converged estimates of the unknown coefficients are achieved. We have an approximation to the BVP (13) by replacing the final quantities in (14).

Satisfying the homogeneous boundary conditions and using degree of polynomial 3, we may obtain the approximate solutions of $p\left(x\right)$ and $q\left(x\right)$ are

$\tilde{p}\left(x\right)=0.00054548x{\left(1-x\right)}^2+2.99843577x^2\left(1-x\right)$ and

$\tilde{q}\left(x\right)=0.39311569x{\left(1-x\right)}^2-2.07669616x^2\left(1-x\right)$ , respectively.

Table 1 shows the numerical outcomes for the given problem. Here we have used 2, 3 and 4 Bernstein polynomials in column 2, 3 and 4 respectively. Consider the fact that we iterate four times in order to achieve the approximations.

Here, we see that the maximum absolute error of $p\left(x\right)$ using 3, 4 and 5 degree polynomials are 1.955756 × 10⁻⁴, 3.083190 × 10⁻⁷ and 4.025768 × 10⁻⁸, respectively whereas the maximum absolute errors using Homotopy Perturbation Method (HPM) in [21] and Block method in [22] are 6.24822 × 10⁻¹ and 6.25 × 10⁻⁴, respectively. On the other hand, Table 2 shows that the maximum absolute error of $q\left(x\right)$ using 3, 4 and 5 degree polynomials are 3.090514 × 10⁻², 3.613628 × 10⁻⁵ and 8.739233 × 10⁻⁶, respectively. The maximum absolute errors using HPM in [21] and Block method in [22] are 2.00026 × 10⁻¹ and 1.48 × 10⁻², respectively. Thus, our proposed method reflects the better results than the previous results. In Figure 1 and Figure 2, we compare the exact solutions with the approximate solutions for $p\left(x\right)$ and $q\left(x\right)$ , respectively. What we’ve noticed is that when we use polynomials of degree 5, the results show very good agreement between exact and approximate solutions.

Example 2. Consider the following nonlinear system of third-order boundary value problem [21] [22]

Figure 2. Exact and approximate solutions of $q\left(x\right)$ for Example 1.

$\{\begin{array}{l}{p}^{‴}\left(x\right)-4q^{″}\left(x\right)+p^{″}\left(x\right)q^{\prime }\left(x\right)=f\left(x\right),0<x<1\\ q^{‴}\left(x\right)+4q^{\prime }\left(x\right)-p^{″}\left(x\right)+p^{\prime }\left(x\right)q^{″}\left(x\right)=g\left(x\right),0<x<1\\ p\left(0\right)=q\left(0\right)=0,\mathrm{<mi>\; </mi><mi>\; </mi>}p\left(1\right)=q\left(1\right)=1,\mathrm{<mi>\; </mi><mi>\; </mi>}{p}^{\prime }\left(0\right)=q^{\prime }\left(0\right)=0\end{array}$ (19)

where

$f\left(x\right)=36x^4\text{and}g\left(x\right)=24x^4+6$

The exact solutions are $p\left(x\right)=x^4$ and $q\left(x\right)=x^3$ .

Let us consider the trial solution of the form

$\{\begin{array}{l}\tilde{p}\left(x\right)={\theta }_0\left(x\right)+{\displaystyle \sum _{i=1}^n}{a}_i{\varphi }_i\left(x\right)\hfill \\ \tilde{q}\left(x\right)={\varphi }_0\left(x\right)+{\displaystyle \sum _{i=1}^n}{b}_i{\varphi }_i\left(x\right)\hfill \end{array}$(20)

where $a_i$ and $b_i$ are parameters and ${\varphi }_i\left(x\right)$ are trial functions (Bernstein polynomials) which satisfy the boundary conditions. Here, we can choose ${\theta }_0\left(x\right)=x^2$ and ${\varphi }_0\left(x\right)=x^2$ in order to satisfy the boundary conditions.

Now applying the Galerkin method, we get weighted residual system of equations

$\{\begin{array}{l}{\displaystyle {\int }_0^1}\left({\tilde{p}}^{‴}\left(x\right)-4{\tilde{q}}^{″}\left(x\right)+{\tilde{p}}^{″}\left(x\right){\tilde{q}}^{\prime }\left(x\right)\right){\varphi }_i\left(x\right)\text{d}x={\displaystyle {\int }_0^1}f\left(x\right){\varphi }_i\left(x\right)\text{d}x\hfill \\ {\displaystyle {\int }_0^1}\left({\tilde{q}}^{‴}\left(x\right)+4{\tilde{q}}^{\prime }\left(x\right)-{\tilde{p}}^{″}\left(x\right)+{\tilde{p}}^{\prime }\left(x\right){\tilde{q}}^{″}\left(x\right)\right){\varphi }_i\left(x\right)\text{d}x={\displaystyle {\int }_0^1}g\left(x\right){\varphi }_i\left(x\right)\text{d}x\hfill \end{array}$(21)

This system can be converted to modified Galerkin form in the usual way and using (20) we finally obtain the matrix form

${\displaystyle \sum _{j=1}^n}\left(A_{j,i}{a}_j+H_{j,i}{b}_j\right)=F_i+N_i,i=1,2,\cdots ,n$ (22a)

${\displaystyle \sum _{j=1}^n}\left(C_{j,i}{b}_j+D_{j,i}{a}_j\right)=G_i+M_i,i=1,2,\cdots ,n$ (22b)

where

$A_{j,i}={\displaystyle {\int }_0^1}\left[{{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{″}}_i\left(x\right)+2x{{\varphi }^{″}}_j\left(x\right){\varphi }_i\left(x\right)\right]\text{d}x-{\left[{{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{\prime }}_i\left(x\right)\right]}_{x=1}$

$H_{j,i}={\displaystyle {\int }_0^1}\left(4{{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{\prime }}_i\left(x\right)+2{{\varphi }^{\prime }}_j\left(x\right){\varphi }_i\left(x\right)\right)\text{d}x$

$F_i={\displaystyle {\int }_0^1}\left(f\left(x\right){\varphi }_i\left(x\right)-4x{\varphi }_i\left(x\right)-8x{{\varphi }^{\prime }}_i\left(x\right)-2x{{\varphi }^{″}}_i\left(x\right)\right)\text{d}x+{\left[2x{{\varphi }^{\prime }}_i\left(x\right)\right]}_{x=1}$

$N_i=-{\displaystyle {\int }_0^1}\left({\displaystyle \sum _{k=1}^n}{a}_k{{\varphi }^{″}}_k\left(x\right)\right)\left({\displaystyle \sum _{k=1}^n}{b}_k{{\varphi }^{\prime }}_k\left(x\right)\right){\varphi }_i\left(x\right)\text{d}x$

$C_{j,i}={\displaystyle {\int }_0^1}\left({{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{″}}_i\left(x\right)+4{{\varphi }^{\prime }}_j\left(x\right){\varphi }_i\left(x\right)+2x{{\varphi }^{″}}_j\left(x\right){\varphi }_i\left(x\right)\right)\text{d}x-{\left[{{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{\prime }}_i\left(x\right)\right]}_{x=1}$

$D_{j,i}={\displaystyle {\int }_0^1}\left(4{{\varphi }^{\prime }}_j\left(x\right){{\varphi }^{\prime }}_i\left(x\right)+2{{\varphi }^{\prime }}_j\left(x\right){\varphi }_i\left(x\right)\right)\text{d}x$

$G_i={\displaystyle {\int }_0^1}\left(g\left(x\right){\varphi }_i\left(x\right)-12x{\varphi }_i\left(x\right)-2x{{\varphi }^{\prime }}_i\left(x\right)-2x{{\varphi }^{″}}_i\left(x\right)\right)\text{d}x+{\left[2x{{\varphi }^{\prime }}_i\left(x\right)\right]}_{x=1}$

$M_i=-{\displaystyle {\int }_0^1}\left({\displaystyle \sum _{k=1}^n}{a}_k{{\varphi }^{\prime }}_k\left(x\right)\right)\left({\displaystyle \sum _{k=1}^n}{b}_k{{\varphi }^{″}}_k\left(x\right)\right){\varphi }_i\left(x\right)\text{d}x$

$i\mathrm{,}j=\mathrm{1,2,}\cdots \mathrm{,}n$

Neglecting the nonlinear terms $N_i$ and $M_i$ in (22a) and (22b) respectively, the initial values of $a_j$ and $b_j$ are obtained. Therefore, we obtain the initial values solving the following system

$\{\begin{array}{l}{\displaystyle \sum _{j=1}^n}\left(A_{j,i}{a}_j+H_{j,i}{b}_j\right)=F_i,i=1,2,\cdots ,n\hfill \\ {\displaystyle \sum _{j=1}^n}\left(C_{j,i}{b}_j+D_{j,i}{a}_j\right)=G_i,i=1,2,\cdots ,n\hfill \end{array}$ (23)

Then we substitute these values in (22a) and (22b) and obtain the new values of $a_j$ and $b_j$ . The iterative process keeps going till the converged estimates of the unknown coefficients are achieved. We have an approximation to the BVP (19) by replacing the final quantities in (20). The approximate solutions of $p\left(x\right)$ and $q\left(x\right)$ using two parameters with polynomial degree 3 are:

$\tilde{p}\left(x\right)=2.36893939x^3-1.77969288x^2+0.41075348x$ , and

$\tilde{q}\left(x\right)=1.03113495x^3-0.05887466x^2+0.02773971x$ , respectively.

We’ve compared exact and approximate solutions for $p\left(x\right)$ and $q\left(x\right)$ in Figure 3 and Figure 4, respectively for Example 2. In Figure 3, it’s evident that employing fifth-degree polynomials yields better outcomes for the function $p\left(x\right)$ when compared to third and fourth-degree polynomials. Table 3 and Table 4 show the numerical outcomes for the given problem. Here we have used 2, 3 and 4 Bernstein polynomials in column 2, 3 and 4 respectively. Consider the

Figure 3. Exact and approximate solutions of $p\left(x\right)$ for Example 2.

Figure 4. Exact and approximate solutions of $q\left(x\right)$ for Example 2.

fact that we iterate four times in order to achieve the approximations. Here, we see that the maximum absolute error of $p\left(x\right)$ in Table 3 using 2, 3 and 4 polynomials are 2.831450 × 10⁻², 1.616954 × 10⁻⁷ and 8.469325 × 10⁻⁸, respectively

Table 3. Absolute errors for $\left|p\left(x\right)-\tilde{p}\left(x\right)\right|$ for Example 2.

x	Present Method (GWRM)			HPM [21]
	Degree 3	Degree 4	Degree 5
0.1	2.554736 × 10−2	1.060524 × 10−7	7.896618 × 10−8	8.74250 × 10−8
0.2	2.831450 × 10−2	6.744158 × 10−8	1.381743 × 10−8	1.08399 × 10−3
0.3	1.891505 × 10−2	2.037354 × 10−8	2.377000 × 10−8	4.71442 × 10−3
0.4	5.562655 × 10−3	9.053284 × 10−8	3.918424 × 10−8	1.14828 × 10−2
0.5	5.929053 × 10−3	1.047750 × 10−7	5.790583 × 10−8	2.18865 × 10−2
0.6	1.214644 × 10−2	5.343750 × 10−8	7.518965 × 10−8	3.5986 × 10−2
0.7	1.207586 × 10−2	4.454330 × 10−8	8.469325 × 10−8	6.87692 × 10−2
0.8	7.103685 × 10−3	1.416323 × 10−7	7.947025 × 10−8	7.66527 × 10−2
0.9	1.016277 × 10−3	1.616954 × 10−7	5.296380 × 10−8	6.20685 × 10−2

Table 4. Absolute errors for $\left|q\left(x\right)-\tilde{q}\left(x\right)\right|$ for Example 2.

x	Present Method (GWRM)			HPM [21]
	Degree 3	Degree 4	Degree 5
0.1	2.216359 × 10−3	8.201831 × 10−8	5.335278 × 10−9	6.73806 × 10−4
0.2	3.442034 × 10−3	7.169045 × 10−8	1.268373 × 10−9	2.82326 × 10−3
0.3	3.863836 × 10−3	3.608265 × 10−8	7.607422 × 10−9	6.61192 × 10−3
0.4	3.668574 × 10−3	1.982045 × 10−8	1.026833 × 10−8	1.21045 × 10−3
0.5	3.043057 × 10−3	4.508859 × 10−8	1.119354 × 10−8	1.91209 × 10−2
0.6	2.174096 × 10−3	1.116311 × 10−7	1.424742 × 10−8	2.69597 × 10−2
0.7	1.248500 × 10−3	1.967514 × 10−7	2.178228 × 10−8	3.39306 × 10−2
0.8	4.530788 × 10−4	2.553118 × 10−7	3.120430 × 10−8	6.57775 × 10−2
0.9	2.535774 × 10−5	2.197343 × 10−7	3.153957 × 10−8	2.88236 × 10−2

whereas the maximum absolute errors using Homotopy Perturbation Method (HPM) in [21] and Block method in [22] are 7.66527 × 10⁻² and 5.50 × 10⁻³, respectively.

On the other hand, Table 4 shows that the maximum absolute error of $q\left(x\right)$ using 2, 3 and 4 polynomials are 3.863836 × 10⁻³, 2.553118 × 10⁻⁷ and 3.153957 × 10⁻⁸, respectively. The maximum absolute errors using HPM in [21] and Block method in [22] are 6.57775 × 10⁻² and 2.39 × 10⁻², respectively. In this example, we show that the proposed method is far better than the existing methods in [21] and in [22].

5. Applications

Although in [24], Agarwal has extensively covered the theorems of the criteria for the existence and uniqueness of solutions of the sixth-order BVPs in a book, it does not include any numerical techniques. Islam and Hossain [23] solved the sixth-order BVPs using Galerkin method. Modified decomposition method was used in [25] to find the solution of the sixth-order BVPs by Wazwaz. Aasma Khalid et al. used Cubic B-spline in [26] and [27] in order to solve the linear and nonlinear sixth-order BVPs, respectively. Cubic-nonpolynomial spline (CNPS) and cubic-polynomial spline (CPS) were used to obtain the solutions of BVPs in [28] arising in hydrodynamic and magnetohydro-dynamic stability theory. Noor and Mohyud-Din [29] solved the sixth-order BVP using homotopy perturbation method. However, in this section, we show that higher even order BVP may be solved, in the alternative way, by the method of reduction order into system of lower order BVPs. For this, we experiment the proposed method to sixth order BVP.

Example 3 Consider the linear sixth-order boundary value problem [23] [26] [29]

$\frac{{\text{d}}^{6}p}{\text{d}{x}^6}-p=-6e^x,\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi>0}\le x\le 1$ (24)

Subject to the boundary conditions

$p\left(0\right)=1,\mathrm{<mi>\; </mi><mi>\; </mi>}p\left(1\right)=0,\mathrm{<mi>\; </mi><mi>\; </mi>}{p}^{\prime }\left(0\right)=0,\mathrm{<mi>\; </mi><mi>\; </mi>}{p}^{\prime }\left(1\right)=-e,\mathrm{<mi>\; </mi><mi>\; </mi>}{p}^{″}\left(0\right)=-1,\mathrm{<mi>\; </mi><mi>\; </mi>}{p}^{″}\left(1\right)=-2e$

The analytic solution of the above problem is, $p\left(x\right)=\left(1-x\right)e^x$ .

If we introduce a new function, say, $q\left(x\right)$ such that

$q=\frac{{\text{d}}^{3}p}{\text{d}{x}^3}$ (25)

then (28) is clearly equivalent to two third-order differential equations

$\{\begin{array}{l}\frac{{\text{d}}^{3}q}{\text{d}{x}^3}-p=-6e^x\\ \frac{{\text{d}}^{3}p}{\text{d}{x}^3}-q=0\end{array}$ (26)

with boundary conditions

$p\left(0\right)=1,\mathrm{<mi>\; </mi><mi>\; </mi>}p\left(1\right)=0,\mathrm{<mi>\; </mi><mi>\; </mi>}{p}^{\prime }\left(0\right)=0,\mathrm{<mi>\; </mi><mi>\; </mi>}q\left(0\right)=-2,\mathrm{<mi>\; </mi><mi>\; </mi>}q\left(1\right)=-3e,\mathrm{<mi>\; </mi><mi>\; </mi>}{q}^{\prime }\left(0\right)=-3$ (27)

Table 5 shows the numerical outcomes for the given problem. We have used 4 Bernstein polynomials of degree 5. Here, we see that the maximum absolute error of $p\left(x\right)$ in Table 5 is 1.84 × 10⁻⁵ whereas the maximum absolute errors using Cubic B-Spline method in [26] and Homotopy Perturbation Method (HPM) in [29] are 1.67 × 10⁻⁴ and 1.32 × 10⁻³, respectively.

Using the method illustrated in the previous section and for different number of polynomials, the maximum absolute errors and the previous results available in the literature are summarized in Table 6. The accuracy of the present method is remarkable.

Example 4. Consider the nonlinear sixth-order boundary value problem [25] [27] [28] [29]

Table 5. Exact, approximate and absolute errors of Example 3.

x	Present Method (GWRM)			Cubic B-Spline [26]	HPM [29]
	Exact	Approximate	Abs Error
0.1	0.99465383	0.99464299	1.08 × 10−5	1.18 × 10−5	4.09 × 10−4
0.2	0.97712221	0.97713362	1.14 × 10−5	4.29 × 10−5	7.78 × 10−4
0.3	0.94490117	0.94491960	1.84 × 10−5	8.53 × 10−5	1.07 × 10−3
0.4	0.89509482	0.89510135	6.53 × 10−6	1.28 × 10−4	1.26 × 10−3
0.5	0.82436064	0.82435116	9.48 × 10−6	1.59 × 10−4	1.32 × 10−3
0.6	0.72884752	0.72883236	1.52 × 10−5	1.67 × 10−4	1.26 × 10−3
0.7	0.60412581	0.60411856	7.26 × 10−6	1.45 × 10−4	1.07 × 10−3
0.8	0.44510819	0.44511280	4.61 × 10−6	9.47 × 10−5	7.78 × 10−4
0.9	0.24596031	0.24596678	6.46 × 10−6	3.33 × 10−5	4.09 × 10−4

Table 6. Absolute errors for $\left|p\left(x\right)-\tilde{p}\left(x\right)\right|$ for Example 3.

Number of Polynomial used	Present Method (GWRM)	In [23]
10	1.165 × 10−14	1.126 × 10−13
11	1.986 × 10−16	2.311 × 10−15
12	4.771 × 10−18	2.220 × 10−16
13	7.228 × 10−20	2.220 × 10−16

$\frac{{\text{d}}^{6}p}{\text{d}{x}^6}=e^{-x}{p}^2\left(x\right)\mathrm{,<mi>\; </mi><mi>\; </mi><mi>\; </mi>0}\le x\le 1$ (28)

Subject to the boundary conditions

$p\left(0\right)=p^{″}\left(0\right)=p^{iv}\left(0\right)=1,\mathrm{<mi>\; </mi><mi>\; </mi>}p\left(1\right)=p^{″}\left(1\right)=p^{iv}\left(1\right)=e$

The analytic solution of the above problem is, $p\left(x\right)=e^x$ .

Introducing a new function, $q\left(x\right)$ , Equation (28) is equivalent to two third-order differential equations

$\{\begin{array}{l}\frac{{\text{d}}^{3}q}{\text{d}{x}^3}=e^{-x}{p}^2\hfill \\ \frac{{\text{d}}^{3}p}{\text{d}{x}^3}-q=0\hfill \end{array}$ (29)

with boundary conditions

$p\left(0\right)=1,\mathrm{<mi>\; </mi><mi>\; </mi>}p\left(1\right)=e,\mathrm{<mi>\; </mi><mi>\; </mi>}{p}^{\prime }\left(0\right)=1,\mathrm{<mi>\; </mi><mi>\; </mi>}q\left(0\right)=1,\mathrm{<mi>\; </mi><mi>\; </mi>}q\left(1\right)=e,\mathrm{<mi>\; </mi><mi>\; </mi>}{q}^{\prime }\left(0\right)=1$ (30)

Table 7 shows the numerical outcomes for the given problem. We have used 4 Bernstein polynomials of degree 5. Here, we see that the maximum absolute error of $p\left(x\right)$ in Table 7 is 5.04 × 10⁻⁵ whereas the maximum absolute errors using Modified Decomposition Method (MPM) in [25], Cubic B-Spline method in [27], Cubic Polynomial Spline (CPS) in [28] and Homotopy Perturbation

Table 7. Exact, approximate and absolute errors of Example 4.

x	GWRM	MDM [25]	Cubic B-Spline [27]	CPS [28]	HPM [29]
0.1	1.00 × 10−5	1.23 × 10−4	3.95 × 10−6	5.05 × 10−5	1.23 × 10−4
0.2	1.41 × 10−5	2.35 × 10−4	1.43 × 10−5	9.23 × 10−5	2.35 × 10−4
0.3	3.97 × 10−5	3.25 × 10−4	2.79 × 10−5	1.25 × 10−4	3.25 × 10−4
0.4	5.04 × 10−5	3.85 × 10−4	4.07 × 10−5	1.47 × 10−4	3.85 × 10−4
0.5	4.35 × 10−5	4.08 × 10−4	4.88 × 10−5	1.59 × 10−4	4.08 × 10−4
0.6	2.63 × 10−5	3.91 × 10−4	4.92 × 10−5	1.58 × 10−4	3.92 × 10−4
0.7	1.04 × 10−5	3.36 × 10−4	4.09 × 10−5	1.44 × 10−4	3.36 × 10−4
0.8	5.42 × 10−6	2.45 × 10−4	2.56 × 10−5	1.14 × 10−4	2.46 × 10−4
0.9	9.77 × 10−6	1.29 × 10−4	8.63 × 10−6	6.70 × 10−5	1.29 × 10−4

Method in [29] are 4.08 × 10⁻⁴, 4.92 × 10⁻⁵, 1.59 × 10⁻⁴ and 4.08 × 10⁻⁴, respectively. We may observe that the solutions converge monotonically from below. This concludes that the present method may be exploited with considerable accuracy.

6. Conclusion

We have deduced the formulation of the Galerkin weighted residual method for system of third-order boundary value problems in detail. We can determine the solutions at each point within the problem’s domain by using this method. Bernstein polynomials have been used in this method as the trial functions in the approximation. The focus has been on the formulations as well as on the performance of the accuracy. Some examples are tested to verify the effectiveness of the desired formulations. A good agreement has been established when comparing the approximate solutions with the exact solutions. Unlike finite difference or finite element methods, which rely on specific mesh structures, the Galerkin method operates directly in the function space, making it applicable to problems defined on irregular domains or with discontinuous coefficients. We have shown that this method may be applied to higher-order systems and/or single BVPs to get the desired accuracy.

Acknowledgements

The authors are thankful for the reviewer’s comments and suggestions for improving the quality of the first version of the manuscript. The first author is also grateful to the Ministry of Science and Technology, Bangladesh, for providing partial financial support (NST Fellowship) during the research period.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Lewis, P.E. and Ward, J.P. (1991) The Finite Element Method: Principles and Applications. Addison-Wesley Publishing, Boston.
[2]	Bhatti, M.I. and Bracken, P. (2007) Solutions of Differential Equations in a Bernstein Polynomial Basis. Journal of Computational and Applied Mathematics, 205, 272-280. https://doi.org/10.1016/j.cam.2006.05.002
[3]	Islam, M.S. and Shirin, A. (2011) Numerical Solutions of a Class of Second Order Boundary Value Problems on Using Bernoulli Polynomials. Applied Mathematics, 2, 1059-1067. https://doi.org/10.4236/am.2011.29147
[4]	Rashidinia, J., Mohammadi, R., Jalilian, R. and Ghasemi, M. (2007) Convergence of Cubic-Spline Approach to the Solution of a System of Boundary-Value Problems. Applied Mathematics and Computation, 192, 319-331. https://doi.org/10.1016/j.amc.2007.03.008
[5]	Siddiqi, S.S. and Akram, G. (2007) Solution of the System of Fourth-Order Boundary Value Problems Using Non-Polynomial Spline Technique. Applied Mathematics and Computation, 185, 128-135. https://doi.org/10.1016/j.amc.2006.07.014
[6]	Rupa, M.J. and Islam, M.S. (2018) Numerical Solutions of System of Second Order Boundary Value Problems Using Galerkin Method. GANIT: Journal of Bangladesh Mathematical Society, 37, 161-174. https://doi.org/10.3329/ganit.v37i0.35734
[7]	Geng, F. and Cui, M. (2007) Solving a Nonlinear System of Second Order Boundary Value Problems. Journal of Mathematical Analysis and Applications, 327, 1167-1181. https://doi.org/10.1016/j.jmaa.2006.05.011
[8]	Lu, J. (2007) Variational Iteration Method for Solving a Nonlinear System of Second-Order Boundary Value Problems. Computers & Mathematics with Applications, 54, 1133-1138. https://doi.org/10.1016/j.camwa.2006.12.060
[9]	Dehghan, M. and Saadatmandi, A. (2007) The Numerical Solution of a Nonlinear System of Second-Order Boundary Value Problems Using the Sinc-Collocation Method. Mathematical and Computer Modelling, 46, 1434-1441. https://doi.org/10.1016/j.mcm.2007.02.002
[10]	Du, J. and Cui, M. (2009) Solving a Nonlinear System of Second Order Two-Point Boundary Value Problem. Applied Mathematical Sciences, 3, 1141-1151.
[11]	Geng, F. and Cui, M. (2011) Homotopy Perturbation-Reproducing Kernel Method for Nonlinear Systems of Second Order Boundary Value Problems. Journal of Computational and Applied Mathematics, 235, 2405-2411. https://doi.org/10.1016/j.cam.2010.10.040
[12]	Biala, T.A. and Jator, S.N. (2017) A Family of Boundary Value Methods for Systems of Second-Order Boundary Value Problems. International Journal of Differential Equations, 2017, Article ID: 2464759. https://doi.org/10.1155/2017/2464759
[13]	Abo-Hammour, Z., Abu Arqub, O., Alsmadi, O., Momani, S. and Alsaedi, A. (2014) An Optimization Algorithm for Solving Systems of Singular Boundary Value Problems. Applied Mathematics & Information Sciences, 8, 2809-2821. https://doi.org/10.12785/amis/080617
[14]	Arqub, O.A. and Abo-Hammour, Z. (2014) Numerical Solution of Systems of Second-Order Boundary Value Problems Using Continuous Genetic Algorithm. Information Sciences, 279, 396-415. https://doi.org/10.1016/j.ins.2014.03.128
[15]	Lima, S.A., Islam, M.S., Ali, H. and Kamrujjaman, M. (2023) Numerical Method to Solve Generalized Nonlinear System of Second Order Boundary Value Problems: Galerkin Approach. Advances in the Theory of Nonlinear Analysis and Its Application, 7, 280-291. https://doi.org/10.31197/atnaa.1141150
[16]	Hao, T., Cong, F. and Shang, Y. (2018) An Efficient Method for Solving Coupled Lane-Emden Boundary Value Problems in Catalytic Diffusion Reactions and Error Estimate. Journal of Mathematical Chemistry, 56, 2691-2706. https://doi.org/10.1007/s10910-018-0912-7
[17]	Singh, R., Singh, G. and Singh, M. (2021) Numerical Algorithm for Solution of the System of Emden-Fowler Type Equations. International Journal of Applied and Computational Mathematics, 7, Article No. 136. https://doi.org/10.1007/s40819-021-01066-7
[18]	Ala'yed, O., Saadeh, R. and Qazza, A. (2023) Numerical Solution for the System of Lane-Emden Type Equations Using Cubic B-Spline Method Arising in Engineering. AIMS Mathematics, 8, 14747-14766. https://doi.org/10.3934/math.2023754
[19]	Iyase, S.A. (2018) On a Third-Order Boundary Value Problem at Resonance on the Half-Line. Arabian Journal of Mathematics, 8, 43-53. https://doi.org/10.1007/s40065-018-0209-5
[20]	Talay Akyildiz, F., Bellout, H., Vajravelu, K. and Van Gorder, R.A. (2011) Existence Results for Third Order Nonlinear Boundary Value Problems Arising in Nano Boundary Layer Fluid Flows over Stretching Surfaces. Nonlinear Analysis: Real World Applications, 12, 2919-2930. https://doi.org/10.1016/j.nonrwa.2011.02.017
[21]	Ezzati, R. and Aqhamohamadi, M. (2009) Solving the Non-Linear System of Third-Order Boundary Value Problems by Using He’s Homotopy Perturbation Method. International Journal of Industrial Mathematics, 1, 351-363.
[22]	See, P.P., Majid, Z.A., Suleiman, M., Ismail, F.B. and Othman, K.I. (2015) Solving Nonlinear System of Third-Order Boundary Value Problems Using Block Method. AIP Conference Proceedings, Penang, 28-30 May 2014, 1-7. https://doi.org/10.1063/1.4915691
[23]	Hossain, M.B. and Islam, M.S. (2014) Numerical Solutions of Sixth Order Linear and Nonlinear Boundary Value Problems. Journal of Advances in Mathematics, 7, 1180-1190.
[24]	Agarwal, R.P. (1986) Boundary Value Problems from Higher Order Differential Equations. World Scientific, Singapore, 310. https://doi.org/10.1142/0266
[25]	Wazwaz, A. (2001) The Numerical Solution of Sixth-Order Boundary Value Problems by the Modified Decomposition Method. Applied Mathematics and Computation, 118, 311-325. https://doi.org/10.1016/s0096-3003(99)00224-6
[26]	Khalid, A., Naeem, M.N., Agarwal, P., Ghaffar, A., Ullah, Z. and Jain, S. (2019) Numerical Approximation for the Solution of Linear Sixth Order Boundary Value Problems by Cubic B-Spline. Advances in Difference Equations, 2019, Article No. 492. https://doi.org/10.1186/s13662-019-2385-9
[27]	Khalid, A.A. and Nawaz Naeem, M. (2018) Cubic B-Spline Solution of Nonlinear Sixth Order Boundary Value Problems. Punjab University Journal of Mathematics, 50, 91-103.
[28]	Khalid, A., Ghaffar, A., Naeem, M.N., Nisar, K.S. and Baleanu, D. (2021) Solutions of BVPs Arising in Hydrodynamic and Magnetohydro-Dynamic Stability Theory Using Polynomial and Non-Polynomial Splines. Alexandria Engineering Journal, 60, 941-953. https://doi.org/10.1016/j.aej.2020.10.022
[29]	Noor, M.A. and Mohyud-Din, S.T. (2008) Homotopy Perturbation Method for Solving Sixth-Order Boundary Value Problems. Computers & Mathematics with Applications, 55, 2953-2972. https://doi.org/10.1016/j.camwa.2007.11.026