Radial basis function methods for the Rosenau equation and other higher order PDEs

TL;DR

This paper demonstrates the effectiveness of radial basis function (RBF) meshfree methods, including fictitious point and resampling approaches, for solving the Rosenau equation and other higher-order PDEs in complex geometries with high accuracy.

Contribution

It extends RBF collocation methods with fictitious point and resampling techniques to efficiently solve the Rosenau equation in multiple dimensions and geometries, with theoretical and numerical accuracy analysis.

Findings

High order convergence achieved for Rosenau equation

RBF methods perform comparably or better than pseudospectral methods

Effective in irregular geometries like starfish-shaped domains

Abstract

Meshfree methods based on radial basis function (RBF) approximation are of interest for numerical solution of partial differential equations (PDEs) because they are flexible with respect to the geometry of the computational domain, they can provide high order convergence, they are not more complicated for problems with many space dimensions and they allow for local refinement. The aim of this paper is to show that the solution of the Rosenau equation, as an example of an initial-boundary value problem with multiple boundary conditions, can be implemented using RBF approximation methods. We extend the fictitious point method and the resampling method to work in combination with an RBF collocation method. Both approaches are implemented in one and two space dimensions. The accuracy of the RBF fictitious point method is analysed partly theoretically and partly numerically. The error estimates indicate that a high order of convergence can be achieved for the Rosenau equation. The numerical experiments show that both methods perform well. In the one-dimensional case, the accuracy of the RBF approaches is compared with that of a pseudospectral resampling method, showing similar or slightly better accuracy for the RBF methods. In the two-dimensional case, the Rosenau problem is solved both in a square domain and in a starfish-shaped domain, to illustrate the capability of the RBF-based methods to handle irregular geometries.

Tables

Table 1. Table 1: Values of the different powers involved in the error estimates.

$q$	$\left(1-q/5\right)$	$q\left(1-q/5\right)$	$\tilde{q}$	$\tilde{q}\left(1-q/5\right)$
1	0.8	0.8	0.8	0.64
2	0.6	1.2	0.6	0.36
3	0.4	1.2	0.8	0.32
4	0.2	0.8	0.6	0.12
5	0	0	0	0
6	-0.2	-1.2	-0.2	0.04
7	-0.4	-2.8	-0.4	0.16
8	-0.6	-4.8	-0.6	0.36
9	-0.8	-7.2	-0.8	0.64
10	-1	-10	-1	1

Equations

$$u_t+\alpha (\mathbf{x},t){\Delta }^2{u}_t=\nabla \cdot g(u),(\mathbf{x},t)\in \Omega \times (0,T],$$

$$u(\mathbf{x},t)$$

$$\frac{\partial u}{\partial n}(\mathbf{x},t)$$

$$u(\mathbf{x},0)=f_0(\mathbf{x}),x\in \Omega .$$

$$s(x)=\sum _{j=1}^N{\lambda }_j\varphi (\Vert x-x_j\Vert ),$$

$$A\bar{\lambda }=\bar{u},$$

$$s(x)=\bar{\varphi }(x)\bar{\lambda }=\bar{\varphi }(x)A^{-1}\bar{u},$$

$$s(x)=\bar{\psi }(x)\bar{u}=\sum _{j=1}^N{\psi }_j(x)u(x_j),$$

$$\mathcal{L}s(x)=\mathcal{L}\bar{\psi }(x)\bar{u}=\sum _{j=1}^N\mathcal{L}{\psi }_j(x)u(x_j).$$

$${\Psi }_{\mathcal{L}}={\Phi }_{\mathcal{L}}{A}^{-1},$$

$$s(x,t)=\sum _{j=1}^N{\psi }_j(x)u_j(t),$$

$$u_t+\alpha (x,t)u_{xxxxt}=g_u(u)u_x,(x,t)\in [-L,L]\times (0,T],$$

$$u(\pm L,t)$$

$$u_x(\pm L,t)$$

$$u_t+\alpha (x,t)w_{xxxt}$$

$$w_t-u_{xt}$$

$$s(x,t)=\sum _{j=0}^{N+1}{\psi }_j(x)u_j(t).$$

$$\underbrace{\left(\begin{array}[]{cc}\psi_{0}^{\prime}(x_{1})&\psi_{N+1}^{\prime}(x_{1})\\ \psi_{0}^{\prime}(x_{N})&\psi_{N+1}^{\prime}(x_{N})\end{array}\right)}_{B_{f}}S_{f}+\underbrace{\left(\begin{array}[]{ccc}\psi_{2}^{\prime}(x_{1})&\cdots&\psi_{N-1}^{\prime}(x_{1})\\ \psi_{2}^{\prime}(x_{N})&\cdots&\psi_{N-1}^{\prime}(x_{N})\end{array}\right)}_{B_{d}}S_{d}+\underbrace{\left(\begin{array}[]{cc}\psi_{1}^{\prime}(x_{1})&\psi_{N}^{\prime}(x_{1})\\ \psi_{1}^{\prime}(x_{N})&\psi_{N}^{\prime}(x_{N})\\ \end{array}\right)}_{B_{b}}F_{1}(t)=F_{2}(t),$$

$$s(x,t)$$

$$+\left([{\psi }_1(x),{\psi }_N(x)]-[{\psi }_0(x),{\psi }_{N+1}(x)]B_f^{-1}{B}_b\right)F_1(t)+[{\psi }_0(x),{\psi }_{N+1}(x)]B_f^{-1}{F}_2(t).$$

$$s(x,t)=\sum _{j=2}^{N-1}{\tilde{\psi }}_j(x)u_j(t)+F(x,t),$$

$${\tilde{\psi }}_j(x)={\psi }_j(x)$$

$$+\frac{{\psi }_0^{\prime }(x_N){\psi }_j^{\prime }(x_1)-{\psi }_0^{\prime }(x_1){\psi }_j^{\prime }(x_N)}{{\psi }_{N+1}^{\prime }(x_N){\psi }_0^{\prime }(x_1)-{\psi }_{N+1}^{\prime }(x_1){\psi }_0^{\prime }(x_N)}{\psi }_{N+1}(x).$$

$${\tilde{\Psi }}_{\mathcal{L}}={\Psi }_{\mathcal{L},d}-{\Psi }_{\mathcal{L},f}{B}_f^{-1}{B}_d,$$

$$[F_{\mathcal{L}}(x_2,t),\dots ,F_{\mathcal{L}}(x_{N-1},t){]}^T=({\Psi }_{\mathcal{L},b}-{\Psi }_{\mathcal{L},f}{B}_f^{-1}{B}_b)F_1(t)+{\Psi }_{\mathcal{L},f}{B}_f^{-1}{F}_2(t).$$

$$u_i^{\prime }(t)$$

$$-\alpha (x_i,t)F_{xxxxt}(x_i,t)+g_u(u_i(t))F_x(x_i,t),i=2,\dots ,N-1.$$

$$\underset{Q(t)}{\underbrace{(I+A_{\alpha }(t){\tilde{\Psi }}_{xxxx})}}{S}_d^{\prime }=\underset{D(S_d)}{\underbrace{G_u(S_d){\tilde{\Psi }}_x}}{S}_d+\underset{F(S_d,t)}{\underbrace{G_u(S_d)F_x(t)-A_{\alpha }(t)F_{xxxx}^{\prime }(t)}},$$

$$A_{\alpha }(t)$$

$$G_u(S_d)$$

$$\left(\begin{array}[]{c}I_{d}+A_{\alpha}\Psi_{xxxx}\\ 0\\ 0\end{array}\right)\left(\begin{array}[]{c}S_{d}^{\prime}\\ S_{b}^{\prime}\\ S_{f}^{\prime}\end{array}\right)=\left(\begin{array}[]{c}G_{u}(S_{d})\Psi_{x}\\ I_{b}\\ B_{d}\ \ B_{b}\ \ B_{f}\end{array}\right)\left(\begin{array}[]{c}S_{d}\\ S_{b}\\ S_{f}\end{array}\right)-\left(\begin{array}[]{c}0\\ F_{1}(t)\\ F_{2}(t)\end{array}\right).$$

$$s(x,t)=\sum _{j=1}^N{\psi }_j(x)u_j(t),$$

$$u_1(t)=f_1(x_1,t),$$

$$\sum _{j=1}^N\frac{d{\psi }_j}{dx}(x_1)u_j(t)=f_2(x_1,t),$$

Full text

††journal: Journal …

1 Introduction

The Rosenau equation has become an established research subject in the field of mathematical physics since its introduction in the late 80s by Philip Rosenau rosenau . The equation is intended to overcome shortcomings of the already famous Korteweg–de Vries (KdV) equation forwhi in describing phenomena of solitary wave interaction. Knowledge about this interaction, particularly when two or more wave packets called solitons are colliding with one another, is indispensable in digital transmission through optical fibers. As data carriers, we need solitons that interact “cleanly” in the sense that none of the solitons loose any information, shape, or other conserved quantities, when they pass through each other. One may consult cipra for a fascinating history behind this subject. The Rosenau equation in its general form is given by

[TABLE]

where $\Omega$ is a bounded domain in $\mathbb{R}^{d}$ ($d\leq 3$), the coefficient $\alpha(x,t)\geq\alpha_{0}>0$ is bounded in the domain $\Omega\times[0,T]$, and the nonlinear function $g(u)$ is polynomial of degree $q+1$, $q\geq 1$. Multiple boundary conditions are required at the boundary $\partial\Omega$, such as

[TABLE]

where ${n}$ is the outward normal direction from the boundary, and we need an initial condition

[TABLE]

The well-posedness of the Rosenau equation has been studied theoretically by Park park ; phdpark . Yet in practice, the equation poses difficulties to solve numerically due to the presence of non-linearity, high spatial derivatives, multiple boundary conditions, and mixed time and space derivatives. Numerical studies based on Galerkin formulations can be found in ChuHa94 ; kimlee ; leeahn ; ChuPa01 , and numerical studies based on finite difference methods are found in Chung98 ; OAAK08 ; HuZhe10 ; AtOm15 .

The objective of this paper is to derive numerical methods based on radial basis function (RBF) collocation methods LaFo03 ; FoFly15 for the Rosenau equation, that can be applied to problems in one, two, and three space dimensions, for non-trivial geometries. These methods will also be applicable to other higher order partial differential equations. We derive and implement a fictitious point RBF collocation method and a resampling RBF collocation method, and perform experiments in one and two space dimensions. We investigate the accuracy and behavior of the derived methods theoretically and numerically. We also compare the RBF methods with a pseudospectral (PS) method Fornberg96 ; Tref00 with respect to accuracy in one space dimension.

The outline of the paper is as follows: In Section 2, a basic RBF collocation scheme is introduced. Section 3 describes different approaches to handle the multiple boundary conditions. Then in Section 4 the approximation properties of the RBF method for the one-dimensional Rosenau problem are analyzed theoretically. Implementation aspects are discussed in Section 5, followed by numerical results in Section 6. Finally, Section 7 contains conclusions and discussion.

2 The basic RBF collocation scheme

The approaches for handling multiple boundary conditions implemented in this paper are combined with an RBF collocation method. In this section, we introduce the general notation and quantities we need for RBF approximation of (time-dependent) PDEs. We start from given scalar function values $u(x_{j})$ at scattered distinct node locations $x_{j}\in\mathbb{R}^{d}$, $j=1,\ldots,N$. We assume that the function is approximated by a standard RBF interpolant

[TABLE]

where $\|\cdot\|$ is the Euclidean norm, $\phi$ is a real-valued function such as the inverse multiquadric $\phi(r)=\frac{1}{\sqrt{\varepsilon^{2}r^{2}+1}}$. The coefficients $\lambda_{j}\in\mathbb{R}$ are determined by the interpolation conditions $s(x_{j})=u(x_{j}),$ $j=1,\ldots,N$. On matrix form we have

[TABLE]

where the matrix elements $A_{ij}=\phi(\|x_{i}-x_{j}\|),~{}i,j=1,\ldots,N$, the vector $\bar{\lambda}=[\lambda_{1},\ldots,\lambda_{N}]^{T}$, and $\bar{u}=[u(x_{1}),\ldots,u(x_{N})]^{T}$. When solving a PDE, we prefer to work with the discrete function values instead of the coefficients. Using (2.1) and (2.2) together, we see that the interpolant can be written

[TABLE]

where $\bar{\phi}(x)=[\phi(\|x-x_{1}\|),\ldots,\phi(\|x-x_{N}\|)]$, assuming that $A$ is non-singular. This holds for commonly used RBFs such as Gaussians, inverse multiquadrics and multiquadrics Scho38 ; micch for distinct node points $x_{j}$. We can furthermore, use (2.3) to identify cardinal basis functions such that we can write the approximation on the FEM like form

[TABLE]

where $\bar{\phi}(x)A^{-1}=\bar{\psi}(x)=[\psi_{1}(x),\ldots,\psi_{N}(x)]$. Because our final target is to solve PDEs, we need to look at how to apply a linear operator $\mathcal{L}$ to the RBF approximation to compute $s_{\mathcal{L}}=[\mathcal{L}s(x_{1}),\ldots,\mathcal{L}s(x_{N})]^{T}$. In cardinal form, we get

[TABLE]

Using relation (2.3), the differentiation matrix $\Psi_{\mathcal{L}}=\{\mathcal{L}\psi_{j}(x_{i})\}_{i,j=1,\ldots,N}$ under operator $\mathcal{L}$ is given by

[TABLE]

where $\Phi_{\mathcal{L}}=\{\mathcal{L}\phi(\|x_{i}-x_{j}\|)\}_{i,j=1,\ldots,N}$.

For time-dependent PDE problems, we use the RBF approximation in space and then discretize the time interval. The solution $u(x,t)$ is approximated by

[TABLE]

where $u_{j}(t)\approx u(x_{j},t)$ are the unknown functions to determine.

3 Dealing with multiple boundary conditions

If we consider the one-dimensional version of equations (1.1)–(1.3) on a symmetric interval $x\in[-L,\,L]$ we have

[TABLE]

with boundary conditions

[TABLE]

Even for the one-dimensional case, how to implement multiple boundary conditions for a time-dependent global collocation problem is not obvious. In our case, we need to enforce two boundary conditions at each end point resulting in a total of four boundary conditions at the two boundary points. Collocating the PDE at all interior node points leads to a situation where we have more equations than node points. That is, unless we accept an overdetermined system, we need to either increase the number of degrees of freedom or discard some of the equations.

In fact, the subtleties of implementing multiple BCs are not tied to RBF methods. They have been actively researched in conjunction with other global collocations methods, particularly pseudospectral methods, since the 70s. We list the following five popular methods:

Mixed hard and weak BCs hesgot 2. 2.

Spectral penalty method hesta 3. 3.

Transforming to lower order system singer 4. 4.

Fictitious/ghost point method fornberg 5. 5.

Resampling method Drishal

In this paper, we only consider methods (3)–(5) as we currently do not have a way to find penalty parameters for methods (1) and (2) that give numerically stable solutions.

3.1 Transforming to lower order system

A common approach when solving PDEs containing high order derivatives is to transform them into a system with lower order derivatives. By letting $w=u_{x}$, the Rosenau equation (3.1) is transformed into

[TABLE]

with boundary conditions $u(\pm L,t)=f_{1}(\pm L,t)$ and $w(\pm L,t)=f_{2}(\pm L,t)$. The advantage of this method is that the Neumann conditions for $u$ at the boundaries become Dirichlet conditions for $w$. However, the system to solve becomes twice as large, as we need a total of $2N$ degrees of freedom for $u$ and $w$. Due to this reason, and especially for global RBFs where differentiation matrices are dense, we are not pursuing this method. However, for RBF methods in finite difference mode where differentiation matrices are sparse, this method may still be worth trying.

3.2 Fictitious point method

Fictitious or ghost point methods have been commonly used as a way to enforce multiple boundary conditions in finite difference methods. The implementation for global collocation methods such as pseudospectral methods is due to Fornberg fornberg .

The Dirichlet conditions (1.2) can be imposed by fixing the values for $u(x_{1})$ and $u(x_{N})$, but for the Neumann conditions (1.3) we use the fictitious point approach proposed by Fornberg fornberg , and introduce two additional points at some arbitrary locations such as $x_{0}$ and $x_{N+1}$.

We introduce an RBF approximation $s(x,t)$ as in (2.7), extended to include the fictitious points, for the spatial approximation of the solution $u(x,t)$,

[TABLE]

Loosely following the fictitious point approach, we will modify this ansatz so that the boundary conditions are fulfilled. Conditions (1.2) are easily fulfilled by replacing $u_{j}(t)$ with $f_{1}(x_{j},t)$ for $j=1,N$. For the conditions (1.3), we need to formally solve a linear system. Define the vectors $S_{f}=[u_{0}(t),\,u_{N+1}(t)]^{T}$ with values at the two fictitious points, and $S_{d}=[u_{2}(t),\ldots,u_{N-1}(t)]^{T}$ containing the approximate solution values at points in the interior of the domain, then we have

[TABLE]

where $F_{j}(t)=[f_{j}(x_{1},t),\,f_{j}(x_{N},t)]^{T}$. Inserting the boundary values $F_{1}(t)$ and the expression we get for $S_{f}$ by solving (3.5) into (3.4) leads to

[TABLE]

This expression is awkward to work with directly. We introduce the shorthand notation

[TABLE]

where $\tilde{\psi}_{j}(x)$ and $F(x,t)$ can be directly identified from (3.6). In this simple two point boundary case, we can actually derive the explicit form of the modified basis for illustration. This yields

[TABLE]

In order to use the RBF approximation (3.7) for a PDE problem, we need to compute the effect of applying a spatial differential operator $\mathcal{L}$ at the interior node points. That is, we need a method to evaluate $\mathcal{L}\tilde{\psi}_{j}(x_{i})$, $i,j=2,\ldots,N-1$, and $\mathcal{L}F(x_{i},t)$, $i=2,\ldots,N-1$. This is done in two steps. First, we use (2.6) to compute $\Psi_{\mathcal{L}}$ for interior node points $x_{i}$, $i=2,\ldots,N-1$. Note however that we include all basis functions $\psi_{j}(x)$, $j=0,\ldots,N+1$. Then we extract the columns pertaining to the fictitious points into $\Psi_{\mathcal{L},f}$, the columns pertaining to the boundary points into $\Psi_{\mathcal{L},b}$, and the remaining columns into $\Psi_{\mathcal{L},d}$. Then the modified differentiation matrix and the contribution in the forcing function can be computed as

[TABLE]

[TABLE]

Note that from (2.6), if no operator is applied, we have $\Psi_{d}=I$ and $\Psi_{f}=\Psi_{b}=0$ leading to $F(x_{i},t)=F_{t}(x_{i},t)=0$.

Collocating the modified RBF approximation (3.7) with the PDE (1.1) at the node points leads to the system of ODEs

[TABLE]

In matrix form, we get the method of lines formulation

[TABLE]

where the diagonal coefficient matrices are

[TABLE]

and the vectors in the right hand side are defined as $F_{\mathcal{L}}(t)=[F_{\mathcal{L}}(x_{2},t),\ldots,F_{\mathcal{L}}(x_{N-1},t)]^{T}$. The problem (3.12) can be solved by employing a solution method for nonlinear ODE systems.

The coefficient matrix $Q(t)$ is in general invertible but non-singularity cannot be guaranteed. Kansa kansa argued that if the centers of the RBFs are distinct and the PDE problem is well-posed, the coefficient matrix is generally found to be non-singular. Hon and Schaback honscha showed that occurrences of singular coefficient matrix are very rare, but do exist.

When $\alpha(x,t)$ is constant, the coefficient matrix is constant over time. Then we can LU-factorize the coefficient matrix once and use this factorization throughout the time stepping algorithm.

An alternative to the derivation above is to use the original cardinal basis functions, and include the boundary condition equations in the final system. Define rectangular identity matrices $I_{k}$ such that $I_{k}(S_{d},\,S_{b},\,S_{f})^{T}=S_{k}$, for $k=d,b,f$. Also, we order the columns in differentiation matrices in accordance with the order of the unknowns such that $\Psi_{\mathcal{L}}=[\Psi_{\mathcal{L},d},\,\Psi_{\mathcal{L},b},\Psi_{\mathcal{L},f}]$. Then we can express (3.12) as

[TABLE]

3.3 Resampling method

In the resampling method, we do not add any points as for the fictitious point method of the previous section. The four boundary conditions are still enforced at the boundary points as algebraic equations, but the PDE is instead collocated at $N-4$ auxiliary interior points. Let the solution $u(x,t)$ be approximated in Lagrange form by

[TABLE]

where $x_{1}$ and $x_{N}$ are boundary points and $x_{2},\ldots,x_{N-1}$ are interior points. The four algebraic equations arising from the boundary conditions are

[TABLE]

To write the equations on matrix form, we again divide the unknown functions into parts, $S_{d}$ for interior node points and $S_{b}$ for boundary node points. Then the boundary conditions can be expressed as

[TABLE]

where the boundary condition matrices $\tilde{B}_{d}$ and $\tilde{B}_{b}$ are analogous to those in the previous section, but with slightly different basis functions as there are no fictitious points. We have $N$ unknown functions, and we have used four equations for the boundary conditions. This means that we need $N-4$ additional equations. The idea in the resampling method is to collocate the PDE at $N-4$ auxiliary interior points, instead of collocating at the node points. We define the auxiliary points $\tilde{x}_{i}$, $i=1,\ldots,N-4$ and collocate the PDE to get the equations

[TABLE]

Define the resampling matrices $\Psi_{\mathcal{L}}^{R}=\{\mathcal{L}\psi_{j}(\tilde{x}_{i})\}_{i=1,\ldots,N-4,\,j=2,\ldots,N-1,1,N}$ (columns ordered according to the unknowns). The resampled equations (3.18), together with the algebraic equations (3.17), yield an $N\times N$ differential algebraic system of equations (DAE) of index 1,

[TABLE]

where $\tilde{A}_{\alpha}=\mathrm{diag}(\alpha(\tilde{x}_{1},t),\ldots,\alpha(\tilde{x}_{N-4},t))$ and \tilde{S}=\Psi^{R}\left(\begin{array}[]{c}S_{d}\\ S_{b}\\ \end{array}\right).

The system of equations (3.19) can be solved using a differential algebraic solver. See DASPK brohinpet ; petzold . An example of how this can be implemented in MATLAB is given in section 5.

3.4 Generalization to more space dimensions

The main differences when moving to more than one space dimensions is that we have a boundary curve or a boundary surface that is discretized in the same way as the interior of the domain instead of just two boundary points. The formulation of the two methods is in all essential parts the same, and the formulations (3.13) and (3.19) are valid in the same form, but when we before had two boundary points and two fictitious points, we instead have $N_{b}$ boundary points and $N_{b}$ fictitious points. Similarly, for the resampling method, we have $2N_{b}$ boundary conditions, and therefore we collocate the PDE at $N-2N_{b}$ auxiliary points. Experiments for problems in two space dimensions are presented in Section 6.

4 Error estimates

In this section, we derive semi-discrete error estimates for the one-dimensional problem. The analysis is carried out for the fictitious point approach. For the analysis, we assume that $\alpha>0$ is constant and that the function $g$ is a polynomial of degree $q+1$ with $q\geq 1$.

4.1 ODE system for the semi-discrete approximation error

We work with spatially discrete solution and approximation values evaluated at the interior node points $x_{i}$, $i=2,\ldots,N-1$. We denote the exact solution vector and its derivatives by $U_{\mathcal{L}}(t)=(\mathcal{L}u(x_{2},t),\ldots,\mathcal{L}u(x_{N-1},t))^{T}$. From the Rosenau equation (3.1), we have

[TABLE]

For the RBF approximation, we use (3.12), but to simplify notation we replace $S_{d}$ with $S$ and the matrix $A_{\alpha}$ with the constant $\alpha$.

[TABLE]

We introduce the auxiliary function $z(x,t)$, which interpolates the exact solution at each time and also satisfies the boundary conditions (1.2) and (1.3),

[TABLE]

We denote the auxiliary function vector by $Z(t)$ and note that

[TABLE]

Furthermore, $Z(t)=U(t)$, $Z^{\prime}(t)=U^{\prime}(t)$, and $Z(0)=S(0)=U(0)$, while $Z_{\mathcal{L}}\neq U_{\mathcal{L}}$.

We define the discrete interpolation error vector $\epsilon(t)=U(t)-Z(t)$. The interpolation error itself is zero at the node points, but its derivatives $\epsilon_{\mathcal{L}}$ are non-zero. By noting that $U(t)=\epsilon(t)+Z(t)$ and by using (4.4) for the derivatives of $Z(t)$, we can rewrite (4.1) as

[TABLE]

Finally, we introduce the discrete error $E(t)=U(t)-S(t)$, and subtract (4.2) from (4.5) to get

[TABLE]

This equation can be seen as an ODE-system for the error $E(t)$ by writing it as

[TABLE]

where $H(t)=-\alpha\epsilon^{\prime}_{xxxx}(t)+G_{u}(U)\tilde{\Psi}_{x}U-G_{u}(S)\tilde{\Psi}_{x}S+G_{u}(U)\epsilon_{x}(t)+(G_{u}(U)-G_{u}(S))F_{x}(t)$. In the following, we will consider $H(t)$ as a forcing function. For a discussion of the non-singularity of $Q$, see Section 3.2. The system of ODEs (4.7) can be formally integrated to yield

[TABLE]

In the following subsections, we will look into each part of the forcing function.

4.2 Estimates for the non-linear term

In order to determine the influence of the non-linear term on convergence and stability, we consider the particular form $g(u)=u^{q+1}$, $q\geq 1$. This matches the functions typically used in the literature. Furthermore, we will use an estimate by Park from park showing that with this form of $g(u)$,

[TABLE]

We can then form the following estimate

[TABLE]

Note that the exponent can never become larger than $1.25$, which occurs at $q=2.5$. For the second derivative, we have

[TABLE]

If we instead consider a point $s$, close to $u$ we have

[TABLE]

which allows the following estimate

[TABLE]

where $\tilde{q}=\max((q-1)(1-\frac{q}{5}),(1-\frac{q}{5}))$.

4.3 Term by term estimates for the error contributions

In this section, we are going to derive an estimate for the discrete approximation error. We first note that $H(t)$ is a sum over number of terms $H_{j}(t)$ and split the integral in (4.8) to get

[TABLE]

For the first error term, we have $H_{1}(t)=-\alpha\epsilon_{xxxx}^{\prime}(t)$. Integration leads to

[TABLE]

then we can get the following estimate

[TABLE]

The second error term that we consider is generated by $H_{2}(t)=G_{u}(U)\epsilon_{x}(t)$, leading to

[TABLE]

where we used (4.10) for the term involving $G_{u}(U)$.

The final part focuses on the non-linear term and is the most complicated. We start by rewriting the generating term

[TABLE]

If we rewrite equation (3.10) in matrix form, the function $F_{x}(t)$ can be expressed as

[TABLE]

Using this, we can provide the first estimate for contribution to the error from the term $H_{3}(t)$ given by (4.22)

[TABLE]

Using equations (4.9), (4.11), and (4.18) together with the inequality $|U+(S-U)|\leq|U|+|S-U|$ yields

[TABLE]

Because this term contains the error in the right hand side, it is the most difficult one to include in the error estimate. We simplify it as far as possible. First we write

[TABLE]

where

[TABLE]

Then we make the observation that either the error is small and $\|E(\tau)\|\geq\|E(\tau)\|^{p}$ for $p>1$ or the error is larger than one, in which case $\|E(\tau)\|^{q+1}\geq\|E(\tau)\|^{p}$ for $p\leq q$. We let $\nu=1$ or $\nu=q+1$, sum up the coefficients and pick the highest power of $(1+\tau)$ to get

[TABLE]

where

[TABLE]

Table 1 shows the different powers that are involved as a function of $q$. The final power in the estimate grows with $q$.

4.4 Global error bound

By combining the error terms (4.20), (4.21) and (4.30), we get the following relation for the error due to the spatial discretization

[TABLE]

where

[TABLE]

To convert this into an error estimate, we need the following Gronwall inequality drago :

Lemma 1 (Gronwall inequality)

Let the functions $E(t)$, $a(t)$ and $k(t)\geq 0$ be continuous functions defined for $t\in[0,b]$. We assume that for $t\in[0,b]$ we have the inequality

[TABLE]

Then for the case $n=1$ it holds

[TABLE]

If the function $a(t)$ is also non-decreasing, then

[TABLE]

For the case $n\geq 2$

[TABLE]

where $\beta_{n}=\sup\{t:(n-1)\int_{0}^{t}k(\tau)a^{n-1}(\tau)\,\mathrm{d}\tau<1\}$.

In our case, it can easily be verified that the function

[TABLE]

is non-decreasing in time. For the case of small errors, $n=\nu=1$, and

[TABLE]

the Gronwall inequality leads to

[TABLE]

For the case of errors larger than one, we do not carry out the full derivation, but note that the limit on the time interval becomes less severe when $a(t)$ is small enough, which is the case when the spatial resolution is high enough.

It remains to insert estimates for the derivatives of the RBF interpolation errors. These are typically measured in terms of the fill distance $h$, which in one space dimension with uniform nodes becomes $h=\frac{1}{2}|x_{j+1}-x_{j}|$, and more generally for a domain $\Omega$ in $d$ dimensions and a node set $X$ is defined as

[TABLE]

Exponential estimates for inverse multiquadric interpolants are given in rigzwi provided that $\Omega$ is a bounded domain with Lipschitz boundary that satisfies an interior cone condition.

[TABLE]

where the constant $\gamma$ depends on the order of the differential operator $\mathcal{L}$, the dimension $d$, and the shape parameter of the RBF $\varepsilon$, and $c$ is a positive constant. The norm in the right hand side denoted by $\|.\|_{\mathcal{N}}(.)$ is the native space norm. For more details about this norm see Fass ; rigzwi .

Now, by inserting the interpolation error estimate (4.44) into the global error estimate (4.42) we get

[TABLE]

where $C(t)=c(C_{1}+C_{2}(t))$. The final estimate tells us that the spatial RBF discretization does allow for exponential convergence in $h$, but we should expect the error to grow in time. For any finite interval $t\in[0,T]$ the growth in time is however bounded.

4.5 Numerical investigation of the matrix norms in the estimates

There are three different matrix norms that appear in the error estimates. We do not have theoretical bounds for these, and therefore we perform a numerical investigation of their behaviors. Based on previous experience of RBF approximation, we have selected the following representation of the method parameters, the fill distance $h$, the relative shape parameter value $\varepsilon h$, and the (half) domain size $L$. In Figure 1, the norms are plotted as a function of $h$ for different combinations of the parameters. The chosen values are such that they are reasonable for approximation. By choosing extreme values, different behaviors can be observed. We see that the norm $\|Q^{-1}\|$ does not change much with $h$ or $\varepsilon h$, but slightly with $L$. Both of the norms $\|B_{x}\|$ and $\|\tilde{\Psi}_{x}\|$ grow as $h^{-1}$, and we can observe a little bit of instability in the value for small $\varepsilon$. This rate is what would be expected for a first derivative such as is represented by these matrices. Changing $L$ has a very small effect also in this case.

5 MATLAB Implementation

In this section, sample MATLAB implementations of the fictitious point and resampling RBF methods for the one-dimensional Rosenau equation (3.1)–(3.3) are presented and discussed. We use an example with a known solution. For $g(u)=10u^{3}-12u^{5}-\frac{3}{2}u$ and $\alpha(x,t)=0.5$ it holds that $u(x,t)=\textnormal{sech}(x-t)$ is a solution park . For both methods, equally spaced nodes are used, and the spatial domain is $[-L,L]$.

5.1 Implementation of the fictitious point method

Following the idea of the fictitious point method in subsection 3.2, we complement the interior and boundary RBF nodes with two (the number of boundary nodes) fictitious points outside the computational domain, see Figure 2 for an illustration.

We generate the differentiation matrices using the modified basis functions according to Equation (3.9). Collocating the Rosenau equation by applying the modified differentiation matrices leads to ODE system (3.12), which we here solve by using the built-in MATLAB ODE solver ode15s. The two functions below constitute a complete MATLAB implementation of the problem.

function [S,x,t]=fictitious(N,L,T) % N : The number of node points % L : [-L,L] is the domain % T : Final time % Exact solution and derivatives for test case u = @(x,t) sech(x-t); ux = @(x,t) sech(t-x).tanh(t-x); ut = @(x,t) -sech(t-x).tanh(t-x); uxt = @(x,t) sech(t - x).^3 - sech(t - x).tanh(t - x).^2; % Generate nodes x. Map x to [-L,L] such that x(2) = -L, x(N-1) = L % and x(1), x(N) are the left and right fictitious points respectively. x = linspace(-L,L,N); linmap = @(x,x1,x2,y1,y2) (y2-y1)(x-x1)/(x2-x1) + y1; x = linmap(x,x(2),x(N-1),-L,L); x = x(:); % Differentiation matrices for inverse multiquadric RBF phi = @(ep,r) 1./sqrt(1+(epr).^2); ep = 0.08/min(diff(x)); % Shape parameter dx = bsxfun(@minus,x,x.’); A = phi(ep,dx); Dx = (-ep^2dx.A.^3)/A; % 1st Derivative matrix D4x = (3ep^4*(3-24*(epdx).^2+8(ep*dx).^4).*A.^9)/A; % 4th Derivative Bf = Dx([2 N-1],[1 N]); Bd = Dx([2 N-1],3:N-2); Bb = Dx([2 N-1],[2,N-1]); % Modify differentiation matrices Dxd = Dx(3:N-2,3:N-2) - (Dx(3:N-2,[1 N])/Bf)*Bd; Dxb = [Dx(3:N-2,[2 N-1]) - (Dx(3:N-2,[1 N])/Bf)*Bb Dx(3:N-2,[1 N])/Bf]; D4xd = D4x(3:N-2,3:N-2) - (D4x(3:N-2,[1 N])/Bf)*Bd; D4xb = [D4x(3:N-2,[2 N-1]) - (D4x(3:N-2,[1 N])/Bf)*Bb D4x(3:N-2,[1 N])/Bf]; % Initial condition S0 = u(x(3:N-2),0); opt = odeset(’RelTol’,1e-10); % Solve the ODE for the approximate solution S fun = @(t,S) odefun(t,S,x,N,Dxd,Dxb,D4xd,D4xb,u,ux,ut,uxt); [t,S] = ode15s(fun,[0 T],S0,opt); % Plot the solution for all times figure(1),clf,plot(x(3:N-2),S) % Plot the error at the final time E = abs(S(end,:)-u(x(3:N-2),T)’); figure(1),clf,plot(x(3:N-2),E)

function Sprime = odefun(t,S,x,N,Dxd,Dxb,D4xd,D4xb,u,ux,ut,uxt) Fx = [u(x([2 N-1]),t); ux(x([2 N-1]),t)]; F4xt = [ut(x([2 N-1]),t); uxt(x([2 N-1]),t)]; Gu = diag(-1.5 - 60S.^4 + 30S.^2); Sprime = (eye(N-4) + 0.5D4xd)(GuDxdS + GuDxbFx - 0.5D4xb*F4xt);

It can be noted that when we use the modified basis functions, we need to provide the time derivatives of the boundary conditions as well as the boundary conditions themselves. This is not needed with the alternative formulation (3.13), but instead the system is stated in DAE form.

5.2 Implementation of the resampling RBF method

For the resampling method, we start directly from the DAE form derived in Section 3.3, where the four boundary conditions enforced at the boundary points constitute the algebraic part. The $N-4$ auxiliary points where the PDE is enforced are uniformly distributed in the interior of the computational domain and do not in general coincide with the RBF node points where the solution is sought. We organize the solution vector as $S=[S_{d}~{}S_{b}]^{T}$, where as before, $S_{d}$ contains solution values at the interior RBF nodes, and $S_{b}$ contains the two boundary values. Then the DAE discretization scheme can be schematically be displayed as

[TABLE]

where $\Psi^{R}$ is an $N-4\times N$ resampling matrix that provides values at the auxiliary points given values at the node points, $\tilde{S}=\Psi^{R}S$, $G_{u}(\tilde{S})$ is an $(N-4)\times(N-4)$ diagonal matrix, $\Psi^{R}_{x}$ and $\Psi^{R}_{xxxx}$ are $N-4\times N$ resampled first and fourth derivative matrices respectively. The derivative boundary condition matrix $\tilde{B}=(\tilde{B}_{d}~{}~{}\tilde{B}_{b})$ is a $2\times N$ matrix, see Equation (3.17) for details.

Both ODEs and DAEs of index 1 can be solved in MATLAB using ode15s, previously employed for the fictitious point method. One may also use the syntactically similar open source software OCTAVE and there use dasspk as DAE solver. The following two functions show the MATLAB implementation of the resampling RBF method.

function [S,x,t]=resampling(N,L,T) % N : The number of node points % L : [-L,L] is the domain % T : Final time % Exact solution and derivatives for test case u = @(x,t) sech(x-t); ux = @(x,t) sech(t-x).tanh(t-x); % Generate N uniform RBF nodes with boundary pts last x = linspace(-L,L,N).’; x = [x(2:end-1); x([1 end])]; % Generate N-4 uniform auxiliary interior points xt = linspace(-L,L,N-3).’; xt(end) = []; xt = xt + 0.5min(diff(xt)); % Differentiation matrices for inverse multiquadric RBF phi = @(ep,r) 1./sqrt(1+(epr).^2); ep = 0.08/(x(2)-x(1)); % Shape parameter r = bsxfun(@minus,x,x.’); A = phi(ep,r); % First derivative matrix at x = -L and x = L r = bsxfun(@minus,x([N-1 N]),x.’); Bt = (-ep^2r.phi(ep,r).^3)/A; % Rectangular projection from x to xt r = bsxfun(@minus,xt,x.’); R = phi(ep,r); PRx = (-ep^2r.R.^3)/A; PR = R/A; PR4x = (3ep^4*(3-24*(epr).^2+8(ep*r).^4).R.^9)/A; % Initial condition S0 = u(x,0); M = [(PR + 0.5PR4x); zeros(4,N)]; opt = odeset(’mass’,M,’masssing’,’yes’,’RelTol’,1e-9); [t,S] = ode15s(@(t,S) daefun(t,S,x,N,PR,PRx,Bt,u,ux),[0 T],S0,opt); % Plot the solution for all times ind = [N-1 2:N-2 N]; figure(1),clf,plot(x(ind),S(:,ind)) % Plot the error at the final time E = abs(S(end,:)-u(x,T)’); figure(2),clf,plot(x(ind),E(ind))

function Sprime = daefun(t,S,x,N,PR,PRx,Bt,u,ux) F = [zeros(N-4,1); u(x([N-1 N]),t); ux(x([N-1 N]),t)]; Gu = diag(- 60*(PRS).^4 + 30(PRS).^2 - 1.5); Sprime = [GuPRx; [zeros(2,N-2) eye(2,2)]; Bt]*S - F;

6 Numerical results

In this section, we perform numerical experiments to investigate the accuracy and convergence of the proposed schemes. Both one-dimensional and two-dimensional test cases are considered. In all tests, the inverse multiquadric RBF is used. The shape parameter has not been optimized for accuracy. Instead, the parameter range has been chosen such that ill-conditioning is avoided. For the one-dimensional test case, we compare the results with those of a pseudo-spectral resampling method. We have not included the code here, but it can be downloaded from the authors’ web pages.

6.1 The one-dimensional case

We consider the same test problem park for the Rosenau equation (3.1) as in the previous section, with the known exact solution $u(x,t)=\textrm{sech}(x-t)$, obtained for $g(u)=10u^{3}-12u^{5}-\frac{3}{2}u$ and $\alpha(t)=0.5$. The initial and boundary conditions are taken from the exact solution.

The exact solution over the interval $[-L,\,L]$ for $L=1$ and $L=10$ is shown at different times $t$ together with the numerical solution using the fictitious point method in Figure 3. The solution is a pulse that travels to the right with time.

In Figure 3, a shape parameter $\varepsilon=\frac{0.08}{h}=\frac{0.04(N-1)}{L}$ is used. This relation is experimentally determined to ensure stable computations and highest possible solution accuracy. Figure 4 shows how the errors of the two RBF methods vary with $\varepsilon$. Using the formula leads to $\varepsilon=1.16$ and $\varepsilon=0.4$ for the two cases, which is within the best region for each method. It can be noted that the stable range of $\varepsilon$ is narrower for the resampling method and that both methods rapidly become ill-conditioned when $\varepsilon$ is too small.

To illustrate the capability of the proposed methods, we start with comparing the approximation accuracy with that of a pseudo-spectral resampling method. The pseudo-spectral method employs the same number of Chebyshev nodes as the number of uniform nodes used by the RBF methods. For a description of its implementation, see Drishal . The absolute errors for two different values of $L$ are plotted in Figure (5). As shown in the figure, the errors of both the RBF based methods and the pseudospectral resampling method are similar in magnitude. For the shorter interval, the pseudo-spectral method has smaller errors near the boundaries, which is consistent with the clustering of the Chebyshev nodes. However, for the larger interval, where the solution is small at the boundary, this effect is not visible.

The $L_{\infty}$ errors over time for the approximated solutions are illustrated in Figure 6. If we go back to the global error estimate (4.45), and insert $q=4$ (for this test case), the exponential time-dependent growth factor becomes $\exp(C_{3}((1+t)^{1.12}-1)/1.12)$. We do not know the precise value of $C_{3}$, but based on our experiments a value larger than one should be expected, in which case the predicted growth would be at least two orders of magnitudes larger than what is observed. However, the growth factor in the estimate arises from the growth estimate (4.9) for the solution $u(x,t)$, which in our particular case does not grow at all. Hence, it would be possible to make a more favorable estimate for our special case. Both the accuracy and the growth rate of the errors of the three different solutions are similar. For the shorter interval $L=1$, the resampling RBF method is slightly worse than the other two methods, while for the longer interval $L=10$, both RBF methods are slightly more accurate than the pseudo-spectral method.

Figure 7 displays the convergence behavior as a function of $N$ for the two RBF methods compared with the PS method. For all three methods, the highest attainable accuracy is the same. When $\varepsilon h$ is constant, as in this experiment, we would expect the error to reach a saturation level, but accuracy is also limited by conditioning, and this may be the effect that we see here. In both cases, the fictitious point RBF method reaches the highest accuracy faster than the resampling RBF method. The PS method performs best for the short interval, and performs worse than the RBF methods for the longer interval. One explanation for this can be that the node density for the Chebyshev nodes compared with the uniform nodes is lower in the interesting region (middle of the domain) in this case.

6.2 Two-dimensional square domain

In this section, we demonstrate how the flexibility of the RBF approximations allows us to implement the fictitious point method and resampling RBF method in a two-dimensional domain with a similar effort as for the one-dimensional problem. Here, we do not compare with the resampling PS method, which is less straightforward to implement. We consider the square domain $\Omega=[-L,\,L]\times[-L,\,L]$ and the Rosenau equation (1.1) with $\alpha=1$, $g(u)=u^{3}+u^{2}$ and initial condition $f_{0}(x,y)=\textrm{sech}(x+y)$ and boundary conditions

[TABLE]

where the derivative in the second condition is taken as either $u_{x}$ or $u_{y}$ depending on which part of the boundary is involved. For the two-dimensional test cases, we do not have any analytical solutions. The approximate solution at two different times is displayed in Figure 8.

We start from a uniform discretization of the domain $\Omega$ with $N=n^{2}$ points. We denote the number of interior points by $N_{d}$ and the number of boundary points by $N_{b}$. For the fictitious point method we need to add $N_{b}$ fictitious points outside the domain to enforce the Neumann boundary conditions. Note that if we simply choose an extension of the uniform grid, the resulting number of fictitious points is too large. For the resampling RBF method we generate $N_{d}-2N_{b}$ auxiliary points inside of the domain. Note that here the number of boundary points is modified (and they are therefore not on the uniform grid) to make the interior and auxiliary node numbers compatible. Sample node distributions for the two methods are displayed in Figure 9.

According to the error estimate (4.45) for the fictitious point method, we expect exponential convergence in $1/\sqrt{h}$ for fixed $\varepsilon$. In practice, we often observe exponential convergence in $1/h$. In Figure 10, we plot the error as a function of $n\propto 1/h$. For this range of $n$-values, the conditioning is low enough to not influence the error, and exponential convergence can be observed for both RBF methods. We see that the estimated slopes in the right subfigure are precisely double those in the left subfigures. If we take into account that $h$ is also twice as large for the case $L=2$, we can conclude that the rate of convergence in terms of $1/h$ is the same in both cases. Again, the fictitious point is more accurate, while the resampling method converges faster. The fictitious point uses a larger number of basis functions, which could account for the smaller error, but it is not clear why the rates differ in the way they do.

6.3 Two-dimensional starfish domain

We now take a step further in demonstrating the flexibility of the RBF based methods by considering an irregular two-dimensional domain. As a test problem, we consider the starfish domain with boundary defined by the parametric equation

[TABLE]

We also need the derivative of the boundary equation in order to compute the outward normal direction $n=(n_{x},n_{y})$, which is needed for the boundary conditions. We have

[TABLE]

We use a similar test problem as for the square domain with boundary Dirichlet data

[TABLE]

For the normal derivative condition, we impose

[TABLE]

The approximate solution at three different times is shown in Figure 11.

Sample node distributions for the fictitious point and resampling RBF methods are illustrated in Figure 12.

Just as for the square domain, the number of fictitious points outside the domain is $N_{b}$ and the number of auxiliary points inside the domain in the resampling method is $N_{d}-2N_{b}$.

The max error as function of $N$ for a fixed shape parameter value, $\varepsilon=2$, is illustrated in Figure 13. The reference solution is computed using the fictitious point method with $N=540$ nodes. The max error is estimated from evaluation at $600$ radially uniformly distributed points in the domain.

The error trends are similar to those for the square domain, showing that the RBF methods provide a well functioning generalization of both the fictitious point method and the resampling method to general domains.

7 Conclusion

The Rosenau equation, which is used as an application throughout this paper, is an example of a non-linear PDE with multiple boundary conditions as well as mixed space-time derivatives. Multiple boundary conditions provides an extra challenge when solving PDE problems. The standard form of a typical collocation method assumes that one condition is imposed at each node point/grid point. Hence, the additional conditions at the boundary nodes lead to a mismatch between the number of conditions and the number of unknowns.

Two approaches to manage multiple boundary conditions that have introduced for spectral methods are fictitious point methods and resampling methods. In this paper we have shown how to implement these approaches in the context of RBF collocation methods. From numerical experiments for a one-dimensional test problem, we could see that the behavior of the method with respect to accuracy in space and time is very similar to that of the corresponding pseudo-spectral method.

For two-dimensional problems, already in a regular geometry such as the square, the application of spectral methods becomes more complicated. Approximations are typically based on tensor product grids, but if we use the one-dimensional extension techniques for each grid line, again the number of extra conditions and extra points do not naturally match. The problem can for example be resolved by choosing one of the directions for the corner points, but then the approximations in the other direction needs to be of lower order.

We show that with the two RBF methods, due to the freedom of node placement, we can distribute the fictitious points or the resampled nodes uniformly and symmetrically with respect to the domain. Furthermore, we show that the concept can be transferred also to irregularly shaped domains.

We have also analyzed the theoretical properties of the fictitious point RBF approximation for the one-dimensional Rosenau equation. We could show that the spectral convergence of the spatial approximation carries over to the PDE solution, while the growth of the error in time in our estimate strongly depends on the growth estimates for the exact solution of the problem. Using a generally valid growth estimate, as we did, led to an overestimate of the error growth in time for the test problem we used, where the norm of the exact solution does not grow.

To conclude, both the implemented approaches are promising for problems with multiple boundary conditions, especially for geometries where spectral methods cannot easily be applied. Global RBF approximations as the ones used here are competitive for problems in one or two space dimensions, but the computational cost can become prohibitive for higher-dimensional problems due to the need to solve dense linear systems. Therefore, an interesting future direction is to see how resampling and fictitious point methods can be combined with localized (stencil or partition based) RBF methods.