Solution of a polylocal problem with a pseudospectral method

Consider the problem (PVP):

where $a,b,\alpha ,\beta \in \mathbb{R}$$a,b,\alpha ,\beta \in \mathbb{R}$a,b,alpha,beta inRa, b, \alpha, \beta \in \mathbb{R}$a,b,\alpha ,\beta \in \mathbb{R}$. This is not a two-point boundary value problem, since $a,b\in (0,1)$$a,b\in (0,1)$a,b in(0,1)a, b \in(0,1)$a,b\in (0,1)$.

We try to solve the problem using a pseudospectral collocation method with Chebyshev extrema combined with a Runge Kutta method.

Then, we compare them in terms of error and cost.
Our choice to use this method is based on the following reasons :

As drawbacks, we mention:

We also consider the BVP :

To apply the collocation theory we need to have an isolated solution $y(x)$$y(x)$y(x)y(x)$y(x)$ of the problem $(1.4)+(1.5)+(1.6)$$(1.4)+(1.5)+(1.6)$(1.4)+(1.5)+(1.6)(1.4)+(1.5)+(1.6)$(1.4)+(1.5)+(1.6)$, and this occurs if the above linearized problem for $y(x)$$y(x)$y(x)y(x)$y(x)$ is uniquely solvable. R.D Russel and L.F.Shampine [3] study the existence and the uniqueness of the isolated solution.

Theorem 1.1 [3] Suppose that $y(x)$$y(x)$y(x)y(x)$y(x)$ is a solution of the boundary value problem $(1.4)+(1.5)+(1.6)$$(1.4)+(1.5)+(1.6)$(1.4)+(1.5)+(1.6)(1.4)+(1.5)+(1.6)$(1.4)+(1.5)+(1.6)$, that the functions :

are defined and continuous for $a\le x\le b$$a\le x\le b$a <= x <= ba \leq x \leq b$a\le x\le b$, and $|y_k-y^{(k)}(x)|\le \delta ,k=0,1$$\left|y_k-y^{(k)}(x)\right|\le \delta ,k=0,1$|y_(k)-y^((k))(x)| <= delta,k=0,1\left|y_{k}-y^{(k)}(x)\right| \leq \delta, k=0,1$|y_k-y^{(k)}(x)|\le \delta ,k=0,1$; $\delta >0$$\delta >0$delta > 0\delta>0$\delta >0$, and the homogeneous equations $y^{\prime \prime }(x)=0$$y^{\prime \prime }(x)=0$y^('')(x)=0y^{\prime \prime}(x)=0$y^{\prime \prime }(x)=0$ subject to the homogeneous boundary conditions (1.5)+(1.6) has only the trivial solution. If the linear homogeneous equations:

has only trivial solution, then this is sufficient to guarantee that there exists a $\sigma >0$$\sigma >0$sigma > 0\sigma>0$\sigma >0$ such $y(x)$$y(x)$y(x)y(x)$y(x)$ is the unique solution of problem BVP in the sphere:

For the existence and uniqueness of an IVP, we recall the following:
Theorem 1.2 [4, pp. 112-113] Suppose that $D=\{a\le x\le b,-\mathrm{\infty }<y<\mathrm{\infty }\}$$D=\{a\le x\le b,-\mathrm{\infty }<y<\mathrm{\infty }\}$D={a <= x <= b,-oo < y < oo}D=\{a \leq x \leq b,-\infty<y< \infty\}$D=\{a\le x\le b,-\mathrm{\infty }<y<\mathrm{\infty }\}$ and $f(x,y)$$f(x,y)$f(x,y)f(x, y)$f(x,y)$ is continuous on $D$$D$DD$D$. If $f$$f$ff$f$ satisfies a Lipschitz condition on $D$$D$DD$D$ in the variable $y$$y$yy$y$, then the initial value problem (IVP)

has a unique solution $y(x)$$y(x)$y(x)y(x)$y(x)$ for $a\le x\le b$$a\le x\le b$a <= x <= ba \leq x \leq b$a\le x\le b$.
If the problem BVP has the unique solution, the requirement $y(x)\in C^2[0,1]$$y(x)\in C^2[0,1]$y(x)inC^(2)[0,1]y(x) \in C^{2}[0,1]$y(x)\in C^2[0,1]$ ensure the existence and the uniqueness of the solution of PVP.

Our second method consists into decomposition of the the problem $(1.1)+(1.2)+(1.3)$$(1.1)+(1.2)+(1.3)$(1.1)+(1.2)+(1.3)(1.1)+(1.2)+(1.3)$(1.1)+(1.2)+(1.3)$ into three problems:

Also we suppose that the problem (1.4)+(1.5)+(1.6) satisfy Theorem 1.1, which assures a sufficient condition to guarantee that there exists a $\sigma >0$$\sigma >0$sigma > 0\sigma>0$\sigma >0$ such $y(x)$$y(x)$y(x)y(x)$y(x)$ is the unique solution of problem BVP in the sphere:

For $0\le x\le a$$0\le x\le a$0 <= x <= a0 \leq x \leq a$0\le x\le a$ and $b\le x\le 1$$b\le x\le 1$b <= x <= 1b \leq x \leq 1$b\le x\le 1$ we have two initial value problems.
Due to conditions in Theorems 1.1 and 1.2, the problem (1.1)+ $(1.2)+(1.3)$$(1.2)+(1.3)$(1.2)+(1.3)(1.2)+(1.3)$(1.2)+(1.3)$ has a unique solution.

Consider the grid

We shall use a pseudospectral method for the solution of $(1.4)+(1.5)+$$(1.4)+(1.5)+$(1.4)+(1.5)+(1.4)+(1.5)+$(1.4)+(1.5)+$ (1.6) and a Runge-Kutta method for the two IVPs.

The solution of the two-point boundary value problems follow the ideas in [6] Let $y(x)$$y(x)$y(x)y(x)$y(x)$ be the solution of (1.4)+(1.5)+(1.6). Considering the Lagrange basis $\left({\ell }_k\right)$$\left({\ell }_k\right)$(ℓ_(k))\left(\ell_{k}\right)$\left({\ell }_k\right)$, we have

where

is the remainder of Lagrange interpolation and

Since $y(x)$$y(x)$y(x)y(x)$y(x)$ fulfills the differential equation (1.4), we obtain

Setting $y\left(x_k\right)=y_k$$y\left(x_k\right)=y_k$y(x_(k))=y_(k)y\left(x_{k}\right)=y_{k}$y\left(x_k\right)=y_k$ and ignoring the rest, one obtains the nonlinear system

with unknowns $y_k,k=1,\dots ,N-1$$y_k,k=1,\dots ,N-1$y_(k),k=1,dots,N-1y_{k}, k=1, \ldots, N-1$y_k,k=1,\dots ,N-1$; here $y_0=y\left(x_0\right)=\alpha$$y_0=y\left(x_0\right)=\alpha$y_(0)=y(x_(0))=alphay_{0}=y\left(x_{0}\right)=\alpha$y_0=y\left(x_0\right)=\alpha$ and $y_N=y\left(x_N\right)=\beta$$y_N=y\left(x_N\right)=\beta$y_(N)=y(x_(N))=betay_{N}= y\left(x_{N}\right)=\beta$y_N=y\left(x_N\right)=\beta$. The approximate solution (that is the collocation polynomial for $(1.4)+(1.5)+(1.6))$$(1.4)+(1.5)+(1.6))$(1.4)+(1.5)+(1.6))(1.4)+(1.5)+(1.6))$(1.4)+(1.5)+(1.6))$ is the Lagrange interpolation polynomial at nodes $\left\{x_k\right\}$$\left\{x_k\right\}${x_(k)}\left\{x_{k}\right\}$\left\{x_k\right\}$

The nonlinear system (2.4) can be rewritten as

where

If the nodes $\left(x_k\right),k=0,\dots ,N$$\left(x_k\right),k=0,\dots ,N$(x_(k)),k=0,dots,N\left(x_{k}\right), k=0, \ldots, N$\left(x_k\right),k=0,\dots ,N$, are symmetric with respect of $(c+d)/2$$(c+d)/2$(c+d)//2(c+d) / 2$(c+d)/2$, then $A$$A$AA$A$ is centrosymmetric (see [6] for the proof), so nonsingular. We recall the definition [7]: an $m\times m$$m\times m$m xx mm \times m$m\times m$ matrix $A$$A$AA$A$ is centrosymmetric if

Examples of such nodes are given by

i.e. the equispaced and the Chebyshev-Lobatto nodes.
We introduce

To solve numerically (1.1) $+(1.2)+(1.3)$$+(1.2)+(1.3)$+(1.2)+(1.3)+(1.2)+(1.3)$+(1.2)+(1.3)$ on $\mathrm{\Delta }$$\mathrm{\Delta }$Delta\Delta$\mathrm{\Delta }$ given by (2.1), we apply pseodospectral collocation at points in $[c,d]$$[c,d]$[c,d][c, d]$[c,d]$ and then a Runge-Kutta to the other points. To apply the Runge-Kutta method for the solution of two IVP on $[a,c]$$[a,c]$[a,c][a, c]$[a,c]$ and $[d,b]$$[d,b]$[d,b][d, b]$[d,b]$ we need the derivatives $y^{\prime }(c)$$y^{\prime }(c)$y^(')(c)y^{\prime}(c)$y^{\prime }(c)$ and $y^{\prime }(d)$$y^{\prime }(d)$y^(')(d)y^{\prime}(d)$y^{\prime }(d)$. This can be computed by deriving the formula (2.5).

In [5] the authors prove the following theorems.

Theorem 2.1 If $f$$f$ff$f$ verifies a Lipschitz condition with respect to the second variable

and $\Vert A^{-1}\Vert L<1$$‖A^{-1}‖L<1$||A^(-1)||L < 1\left\|A^{-1}\right\| L<1$\Vert A^{-1}\Vert L<1$, then the system (2.4) has a unique solution which can be calculated by the successive approximation method

with $Y^{(0)}$$Y^{(0)}$Y^((0))Y^{(0)}$Y^{(0)}$ fixed and $G$$G$GG$G$ given by (2.6).
Theorem 2.2 If $Y={[y\left(x_1\right),\dots ,y\left(x_{N-1}\right)]}^T$$Y={\left[y\left(x_1\right),\dots ,y\left(x_{N-1}\right)\right]}^T$Y=[y(x_(1)),dots,y(x_(N-1))]^(T)Y=\left[y\left(x_{1}\right), \ldots, y\left(x_{N-1}\right)\right]^{T}$Y={[y\left(x_1\right),\dots ,y\left(x_{N-1}\right)]}^T$ where $y(x)$$y(x)$y(x)y(x)$y(x)$ is the solution of problem $(1.4)+(1.5)+(1.6);Y_N={[y_1,\dots ,y_{N-1}]}^T$$(1.4)+(1.5)+(1.6);Y_N={\left[y_1,\dots ,y_{N-1}\right]}^T$(1.4)+(1.5)+(1.6);Y_(N)=[y_(1),dots,y_(N-1)]^(T)(1.4)+(1.5)+(1.6) ; Y_{N}=\left[y_{1}, \ldots, y_{N-1}\right]^{T}$(1.4)+(1.5)+(1.6);Y_N={[y_1,\dots ,y_{N-1}]}^T$ where $y_i$$y_i$y_(i)y_{i}$y_i$ are the values of approximated solution at ${\dot{x}}_i$${\dot{x}}_i$x^(˙)_(i)\dot{x}_{i}${\dot{x}}_i$ computed by (2.4) and $R=$$R=$R=R=$R=$
${[-{\left(R_{N}y\right)}^{\prime \prime }\left(x_1\right),\dots ,-{\left(R_{N}y\right)}^{\prime \prime }\left(x_{N-1}\right)]}^T$${\left[-{\left(R_{N}y\right)}^{\prime \prime }\left(x_1\right),\dots ,-{\left(R_{N}y\right)}^{\prime \prime }\left(x_{N-1}\right)\right]}^T$[-(R_(N)y)^('')(x_(1)),dots,-(R_(N)y)^('')(x_(N-1))]^(T)\left[-\left(R_{N} y\right)^{\prime \prime}\left(x_{1}\right), \ldots,-\left(R_{N} y\right)^{\prime \prime}\left(x_{N-1}\right)\right]^{T}${[-{\left(R_{N}y\right)}^{\prime \prime }\left(x_1\right),\dots ,-{\left(R_{N}y\right)}^{\prime \prime }\left(x_{N-1}\right)]}^T$ then for the error $\Vert Y-Y_N\Vert$$‖Y-Y_N‖$||Y-Y_(N)||\left\|Y-Y_{N}\right\|$\Vert Y-Y_N\Vert$ it holds

Combining the results of Theorems 2.1 and 2.2 with the stability and convergence of Runge-Kutta methods we have:

Theorem 2.3 If

and

then for each point $x_i$$x_i$x_(i)x_{i}$x_i$ in $\mathrm{\Delta },i=-q,\dots ,N+p,|y\left(x_i\right)-y_i|=O\left(h^k\right)$$\mathrm{\Delta },i=-q,\dots ,N+p,\left|y\left(x_i\right)-y_i\right|=O\left(h^k\right)$Delta,i=-q,dots,N+p,|y(x_(i))-y_(i)|=O(h^(k))\Delta, i=-q, \ldots, N+p,\left|y\left(x_{i}\right)-y_{i}\right|=O\left(h^{k}\right)$\mathrm{\Delta },i=-q,\dots ,N+p,|y\left(x_i\right)-y_i|=O\left(h^k\right)$.
Proof. The condition $\Vert A^{-1}\Vert L<1$$‖A^{-1}‖L<1$||A^(-1)||L < 1\left\|A^{-1}\right\| L<1$\Vert A^{-1}\Vert L<1$ assures us that $G$$G$GG$G$ is a contraction. From Banach's fixpoint theorem it follows that $\left(Y^{(n)}\right)$$\left(Y^{(n)}\right)$(Y^((n)))\left(Y^{(n)}\right)$\left(Y^{(n)}\right)$ given by (2.7) is convergent to the exact solution $\overline{Y}$$\overline{Y}$bar(Y)\bar{Y}$\overline{Y}$ of (2.4) and the following estimation holds

If the accuracy of the collocation method for the BVP is $O\left(h^k\right)$$O\left(h^k\right)$O(h^(k))O\left(h^{k}\right)$O\left(h^k\right)$ (that is, the approximate solution $y$$y$yy$y$ and its derivative $y^{\prime }$$y^{\prime }$y^(')y^{\prime}$y^{\prime }$ at $c$$c$cc$c$ and $d$$d$dd$d$ are within this accuracy limit), and if the Runge-Kutta method is stable and has the order $k$$k$kk$k$, then the final solution has the same accuracy. The stability and convergence of Runge-Kutta method are guaranteed by Theorems 5.3.1, page 285 and 5.3.2, page 288 in [8]. $◻$$◻$◻\square$◻$

Our combined method was implemented in MATLAB, using the functions cebdif, cebint and cebdifft contained in dmsuite and described
in [1]. We chose $\{x_k,k=0,N\}$$\left\{x_k,k=0,N\right\}${x_(k),k=0,N}\left\{x_{k}, k=0, N\right\}$\{x_k,k=0,N\}$ as extreme Chebysev points and the other of $\mathrm{\Delta }$$\mathrm{\Delta }$Delta\Delta$\mathrm{\Delta }$ were computed using the MATLAB solver ode45. Since the successive approximation method is slow, we solve the nonlinear system (2.4) by Newton's method. The MATLAB function solvepolylocalceb solve the nonlinear system and call the Runge-Kutta solver. The derivatives at $c$$c$cc$c$ and $d$$d$dd$d$ were computed by calling cebdifft.

Let us consider two nonlinear examples.

Example 3.1 [9, page 491]Bratu's equation for $\lambda =1$$\lambda =1$lambda=1\lambda=1$\lambda =1$ :

We took $N=128$$N=128$N=128N=128$N=128$ and the tolerance $\epsilon ={10}^{-10}$$\epsilon ={10}^{-10}$epsi=10^(-10)\varepsilon=10^{-10}$\epsilon ={10}^{-10}$. The starting value is $y^{(0)}(x)=\frac{x(1-x)}{4}$$y^{(0)}(x)=\frac{x(1-x)}{4}$y^((0))(x)=(x(1-x))/(4)y^{(0)}(x)=\frac{x(1-x)}{4}$y^{(0)}(x)=\frac{x(1-x)}{4}$. The solution of nonlinear system is obtained after 4 iterations. The graph of the solution is given in Figure 1.

Example 3.2 Consider the problem

The positive solution of this problem represents the average temperature in a reaction-diffusion process. In [10], the authors proved that, for $p=3$$p=3$p=3p=3$p=3$, the problem (3.1) has a unique positive solution. We consider here the "variant"

In this example, $N=128,\epsilon ={10}^{-8},y^{(0)}=\frac{6\pi }{\sqrt{2}}x(1-x)$$N=128,\epsilon ={10}^{-8},y^{(0)}=\frac{6\pi }{\sqrt{2}}x(1-x)$N=128,epsi=10^(-8),y^((0))=(6pi)/(sqrt2)x(1-x)N=128, \varepsilon=10^{-8}, y^{(0)}=\frac{6 \pi}{\sqrt{2}} x(1-x)$N=128,\epsilon ={10}^{-8},y^{(0)}=\frac{6\pi }{\sqrt{2}}x(1-x)$. The desired tolerance was obtained after 8 iterations. Figure 2 shows the graph of the solution.

[1] J. A. C. Weideman, S. C. Reddy, A MATLAB Differentiation Matrix Suite, ACM Transactions on Mathematical Software, 26(4), 2000, pp. 465-519.
[2] W. Heinrichs, Strong convergence estimates for pseudospectral methods, Applications of Mathematics, 37(6), 1992, pp. 401-417.
[3] R.D. Russel, L.F.Shampine, A Collocation Method for boundary Value Problems, Numer. Math. 19, pp. 1-28, Springer Verlag, 1972.
[4] Henrici, P., Discrete variable methods in ordinary differential equations, John Willey, New York, pp. 345-384, 1962.
[5] F. Costabile, A. Napoli, A Collocation Method for Global Approximation of General Second Order BVPs, Computing Letters, 3(1), 2007, pp. 23-34.
[6] F. A. Costabile, F. Dell'Accio, On the nonsingularity of a special class of centrosymmetric matrices arising in spectral methods in BVPs, Applied Mathematics and Computation, 206, 2008, pp. 991-993.
[7] A. L. Andrew, Centrosymmetric Matrices, SIAM review, 40(3), 1998, pp. 697-698.
[8] W. Gautschi, Numerical Analysis. An Introduction, Birkhäuser, Boston, Basel, Berlin, 1997.
[9] U. Ascher, R.M. Mattheij, and R.D. Russel, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, SIAM, 1997.
[10] C. I. Gheorghiu, D. Trif, The Numerical Approximation to Positive Solution for some Reaction-Diffusion Problems, Pu.M.A., 11, 2000, pp. 243-253.