A new algorithm based on compressed Legendre polynomials for solving boundary value problems

1.   Introduction

Boundary value problems(BVPs) of differential equations originate from many scientific fields such as applied mathematics, physics, chemistry, biology, medical science, economics, engineering science, etc. ^[1,2,3,4]. BVPs have always been the frontier of scientific research due to its important applications and scientific implications. It is difficult to obtain analytic solutions of these equations, which promotes the research of approximate methods to get numerical solutions. For instance, the authors in ^[5] develop a numerical method for second-order three-point BVPs by using the idea of piecewise approximation. A computational approach for multi-point BVPs is discussed based on the least squares approximation method and the Lagrange-multiplier method in ^[6]. The authors in ^[7] investigate the collocation method with linear/linear rational splines for the numerical solution of two-point boundary value problems. Numerical approaches based on B-spline method ^[8,9,10,11] and reproducing kernel techniques ^[12,13,14] are proposed respectively for solving singular linear and nonlinear equations. In ^[15] and ^[16], the authors applied Legendre polynomials to handle numerical solutions. Recently, multiscale algorithms with different multiscale orthogonal basis are introduced in ^[17,18].

Numerical methods for solving linear equations have been studied extensively. This paper is concerned with the the following BVPs:

where $a_i(x), i = 1, \cdots, n, g(x)$ are sufficiently smooth functions and $\mathcal{B}_iy = \alpha_i, i = 1, \cdots, n$ are linear boundary conditions. In this work, we will give a new approach to Eq (1.1) based on compressed Legendre polynomials. The method can reduce computational cost and provide highly accurate approximate solutions. What's more, this method can avoid Runge phenomenon caused by high-order polynomial approximation. Without loss of generality, we will discuss a class of equations with boundary conditions as in Eq (1.2) to show our method conveniently. The whole analysis is also applicable to higher order equations or equations with more complex boundary conditions by proper modifications.

This paper is organized as follows. In section 2, we introduce an orthonormal basis generated from compressed Legendre polynomials. The numerical algorithm is established in section 3, and the convergence analysis is given in section 4. In section 5, we give some numerical examples to testify the effectiveness of our method.

2.   A new basis based on compressed Legendre polynomials

For convenience, we homogenize the boundary conditions in Eq (1.2) so to get Eq (2.1).

We define

$W_2^2[0, 1] = \{ {u(x)|u' \; \mbox{is absolutely continuous on $[0, 1]$, }u''\in L^2[0, 1], u(0) = u(1) = 0} \}$.

It is a reproducing kernel space as introduced in ^[19].

The inner product in $W_2^2[0, 1]$ is given by

As we all know that Legendre polynomials $L_n(x)$ is an orthogonal basis in $L^2[-1, 1]$. Therefore,

is an orthonormal basis of $L^2[0, 1]$.

By using Eq (2.2), we will construct an orthonormal basis of $W_2^2[0, 1]$. First, we compress $P_n(x)$ to obtain a new basis in $L^2[0, 1]$.

Lemma 2.1. Given $k, p\in\mathbb Z^+$, with $1\leq k\leq p$, define

where $p$ is the compression coefficient, then $\{ {\varphi _{nk}}\}_{n = 0, k = 1}^{\infty\quad p}$ is an orthonormal basis of $L^2[0, 1]$.

Proof. $\forall {\varphi _{nk}}, {\varphi _{mj}}, n \ne m$,

i) $k \ne j$, ${\varphi _{nk}}, {\varphi _{mj}}$ are orthogonal according to the definition;

ii) $k{\rm{ = }}j$, $\int_{\rm{0}}^{\rm{1}} {{\varphi _{nk}}{\varphi _{mj}}} dx = \int_{\rm{0}}^{\rm{1}} {{\varphi _{nk}}{\varphi _{mk}}} dx = \int_{\frac{{k - 1}}{p}}^{\frac{k}{p}} {p{P_n}(px - k + 1){P_m}(px - k + 1)} dx$.

Let $s = px - k + 1$,

${{\varphi _{nk}(x)}}$ maintain the orthogonality.

And

It follows that ${{\varphi _{nk}(x)}}$ are orthonormal.

In addition,

$\forall u \in {L^2}[0, 1],$ if $\int_{\rm{0}}^{\rm{1}} {u{\varphi _{nk}}} dx = 0$, that is

then $u = 0$.

So, ${{\varphi _{nk}(x)}}$ are orthonormal and complete. Therefore $\{ {\varphi _{nk}}\}_{n = 0, k = 1}^{\infty\quad p}$ is an orthonormal basis of ${L^2}[0, 1]$.

Next, we construct an orthonormal basis of $W_2^2[0, 1]$.

Define ${J^2}{\varphi _{nk\; }}$ as

Note that ${J^2}{\varphi _{nk\; }}(0) = {J^2}{\varphi _{nk\; }}(1) = 0$. We will prove that $\{ {J^2}{\varphi _{nk}}\}_{n = 0, k = 1}^{\infty\quad p}$ is a basis of $W_2^2[0, 1]$.

Theorem 2.1. $\{ {J^2}{\varphi _{nk}}(x)\}_{n = 0, k = 1}^{\infty\quad p}$ is an orthonormal basis of $W_2^2[0, 1]$.

Proof. ${{J^2}{\varphi _{nk\; }}(x)}$ are orthonormal. In fact,

$\forall u \in W_2^2[0, 1],$ if $\int_{\rm{0}}^{\rm{1}} {u''({J^2}{\varphi _{nk}})''} dx = 0$, then $\int_{\rm{0}}^{\rm{1}} {u''{\varphi _{nk}}} dx = 0$, which implies that $u'' = 0$.

Because $u'$ is absolutely continuous on $[0, 1]$, we have $u' = C$, $C$ is a constant.

Notice that $u(0) = u(1) = 0$, so we have $u = 0$ on $[0, 1]$ which shows the completeness.

Therefore, $\left\{ {{J^2}{\varphi _{nk\; }}(x)} \right\}$ is an orthonormal basis of $W_2^2[0, 1]$.

3.   Numerical method

In this section, we introduce a numerical algorithm for solving Eq (2.1).

Define a linear operator

by

Eq (2.1) can be transformed into the following operator form

It can be proved that the operator is bounded. Next, we demonstrate the implementation steps of our numerical algorithm for solving Eq (3.1). The theoretical background and the error analysis will be discussed in the upcoming section.

Definition 3.1. For any $\varepsilon > 0$, $u$ is called an $\varepsilon$-approximate solution if ${\left\| {\mathcal{L}u - f} \right\|_{L^2[0, 1]}^2} < \varepsilon$.

Let $F({c_{01}}, \cdots \cdots, {c_{np}}) = {\left\| {\sum\limits_{j = 0}^n {\sum\limits_{k = 1}^p {{c_{jk}}\mathcal{L}{J^2}{\varphi _{jk}}} } - f} \right\|_{L^2[0, 1]}^2}$, and $(c_{01}^{\rm{*}}, \cdots \cdots, c_{np}^*)$ is the minimum value point of $F({c_{01}}, \cdots \cdots, {c_{np}})$.

Theorem 3.1. For any $\varepsilon > 0$, $\exists N$, when $n > N$, ${u_{np}} = \sum\limits_{j = {\rm{0}}}^n {\sum\limits_{k = 1}^p {c_{jk}^*{J^2}{\varphi _{jk}}} }$ is an $\varepsilon$-approximate solution of Eq (3.1).

Proof. Assume that $u$ is the exact solution of Eq (3.1), $u = \sum\limits_{j = {\rm{0}}}^\infty {\sum\limits_{k = 1}^p {{a_{jk}}{J^2}{\varphi _{jk}}} }$.

$\forall \varepsilon > 0$, $\exists N$, when $n > N$, we have $w_{np} = \sum\limits_{j = {\rm{0}}}^n {\sum\limits_{k = 1}^p {{a_{jk}}{J^2}{\varphi _{jk}}} }$ satisfies

Because

we get

According to Theorem 3.1, we obtain an $\varepsilon$-approximate solution ${u_{np}}$ of Eq (3.1) by using the new basis. The coefficient is the minimum value point of $F$. Now, we will prove the minimum value point of $F$ is unique.

Take partial derivatives of the function $F$ with respect to ${c_{ml}}, m = 0 \cdots n, l = 1 \cdots p$.

Let $\frac{{\partial F}}{{\partial {c_{ml}}}} = 0$, then

Define a matrix $\textbf{D} = {\left({ \langle \mathcal{L}{J^2}{\varphi _{jk}}, \mathcal{L}{J^2}{\varphi _{ml}} \rangle } \right)_{(n + 1)p \times (n + 1)p}}$, and a vector $\textbf{b} = (\langle \mathcal{L}{J^2}{\varphi _{ml}}, f \rangle)_{(n + 1)p}^T$. Then $(c_{01}^{\rm{*}}, \cdots \cdots c_{np}^*)$ is the solution of $\textbf{Dc} = \textbf{b}$, where $\textbf{c} = (c_{01}, \cdots \cdots c_{np})^T$.

Theorem 3.2. Suppose that $\mathcal{L}$ is invertible, then the solution of (3.2) is unique.

Proof. The homogeneous equation corresponding to (3.2) is

Let's multiply both sides of the above equations by ${c_{ml}}, m = 0 \cdots n, l = 1 \cdots p$, and add them all. We have

which is just ${\left\| {\sum\limits_{j = 0}^n {\sum\limits_{k = 1}^p {{c_{jk}}\mathcal{L}{J^2}{\varphi _{jk}}} } } \right\|_{L^2[0, 1]}^2} = 0$. Then $\sum\limits_{j = 0}^n {\sum\limits_{k = 1}^p {{c_{jk}}\mathcal{L}{J^2}{\varphi _{jk}}} } = 0$.

$\mathcal{L}$ is invertible and $\{ {J^2}{\varphi _{nk}}(x)\}_{n = 0, k = 1}^{\infty\quad p}$ is an orthonormal basis, so

Hence, Eq (3.2) has a unique solution.

4.   Convergence and stability analysis

In this section, we will prove the convergence and stability of our method.

Let $u(x) = \sum\limits_{j = {\rm{0}}}^\infty {\sum\limits_{k = 1}^p {{a_{jk}}{J^2}{\varphi _{jk}}} }$ is the exact solution of Eq (3.1), ${u_{np}}(x) = \sum\limits_{j = {\rm{0}}}^n {\sum\limits_{k = 1}^p {c_{jk}^*{J^2}{\varphi _{jk}}} }$ is the $\varepsilon$-approximate solution as constructed in the previous section.

Theorem 4.1. ${u_{np}}(x)$ uniformly converges to $u(x)$.

Proof. Theorem 3.1 implies that $\left\| {\mathcal{L}u_{np} - f} \right\|_{L^2[0, 1]}^2 \to {\rm{0}}.$ We obtain that

So,

Furthermore,

where $R_x(t)$ is the reproducing kernel of $W_2^2[0, 1]$.

Because $R_x(t)$ is bounded on $[0, 1]$, we get

We now present the error analysis of our method.

Let

By (4.1), we get

By Theorem 3.1, it follows that

So,

In the above formula,

Since $u'' = \sum\limits_{j = {\rm{0}}}^\infty {\sum\limits_{k = 1}^p {{a_{jk}}{\varphi _{jk}}} }, w''_{np} = \sum\limits_{j = 0}^n {\sum\limits_{k = 1}^p {{a_{jk}}{\varphi _{jk}}} }$, we get

Assume $v(s, k) = u''(\frac{{s + k - 1}}{p})$, ${v_n}(s, k) = {w''_{np}}(\frac{{s + k - 1}}{p})$, then

where $b_{jk}$ are the coefficients of the generalized Fourier expansion with respect to the basis $\{P_j(x)\}_{j = 1}^\infty$. So $v(s, k) = \sum\limits_{j = 0}^{\infty } {{b_{jk}}{P_j}(s)}, {v_n}(s, k) = \sum\limits_{j = 0}^n {{b_{jk}}{P_j}(s)}.$

Now Eq (4.3) becomes

According to ^[20], we can get

where $C$ is a constant and $m$ makes $u^{(m+2)}\in L^2[0, 1].$

Since $v(s, k) = u''(\frac{{s{\rm{ + }}k - 1}}{p})$, $\partial _s^lv(s, k) = \frac{1}{{{p^l}}}{u^{(l+2)}}(x)$.

Hence,

where $C_0$ is a constant.

By (4.3) and (4.4), we have that

Plug it into (4.2), we conclude that

where $M_0$ is a constant.

Inequality (4.5) gives the error bound of our new approach. We summarize our result in the following theorem.

Theorem 4.2. If ${u^{(m + 2)}} \in {L^2}[0, 1]$ with $m\geq 0$, and $n\geq m-1$, then there exists a positive constant $M_0$ so that $\left\| {u_{np}-u} \right\|_{{W_2^2}[0, 1]}\le {\sqrt{M_0}}\frac{1}{n^{m}} {\frac{1}{p^{m}}}$.

Theorem 4.2 indicates the convergence of our algorithm is influenced by both of $n$ and $p$. The rate of convergence gets faster with the increasing of the two parameters.

In the end of this section, the stability of our algorithm is discussed by using the condition number of the matrix $\bf{D}$.

Let $\left\{ {{\psi _i}} \right\}_{i = 1}^\infty = {\rm{\{ }}{J^2}{\varphi _{01}}, \cdots, {J^2}{\varphi _{0p}}, \cdots, {J^2}{\varphi _{n1}}, \cdots, {J^2}{\varphi _{np}}, \cdots \}$, then $\bf{D} = {\left({ \langle \mathcal{L}{J^2}{\varphi _{jk}}, \mathcal{L}{J^2}{\varphi _{ml}} \rangle } \right)_{M \times M}}$ could be rewritten as $\bf{D} = {\left({ \langle \mathcal{L}{\psi _i}, \mathcal{L}{\psi _j} \rangle } \right)_{M \times M}}$, $M = (n + 1)p$.

Lemma 4.1 (^[17]) If $u \in W_2^2[0, 1]$ with $\left\| u\right\|_{{W_2^2[0, 1]}} = 1$, then ${\left\| {\mathcal{L}u} \right\|_{{L^2[0, 1]}}} \ge \frac{{\rm{1}}}{{\left\| {{\mathcal{L}^{ - 1}}} \right\|}}$.

Theorem 4.3. The condition number of $\bf{D}$ satisfies $Cond(\bf{D})_2 \le {\left\| \mathcal{L} \right\|^2}\left\| {{\mathcal{L}^{ - 1}}} \right\|^2$.

Proof. Assume $\bf{x} = {({x_1}, \cdots {x_M})^T}$ is a unit eigenvector belonging to $\lambda$, i.e., $\lambda \bf{x} = \bf{D}\bf{x}$, $\left\| \bf{x} \right\| = 1$.

We have

Multiplying Eq (4.6) by ${x_i}, i = 1 \cdots M$, and add all the equations.

On the other hand,

where $u{\rm{ = }}\sum\limits_{i = 1}^M {{x_i}{\psi _i}}$.

From Lemma 4.1, $\lambda \ge \frac{{\rm{1}}}{{\left\| {{\mathcal{L}^{ - 1}}} \right\|^2}}$.

So

The stability is proved.

5.   Numerical examples

In this part, we will test four examples of BVPs which have been discussed by different algorithms in ^{[17,21,22,23,24]}. In the following examples, comparison of numerical results demonstrate the efficiency and stability of our method. The absolute error function is defined as $E = \left|{ {y_{np}}(x)-y(x)} \right|, 0 \le x \le 1$, where $y(x)$ is the exact solution and ${y_{np}}(x)$ is an $\varepsilon$-approximate solution.

Example 1. (^[17,21]) Consider an equation with Robin boundary condition:

$y = x^2$ is the exact solution of Eq (5.1). The numerical results using our method are compared with ^[17,21] in Table 1. It shows that our method converges rapidly with higher accuracy.

Example 2. Consider the following two-point BVP:

where $f(x) = \left\{ \begin{array}{l} x\cos x, {\kern 1pt} \quad \qquad \qquad 0 \le x < \frac{1}{2}, \\ \left({ - x/2 + 1} \right)\sinh x, \quad \frac{1}{2} \le x \le 1. \end{array} \right.$. Since the exact solution of Eq (5.2) is not known, we discuss the absolute residual error function defined as $R = |{\mathcal{L}u_{np} - f}|$. The results are listed in Table 2. It shows that the numerical error is rather smaller by our method.

Example 3. (^[22,23] We consider the following sixth-order BVP:

$y = (1-x){e^x}$ is the exact solution of Eq (5.3). The results are reported in the Table 3. The absolute error of our method is smaller than those obtained in ^[22] by homotopy perturbation method and in ^[23] by Legendre wavelet collocation method.

Example 4. (^[24]) Consider the following nonlinear fourth-order BVP:

where $f(x) = 1+\sin(1+\sinh(x))-(e^x-2)\sinh(x)$. The exact solution of Eq (5.4) is $y(x) = 1+\sinh(x)$. Firstly, Eq (5.4) is linearized by Newton iterative method ^[13]. Then, we use our method to solve the linear equation. The numerical results obtained after three iterations are given in Table 4. It is seen from Table 4 that we have got a better approximation to the exact solution of the nonlinear problem.

6.   Conclusions

In this paper, a new basis $\{ {J^2}{\varphi _{nk}}(x)\}_{n = 0, k = 1}^{\infty\quad p}$ of $W_2^2[0, 1]$ is constructed from compressed Legendre polynomials. This basis is simple, efficient, and practical. Our method can avoid Runge phenomenon caused by high-order polynomial approximation. Using the concept of $\varepsilon$-approximate solution, we describe a method with the new basis to solve boundary value problems. The convergence and stability of our algorithm are proved. Our algorithm is tested by some numerical examples with different boundary conditions. The numerical results show that our method have higher accuracy for solving linear and nonlinear problems.

Acknowledgments

This work has been supported by characteristic innovation project of Guangdong Province (2019KTSCX217) and Big data research center of Zhuhai Campus, Beijing Institute of Technology (XJ-2018-05).

Conflict of interest

The authors declare that they have no competing interests.