- Email: [email protected]
BS2 methods for semi-linear second order boundary value problems
Applied Mathematics and Computation xxx (2014) xxx–xxx
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
BS2 methods for semi-linear second order boundary value problems Carla Manni a, Francesca Mazzia b, Alessandra Sestini c,⇑, Hendrik Speleers a,d a
Dipartimento di Matematica, Università di Roma ‘‘Tor Vergata’’, Italy Dipartimento di Matematica, Università di Bari, Italy c Dipartimento di Matematica e Informatica, Università di Firenze, Italy d Departement Computerwetenschappen, Katholieke Universiteit Leuven, Belgium b
a r t i c l e
i n f o
a b s t r a c t A new class of Linear Multistep Methods based on B-splines for the numerical solution of semi-linear second order Boundary Value Problems is introduced. The presented schemes are called BS2 methods, because they are connected to the BS (B-spline) methods previously introduced in the literature to deal with first order problems. We show that, when using an even number of steps, schemes with good general behavior are obtained. In particular, the absolute stability of the 2-step and 4-step BS2 methods is shown. Like BS methods, BS2 methods are of particular interest, because it is possible to associate with the discrete solution a spline extension which collocates the differential equation at the mesh points. Ó 2014 Elsevier Inc. All rights reserved.
Keywords: Boundary Value Problems Linear Multistep Methods Boundary Value Methods B-splines BS methods
1. Introduction Let L be the following linear second order differential functional,
LðyÞðÞ :¼ y00 ðÞ þ 2ly0 ðÞ;
ð1Þ
with l a given real constant. We are interested in numerical methods for solving the associated second order Boundary Value Problem (BVP),
LðyÞðtÞ ¼ f ðt; yÞ;
t 2 ½a; b;
yðaÞ ¼ ya ;
yðbÞ ¼ yb ;
ð2Þ
where f : ½a; b R ! R is a continuous given function and ya ; yb are used to denote the separated boundary conditions. Note that the differential equation considered in (2) is not a general quasi-linear second order equation, but it is more general than a fully linear one because f is not necessarily linear. Note also that problem (2) could be easily reformulated as a first order BVP with doubled dimension, but we want to avoid such a transformation in order to increase the efficiency of the proposed numerical approach, similarly to the strategy in [2,3]. Finally, we want to point out that we are dealing with scalar-valued problems only for easiness of presentation, but no difficulties are expected when extending the proposed methodology to vector-valued problems. ⇑ Corresponding author. E-mail addresses: [email protected] (C. Manni), [email protected] (F. Mazzia), alessandra.sestini@unifi.it (A. Sestini), speleers@mat. uniroma2.it, [email protected] (H. Speleers). http://dx.doi.org/10.1016/j.amc.2014.08.046 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
2
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
In this paper, we develop the so-called BS2 methods to solve numerically problem (2) without doubling its dimension. They are a specific variation of the B-spline (BS) methods, introduced in [10,11], which are a class of Linear Multistep Methods (LMMs) based on B-splines and suited for addressing first order problems of any dimension. Similarly to several other classes of LMMs, in order to get methods with desirable stability properties, it has been shown in [10,11] that BS methods must be used as Boundary Value Methods (BVMs) [6], i.e., they must be combined with a suitable number of additional left and right auxiliary schemes. The peculiarity of BS methods is that they can be interpreted as spline collocation methods. This has the attractive consequence that the corresponding numerical discrete solution can be extended to a continuous spline which collocates the differential problem at the mesh points [9,12]. We emphasize that for the BS and BS2 methods the mesh points and the collocation points are exactly the same. This makes them different from standard spline collocation methods used for the numerical solution of BVPs [3]. The paper is organized as follows. Section 2 introduces the reference problem which has to be used to investigate the stability properties of a numerical scheme for a second order BVP, and in Section 3 the BVM formulation of an LMM for problem (2) is given. Then, Section 4 presents the BS2 methods and analyzes their convergence and stability behavior. In the next section we discuss the computation of the associated collocating spline extension. Section 6 presents some numerical results, and some concluding remarks are collected in Section 7. 2. The reference problem In order to evaluate the stability behavior of any numerical scheme, the differential problem is locally linearized around a stable solution, and the scheme has to produce a stable numerical solution for such a linearized problem. In our situation, we are interested in the following linear homogeneous reference problem,
LðyÞðÞ ¼ cyðÞ;
yðaÞ ¼ ya ;
yðbÞ ¼ yb ;
ð3Þ
with c a given real positive constant, as proposed in [6, Section 9.4]. With some standard computations, it can be easily verified that our reference problem (3) has a unique solution for any value of T :¼ b a. The solution is
yðsÞ ¼ ya f l ðsÞ þ yb f r ðsÞ;
s :¼
ta 2 ½0; 1; ba
ð4Þ
where
f l ðsÞ :¼ gð1 sÞ expðlT sÞ;
f r ðsÞ :¼ gðsÞ expðlTð1 sÞÞ;
and
gðsÞ :¼
sinhðkT sÞ ; sinhðkTÞ
k :¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi l2 þ c:
The functions f l and f r are monotone and vary in ½0; 1, for any T > 0 and s 2 ½0; 1. Therefore, our reference problem is well conditioned with respect to perturbations on the boundary conditions ya and yb , for any T > 0. The positivity of c is necessary and sufficient to ensure that the characteristic polynomial associated with (3) has one positive and one negative real root. This means that the exponential dichotomy assumption holds [3, Theorem 3.106]. Dichotomy is the key assumption for ensuring the well-conditioning of a BVP, not only in the simple case of perturbations on the boundary conditions (as discussed above) but also in the case of perturbations on the coefficients of the problem. We refer to [3, Section 3.4] for the full analysis. Note that the reference problem (3) has the time-reversal symmetry property: when formulated in the new variable x ¼ a þ b t, the solution is just uðxÞ ¼ yða þ b xÞ, where yðtÞ is the solution of (3). 3. BVMs for semi-linear second order BVPs Let us consider a uniform mesh ti :¼ a þ ih; i ¼ 0; . . . ; N, in the integration interval, with h :¼ ðb aÞ=N denoting the mesh size. As usual, the numerical solution is denoted by fyi gNi¼0 (with yi ’ yðti Þ), and define f i :¼ f ðt i ; yi Þ. A k-step LMM for the second order problem (1) can be formulated as follows on a uniform mesh with N P k, k2 X
2 ð1Þ að2Þ iþk1 þ 2hl aiþk1 yjþi ¼ h
i¼k1
k2 X
biþk1 f jþi ;
j ¼ k1 ; . . . ; N k2 ;
ð5Þ
i¼k1 ð2Þ
ð1Þ
where k1 P 1; k2 P 1; k1 þ k2 ¼ k P 2, and where the coefficients ai ; ai and bi ; i ¼ 0; . . . ; k, characterize the main method together with the choice of k1 . For brevity, in the sequel the following notation will be also used:
al :¼ alð2Þ þ 2hlað1Þ l ¼ 0; . . . ; k: l ;
ð6Þ
Note that the following three polynomials of degree k can be associated with the main method, Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
q2 ðzÞ :¼
k X
k X
k X bi zi :
i¼0
i¼0
i¼0
i að2Þ q1 ðzÞ :¼ i z;
aið1Þ zi ; rðzÞ :¼
The method is called symmetric if the polynomials q2 and ð2Þ l
a
ð2Þ kl ;
¼a
ð1Þ l
a
ð1Þ kl ;
¼ a
ð7Þ
r are symmetric and q1 is skew-symmetric, namely
l ¼ 0; . . . ; k:
bl ¼ bkl ;
3
ð8Þ
Since we have N þ 1 unknowns and only N k þ 1 equations in (5), the above main method needs to be completed with k1 P 1 left and k2 P 1 right linear conditions. Obviously, one left and one right condition are directly given by the boundary conditions in (2). This means that, when k > 2, the main method has to be combined with k1 1 additional left methods and k2 1 right ones. For brevity, we report only the general expression of the left ones, kj X
2 ð2;L;jÞ ð1;L;jÞ yjþi ¼ h aiþj þ 2hlaiþj
i¼j
kj X
bL;j iþj f jþi ;
j ¼ 1; . . . ; k1 1;
i¼j ð2; L; jÞ
ð1; L; jÞ
j which are characterized by the coefficients al ; al and bL; l ; l ¼ 0; . . . ; k; j ¼ 1; . . . ; k1 1. In matrix form, the entire LMM corresponds to the following system,
y0 ; yN given;
2 Að2Þ þ 2hlAð1Þ y h Bf ¼ 0; T
where y :¼ ðy0 ; y1 ; . . . ; yN ÞT , f :¼ ðf 0 ; f 1 ; . . . ; f N Þ , with Að2Þ ; Að1Þ and B banded matrices of dimension ðN 1Þ ðN þ 1Þ containing respectively the coefficients að2Þ ; að1Þ and b of the main method and of the additional methods. ^ y, where y ^ :¼ ðyðt0 Þ; yðt1 Þ; . . . ; yðt N ÞÞT , when the mesh size h tends to We are interested in the behavior of the error e :¼ y zero. When assuming for simplicity that f does not depend on y and l ¼ 0, we can verify that e satisfies the following conditions,
e ð2Þ e ¼ ð0; sT ; 0ÞT ; A where s is the vector whose N 1 entries are the local truncation errors of all the left, the main and the right methods. The e ð2Þ is defined as follows, ðN þ 1Þ ðN þ 1Þ matrix A
T e ð2Þ :¼ i1 ; ðAð2Þ ÞT ; iNþ1 ; A with ij being the j-th canonical basis vector. This is a quasi-Toeplitz matrix, because it can be written as the sum of a Toeplitz e ð2Þ plus a low rank matrix with at most k þ 1 non-zero entries in the first k1 and the last k2 rows. Thus, to study the matrix A T e ð2Þ . The size of A e ð2Þ grows when h tends to zero, and the behavior of the error, we need to analyze the norm of the inverse of A e ð2Þ , see [6]. From the theory of Toeplitz matrices (see [1, norm of its inverse behaves like the norm of the inverse of A T ð2Þ 1
e Þ k depends on the roots of the associated polynomial q . Because of consistency requireTheorem 3]), we know that kð A 2 T e ð2Þ Þ1 k is at least Oðh2 Þ when h tends to ments, we have q2 ðzÞ ¼ ðz 1Þ2 v ðzÞ with v ðzÞ of degree k 2. This implies that kð A T e ð2Þ Þ1 k ¼ Oðh2 Þ for k > 2, we need to require that v has k2 1 roots with modulus greater than 1 and zero. To ensure that kð A T the remaining k1 1 with modulus less or equal than 1 but with multiplicity at most 2 if they are on the unit circle of the complex plane. This is known as the generalized zero-stability requirement [6]. Therefore, in the case of generalized zero-stap
pþ2
bility, the error e behaves like Oðh Þ when s ¼ Oðh Þ; p P 1. Note that the same conclusion also holds in more general situations, i.e., when f also depends on y and/or l – 0, see [6]. In addition, it is important to ensure that the discrete solution yi ; i ¼ 0; . . . ; N preserves the behavior of the solution of the continuous reference problem (3) whenever c > 0, see the discussion in Section 2. This results in a well conditioned discrete solution for all h. To this aim, let us represent the discrete solution of the reference problem as follows (see (4)):
yj ¼ M jþ1;1 ya þ M jþ1;Nþ1 yb ;
j ¼ 0; . . . ; N;
where M jþ1;1 and M jþ1;Nþ1 are the entries of the first and the last column of the matrix
0
T
i1
11
C B ð2Þ 2 ð1Þ C M :¼ B @ A þ 2hl A h cB A :
ð9Þ
T iNþ1
Recall that Að2Þ ; Að1Þ and B are a low rank perturbation of a Toeplitz matrix, and thus M is the inverse of a matrix of this kind. It follows that the discrete problem is well conditioned for any h, if the entries of the first column of M decrease from 1 to 0 and those of the last column have the opposite behavior or at most their absolute values are bounded by a constant not depending on h [6]. From the theory of Toeplitz matrices, see [1, Theorem 3], we know that this is satisfied when we have a precise qualitative distribution of the roots of the following polynomial,
pðz; q; gÞ :¼ q2 ðzÞ þ gq1 ðzÞ qrðzÞ;
ð10Þ
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
4
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx 2
where g :¼ 2hl and q :¼ h c. More precisely, for any g and for any positive q, the discrete solution has the required behavior if p has k1 roots with modulus less than 1 and k2 roots with modulus greater than 1. This is known as the generalized absolute stability requirement [6]. A final reasonable requirement for an LMM is its capability to preserve in the discrete case the time-reversal symmetry when applied to (3). In case k is even and k1 ¼ k2 ¼ k=2, this is true if the main method is symmetric, see (8), provided that the additional k2 1 right auxiliary methods are just the reversed version of the left ones. Note that, when k is odd and l ¼ 0, for a symmetric LMM we have q2 ð1Þ ¼ rð1Þ ¼ 0, see (7) and (8). Therefore, pð1; q; 0Þ ¼ 0 and the generalized absolute stability requirement is not fulfilled in this case. This motivates why we will only address symmetric LMMs with k even. 4. BS2 methods The coefficients characterizing the main k-step BS2 method for the second order BVP in (2) are given by ð1Þ að2Þ :¼ B00 ðk þ 1 iÞ; ai :¼ B0 ðk þ 1 iÞ; bi :¼ Bðk þ 1 iÞ; i ¼ 0; . . . ; k; i
ð11Þ
where B is the cardinal B-spline of degree k þ 1 ðk P 2Þ with integer knots 0; . . . ; k þ 2, see [7]. Note that the symmetry of the cardinal B-spline implies that the corresponding main method is symmetric. By considering the definition (5), it can be easily proved that for all l ¼ 1 k; . . . ; N 1 the h-scaled translated spline Bl;h ðtÞ :¼ Bðta lÞ satisfies the following relations, h k2 X
2 ð2Þ aiþk Bl;h ðtjþi Þ ¼ h 1
i¼k1 k2 X
k2 X
biþk1 B00l;h ðtjþi Þ; j ¼ k1 ; . . . ; N k2 ;
i¼k1 ð1Þ iþk1 B l;h ðt jþi Þ
a
¼h
i¼k1
k2 X
biþk1 B0l;h ðt jþi Þ;
ð12Þ j ¼ k1 ; . . . ; N k2 :
i¼k1
Since any spline of degree k þ 1 and knots at the mesh points can be represented in the form
sh ðÞ ¼
N1 X
cl Bl;h ðÞ;
ð13Þ
l¼1k
we can conclude from the linearity that the relations in (12) are satisfied by any function sh of this kind. As a consequence, we can also state that any spline sh of degree k þ 1 and knots at the mesh points satisfies the following relations which involve the differential operator L, k2 X
aiþk1 sh ðtjþi Þ ¼ h2
i¼k1
k2 X
biþk1 Lðsh Þðt jþi Þ;
j ¼ k1 ; . . . ; N k2 ;
ð14Þ
i¼k1
where the coefficients al are defined according to (6). We can conclude that the order p of the k-step BS2 method is greater than or equal to k, because the relations in (14) imply that the equations in (5) are surely verified by any spline of degree k þ 1 with knots at the mesh points and so in particular by any polynomial of degree less than or equal to k þ 1. Since all the BS2 methods are symmetric methods, it follows that they cannot be absolutely stable when an odd number of steps is used and l ¼ 0 (see the last paragraph in the previous section). In the sequel, we will analyze in detail the BS2 methods for k ¼ 2 and k ¼ 4. For both the cases, we choose k1 ¼ k2 ¼ k=2. 4.1. The 2-step BS2 method When k ¼ 2, we consider the cubic cardinal B-spline in (11), and the corresponding coefficients of the main BS2 method are ð2Þ
ð2Þ
ð2Þ
ð1Þ
ð1Þ
ð1Þ
ða0 ; a1 ; a2 Þ ¼ ð1; 2; 1Þ; ða0 ; a1 ; a2 Þ ¼ 12 ð1; 0; 1Þ; ¼ 16 ð1; 4; 1Þ: ðb0 ; b1 ; b2 Þ Since we have q2 ðzÞ ¼ ðz 1Þ2 , it is easy to see that generalized zero-stability is ensured. Clearly, in this case we do not need additional left and right methods because k1 ¼ k2 ¼ 1. Moreover, the 2-step BS2 method possesses the time-reversal symmetry property because the main method is symmetric. One can also verify that the order of the method is p ¼ k ¼ 2. Finally, the stability polynomial defined in (10) and associated with the 2-step BS2 method is the following quadratic polynomial,
pðz; g; qÞ ¼ ð6 þ 3g qÞz2 4ð3 þ qÞz þ ð6 3g qÞ =6; whose discriminant is D ¼ q2 =3 þ 4q þ g2 > 0. So, its roots z1 and z2 are real and they tend to coalesce into 1 when h ! 0. Since pð1; g; qÞ ¼ q < 0 and pð1; g; qÞ ¼ 4 þ q=3 > 0, one root of this polynomial has modulus less than 1. Actually, by considering that it is quadratic, this also implies that the other root has modulus greater than 1. Hence, the 2-step BS2 method is absolutely stable, without any restriction on the mesh size h. Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
5
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
4.2. The 4-step BS2 method When k ¼ 4, we consider the quintic cardinal B-spline in (11), and the corresponding coefficients of the main BS2 method are ð2Þ
ð2Þ
ð2Þ
ð2Þ
ð2Þ
ð1Þ
ð1Þ
ð1Þ
ð1Þ
ð1Þ
ða0 ; a1 ; a2 ; a3 ; a4 Þ ¼ 16 ð1; 2; 6; 2; 1Þ; 1 ð1; 10; 0; 10; 1Þ; ða0 ; a1 ; a2 ; a3 ; a4 Þ ¼ 24 1 ¼ 120 ð1; 26; 66; 26; 1Þ:
ðb0 ; b1 ; b2 ; b3 ; b4 Þ
Generalized zero-stability is ensured in this case, because q2 ðzÞ ¼ ðz 1Þ2 v ðzÞ with v ðzÞ ¼ 16 ðz2 þ 4z þ 1Þ, which has two real negative roots: one with modulus less and the other with modulus greater than 1. The two additional methods are selected by taking the Generalized Backward Differentiation Formulas with 4 steps [6], centered respectively at t 1 and t N1 . The right additional method is the reversed form of the left one, and so the time-reversal symð2;L;1Þ
metry property of the 4-step BS2 method is ensured. For brevity, we report only the coefficients aj
ð1;L;1Þ
; aj
ðL;1Þ
; bj
; j ¼ 0; . . . ; k,
of the left one: ð2;L;1Þ
ða0
ð2;L;1Þ
; a1
ð2;L;1Þ
; a2
ð2;L;1Þ
; a3
ð1;L;1Þ ð1;L;1Þ ð1;L;1Þ ð1;L;1Þ ; 1 ; 2 ; 3 ; ð 0 L;1 L;1 L;1 L;1 L;1 ðb0 ; b1 ; b2 ; b3 ; b4 Þ
a
a
a
a
ð2;L;1Þ
; a4
1 Þ ¼ 12 ð11; 20; 6; 4; 1Þ;
ð1;L;1Þ Þ 4
a
1 ¼ 12 ð3; 10; 18; 6; 1Þ;
¼ ð0; 1; 0; 0; 0Þ:
One can verify that in this case the order is p ¼ k ¼ 4 for the main method, but only p ¼ 3 for the additional methods. This reduced order at the boundary does not lead to a reduction of the order of the full scheme, as mentioned in [6, Remark 4.4.1]. We also refer to the numerical orders reported in Section 6 for the considered examples. The stability polynomial defined in (10) and associated with this method is now the following quartic polynomial,
pðz; g; qÞ ¼ ð20 þ 5g qÞz4 þ ð40 þ 50g 26qÞz3 ð120 þ 66qÞz2 þ ð40 50g 26qÞz þ ð20 5g qÞ =120: We first show (by contradiction) that this polynomial does not have any root on the unit circle of the complex plane. Assuming there exists such a root, it must satisfy pðexpðihÞ; g; qÞ ¼ 0, and we get
q ¼ qðhÞ ¼ q2 ðhÞ þ gq1 ðhÞ; with q2 ðhÞ :¼ q2 ðexpðihÞÞ=rðexpðihÞÞ and q1 ðhÞ :¼ q1 ðexpðihÞÞ=rðexpðihÞÞ, see (10). For fixed g, one can verify that qðhÞ describes in the complex plane a closed curve passing through the origin (because q2 ð1Þ ¼ q1 ð1Þ ¼ 0 and rð1Þ ¼ 1). Moreover, since
1 3 1 q1 ðexpðihÞÞ ¼ i expði2hÞ½ sinð2hÞ 10 sinðhÞ; 12 1 rðexpðihÞÞ ¼ expði2hÞ½cosð2hÞ þ 26 cosðhÞ þ 33; 60
q2 ðexpðihÞÞ ¼ expði2hÞ½cosð2hÞ þ 2 cosðhÞ 3;
the curve qðhÞ belongs to the half-plane ReðqÞ 6 0. However, from its definition it follows that q is a real positive value (see (10)), so the stability polynomial cannot have a root on the unit circle of the complex plane. Furthermore, when considering the roots of the stability polynomial as continuous functions of h, they tend to the roots of q2 for h ! 0. We know that q2 has two real negative roots, one strictly inside and the other strictly outside the unit circle. Hence, we can conclude that for any h at least one root of p is strictly inside the unit circle and another is strictly outside. On the other hand, because pð1; g; qÞ ¼ 2ð10 þ qÞ=15 < 0 and pð1; g; qÞ ¼ q < 0, we can say that p must have an even number of real roots in ð1; 1Þ, counted with their multiplicity. This means that the stability polynomial must have exactly two complex conjugate roots or two real roots inside the unit circle. In any case, we can conclude that the other two roots are outside the unit circle. This implies that our 4-step BS2 method is absolutely stable. 5. The continuous spline extension Like for BS methods, it is possible to extend the discrete solution generated by the k-step BS2 method with a spline of degree k þ 1 and knots at the mesh points. This spline collocates the differential equation at the mesh points. However, in the BS2 case, in order to guarantee a unique continuous spline extension, it will be shown that special combinations of h and l have to be avoided. Under this assumption, the coefficients of such a spline in the B-spline form can be computed in a similar way as for the BS methods [10]. For brevity reasons we only summarize the main idea, and we refer to [10] for more details. We look for a spline sh of degree k þ 1 and knots at the mesh points such that
sh ðti Þ ¼ yi ;
Lðsh Þðt i Þ ¼ f i ;
i ¼ 0; . . . ; N;
ð15Þ
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
6
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
where we can represent sh in the B-spline form introduced in (13). Then, we have to determine the coefficient vector c :¼ ðc1k ; . . . ; cN1 ÞT such that 2
T
2
Sc ¼ ðy0 ; . . . ; yN ; h f 0 ; . . . ; h f N Þ ;
S :¼
SI SL
;
ð16Þ
where SI and SL are banded Toeplitz matrices of dimension ðN þ 1Þ ðN þ k þ 1Þ with bandwidth k, and where ðSI Þ1;j :¼ bj1 ð2Þ
ð1Þ
and ðSL Þ1;j :¼ aj1 þ 2hlaj1 ; j ¼ 1; . . . ; k þ 1. The collocation system in (16) is overdetermined. The following proposition ensures when it has a unique solution. Proposition 1. The rectangular linear system in (16) has a unique solution if rankðSÞ ¼ N þ k þ 1 and if the entries of 2
T
2
ð2Þ
ð1Þ
ðy0 ; . . . ; yN ; h f 0 ; . . . ; h f N Þ satisfy (5) where the coefficients ai ; ai ; bi ; i ¼ 0; . . . ; k, are given in (11). Proof. Since the proof is similar to the one in [10, Theorem 6], we only give a brief outline of it. The vector 2
T
2
ðy0 ; . . . ; yN ; h f 0 ; . . . ; h f N Þ satisfies the N k þ 1 linear conditions given in (5). On the other hand, from (12)–(14) we know that for the l-th column of S; l ¼ 1; . . . ; N þ k þ 1, we have k2 X
aiþk1 ðSI Þjþi;l ¼
i¼k1
k2 X
biþk1 ðSL Þjþi;l ;
j ¼ k1 ; . . . ; N k2 :
i¼k1
It follows that there are N k þ 1 redundancies in the rectangular system (16), and this implies that there exists a unique solution since S is a full rank matrix of dimension ð2N þ 2Þ ðN þ k þ 1Þ. h 2
2
T
If the vector ðy0 ; . . . ; yN ; h f 0 ; . . . ; h f N Þ is generated by the k-step BS2 method, then (5) is satisfied. Therefore, Proposition 1 implies that the collocation system in (16) has a unique solution whenever the matrix S has full rank. In the next proposition we analyze the rank of S. Proposition 2. Let Dk be the 2k 2k matrix, where ðDk Þi;j ¼ ðSI Þi;j and ðDk Þkþi;j ¼ ðSL Þi;j for i ¼ 1; . . . ; k and j ¼ 1; . . . ; 2k. The rectangular matrix S of the linear system in (16) has full rank (i.e. N þ k þ 1) if and only if Dk is non-singular. Proof. We will show that the N þ k þ 1 columns of the matrix S are linearly independent. To this end, it suffices to prove that Sd ¼ 0 if and only if the vector d is the zero-vector. Let us start with considering the first 2k columns of S. From the sparsity structure of S and the definition of Dk , we can conclude that these columns are linearly independent if and only if Dk is nonsingular. Moreover, because all the entries ðSI Þi;2kþ1 ¼ . . . ¼ ðSI Þi;Nþkþ1 ¼ 0 and ðSL Þi;2kþ1 ¼ . . . ¼ ðSL Þi;Nþkþ1 ¼ 0; i ¼ 1; . . . ; k, it follows that the first 2k entries of d must be zero to obtain Sd ¼ 0. In a similar way, when considering any set of 2k successive columns of S, we obtain that the corresponding 2k entries of d must be zero, which concludes the proof. h ^2 and h R fh ^ ; h ^ þ g respectively, where In the case of k ¼ 2 and k ¼ 4, the condition in Proposition 2 is satisfied when h – h 4 4
^2 ¼ h
pffiffiffi 3 ; jlj
^ ¼ h 4
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi 345 29 105 ; 8jlj
l – 0;
ð17Þ
and there are no restrictions on h if l ¼ 0. This can be verified by a direct computation. More precisely, recalling g :¼ 2hl, the submatrix Dk described in the proposition is given by
2
1=6 2=3 60 1=6 6 D2 ¼ 6 4 1 g=2 2
2
1=5
60 6 6 60 6 1 6 60 D4 ¼ 6 24 6 4 g 6 60 6 6 40 0
26=5
0
2=3 1=6 1 þ g=2 0
1 g=2 2
0 and
1=6
66=5
3
7 7 7; 5
1 þ g=2 26=5
1=5
0
0
0
3
7 7 7 7 0 1=5 26=5 66=5 26=5 1=5 0 7 0 0 1=5 26=5 66=5 26=5 1=5 7 7 7; 7 8 10g 24 8 þ 10g 4 þ g 0 0 0 7 7 4g 8 10g 24 8 þ 10g 4 þ g 0 0 7 7 5 0 4g 8 10g 24 8 þ 10g 4 þ g 0 8 10g 24 8 þ 10g 4 þ g 0 0 4g 1=5
26=5
66=5
26=5
1=5
0
0
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
7
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
with
detðD2 Þ ¼
g2 12
þ 1;
detðD4 Þ ¼
4 1 g 23g2 þ1 : 32400 120 64
Note that this restriction on the allowed mesh sizes is required only to compute the continuous extension and does not affect the stability properties of the BS2 methods. Moreover, when l 0, usually the problem is singularly perturbed. In such a case, the use of a numerical scheme based on a uniform mesh size is not very efficient and a mesh size variation strategy should be considered. To compute the coefficient vector c of the unique spline solution sh of (15), we could solve (16) in the least-square sense, or use a quasi-interpolation approach similar to the one for BS methods described in [9]. 6. Numerical results In this section we present some numerical results for two linear problems with f ðt; yÞ ¼ cy þ gðtÞ, both depending on a parameter > 0 which determines their stiffness level and both taken from the test set for BVP solvers at the site https:// www.dm.uniba.it/bvpsolvers/testsetbvpsolvers. Furthermore, we consider a standard harmonic oscillator problem – already investigated in [10] – to confirm the computational gain which can be obtained for second order problems by using the BS2 methods instead of the BS ones. 6.1. Two linear problems with boundary layers We consider two linear problems depending on a parameter > 0 for which we have selected two different values, say 1 and 2 , in order to deal once with a non-stiff problem and once with a quite stiff one. The BS2 methods with k ¼ 2 and k ¼ 4
Table 1 Problem 1 with ¼ 1 ¼ 4 102 (top) and ¼ 2 ¼ 4 104 (bottom). The maximal values of the errors at the mesh points (err m ), the related numerical approximation orders, and the maximal spline approximation errors (err s ) are given. N
1
2
k¼2
k¼4
err m
Order
err s
err m
Order
err s
10
1:26 102
/
1:31 102
4:58 103
/
8:33 103
20
3:06 103
2:04
3:08 103
6:43 104
2:83
9:55 104
40
4
7:65 10
2:00
4
7:66 10
5
3:94 10
4:03
5:29 105
80
1:91 104
2:00
1:91 104
1:70 106
4:53
2:18 106
50
6:35 102
/
1:01 101
2:79 102
/
7:95 102
100
1:69 102
1:90
1:74 102
7:14 103
1:97
1:31 102
200
3
3:92 10
2:11
3
3:94 10
4
7:17 10
3:32
1:07 103
400
9:63 104
2:03
9:63 104
4:08 105
4:14
5:05 105
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
−0.2
−0.2
−0.4
−0.4
−0.6
−0.6
−0.8
−0.8
−1 −1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
−1 −1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
4
Fig. 1. Problem 1 with ¼ 2 ¼ 4 10 and N ¼ 50. The exact solution (blue) is compared with the spline extension (red) of the numerical solution produced by the BS2 methods with k ¼ 2 (left) and k ¼ 4 (right). The numerical solution is depicted with the symbol. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
8
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
Table 2 Problem 2 with ¼ 1 ¼ 2 (top) and ¼ 2 ¼ 0:02 (bottom). The maximal values of the errors at the mesh points (errm ), the related numerical approximation orders, and the maximal spline approximation errors (errs ) are given. N
k¼2
k¼4 Order
err m
1
2
err s
err m
Order
err s
5
7:54 10
/
7:78 10
1:10 10
/
1:85 103
10
1:85 103
2:03
1:87 103
8:14 105
3:76
1:14 104
20
4:64 104
2:00
4:64 104
3:82 106
4:41
4:99 106
40
4
2:00
4
1:46 10
7
4:71
1:84 107
1
3
3
1:16 10
3
1:16 10
3
50
1:37 10
/
2:52 10
6:65 10
/
6:23 102
100
3:53 102
1:94
3:60 102
4:18 103
0:66
1:23 102
200
8:04 103
2:13
8:07 103
5:71 104
2:87
9:74 104
400
3
2:03
3
5
3:90
4:92 105
1:97 10
1
3:82 10
1:97 10
have been tested by successively halving the mesh size h, starting from a suitably selected mesh size h0 . Note that for all the performed tests we have verified that the infinity norm of both the first and last columns of the matrix M introduced in (9) is equal to 1. The first problem (BVPT14) is given by
y00 ¼ y ð1 þ p2 Þ cosðptÞ;
t 2 ½1; 1;
yð1Þ ¼ yð1Þ ¼ expð2=
pffiffiffi Þ;
pffiffiffi pffiffiffi pffiffiffi whose exact solution is yðtÞ ¼ cosðptÞ þ expððt 1Þ= Þ þ expððt þ 1Þ= Þ. It has a boundary layer of width Oð Þ near t ¼ 1. Note that in this case we have l ¼ 0, and so the spline extension of the solution produced by the BS2 methods pffiffiffi pffiffiffi can be defined for any mesh size. For the related experiments we have set h0 ¼ when ¼ 1 ¼ 4 102 and h0 ¼ 2 4 when ¼ 2 ¼ 4 10 . Table 1 reports the maximal errors at the mesh points, the corresponding numerical approximation orders, and the maximal spline approximation errors (sampled at 2000 uniformly spaced points) obtained with the two studied BS2 methods. The results confirm that these methods have order 2 and 4 respectively. Fig. 1 illustrates the case ¼ 2 with N ¼ 50, and it can be seen that only near the boundary layers the exact solution slightly differs from the spline extension sh of the numerical solution produced by the BS2 methods. The second problem (BVPT4) is
y00 þ y0 ¼ ð1 þ Þy;
t 2 ½1; 1;
yð1Þ ¼ 1 þ e2 ;
yð1Þ ¼ 1 þ e2ð1þÞ= ;
whose exact solution is yðtÞ ¼ et1 þ eðtþ1Þð1þÞ= . It has a boundary layer of width OðÞ near t ¼ 1. In view of our notation, l is now equal to 1=ð2Þ, and the experiments have been done for ¼ i ; i ¼ 1; 2 with 1 ¼ 2 and 2 ¼ 2 102 . The maximal mesh size is h0 ¼ 0:4 (N ¼ 5) in the first case (the non-stiff case) and h0 ¼ 2 (N ¼ 50) in the second case. We note that the critical step sizes given in (17) are avoided, and so the spline extensions can be computed. The numerical results are collected in Table 2, which is structured in the same way as Table 1. Fig. 2 shows that the methods are capable of producing a nice continuous approximation of the solution, even when very few points are used for non-stiff problems.
1.3
1.3
1.2
1.2
1.1
1.1
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5 −1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
0.5 −1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
Fig. 2. Problem 2 with ¼ 1 ¼ 2 and N ¼ 5. The exact solution (blue) is compared with the spline extension (red) of the numerical solution produced by the BS2 methods with k ¼ 2 (left) and k ¼ 4 (right). The numerical solution is depicted with the symbol. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
9
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
Table 3 Comparison between the method implemented in the solver colnew and BS2 methods for Problem 1 (top) and 2 (bottom). The maximal values of the errors at the collocation points (err m ) and the maximal spline approximation errors (errs ) are given for the two methods using the same number of collocation points (N c ). Nc
Problem 1, 8
BS2, k ¼ 4
colnew err m
err s
err m
err s
1 ¼ 0:04 6:99 103
1:41 102
4:74 103
9:96 103
16
4:39 104
2:32 103
1:69 103
2:69 103
32
2:72 105
2:22 104
1:16 104
1:61 104
64
1:69 106
1:72 105
5:15 106
6:18 106
Problem 2, 8
1 ¼ 2 4:98 105
5:60 104
3:39 104
5:10 104
16
3:04 106
4:29 105
1:41 106
1:88 105
32
1:93 107
2:97 106
4:93 107
6:28 107
64
1:20 108
1:95 107
1:62 108
2:02 108
5
5
4
4
3
3
2
2
1
1
0
0
−1
−1
−2
−2
−3
−3
−4
−4
−5 −5
−4
−3
−2
−1
0
1
2
3
4
5
−5 −5
5
5
4
4
3
3
2
2
1
1
0
0
−1
−1
−2
−2
−3
−3
−4
−4
−5 −5
−4
−3
−2
−1
0
1
2
3
4
5
−5 −5
−4
−3
−2
−1
0
1
2
3
4
5
−4
−3
−2
−1
0
1
2
3
4
5
Fig. 3. The oscillator test problem given in (18). The exact solution (green) is compared with the spline extension (red) of the BS2 solution in the phase plane for different values of N : N ¼ 10 (top-left), N ¼ 15 (top-right), N ¼ 20 (bottom-left), N ¼ 25 (bottom-right). The numerical solution is produced by the BS2 method with k ¼ 4. The symbol identifies the mesh points. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
We now illustrate the difference between our method and standard collocation applied to Problems 1 and 2 with ¼ 1 . To this aim, we compare the 4-step BS2 method with the method implemented in the solver colnew [4,5] collocating at two Gaussian points in each subinterval.1 They are both fourth order schemes and do not convert the considered differential problem into a system of first order differential equations. More specifically, for the tests we have employed the function bvpcol
1 The solver colnew generates adaptively the mesh, and so the value ¼ 2 could not be used because we want to deal only with uniform meshes for a fair comparison.
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
10
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
implemented in the R package bvpSolve [8]. The results are reported in Table 3, where the two methods are chosen such that they use the same number of collocation points. The table confirms that the behavior of the two methods is similar. The absolute error at the collocation points is mostly smaller for collocation at Gaussian points, whereas the opposite is true for the error in the associated spline. Note that the BS2 method with k ¼ 4 gives a quintic C 4 spline while colnew provides just a cubic C 1 spline. 6.2. A harmonic oscillator problem In this section we show that the BS2 methods preserve the high quality spline approximations of the solutions of second order BVPs produced by the BS methods, but with a reduced computational cost. We consider the harmonic oscillator test problem,
y00 þ x2 y ¼ 0;
t 2 ½0; p=2;
yð0Þ ¼ 3;
yðp=2Þ ¼ 3;
ð18Þ
with x ¼ 5, whose exact solution is yðtÞ ¼ 3ðcos xt sin xtÞ. Let us denote by sh the quintic spline approximation of y produced by the BS2 method with k ¼ 4. Its qualitative behavior can be analyzed by looking at Fig. 3: for different h steps, the 0 0 curve pffiffiffi Xh ðtÞ :¼ ðsh ðtÞ; sh ðtÞ=xÞ is compared with the curve XðtÞ :¼ ðyðtÞ; y ðtÞ=xÞ which is a circumference with radius equal to 3 2. Note that yðtÞ is also the solution of the oscillator test problem considered in Problem 2 in [10], where the BS methods were used to solve numerically the above second order equation reformulated as a bidimensional first order problem on the integration domain ½0; 2p and with boundary conditions fixing y at t ¼ 2p and its derivative at t ¼ 0. Fig. 3 may be compared with Fig. 6.1 in [10] because in both cases the results are obtained using a numerical scheme of fourth order. The comparison shows that the BS2 method attains the same approximation quality with almost half the computational cost of the BS method. Indeed, for example with N ¼ 20, here we are using 21 scalar degrees of freedom; on the other hand, the same accuracy was obtained in [10] with N ¼ 80 on a four times larger integration interval, that is with 162 scalar degrees of freedom. 7. Conclusions We have proposed a class of k-step LMMs based on B-splines (the so-called BS2 methods) to solve numerically the second order BVP given in (2) without doubling its dimension. These methods are a specific formulation of the BS methods introduced in [10]. When k is even, they are symmetric and of order k. Moreover, we have shown that the generalized zero-stability and absolute stability requirements are satisfied for k ¼ 2 and k ¼ 4. Like for BS methods, it is possible to extend the discrete solution generated by the k-step BS2 method with a spline of degree k þ 1 and knots at the mesh points. This spline collocates the differential problem at the mesh points and it can be always defined except for few particular values of the mesh size. Finally, we want to point out that for stiff problems, in order to avoid very fine uniform meshes, it could be better to consider a non-uniform extension of these methods. This topic is currently under investigation. References [1] P. Amodio, L. Brugnano, The conditioning of Toeplitz band matrices, Math. Comput. Model. 23 (1996) 29–42. [2] P. Amodio, I. Sgura, High-order finite difference schemes for the solution of second-order BVPs, J. Comput. Appl. Math. 176 (2005) 59–76. [3] U.M. Ascher, R.M.M. Mattheij, R.D. Russel, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, Prentice Hall, Englewood Cliffs, New York, 1988. [4] U.M. Ascher, J. Christiansen, R.D. Russell, Collocation software for boundary-value ODEs, ACM Trans. Math. Softw. 7 (1981) 209–222. [5] G. Bader, U.M. Ascher, A new basis implementation for a mixed order boundary value ODE solver, SIAM J. Sci. Stat. Comput. 8 (1987) 483–500. [6] L. Brugnano, D. Trigiante, Solving Differential Problems by Multistep Initial and Boundary Value Methods, Gordon and Breach Science Publishers, Amsterdam, 1998. [7] C. de Boor, A Practical Guide to Splines, Revised Edition., Springer-Verlag, New York, 2001. [8] F. Mazzia, J.R. Cash, K. Soetaert, Solving boundary value problems in the open source software R: package bvpSolve, Opuscula Math. 34 (2014) 387–403. [9] F. Mazzia, A. Sestini, The BS class of Hermite spline quasi-interpolants on nonuniform knot distributions, BIT Numer. Math. 49 (2009) 611–628. [10] F. Mazzia, A. Sestini, D. Trigiante, B-spline multistep methods and their continuous extensions, SIAM J. Numer. Anal. 44 (2006) 1954–1973. [11] F. Mazzia, A. Sestini, D. Trigiante, BS linear multistep methods on non-uniform meshes, J. Numer. Anal. Ind. Appl. Math. 1 (2006) 131–144. [12] F. Mazzia, A. Sestini, D. Trigiante, The continuous extension of the B-spline linear multistep methods for BVPs on non-uniform meshes, Appl. Numer. Math. 59 (2009) 723–738.
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
BS2 methods for semi-linear second order boundary value problems Carla Manni a, Francesca Mazzia b, Alessandra Sestini c,⇑, Hendrik Speleers a,d a
Dipartimento di Matematica, Università di Roma ‘‘Tor Vergata’’, Italy Dipartimento di Matematica, Università di Bari, Italy c Dipartimento di Matematica e Informatica, Università di Firenze, Italy d Departement Computerwetenschappen, Katholieke Universiteit Leuven, Belgium b
a r t i c l e
i n f o
a b s t r a c t A new class of Linear Multistep Methods based on B-splines for the numerical solution of semi-linear second order Boundary Value Problems is introduced. The presented schemes are called BS2 methods, because they are connected to the BS (B-spline) methods previously introduced in the literature to deal with first order problems. We show that, when using an even number of steps, schemes with good general behavior are obtained. In particular, the absolute stability of the 2-step and 4-step BS2 methods is shown. Like BS methods, BS2 methods are of particular interest, because it is possible to associate with the discrete solution a spline extension which collocates the differential equation at the mesh points. Ó 2014 Elsevier Inc. All rights reserved.
Keywords: Boundary Value Problems Linear Multistep Methods Boundary Value Methods B-splines BS methods
1. Introduction Let L be the following linear second order differential functional,
LðyÞðÞ :¼ y00 ðÞ þ 2ly0 ðÞ;
ð1Þ
with l a given real constant. We are interested in numerical methods for solving the associated second order Boundary Value Problem (BVP),
LðyÞðtÞ ¼ f ðt; yÞ;
t 2 ½a; b;
yðaÞ ¼ ya ;
yðbÞ ¼ yb ;
ð2Þ
where f : ½a; b R ! R is a continuous given function and ya ; yb are used to denote the separated boundary conditions. Note that the differential equation considered in (2) is not a general quasi-linear second order equation, but it is more general than a fully linear one because f is not necessarily linear. Note also that problem (2) could be easily reformulated as a first order BVP with doubled dimension, but we want to avoid such a transformation in order to increase the efficiency of the proposed numerical approach, similarly to the strategy in [2,3]. Finally, we want to point out that we are dealing with scalar-valued problems only for easiness of presentation, but no difficulties are expected when extending the proposed methodology to vector-valued problems. ⇑ Corresponding author. E-mail addresses: [email protected] (C. Manni), [email protected] (F. Mazzia), alessandra.sestini@unifi.it (A. Sestini), speleers@mat. uniroma2.it, [email protected] (H. Speleers). http://dx.doi.org/10.1016/j.amc.2014.08.046 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
2
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
In this paper, we develop the so-called BS2 methods to solve numerically problem (2) without doubling its dimension. They are a specific variation of the B-spline (BS) methods, introduced in [10,11], which are a class of Linear Multistep Methods (LMMs) based on B-splines and suited for addressing first order problems of any dimension. Similarly to several other classes of LMMs, in order to get methods with desirable stability properties, it has been shown in [10,11] that BS methods must be used as Boundary Value Methods (BVMs) [6], i.e., they must be combined with a suitable number of additional left and right auxiliary schemes. The peculiarity of BS methods is that they can be interpreted as spline collocation methods. This has the attractive consequence that the corresponding numerical discrete solution can be extended to a continuous spline which collocates the differential problem at the mesh points [9,12]. We emphasize that for the BS and BS2 methods the mesh points and the collocation points are exactly the same. This makes them different from standard spline collocation methods used for the numerical solution of BVPs [3]. The paper is organized as follows. Section 2 introduces the reference problem which has to be used to investigate the stability properties of a numerical scheme for a second order BVP, and in Section 3 the BVM formulation of an LMM for problem (2) is given. Then, Section 4 presents the BS2 methods and analyzes their convergence and stability behavior. In the next section we discuss the computation of the associated collocating spline extension. Section 6 presents some numerical results, and some concluding remarks are collected in Section 7. 2. The reference problem In order to evaluate the stability behavior of any numerical scheme, the differential problem is locally linearized around a stable solution, and the scheme has to produce a stable numerical solution for such a linearized problem. In our situation, we are interested in the following linear homogeneous reference problem,
LðyÞðÞ ¼ cyðÞ;
yðaÞ ¼ ya ;
yðbÞ ¼ yb ;
ð3Þ
with c a given real positive constant, as proposed in [6, Section 9.4]. With some standard computations, it can be easily verified that our reference problem (3) has a unique solution for any value of T :¼ b a. The solution is
yðsÞ ¼ ya f l ðsÞ þ yb f r ðsÞ;
s :¼
ta 2 ½0; 1; ba
ð4Þ
where
f l ðsÞ :¼ gð1 sÞ expðlT sÞ;
f r ðsÞ :¼ gðsÞ expðlTð1 sÞÞ;
and
gðsÞ :¼
sinhðkT sÞ ; sinhðkTÞ
k :¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi l2 þ c:
The functions f l and f r are monotone and vary in ½0; 1, for any T > 0 and s 2 ½0; 1. Therefore, our reference problem is well conditioned with respect to perturbations on the boundary conditions ya and yb , for any T > 0. The positivity of c is necessary and sufficient to ensure that the characteristic polynomial associated with (3) has one positive and one negative real root. This means that the exponential dichotomy assumption holds [3, Theorem 3.106]. Dichotomy is the key assumption for ensuring the well-conditioning of a BVP, not only in the simple case of perturbations on the boundary conditions (as discussed above) but also in the case of perturbations on the coefficients of the problem. We refer to [3, Section 3.4] for the full analysis. Note that the reference problem (3) has the time-reversal symmetry property: when formulated in the new variable x ¼ a þ b t, the solution is just uðxÞ ¼ yða þ b xÞ, where yðtÞ is the solution of (3). 3. BVMs for semi-linear second order BVPs Let us consider a uniform mesh ti :¼ a þ ih; i ¼ 0; . . . ; N, in the integration interval, with h :¼ ðb aÞ=N denoting the mesh size. As usual, the numerical solution is denoted by fyi gNi¼0 (with yi ’ yðti Þ), and define f i :¼ f ðt i ; yi Þ. A k-step LMM for the second order problem (1) can be formulated as follows on a uniform mesh with N P k, k2 X
2 ð1Þ að2Þ iþk1 þ 2hl aiþk1 yjþi ¼ h
i¼k1
k2 X
biþk1 f jþi ;
j ¼ k1 ; . . . ; N k2 ;
ð5Þ
i¼k1 ð2Þ
ð1Þ
where k1 P 1; k2 P 1; k1 þ k2 ¼ k P 2, and where the coefficients ai ; ai and bi ; i ¼ 0; . . . ; k, characterize the main method together with the choice of k1 . For brevity, in the sequel the following notation will be also used:
al :¼ alð2Þ þ 2hlað1Þ l ¼ 0; . . . ; k: l ;
ð6Þ
Note that the following three polynomials of degree k can be associated with the main method, Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
q2 ðzÞ :¼
k X
k X
k X bi zi :
i¼0
i¼0
i¼0
i að2Þ q1 ðzÞ :¼ i z;
aið1Þ zi ; rðzÞ :¼
The method is called symmetric if the polynomials q2 and ð2Þ l
a
ð2Þ kl ;
¼a
ð1Þ l
a
ð1Þ kl ;
¼ a
ð7Þ
r are symmetric and q1 is skew-symmetric, namely
l ¼ 0; . . . ; k:
bl ¼ bkl ;
3
ð8Þ
Since we have N þ 1 unknowns and only N k þ 1 equations in (5), the above main method needs to be completed with k1 P 1 left and k2 P 1 right linear conditions. Obviously, one left and one right condition are directly given by the boundary conditions in (2). This means that, when k > 2, the main method has to be combined with k1 1 additional left methods and k2 1 right ones. For brevity, we report only the general expression of the left ones, kj X
2 ð2;L;jÞ ð1;L;jÞ yjþi ¼ h aiþj þ 2hlaiþj
i¼j
kj X
bL;j iþj f jþi ;
j ¼ 1; . . . ; k1 1;
i¼j ð2; L; jÞ
ð1; L; jÞ
j which are characterized by the coefficients al ; al and bL; l ; l ¼ 0; . . . ; k; j ¼ 1; . . . ; k1 1. In matrix form, the entire LMM corresponds to the following system,
y0 ; yN given;
2 Að2Þ þ 2hlAð1Þ y h Bf ¼ 0; T
where y :¼ ðy0 ; y1 ; . . . ; yN ÞT , f :¼ ðf 0 ; f 1 ; . . . ; f N Þ , with Að2Þ ; Að1Þ and B banded matrices of dimension ðN 1Þ ðN þ 1Þ containing respectively the coefficients að2Þ ; að1Þ and b of the main method and of the additional methods. ^ y, where y ^ :¼ ðyðt0 Þ; yðt1 Þ; . . . ; yðt N ÞÞT , when the mesh size h tends to We are interested in the behavior of the error e :¼ y zero. When assuming for simplicity that f does not depend on y and l ¼ 0, we can verify that e satisfies the following conditions,
e ð2Þ e ¼ ð0; sT ; 0ÞT ; A where s is the vector whose N 1 entries are the local truncation errors of all the left, the main and the right methods. The e ð2Þ is defined as follows, ðN þ 1Þ ðN þ 1Þ matrix A
T e ð2Þ :¼ i1 ; ðAð2Þ ÞT ; iNþ1 ; A with ij being the j-th canonical basis vector. This is a quasi-Toeplitz matrix, because it can be written as the sum of a Toeplitz e ð2Þ plus a low rank matrix with at most k þ 1 non-zero entries in the first k1 and the last k2 rows. Thus, to study the matrix A T e ð2Þ . The size of A e ð2Þ grows when h tends to zero, and the behavior of the error, we need to analyze the norm of the inverse of A e ð2Þ , see [6]. From the theory of Toeplitz matrices (see [1, norm of its inverse behaves like the norm of the inverse of A T ð2Þ 1
e Þ k depends on the roots of the associated polynomial q . Because of consistency requireTheorem 3]), we know that kð A 2 T e ð2Þ Þ1 k is at least Oðh2 Þ when h tends to ments, we have q2 ðzÞ ¼ ðz 1Þ2 v ðzÞ with v ðzÞ of degree k 2. This implies that kð A T e ð2Þ Þ1 k ¼ Oðh2 Þ for k > 2, we need to require that v has k2 1 roots with modulus greater than 1 and zero. To ensure that kð A T the remaining k1 1 with modulus less or equal than 1 but with multiplicity at most 2 if they are on the unit circle of the complex plane. This is known as the generalized zero-stability requirement [6]. Therefore, in the case of generalized zero-stap
pþ2
bility, the error e behaves like Oðh Þ when s ¼ Oðh Þ; p P 1. Note that the same conclusion also holds in more general situations, i.e., when f also depends on y and/or l – 0, see [6]. In addition, it is important to ensure that the discrete solution yi ; i ¼ 0; . . . ; N preserves the behavior of the solution of the continuous reference problem (3) whenever c > 0, see the discussion in Section 2. This results in a well conditioned discrete solution for all h. To this aim, let us represent the discrete solution of the reference problem as follows (see (4)):
yj ¼ M jþ1;1 ya þ M jþ1;Nþ1 yb ;
j ¼ 0; . . . ; N;
where M jþ1;1 and M jþ1;Nþ1 are the entries of the first and the last column of the matrix
0
T
i1
11
C B ð2Þ 2 ð1Þ C M :¼ B @ A þ 2hl A h cB A :
ð9Þ
T iNþ1
Recall that Að2Þ ; Að1Þ and B are a low rank perturbation of a Toeplitz matrix, and thus M is the inverse of a matrix of this kind. It follows that the discrete problem is well conditioned for any h, if the entries of the first column of M decrease from 1 to 0 and those of the last column have the opposite behavior or at most their absolute values are bounded by a constant not depending on h [6]. From the theory of Toeplitz matrices, see [1, Theorem 3], we know that this is satisfied when we have a precise qualitative distribution of the roots of the following polynomial,
pðz; q; gÞ :¼ q2 ðzÞ þ gq1 ðzÞ qrðzÞ;
ð10Þ
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
4
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx 2
where g :¼ 2hl and q :¼ h c. More precisely, for any g and for any positive q, the discrete solution has the required behavior if p has k1 roots with modulus less than 1 and k2 roots with modulus greater than 1. This is known as the generalized absolute stability requirement [6]. A final reasonable requirement for an LMM is its capability to preserve in the discrete case the time-reversal symmetry when applied to (3). In case k is even and k1 ¼ k2 ¼ k=2, this is true if the main method is symmetric, see (8), provided that the additional k2 1 right auxiliary methods are just the reversed version of the left ones. Note that, when k is odd and l ¼ 0, for a symmetric LMM we have q2 ð1Þ ¼ rð1Þ ¼ 0, see (7) and (8). Therefore, pð1; q; 0Þ ¼ 0 and the generalized absolute stability requirement is not fulfilled in this case. This motivates why we will only address symmetric LMMs with k even. 4. BS2 methods The coefficients characterizing the main k-step BS2 method for the second order BVP in (2) are given by ð1Þ að2Þ :¼ B00 ðk þ 1 iÞ; ai :¼ B0 ðk þ 1 iÞ; bi :¼ Bðk þ 1 iÞ; i ¼ 0; . . . ; k; i
ð11Þ
where B is the cardinal B-spline of degree k þ 1 ðk P 2Þ with integer knots 0; . . . ; k þ 2, see [7]. Note that the symmetry of the cardinal B-spline implies that the corresponding main method is symmetric. By considering the definition (5), it can be easily proved that for all l ¼ 1 k; . . . ; N 1 the h-scaled translated spline Bl;h ðtÞ :¼ Bðta lÞ satisfies the following relations, h k2 X
2 ð2Þ aiþk Bl;h ðtjþi Þ ¼ h 1
i¼k1 k2 X
k2 X
biþk1 B00l;h ðtjþi Þ; j ¼ k1 ; . . . ; N k2 ;
i¼k1 ð1Þ iþk1 B l;h ðt jþi Þ
a
¼h
i¼k1
k2 X
biþk1 B0l;h ðt jþi Þ;
ð12Þ j ¼ k1 ; . . . ; N k2 :
i¼k1
Since any spline of degree k þ 1 and knots at the mesh points can be represented in the form
sh ðÞ ¼
N1 X
cl Bl;h ðÞ;
ð13Þ
l¼1k
we can conclude from the linearity that the relations in (12) are satisfied by any function sh of this kind. As a consequence, we can also state that any spline sh of degree k þ 1 and knots at the mesh points satisfies the following relations which involve the differential operator L, k2 X
aiþk1 sh ðtjþi Þ ¼ h2
i¼k1
k2 X
biþk1 Lðsh Þðt jþi Þ;
j ¼ k1 ; . . . ; N k2 ;
ð14Þ
i¼k1
where the coefficients al are defined according to (6). We can conclude that the order p of the k-step BS2 method is greater than or equal to k, because the relations in (14) imply that the equations in (5) are surely verified by any spline of degree k þ 1 with knots at the mesh points and so in particular by any polynomial of degree less than or equal to k þ 1. Since all the BS2 methods are symmetric methods, it follows that they cannot be absolutely stable when an odd number of steps is used and l ¼ 0 (see the last paragraph in the previous section). In the sequel, we will analyze in detail the BS2 methods for k ¼ 2 and k ¼ 4. For both the cases, we choose k1 ¼ k2 ¼ k=2. 4.1. The 2-step BS2 method When k ¼ 2, we consider the cubic cardinal B-spline in (11), and the corresponding coefficients of the main BS2 method are ð2Þ
ð2Þ
ð2Þ
ð1Þ
ð1Þ
ð1Þ
ða0 ; a1 ; a2 Þ ¼ ð1; 2; 1Þ; ða0 ; a1 ; a2 Þ ¼ 12 ð1; 0; 1Þ; ¼ 16 ð1; 4; 1Þ: ðb0 ; b1 ; b2 Þ Since we have q2 ðzÞ ¼ ðz 1Þ2 , it is easy to see that generalized zero-stability is ensured. Clearly, in this case we do not need additional left and right methods because k1 ¼ k2 ¼ 1. Moreover, the 2-step BS2 method possesses the time-reversal symmetry property because the main method is symmetric. One can also verify that the order of the method is p ¼ k ¼ 2. Finally, the stability polynomial defined in (10) and associated with the 2-step BS2 method is the following quadratic polynomial,
pðz; g; qÞ ¼ ð6 þ 3g qÞz2 4ð3 þ qÞz þ ð6 3g qÞ =6; whose discriminant is D ¼ q2 =3 þ 4q þ g2 > 0. So, its roots z1 and z2 are real and they tend to coalesce into 1 when h ! 0. Since pð1; g; qÞ ¼ q < 0 and pð1; g; qÞ ¼ 4 þ q=3 > 0, one root of this polynomial has modulus less than 1. Actually, by considering that it is quadratic, this also implies that the other root has modulus greater than 1. Hence, the 2-step BS2 method is absolutely stable, without any restriction on the mesh size h. Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
5
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
4.2. The 4-step BS2 method When k ¼ 4, we consider the quintic cardinal B-spline in (11), and the corresponding coefficients of the main BS2 method are ð2Þ
ð2Þ
ð2Þ
ð2Þ
ð2Þ
ð1Þ
ð1Þ
ð1Þ
ð1Þ
ð1Þ
ða0 ; a1 ; a2 ; a3 ; a4 Þ ¼ 16 ð1; 2; 6; 2; 1Þ; 1 ð1; 10; 0; 10; 1Þ; ða0 ; a1 ; a2 ; a3 ; a4 Þ ¼ 24 1 ¼ 120 ð1; 26; 66; 26; 1Þ:
ðb0 ; b1 ; b2 ; b3 ; b4 Þ
Generalized zero-stability is ensured in this case, because q2 ðzÞ ¼ ðz 1Þ2 v ðzÞ with v ðzÞ ¼ 16 ðz2 þ 4z þ 1Þ, which has two real negative roots: one with modulus less and the other with modulus greater than 1. The two additional methods are selected by taking the Generalized Backward Differentiation Formulas with 4 steps [6], centered respectively at t 1 and t N1 . The right additional method is the reversed form of the left one, and so the time-reversal symð2;L;1Þ
metry property of the 4-step BS2 method is ensured. For brevity, we report only the coefficients aj
ð1;L;1Þ
; aj
ðL;1Þ
; bj
; j ¼ 0; . . . ; k,
of the left one: ð2;L;1Þ
ða0
ð2;L;1Þ
; a1
ð2;L;1Þ
; a2
ð2;L;1Þ
; a3
ð1;L;1Þ ð1;L;1Þ ð1;L;1Þ ð1;L;1Þ ; 1 ; 2 ; 3 ; ð 0 L;1 L;1 L;1 L;1 L;1 ðb0 ; b1 ; b2 ; b3 ; b4 Þ
a
a
a
a
ð2;L;1Þ
; a4
1 Þ ¼ 12 ð11; 20; 6; 4; 1Þ;
ð1;L;1Þ Þ 4
a
1 ¼ 12 ð3; 10; 18; 6; 1Þ;
¼ ð0; 1; 0; 0; 0Þ:
One can verify that in this case the order is p ¼ k ¼ 4 for the main method, but only p ¼ 3 for the additional methods. This reduced order at the boundary does not lead to a reduction of the order of the full scheme, as mentioned in [6, Remark 4.4.1]. We also refer to the numerical orders reported in Section 6 for the considered examples. The stability polynomial defined in (10) and associated with this method is now the following quartic polynomial,
pðz; g; qÞ ¼ ð20 þ 5g qÞz4 þ ð40 þ 50g 26qÞz3 ð120 þ 66qÞz2 þ ð40 50g 26qÞz þ ð20 5g qÞ =120: We first show (by contradiction) that this polynomial does not have any root on the unit circle of the complex plane. Assuming there exists such a root, it must satisfy pðexpðihÞ; g; qÞ ¼ 0, and we get
q ¼ qðhÞ ¼ q2 ðhÞ þ gq1 ðhÞ; with q2 ðhÞ :¼ q2 ðexpðihÞÞ=rðexpðihÞÞ and q1 ðhÞ :¼ q1 ðexpðihÞÞ=rðexpðihÞÞ, see (10). For fixed g, one can verify that qðhÞ describes in the complex plane a closed curve passing through the origin (because q2 ð1Þ ¼ q1 ð1Þ ¼ 0 and rð1Þ ¼ 1). Moreover, since
1 3 1 q1 ðexpðihÞÞ ¼ i expði2hÞ½ sinð2hÞ 10 sinðhÞ; 12 1 rðexpðihÞÞ ¼ expði2hÞ½cosð2hÞ þ 26 cosðhÞ þ 33; 60
q2 ðexpðihÞÞ ¼ expði2hÞ½cosð2hÞ þ 2 cosðhÞ 3;
the curve qðhÞ belongs to the half-plane ReðqÞ 6 0. However, from its definition it follows that q is a real positive value (see (10)), so the stability polynomial cannot have a root on the unit circle of the complex plane. Furthermore, when considering the roots of the stability polynomial as continuous functions of h, they tend to the roots of q2 for h ! 0. We know that q2 has two real negative roots, one strictly inside and the other strictly outside the unit circle. Hence, we can conclude that for any h at least one root of p is strictly inside the unit circle and another is strictly outside. On the other hand, because pð1; g; qÞ ¼ 2ð10 þ qÞ=15 < 0 and pð1; g; qÞ ¼ q < 0, we can say that p must have an even number of real roots in ð1; 1Þ, counted with their multiplicity. This means that the stability polynomial must have exactly two complex conjugate roots or two real roots inside the unit circle. In any case, we can conclude that the other two roots are outside the unit circle. This implies that our 4-step BS2 method is absolutely stable. 5. The continuous spline extension Like for BS methods, it is possible to extend the discrete solution generated by the k-step BS2 method with a spline of degree k þ 1 and knots at the mesh points. This spline collocates the differential equation at the mesh points. However, in the BS2 case, in order to guarantee a unique continuous spline extension, it will be shown that special combinations of h and l have to be avoided. Under this assumption, the coefficients of such a spline in the B-spline form can be computed in a similar way as for the BS methods [10]. For brevity reasons we only summarize the main idea, and we refer to [10] for more details. We look for a spline sh of degree k þ 1 and knots at the mesh points such that
sh ðti Þ ¼ yi ;
Lðsh Þðt i Þ ¼ f i ;
i ¼ 0; . . . ; N;
ð15Þ
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
6
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
where we can represent sh in the B-spline form introduced in (13). Then, we have to determine the coefficient vector c :¼ ðc1k ; . . . ; cN1 ÞT such that 2
T
2
Sc ¼ ðy0 ; . . . ; yN ; h f 0 ; . . . ; h f N Þ ;
S :¼
SI SL
;
ð16Þ
where SI and SL are banded Toeplitz matrices of dimension ðN þ 1Þ ðN þ k þ 1Þ with bandwidth k, and where ðSI Þ1;j :¼ bj1 ð2Þ
ð1Þ
and ðSL Þ1;j :¼ aj1 þ 2hlaj1 ; j ¼ 1; . . . ; k þ 1. The collocation system in (16) is overdetermined. The following proposition ensures when it has a unique solution. Proposition 1. The rectangular linear system in (16) has a unique solution if rankðSÞ ¼ N þ k þ 1 and if the entries of 2
T
2
ð2Þ
ð1Þ
ðy0 ; . . . ; yN ; h f 0 ; . . . ; h f N Þ satisfy (5) where the coefficients ai ; ai ; bi ; i ¼ 0; . . . ; k, are given in (11). Proof. Since the proof is similar to the one in [10, Theorem 6], we only give a brief outline of it. The vector 2
T
2
ðy0 ; . . . ; yN ; h f 0 ; . . . ; h f N Þ satisfies the N k þ 1 linear conditions given in (5). On the other hand, from (12)–(14) we know that for the l-th column of S; l ¼ 1; . . . ; N þ k þ 1, we have k2 X
aiþk1 ðSI Þjþi;l ¼
i¼k1
k2 X
biþk1 ðSL Þjþi;l ;
j ¼ k1 ; . . . ; N k2 :
i¼k1
It follows that there are N k þ 1 redundancies in the rectangular system (16), and this implies that there exists a unique solution since S is a full rank matrix of dimension ð2N þ 2Þ ðN þ k þ 1Þ. h 2
2
T
If the vector ðy0 ; . . . ; yN ; h f 0 ; . . . ; h f N Þ is generated by the k-step BS2 method, then (5) is satisfied. Therefore, Proposition 1 implies that the collocation system in (16) has a unique solution whenever the matrix S has full rank. In the next proposition we analyze the rank of S. Proposition 2. Let Dk be the 2k 2k matrix, where ðDk Þi;j ¼ ðSI Þi;j and ðDk Þkþi;j ¼ ðSL Þi;j for i ¼ 1; . . . ; k and j ¼ 1; . . . ; 2k. The rectangular matrix S of the linear system in (16) has full rank (i.e. N þ k þ 1) if and only if Dk is non-singular. Proof. We will show that the N þ k þ 1 columns of the matrix S are linearly independent. To this end, it suffices to prove that Sd ¼ 0 if and only if the vector d is the zero-vector. Let us start with considering the first 2k columns of S. From the sparsity structure of S and the definition of Dk , we can conclude that these columns are linearly independent if and only if Dk is nonsingular. Moreover, because all the entries ðSI Þi;2kþ1 ¼ . . . ¼ ðSI Þi;Nþkþ1 ¼ 0 and ðSL Þi;2kþ1 ¼ . . . ¼ ðSL Þi;Nþkþ1 ¼ 0; i ¼ 1; . . . ; k, it follows that the first 2k entries of d must be zero to obtain Sd ¼ 0. In a similar way, when considering any set of 2k successive columns of S, we obtain that the corresponding 2k entries of d must be zero, which concludes the proof. h ^2 and h R fh ^ ; h ^ þ g respectively, where In the case of k ¼ 2 and k ¼ 4, the condition in Proposition 2 is satisfied when h – h 4 4
^2 ¼ h
pffiffiffi 3 ; jlj
^ ¼ h 4
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi 345 29 105 ; 8jlj
l – 0;
ð17Þ
and there are no restrictions on h if l ¼ 0. This can be verified by a direct computation. More precisely, recalling g :¼ 2hl, the submatrix Dk described in the proposition is given by
2
1=6 2=3 60 1=6 6 D2 ¼ 6 4 1 g=2 2
2
1=5
60 6 6 60 6 1 6 60 D4 ¼ 6 24 6 4 g 6 60 6 6 40 0
26=5
0
2=3 1=6 1 þ g=2 0
1 g=2 2
0 and
1=6
66=5
3
7 7 7; 5
1 þ g=2 26=5
1=5
0
0
0
3
7 7 7 7 0 1=5 26=5 66=5 26=5 1=5 0 7 0 0 1=5 26=5 66=5 26=5 1=5 7 7 7; 7 8 10g 24 8 þ 10g 4 þ g 0 0 0 7 7 4g 8 10g 24 8 þ 10g 4 þ g 0 0 7 7 5 0 4g 8 10g 24 8 þ 10g 4 þ g 0 8 10g 24 8 þ 10g 4 þ g 0 0 4g 1=5
26=5
66=5
26=5
1=5
0
0
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
7
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
with
detðD2 Þ ¼
g2 12
þ 1;
detðD4 Þ ¼
4 1 g 23g2 þ1 : 32400 120 64
Note that this restriction on the allowed mesh sizes is required only to compute the continuous extension and does not affect the stability properties of the BS2 methods. Moreover, when l 0, usually the problem is singularly perturbed. In such a case, the use of a numerical scheme based on a uniform mesh size is not very efficient and a mesh size variation strategy should be considered. To compute the coefficient vector c of the unique spline solution sh of (15), we could solve (16) in the least-square sense, or use a quasi-interpolation approach similar to the one for BS methods described in [9]. 6. Numerical results In this section we present some numerical results for two linear problems with f ðt; yÞ ¼ cy þ gðtÞ, both depending on a parameter > 0 which determines their stiffness level and both taken from the test set for BVP solvers at the site https:// www.dm.uniba.it/bvpsolvers/testsetbvpsolvers. Furthermore, we consider a standard harmonic oscillator problem – already investigated in [10] – to confirm the computational gain which can be obtained for second order problems by using the BS2 methods instead of the BS ones. 6.1. Two linear problems with boundary layers We consider two linear problems depending on a parameter > 0 for which we have selected two different values, say 1 and 2 , in order to deal once with a non-stiff problem and once with a quite stiff one. The BS2 methods with k ¼ 2 and k ¼ 4
Table 1 Problem 1 with ¼ 1 ¼ 4 102 (top) and ¼ 2 ¼ 4 104 (bottom). The maximal values of the errors at the mesh points (err m ), the related numerical approximation orders, and the maximal spline approximation errors (err s ) are given. N
1
2
k¼2
k¼4
err m
Order
err s
err m
Order
err s
10
1:26 102
/
1:31 102
4:58 103
/
8:33 103
20
3:06 103
2:04
3:08 103
6:43 104
2:83
9:55 104
40
4
7:65 10
2:00
4
7:66 10
5
3:94 10
4:03
5:29 105
80
1:91 104
2:00
1:91 104
1:70 106
4:53
2:18 106
50
6:35 102
/
1:01 101
2:79 102
/
7:95 102
100
1:69 102
1:90
1:74 102
7:14 103
1:97
1:31 102
200
3
3:92 10
2:11
3
3:94 10
4
7:17 10
3:32
1:07 103
400
9:63 104
2:03
9:63 104
4:08 105
4:14
5:05 105
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
−0.2
−0.2
−0.4
−0.4
−0.6
−0.6
−0.8
−0.8
−1 −1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
−1 −1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
4
Fig. 1. Problem 1 with ¼ 2 ¼ 4 10 and N ¼ 50. The exact solution (blue) is compared with the spline extension (red) of the numerical solution produced by the BS2 methods with k ¼ 2 (left) and k ¼ 4 (right). The numerical solution is depicted with the symbol. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
8
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
Table 2 Problem 2 with ¼ 1 ¼ 2 (top) and ¼ 2 ¼ 0:02 (bottom). The maximal values of the errors at the mesh points (errm ), the related numerical approximation orders, and the maximal spline approximation errors (errs ) are given. N
k¼2
k¼4 Order
err m
1
2
err s
err m
Order
err s
5
7:54 10
/
7:78 10
1:10 10
/
1:85 103
10
1:85 103
2:03
1:87 103
8:14 105
3:76
1:14 104
20
4:64 104
2:00
4:64 104
3:82 106
4:41
4:99 106
40
4
2:00
4
1:46 10
7
4:71
1:84 107
1
3
3
1:16 10
3
1:16 10
3
50
1:37 10
/
2:52 10
6:65 10
/
6:23 102
100
3:53 102
1:94
3:60 102
4:18 103
0:66
1:23 102
200
8:04 103
2:13
8:07 103
5:71 104
2:87
9:74 104
400
3
2:03
3
5
3:90
4:92 105
1:97 10
1
3:82 10
1:97 10
have been tested by successively halving the mesh size h, starting from a suitably selected mesh size h0 . Note that for all the performed tests we have verified that the infinity norm of both the first and last columns of the matrix M introduced in (9) is equal to 1. The first problem (BVPT14) is given by
y00 ¼ y ð1 þ p2 Þ cosðptÞ;
t 2 ½1; 1;
yð1Þ ¼ yð1Þ ¼ expð2=
pffiffiffi Þ;
pffiffiffi pffiffiffi pffiffiffi whose exact solution is yðtÞ ¼ cosðptÞ þ expððt 1Þ= Þ þ expððt þ 1Þ= Þ. It has a boundary layer of width Oð Þ near t ¼ 1. Note that in this case we have l ¼ 0, and so the spline extension of the solution produced by the BS2 methods pffiffiffi pffiffiffi can be defined for any mesh size. For the related experiments we have set h0 ¼ when ¼ 1 ¼ 4 102 and h0 ¼ 2 4 when ¼ 2 ¼ 4 10 . Table 1 reports the maximal errors at the mesh points, the corresponding numerical approximation orders, and the maximal spline approximation errors (sampled at 2000 uniformly spaced points) obtained with the two studied BS2 methods. The results confirm that these methods have order 2 and 4 respectively. Fig. 1 illustrates the case ¼ 2 with N ¼ 50, and it can be seen that only near the boundary layers the exact solution slightly differs from the spline extension sh of the numerical solution produced by the BS2 methods. The second problem (BVPT4) is
y00 þ y0 ¼ ð1 þ Þy;
t 2 ½1; 1;
yð1Þ ¼ 1 þ e2 ;
yð1Þ ¼ 1 þ e2ð1þÞ= ;
whose exact solution is yðtÞ ¼ et1 þ eðtþ1Þð1þÞ= . It has a boundary layer of width OðÞ near t ¼ 1. In view of our notation, l is now equal to 1=ð2Þ, and the experiments have been done for ¼ i ; i ¼ 1; 2 with 1 ¼ 2 and 2 ¼ 2 102 . The maximal mesh size is h0 ¼ 0:4 (N ¼ 5) in the first case (the non-stiff case) and h0 ¼ 2 (N ¼ 50) in the second case. We note that the critical step sizes given in (17) are avoided, and so the spline extensions can be computed. The numerical results are collected in Table 2, which is structured in the same way as Table 1. Fig. 2 shows that the methods are capable of producing a nice continuous approximation of the solution, even when very few points are used for non-stiff problems.
1.3
1.3
1.2
1.2
1.1
1.1
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5 −1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
0.5 −1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
Fig. 2. Problem 2 with ¼ 1 ¼ 2 and N ¼ 5. The exact solution (blue) is compared with the spline extension (red) of the numerical solution produced by the BS2 methods with k ¼ 2 (left) and k ¼ 4 (right). The numerical solution is depicted with the symbol. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
9
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
Table 3 Comparison between the method implemented in the solver colnew and BS2 methods for Problem 1 (top) and 2 (bottom). The maximal values of the errors at the collocation points (err m ) and the maximal spline approximation errors (errs ) are given for the two methods using the same number of collocation points (N c ). Nc
Problem 1, 8
BS2, k ¼ 4
colnew err m
err s
err m
err s
1 ¼ 0:04 6:99 103
1:41 102
4:74 103
9:96 103
16
4:39 104
2:32 103
1:69 103
2:69 103
32
2:72 105
2:22 104
1:16 104
1:61 104
64
1:69 106
1:72 105
5:15 106
6:18 106
Problem 2, 8
1 ¼ 2 4:98 105
5:60 104
3:39 104
5:10 104
16
3:04 106
4:29 105
1:41 106
1:88 105
32
1:93 107
2:97 106
4:93 107
6:28 107
64
1:20 108
1:95 107
1:62 108
2:02 108
5
5
4
4
3
3
2
2
1
1
0
0
−1
−1
−2
−2
−3
−3
−4
−4
−5 −5
−4
−3
−2
−1
0
1
2
3
4
5
−5 −5
5
5
4
4
3
3
2
2
1
1
0
0
−1
−1
−2
−2
−3
−3
−4
−4
−5 −5
−4
−3
−2
−1
0
1
2
3
4
5
−5 −5
−4
−3
−2
−1
0
1
2
3
4
5
−4
−3
−2
−1
0
1
2
3
4
5
Fig. 3. The oscillator test problem given in (18). The exact solution (green) is compared with the spline extension (red) of the BS2 solution in the phase plane for different values of N : N ¼ 10 (top-left), N ¼ 15 (top-right), N ¼ 20 (bottom-left), N ¼ 25 (bottom-right). The numerical solution is produced by the BS2 method with k ¼ 4. The symbol identifies the mesh points. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
We now illustrate the difference between our method and standard collocation applied to Problems 1 and 2 with ¼ 1 . To this aim, we compare the 4-step BS2 method with the method implemented in the solver colnew [4,5] collocating at two Gaussian points in each subinterval.1 They are both fourth order schemes and do not convert the considered differential problem into a system of first order differential equations. More specifically, for the tests we have employed the function bvpcol
1 The solver colnew generates adaptively the mesh, and so the value ¼ 2 could not be used because we want to deal only with uniform meshes for a fair comparison.
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046
10
C. Manni et al. / Applied Mathematics and Computation xxx (2014) xxx–xxx
implemented in the R package bvpSolve [8]. The results are reported in Table 3, where the two methods are chosen such that they use the same number of collocation points. The table confirms that the behavior of the two methods is similar. The absolute error at the collocation points is mostly smaller for collocation at Gaussian points, whereas the opposite is true for the error in the associated spline. Note that the BS2 method with k ¼ 4 gives a quintic C 4 spline while colnew provides just a cubic C 1 spline. 6.2. A harmonic oscillator problem In this section we show that the BS2 methods preserve the high quality spline approximations of the solutions of second order BVPs produced by the BS methods, but with a reduced computational cost. We consider the harmonic oscillator test problem,
y00 þ x2 y ¼ 0;
t 2 ½0; p=2;
yð0Þ ¼ 3;
yðp=2Þ ¼ 3;
ð18Þ
with x ¼ 5, whose exact solution is yðtÞ ¼ 3ðcos xt sin xtÞ. Let us denote by sh the quintic spline approximation of y produced by the BS2 method with k ¼ 4. Its qualitative behavior can be analyzed by looking at Fig. 3: for different h steps, the 0 0 curve pffiffiffi Xh ðtÞ :¼ ðsh ðtÞ; sh ðtÞ=xÞ is compared with the curve XðtÞ :¼ ðyðtÞ; y ðtÞ=xÞ which is a circumference with radius equal to 3 2. Note that yðtÞ is also the solution of the oscillator test problem considered in Problem 2 in [10], where the BS methods were used to solve numerically the above second order equation reformulated as a bidimensional first order problem on the integration domain ½0; 2p and with boundary conditions fixing y at t ¼ 2p and its derivative at t ¼ 0. Fig. 3 may be compared with Fig. 6.1 in [10] because in both cases the results are obtained using a numerical scheme of fourth order. The comparison shows that the BS2 method attains the same approximation quality with almost half the computational cost of the BS method. Indeed, for example with N ¼ 20, here we are using 21 scalar degrees of freedom; on the other hand, the same accuracy was obtained in [10] with N ¼ 80 on a four times larger integration interval, that is with 162 scalar degrees of freedom. 7. Conclusions We have proposed a class of k-step LMMs based on B-splines (the so-called BS2 methods) to solve numerically the second order BVP given in (2) without doubling its dimension. These methods are a specific formulation of the BS methods introduced in [10]. When k is even, they are symmetric and of order k. Moreover, we have shown that the generalized zero-stability and absolute stability requirements are satisfied for k ¼ 2 and k ¼ 4. Like for BS methods, it is possible to extend the discrete solution generated by the k-step BS2 method with a spline of degree k þ 1 and knots at the mesh points. This spline collocates the differential problem at the mesh points and it can be always defined except for few particular values of the mesh size. Finally, we want to point out that for stiff problems, in order to avoid very fine uniform meshes, it could be better to consider a non-uniform extension of these methods. This topic is currently under investigation. References [1] P. Amodio, L. Brugnano, The conditioning of Toeplitz band matrices, Math. Comput. Model. 23 (1996) 29–42. [2] P. Amodio, I. Sgura, High-order finite difference schemes for the solution of second-order BVPs, J. Comput. Appl. Math. 176 (2005) 59–76. [3] U.M. Ascher, R.M.M. Mattheij, R.D. Russel, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, Prentice Hall, Englewood Cliffs, New York, 1988. [4] U.M. Ascher, J. Christiansen, R.D. Russell, Collocation software for boundary-value ODEs, ACM Trans. Math. Softw. 7 (1981) 209–222. [5] G. Bader, U.M. Ascher, A new basis implementation for a mixed order boundary value ODE solver, SIAM J. Sci. Stat. Comput. 8 (1987) 483–500. [6] L. Brugnano, D. Trigiante, Solving Differential Problems by Multistep Initial and Boundary Value Methods, Gordon and Breach Science Publishers, Amsterdam, 1998. [7] C. de Boor, A Practical Guide to Splines, Revised Edition., Springer-Verlag, New York, 2001. [8] F. Mazzia, J.R. Cash, K. Soetaert, Solving boundary value problems in the open source software R: package bvpSolve, Opuscula Math. 34 (2014) 387–403. [9] F. Mazzia, A. Sestini, The BS class of Hermite spline quasi-interpolants on nonuniform knot distributions, BIT Numer. Math. 49 (2009) 611–628. [10] F. Mazzia, A. Sestini, D. Trigiante, B-spline multistep methods and their continuous extensions, SIAM J. Numer. Anal. 44 (2006) 1954–1973. [11] F. Mazzia, A. Sestini, D. Trigiante, BS linear multistep methods on non-uniform meshes, J. Numer. Anal. Ind. Appl. Math. 1 (2006) 131–144. [12] F. Mazzia, A. Sestini, D. Trigiante, The continuous extension of the B-spline linear multistep methods for BVPs on non-uniform meshes, Appl. Numer. Math. 59 (2009) 723–738.
Please cite this article in press as: C. Manni et al., BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. (2014), http://dx.doi.org/10.1016/j.amc.2014.08.046