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

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.


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 [28]. The equation is intended to overcome shortcomings of the already famous Korteweg-de Vries (KdV) equation [13] 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 [6] for a fascinating history behind this subject. The Rosenau equation in its general form is given by where Ω is a bounded domain in R d (d ≤ 3), the coefficient α(x,t) ≥ α 0 > 0 is bounded in the domain Ω × [0, T ], and the nonlinear function g(u) is polynomial of degree q + 1, q ≥ 1. Multiple boundary conditions are required at the boundary ∂ Ω , such as where n is the outward normal direction from the boundary, and we need an initial condition The well-posedness of the Rosenau equation has been studied theoretically by Park [25,24]. 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 [4,19,21,5], and numerical studies based on finite difference methods are found in [3,23,17,1].
The objective of this paper is to derive numerical methods based on radial basis function (RBF) collocation methods [20,12] 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 [10,31] 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 onedimensional 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.

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 ∈ R d , j = 1, . . . , N. We assume that the function is approximated by a standard RBF interpolant where · is the Euclidean norm, φ is a real-valued function such as the inverse multiquadric φ (r) = 1 √ ε 2 r 2 +1 . The coefficients λ j ∈ R are determined by the interpolation conditions s(x j ) = u(x j ), j = 1, . . . , N. On matrix form we have Aλ =ū, (2.2) where the matrix elements A i j = φ ( x i − x j ), i, j = 1, . . ., N, the vectorλ = [λ 1 , . . . , λ N ] T , andū = [u(x 1 ), . . ., 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 , assuming that A is non-singular. This holds for commonly used RBFs such as Gaussians, inverse multiquadrics and multiquadrics [29,22] 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 Because our final target is to solve PDEs, we need to look at how to apply a linear operator L to the RBF approximation to compute s L = [L s(x 1 ), . . ., L s(x N )] T . In cardinal form, we get Using relation (2.3), the differentiation matrix Ψ L = {L ψ j (x i )} i, j=1,...,N under operator L is given by 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 where u j (t) ≈ u(x j ,t) are the unknown functions to determine.

Dealing with multiple boundary conditions
If we consider the one-dimensional version of equations (1.1)-(1.3) on a symmetric interval x ∈ [−L, L] we have with boundary conditions 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: 1. Mixed hard and weak BCs [15] 2. Spectral penalty method [14] 3. Transforming to lower order system [30] 4. Fictitious/ghost point method [11] 5. Resampling method [8] 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 with boundary conditions u(±L,t) = f 1 (±L,t) and w(±L,t) = f 2 (±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.

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 [11].
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 [11], 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), (3.4) 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), . . ., u N−1 (t)] T containing the approximate solution values at points in the interior of the domain, then we have Inserting the boundary values F 1 (t) and the expression we get for S f by solving (3.5) into (3.4) leads to This expression is awkward to work with directly. We introduce the shorthand notation whereψ 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ψ 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 L at the interior node points. That is, we need a method to evaluate Lψ j (x i ), i, j = 2, . . ., N − 1, and L F(x i ,t), i = 2, . . . , N − 1. This is done in two steps. First, we use (2.6) to compute Ψ L for interior node points x i , i = 2, . . . , N − 1. Note however that we include all basis functions ψ j (x), j = 0, . . . , N + 1. Then we extract the columns pertaining to the fictitious points into Ψ L , f , the columns pertaining to the boundary points into Ψ L ,b , and the remaining columns into Ψ L ,d . Then the modified differentiation matrix and the contribution in the forcing function can be computed as (3.10) Note that from (2.6), if no operator is applied, we have Ψ d = I and Collocating the modified RBF approximation (3.7) with the PDE (1.1) at the node points leads to the system of ODEs In matrix form, we get the method of lines formulation , (3.12) where the diagonal coefficient matrices are G u (S d ) = diag(g u (u 2 (t)), . . ., g u (u N−1 (t))), and the vectors in the right hand side are defined as F L (t) = [F L (x 2 ,t), . . ., F 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 [18] 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 [16] showed that occurrences of singular coefficient matrix are very rare, but do exist.
When α(x,t) is constant, the coefficient matrix is constant over time. Then we can LUfactorize 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 f . Also, we order the columns in differentiation matrices in accordance with the order of the unknowns such that (3.13)

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 where x 1 and x N are boundary points and x 2 , . . . , x N−1 are interior points. The four algebraic equations arising from the boundary conditions are 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 where the boundary condition matricesB d andB 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 pointsx i , i = 1, . . . , N − 4 and collocate the PDE to get the equations ..,N−4, j=2,...,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 ×N differential algebraic system of equations (DAE) of index 1, The system of equations (3.19) can be solved using a differential algebraic solver. See DASPK [2,26]. An example of how this can be implemented in MATLAB is given in section 5.

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.

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 α > 0 is constant and that the function g is a polynomial of degree q + 1 with q ≥ 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, . . . , N − 1. We denote the exact solution vector and its derivatives by For the RBF approximation, we use (3.12), but to simplify notation we replace S d with S and the matrix A α with the constant α.
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), We denote the auxiliary function vector by Z(t) and note that The interpolation error itself is zero at the node points, but its derivatives ε L are non-zero. By noting that U(t) = ε(t) + Z(t) and by using (4.4) for the derivatives of Z(t), we can rewrite (4.1) as Finally, we introduce the discrete error E(t) = U(t) − S(t), and subtract (4.2) from (4.5) to get This equation can be seen as an ODE-system for the error E(t) by writing it as where In the following, we will consider H(t) as a forcing function. For a discussion of the nonsingularity of Q, see Section 3.2. The system of ODEs (4.7) can be formally integrated to yield In the following subsections, we will look into each part of the forcing function.

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 ≥ 1. This matches the functions typically used in the literature. Furthermore, we will use an estimate by Park from [25] showing that with this form of g(u), We can then form the following estimate Note that the exponent can never become larger than 1.25, which occurs at q = 2.5. For the second derivative, we have If we instead consider a point s, close to u we have which allows the following estimate

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 For the first error term, we have H 1 (t) = −αε ′ xxxx (t). Integration leads to then we can get the following estimate The second error term that we consider is generated by 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 (4.18) If we rewrite equation (3.10) in matrix form, the function F x (t) can be expressed as Using this, we can provide the first estimate for contribution to the error from the term H 3 (t) given by (4.18) Using equations (4.9), (4.11), and (4.14) together with the inequality |U + (S −U)| ≤ |U| + |S −U| yields (4.21) 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 Then we make the observation that either the error is small and E(τ) ≥ E(τ) p for p > 1 or the error is larger than one, in which case E(τ) q+1 ≥ E(τ) p for p ≤ q. We let ν = 1 or ν = q + 1, sum up the coefficients and pick the highest power of (1 + τ) to get whereC q = (q + 1)(2 + max 0≤τ≤t F(τ) ) max(C,C q ,C q C).
(4.27) Table 1 shows the different powers that are involved as a function of q. The final power in the estimate grows with q.

Global error bound
By combining the error terms (4.16), (4.17) and (4.26), we get the following relation for the error due to the spatial discretization where To convert this into an error estimate, we need the following Gronwall inequality [7]: Lemma 1 (Gronwall inequality) Let the functions E(t), a(t) and k(t) ≥ 0 be continuous functions defined for t ∈ [0, b]. We assume that for t ∈ [0, b] we have the inequality Then for the case n = 1 it holds If the function a(t) is also non-decreasing, then For the case n ≥ 2 where β n = sup{t : (n − 1) t 0 k(τ)a n−1 (τ) dτ < 1}. In our case, it can easily be verified that the function is non-decreasing in time. For the case of small errors, n = ν = 1, and the Gronwall inequality leads to . (4.38) 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 = 1 2 |x j+1 − x j |, and more generally for a domain Ω in d dimensions and a node set X is defined as Exponential estimates for inverse multiquadric interpolants are given in [27] provided that Ω is a bounded domain with Lipschitz boundary that satisfies an interior cone condition.
where the constant γ depends on the order of the differential operator L , the dimension d, and the shape parameter of the RBF ε , and c is a positive constant. The norm in the right hand side denoted by . N (.) is the native space norm. For more details about this norm see [9,27]. Now, by inserting the interpolation error estimate (4.40) into the global error estimate (4.38) we get 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 ∈ [0, T ] the growth in time is however bounded.

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 ε 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 ε h, but slightly with L. Both of the norms B x and Ψ x grow as h −1 , and we can observe a little bit of instability in the value for small ε . 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.

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 − 3 2 u and α(x,t) = 0.5 it holds that u(x,t) = sech(x − t) is a solution [25]. For both methods, equally spaced nodes are used, and the spatial domain is [−L, L].

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 differen-  tiation 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. 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.

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

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 illconditioning 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.

The one-dimensional case
We consider the same test problem [25] for the Rosenau equation (3.1) as in the previous section, with the known exact solution u(x,t) = sech(x − t), obtained for g(u) = 10u 3 − 12u 5 − 3 2 u and α(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 ε = 0.08 h = 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 ε . Using the formula leads to ε = 1.16 and ε = 0.4 for the two cases, which is within the best region for each method. It can be noted that the stable range of ε is narrower for the resampling method and that both methods rapidly become ill-conditioned when ε 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 [8]. 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 ∞ errors over time for the approximated solutions are illustrated in Figure 6. If we go back to the global error estimate (4.41), 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. Both RBF methods use ε = 0.08 h and a uniform node distribution, while the PS method employs a Chebyshev node distribution. 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 ε 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.  Fig. 7 The L ∞ error at time t = 1 as a function of the number of node points N for the fictitious point RBF method, the resampling RBF method, and a resampling PS method. Results are shown for L = 1 (left) and L = 10 (right). Both RBF methods use ε = 0.08 h and uniform node distributions, while the PS method employs Chebyshev node distributions.

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 Ω = [−L, L] × [−L, L] and the Rosenau equation (1.1) with α = 1, g(u) = u 3 + u 2 and initial condition f 0 (x, y) = sech(x + y) and boundary conditions 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 Ω 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.41) for the fictitious point method, we expect exponential convergence in 1/ √ h for fixed ε . In practice, we often observe exponential convergence in 1/h. In Figure 10, we plot the error as a function of n ∝ 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   Fig. 10 The error at time t = 1 against the number of points per dimension n for L = 1 and ε = 1.6 (right) and L = 2 and ε = 0.8 (left). In both cases, errors are computed against a reference solution computed with the fictitious point RBF method for n = 28, and the error is evaluated at 25 × 25 interior points.

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 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 (n x , n y ) = r ′ (θ )(sin(θ ), − cos(θ )) + r(θ )(cos(θ ), sin(θ )) We use a similar test problem as for the square domain with boundary Dirichlet data For the normal derivative condition, we impose The approximate solution at three different times is shown in Figure 11. The max error as function of N for a fixed shape parameter value, ε = 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.

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.