- Email: [email protected]
Approximation by integro cubic splines
Applied Mathematics and Computation 175 (2006) 8–15 www.elsevier.com/locate/amc
Approximation by integro cubic splines Hossein Behforooz Department of Mathematics, Utica College, Utica, NY 13502, USA
Abstract Years ago, after I presented a paper on spline functions at Oxford University, Professor Powell criticized us for using, most of the time, the function values rather than the integral values on constructing of the spline functions. His comments and his request became the main motivation for this work. In this paper, we assume that, on each subinterval of the spline interval [a, b], the integral value of the function y = y(x) is known. By using these values, rather than the function values at the knots, we introduce a class of new types of interpolatory cubic splines to approximate the function y = y(x). The selection of the end conditions for our integro cubic splines will be discussed. The numerical examples and computational results, illustrate and guarantee a higher accuracy for this approximation. Ó 2005 Elsevier Inc. All rights reserved. Keywords: Cubic spline; End conditions; Interpolation; Knot; mc-Spline
1. Introduction In practice and application, sometimes we deal with phenomena which involve the function y = y(x) and also its integral. For example, in mechanics, the velocity v(t) and the displacement s(t), or the acceleration a(t) and the E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.07.066
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
9
velocity v(t); in statistics, the probability density function and the cumulative distribution function; and, in electricity, the current function I(t) and charge function q(t) are some real examples about our consideration. In this paper, we construct a piecewise cubic spline S(x) to approximate the function y = y(x) when the integrals of S(x) and y(x) in each of the subinterval match each other. Then numerical examples will be presented in the last section to illustrate the accuracy of the results of this paper. Recently, Behforooz [1,2] has obtained pairs of end conditions by using the integral values at two end subintervals, but not in a global fashion over the whole interval [a, b].
2. Integro cubic splines with first derivative representation Suppose that the interval [a, b] is partitioned by the following k + 1 equally spaced points: a ¼ x0 < x1 < x2 < < xk1 < xk ¼ b;
ð1Þ
such that xi = a + ih, for i = 0, 1, 2, . . . , k, with h ¼ ba . k In traditional interpolation problems by cubic spline functions, we assume that the function values at knots are given and y(xi); i = 0(1)k, are known and we use them to construct the cubic spline S(x) such that S(xi) = y(xi); i = 0(1)k. In this paper, these function values are not given and we assume that the values of all k integrals of y = y(x) on k intervals [xi1, xi] are known and they are equal to Z xi yðxÞ dx ¼ I i ; i ¼ 1ð1Þk. ð2Þ xi1
Now, we are going to use these numbers Ii and construct a class of integro cubic spline functions S(x). These cubic splines S(x) are piecewise cubic polynomials such that S, S 0 and S00 are continuous on [a, b], i.e., they are pieced together at the interior knots in such a way that their left side and right side values coincide at each knot xi; i = 1(1)k 1 and in each subinterval [xi1, xi] we have Z xi Z xi yðxÞ dx ¼ SðxÞ dx ¼ I i ; i ¼ 1ð1Þk. ð3Þ xi1
xi1
For simplicity, we will use the following notations: yi = y(xi), Si = S(xi), mi = S 0 (xi), and Mi = S00 (xi). By using the cubic Hermite interpolation polynomial formula, S(x) can be written as SðxÞ ¼ /1 ðxÞS i1 þ /2 ðxÞS i þ /3 ðxÞmi1 þ /4 mi ;
x 2 ½xi1 ; xi ;
ð4Þ
10
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
where, 8 > > > /1 > > /3 > > > > :/ 4
2
3
¼ h13 f3hðx xi Þ þ 2ðx xi Þ g; ¼ h13 f3hðx xi1 Þ2 2ðx xi1 Þ3 g; 2
ð5Þ
3
¼ h12 fhðx xi Þ þ ðx xi Þ g; 2
3
¼ h12 fhðx xi1 Þ þ ðx xi1 Þ g.
By considering two pieces of S(x) on [xi1, xi] and [xi, xi+1] we can easily verify that S and S 0 are continuous at the interior knots xi; i = 1(1)k 1. The continuity of the second derivatives S00 (x) at xi, i = 1(1)k 1, gives us the following consistency relations: 3 mi1 þ 4mi þ miþ1 ¼ ðS iþ1 S i1 Þ; i ¼ 1ð1Þk 1. ð6Þ h To construct a unique integro cubic spline S(x) we need to have the values of 2k + 2 unknowns m0, m1, . . . , mk and S0, S1, . . . , Sk. But, in the traditional construction, we need only k + 1 unknowns m0, m1, . . . , mk to construct S(x) to interpolate the function at k + 1 knots and satisfy y(xi) = S(xi); i = 0(1)k. The integration of S(x) from (4), together with (5) and (3), gives the following relations between the values of mÕs, SÕs and IÕs: h2 ðmi1 mi Þ þ 6hðS i1 þ S i Þ ¼ 12I i ; 2
h ðmi miþ1 Þ þ 6hðS i þ S iþ1 Þ ¼ 12I iþ1 ;
i ¼ 1ð1Þk; i ¼ 0ð1Þk 1.
ð7Þ ð8Þ
By adding and subtracting (7) and (8) we get 6 12 ð9Þ mi1 miþ1 ¼ ðS i1 þ 2S i þ S iþ1 Þ þ 2 ðI i þ I iþ1 Þ; h h 6 12 mi1 2mi þ miþ1 ¼ ðS i1 S iþ1 Þ þ 2 ðI i I iþ1 Þ. ð10Þ h h So, to construct S(x), we need 2k + 2 unknowns and we have enough linear equations to solve them simultaneously and obtain those unknowns. But fortunately, the above relations are such that we can easily eliminate the parameters S between Eqs. (6) and (10) and obtain a new consistency relations only between unknown numbers mi and known values Ii, 12 mi1 þ 10mi þ miþ1 ¼ 2 ðI iþ1 I i Þ; i ¼ 1ð1Þk 1. ð11Þ h These equations reduce the system of linear equations from (2k + 2) by (2k + 2) to (k + 1) by (k + 1). Obviously, in order to use the Eq. (11) and solve them for k + 1 unknowns m0, m1, . . . , mk, we need (as usual) two additional linear equations. In the literature, these equations are called the end conditions. The best set of end conditions for our purpose are the first derivative end conditions at two end points a and b. Suppose that y 0 (a) = a and y 0 (b) = b are given. Then by setting m0 = a and mk = b, we can easily solve the following
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
11
(k 1) by (k 1) linear tridiagonal equations to obtain a unique set of solutions for m1, m2, . . . , mk1: 8 < 10m1 þ m2 ¼ b1 a; ð12Þ mi1 þ 10mi þ miþ1 ¼ bi ; i ¼ 2ð1Þk 2; : mk1 þ 10mk ¼ bk1 b; ðI iþ1 I i Þ. where bi ¼ 12 h2 When m0, m1, . . . , mk are known, then we can use (7) or (8) to compute S0, S1, . . . , Sk. But we need another additional given value for y(a) or y(b). If y0 = y(a) is given, then by considering S0 = y0 = y(a) and using the Eq. (7), we can obtain the values of S1, S2, . . . , Sk. With these two sets of values of mÕs and SÕs, and using S(x) from (4) with (5), we can approximate y(x) at any x 2 [a, b]. If the third end condition y0 = y(a) is not given and y(a) is not available, then we can use the simple difference forward or backward formula for the first derivative approximation and use one of the following relations: S 1 S 0 ¼ hm0
or
S 1 S 0 ¼ hm1
ð13Þ
as an additional equation to mix with the Eq. (7) to find all k + 1 unknowns S0, S1, . . . , Sk. 3. Integro cubic splines with second derivative representation Instead of the cubic Hermite polynomial (4), we can construct a cubic spline S(x) in terms of the second derivatives. The following is another representation of S(x), for x 2 [xi1, xi]; i = 1(1)k, in terms of Mi1, Mi, Si1 and Si: 1 SðxÞ ¼ fðx xi1 Þ3 M i ðx xi Þ3 M i1 g 6h 1 h2 h2 ðM i1 M i Þ þ ðS i S i1 Þ ðx xi Þ þ S i M i . ð14Þ þ h 6 6 It is easy to see that S(x) and S00 (x) are continuous at the interior points xi; i = 1(1)k 1. In order to have a continuous S 0 (x) at xi, we must have the next consistency relations which is in terms of Mi1, Mi, Mi+1, and Si1, Si and Si+1. 6 ð15Þ M i1 þ 4M i þ M iþ1 ¼ 2 ðS i1 2S i þ S iþ1 Þ; i ¼ 1ð1Þk 1. h By using the integrals of (14) in (3), we obtain the following relations between MÕs, SÕs and IÕs: h3 h ðM i1 þ M i Þ þ ðS i1 þ S i Þ ¼ I i ; i ¼ 1ð1Þk; 2 24 h3 h ðM i þ M iþ1 Þ þ ðS i þ S iþ1 Þ ¼ I iþ1 ; i ¼ 0ð1Þk 1. 2 24
ð16Þ ð17Þ
12
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
By adding and subtracting (16) and (17) we obtain h3 h ðM i1 þ 2M i þ M iþ1 Þ þ ðS i1 þ 2S i þ S iþ1 Þ ¼ I iþ1 þ I i ; ð18Þ 2 24 h3 h ð19Þ ðM iþ1 M i1 Þ þ ðS iþ1 S i1 Þ ¼ I iþ1 I i . 2 24 Unlike the first derivative representation, here, we cannot eliminate SÕs between (15)–(19) to obtain a relation similar to (11) without SÕs. So, the advantage of the Hermite representation or using the relations (11) is that, to find the parameters, we need only to solve a simple (k 1) by (k 1) tridiagonal equations, not a system with a full matrix of order (2k + 2) by (2k + 2). Even using the natural cubic spline end conditions M0 = Mk = 0, or knowing the second derivative values at x = a and x = b and using them as two additional end conditions M0 = y00 (a) and Mk = y00 (b), do not lead us to solving a system of (k 1) by (k 1) tridiagonal equations. So, to construct S(x) by (14) with any end conditions, we have to solve a system of linear equations with a full matrix of order (2k + 2) by (2k + 2).
4. Other end conditions In practice we do not expect to know the values of the first derivative of the function at two end points. But these two end conditions m0 = a and mk = b reduce the (2k + 2) by (2k + 2) system to an easy system of (k 1) by (k 1) tridiagonal equations. Furthermore, if y 2 C4[a, b], then the cubic spline with these two end conditions has the highest order of convergence O(h4). In the literature, we have so many other end conditions for this purpose. Even with the same accuracy and with the same order of convergence O(h4). See for example, the end conditions that are mentioned in Behforooz [1] and the references therein. Most of those end conditions are in terms of the derivatives of the spline functions and the function values. But here, the problem is that we do not have the function values at the knots and when we use those types of end conditions, we have to solve a system of equations with full matrix of order (2k + 2) by (2k + 2). 5. Minimum property The traditional cubic spline S(x) with interpolation conditions S(xi) = y(xi); i = 0(1)k with either end conditions M0 = Mk = 0 (for natural cubic spline) or the first derivative end conditions m0 = y 0 (a), mk = y 0 (b), minimizes the integral Z b ½f 00 ðxÞ2 dx ð20Þ a
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
13
over all the functions y = f(x) 2 C2[a, b]. To investigate this, we start as follows: Z b 2 ½f 00 ðxÞ S 00 ðxÞ dx a Z b Z b Z b 2 2 ½f 00 ðxÞ dx þ ½S 00 ðxÞ dx 2 f 00 ðxÞS 00 ðxÞ dx ¼ a a a Z b Z b Z b ½f 00 ðxÞ2 dx ½S 00 ðxÞ2 dx 2 ½f 00 ðxÞ S 00 ðxÞS 00 ðxÞ dx. ¼ a
a
a
Applying the integration by parts on the last integral we obtain Z b ½f 00 ðxÞ S 00 ðxÞS 00 ðxÞ dx a Z b b 0 00 0 ½f 0 ðxÞ S 0 ðxÞS 000 ðxÞ dx ¼ ½f ðxÞ S ðxÞS ðxÞja a b
¼ ½f 0 ðxÞ S 0 ðxÞS 00 ðxÞja
k X i¼1
0
0
¼ ½f ðxÞ S ðxÞS
00
b ðxÞja
k X i¼1
Z
xi
Ci
½f 0 ðxÞ S 0 ðxÞ dx
xi1
xi C i ½f ðxÞ SðxÞ ;
ð21Þ
xi1
where, Ci = S000 (x) is a constant in each interval [xi1, xi]. So, for the traditional cubic spline S(x) with interpolation conditions S(xi) = y(xi); i = 0(1)k with either end conditions {M0 = Mk = 0} or {m0 = y 0 (a), mk = y 0 (b)} we have Z b ½f 00 ðxÞ S 00 ðxÞS 00 ðxÞ dx ¼ 0 a
and Z
b 2
½S 00 ðxÞ dx ¼ a
Z
b 2
½f 00 ðxÞ dx
a
Therefore for all f 2 C2[a, b] Z b Z b ½S 00 ðxÞ2 dx 6 ½f 00 ðxÞ2 dx; b
Z
b 2
½f 00 ðxÞ S 00 ðxÞ dx 6 a
Z
b 2
½f 00 ðxÞ dx. a
ð22Þ
a
which is called the minimum property of S(x) with either one of those two sets of end conditions. In other words, S(x) is the smoothest function which has the minimum total curvature. In our integro cubic spline with end conditions m0 = y 0 (a) and mk = y 0 (b), the first bracket in (21) is equal to zero. But since our interpolating and fitting conditions are Z xi Z xi yðxÞ dx ¼ SðxÞ dx ¼ I i ; i ¼ 1ð1Þk xi1
xi1
14
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
not S(xi) = y(xi), then we do not make the second bracket equal to zero. But the cubic spline with end conditions m0 = y 0 (a) and mk = y 0 (b) has the highest order of convergence and for f 2 C4[a, b], jf ðxÞ SðxÞj ¼ Oðh4 Þ. The numerical examples below, show that the order of convergence is high enough to see that jf(x) S(x)j is almost zero at xi; i = 0(1)k and this implies that the second bracket in Eq. (21) is almost zero and our integro cubic spline has minimum total curvature.
6. Numerical examples We examined this method on several functions and the accuracy of the approximation is very high. I compared the results of this integro cubic spline with my other ordinary higher degree splines and the results are reasonably comparable, see the numerical results in Behforooz [3,4]. Here we show the results of one of those functions. Let S be the integro cubic spline to interpolate the function y(x) = exp(x) at the knots xi = 0.05i, i = 0(1)20, and satisfies the end conditions m0 = exp(0) = 1 and m20 = exp(1). In Table 1 we have listed the values of the error E(x) = jS(x) exp(x)j computed at some points in the interval [0, 1]. For comparison, we also listed the errors which are computed by using the ordinary natural cubic spline (NCS) and the famous not-a-knot cubic spline. This example shows that our cubic spline S is more accurate than the other two cubic splines (specially at the two end subintervals). It is worth mentioning that, to find the inverse of the coefficient matrix of (11), one may use an algorithm similar to the algorithm of Moon [6], which reduces the number of multiplication from O(n2) to O(n). Also, Delhez in [5], has introduced the mc-spline of degree four, but to construct this spline we have to solve a system Table 1 Comparison of jexp(x) s(x)j for three cubic splines x
N.C.S.
Not-a-knot
E(x)
0.01 0.02 0.07 0.09 0.22 0.36 0.62 0.93 0.96 0.98 0.99
9.95E5 1.23E4 3.36E5 1.34E5 6.17E7 2.04E8 5.43E8 8.92E5 1.37E4 3.32E4 2.74E4
1.54E6 1.84E6 7.02E7 2.82E7 1.83E7 9.95E8 2.84E7 1.74E6 1.76E6 4.54E6 3.87E7
9.00E9 1.62E8 1.21E9 1.04E8 2.65E9 1.07E8 5.61E9 2.30E8 3.61E8 2.35E8 4.32E8
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
15
of (2k + 2) by (2k + 2) blocked matrix linear equations which is more difficult than our (k + 1) by (k + 1) diagonal equations (11).
References [1] H. Behforooz, End conditions for cubic spline interpolation derived from integration, Appl. Math. Comput. 29 (1989) 231–246. [2] H. Behforooz, A new approach to spline functions, Appl. Numer. Math. 13 (1993) 271–276. [3] H. Behforooz, The not-a-knot interpolatory cubic polynomial, Appl. Math. Comput. 52 (1992) 29–35. [4] H. Behforooz, A comparison of the E(3) and not-a-knot cubic spline, Appl. Math. Comput. 72 (1992) 219–223. [5] E.J.M. Delhez, A spline interpolation technique that preserve mass budget, Appl. Math. Lett. 16 (2003) 17–26. [6] B.S. Moon, An explicit solution for the cubic spline interpolation for functions of single variable, J. Appl. Math. Comput. 117 (2001) 251–255.
Approximation by integro cubic splines Hossein Behforooz Department of Mathematics, Utica College, Utica, NY 13502, USA
Abstract Years ago, after I presented a paper on spline functions at Oxford University, Professor Powell criticized us for using, most of the time, the function values rather than the integral values on constructing of the spline functions. His comments and his request became the main motivation for this work. In this paper, we assume that, on each subinterval of the spline interval [a, b], the integral value of the function y = y(x) is known. By using these values, rather than the function values at the knots, we introduce a class of new types of interpolatory cubic splines to approximate the function y = y(x). The selection of the end conditions for our integro cubic splines will be discussed. The numerical examples and computational results, illustrate and guarantee a higher accuracy for this approximation. Ó 2005 Elsevier Inc. All rights reserved. Keywords: Cubic spline; End conditions; Interpolation; Knot; mc-Spline
1. Introduction In practice and application, sometimes we deal with phenomena which involve the function y = y(x) and also its integral. For example, in mechanics, the velocity v(t) and the displacement s(t), or the acceleration a(t) and the E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.07.066
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
9
velocity v(t); in statistics, the probability density function and the cumulative distribution function; and, in electricity, the current function I(t) and charge function q(t) are some real examples about our consideration. In this paper, we construct a piecewise cubic spline S(x) to approximate the function y = y(x) when the integrals of S(x) and y(x) in each of the subinterval match each other. Then numerical examples will be presented in the last section to illustrate the accuracy of the results of this paper. Recently, Behforooz [1,2] has obtained pairs of end conditions by using the integral values at two end subintervals, but not in a global fashion over the whole interval [a, b].
2. Integro cubic splines with first derivative representation Suppose that the interval [a, b] is partitioned by the following k + 1 equally spaced points: a ¼ x0 < x1 < x2 < < xk1 < xk ¼ b;
ð1Þ
such that xi = a + ih, for i = 0, 1, 2, . . . , k, with h ¼ ba . k In traditional interpolation problems by cubic spline functions, we assume that the function values at knots are given and y(xi); i = 0(1)k, are known and we use them to construct the cubic spline S(x) such that S(xi) = y(xi); i = 0(1)k. In this paper, these function values are not given and we assume that the values of all k integrals of y = y(x) on k intervals [xi1, xi] are known and they are equal to Z xi yðxÞ dx ¼ I i ; i ¼ 1ð1Þk. ð2Þ xi1
Now, we are going to use these numbers Ii and construct a class of integro cubic spline functions S(x). These cubic splines S(x) are piecewise cubic polynomials such that S, S 0 and S00 are continuous on [a, b], i.e., they are pieced together at the interior knots in such a way that their left side and right side values coincide at each knot xi; i = 1(1)k 1 and in each subinterval [xi1, xi] we have Z xi Z xi yðxÞ dx ¼ SðxÞ dx ¼ I i ; i ¼ 1ð1Þk. ð3Þ xi1
xi1
For simplicity, we will use the following notations: yi = y(xi), Si = S(xi), mi = S 0 (xi), and Mi = S00 (xi). By using the cubic Hermite interpolation polynomial formula, S(x) can be written as SðxÞ ¼ /1 ðxÞS i1 þ /2 ðxÞS i þ /3 ðxÞmi1 þ /4 mi ;
x 2 ½xi1 ; xi ;
ð4Þ
10
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
where, 8 > > > /1 > > /3 > > > > :/ 4
2
3
¼ h13 f3hðx xi Þ þ 2ðx xi Þ g; ¼ h13 f3hðx xi1 Þ2 2ðx xi1 Þ3 g; 2
ð5Þ
3
¼ h12 fhðx xi Þ þ ðx xi Þ g; 2
3
¼ h12 fhðx xi1 Þ þ ðx xi1 Þ g.
By considering two pieces of S(x) on [xi1, xi] and [xi, xi+1] we can easily verify that S and S 0 are continuous at the interior knots xi; i = 1(1)k 1. The continuity of the second derivatives S00 (x) at xi, i = 1(1)k 1, gives us the following consistency relations: 3 mi1 þ 4mi þ miþ1 ¼ ðS iþ1 S i1 Þ; i ¼ 1ð1Þk 1. ð6Þ h To construct a unique integro cubic spline S(x) we need to have the values of 2k + 2 unknowns m0, m1, . . . , mk and S0, S1, . . . , Sk. But, in the traditional construction, we need only k + 1 unknowns m0, m1, . . . , mk to construct S(x) to interpolate the function at k + 1 knots and satisfy y(xi) = S(xi); i = 0(1)k. The integration of S(x) from (4), together with (5) and (3), gives the following relations between the values of mÕs, SÕs and IÕs: h2 ðmi1 mi Þ þ 6hðS i1 þ S i Þ ¼ 12I i ; 2
h ðmi miþ1 Þ þ 6hðS i þ S iþ1 Þ ¼ 12I iþ1 ;
i ¼ 1ð1Þk; i ¼ 0ð1Þk 1.
ð7Þ ð8Þ
By adding and subtracting (7) and (8) we get 6 12 ð9Þ mi1 miþ1 ¼ ðS i1 þ 2S i þ S iþ1 Þ þ 2 ðI i þ I iþ1 Þ; h h 6 12 mi1 2mi þ miþ1 ¼ ðS i1 S iþ1 Þ þ 2 ðI i I iþ1 Þ. ð10Þ h h So, to construct S(x), we need 2k + 2 unknowns and we have enough linear equations to solve them simultaneously and obtain those unknowns. But fortunately, the above relations are such that we can easily eliminate the parameters S between Eqs. (6) and (10) and obtain a new consistency relations only between unknown numbers mi and known values Ii, 12 mi1 þ 10mi þ miþ1 ¼ 2 ðI iþ1 I i Þ; i ¼ 1ð1Þk 1. ð11Þ h These equations reduce the system of linear equations from (2k + 2) by (2k + 2) to (k + 1) by (k + 1). Obviously, in order to use the Eq. (11) and solve them for k + 1 unknowns m0, m1, . . . , mk, we need (as usual) two additional linear equations. In the literature, these equations are called the end conditions. The best set of end conditions for our purpose are the first derivative end conditions at two end points a and b. Suppose that y 0 (a) = a and y 0 (b) = b are given. Then by setting m0 = a and mk = b, we can easily solve the following
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
11
(k 1) by (k 1) linear tridiagonal equations to obtain a unique set of solutions for m1, m2, . . . , mk1: 8 < 10m1 þ m2 ¼ b1 a; ð12Þ mi1 þ 10mi þ miþ1 ¼ bi ; i ¼ 2ð1Þk 2; : mk1 þ 10mk ¼ bk1 b; ðI iþ1 I i Þ. where bi ¼ 12 h2 When m0, m1, . . . , mk are known, then we can use (7) or (8) to compute S0, S1, . . . , Sk. But we need another additional given value for y(a) or y(b). If y0 = y(a) is given, then by considering S0 = y0 = y(a) and using the Eq. (7), we can obtain the values of S1, S2, . . . , Sk. With these two sets of values of mÕs and SÕs, and using S(x) from (4) with (5), we can approximate y(x) at any x 2 [a, b]. If the third end condition y0 = y(a) is not given and y(a) is not available, then we can use the simple difference forward or backward formula for the first derivative approximation and use one of the following relations: S 1 S 0 ¼ hm0
or
S 1 S 0 ¼ hm1
ð13Þ
as an additional equation to mix with the Eq. (7) to find all k + 1 unknowns S0, S1, . . . , Sk. 3. Integro cubic splines with second derivative representation Instead of the cubic Hermite polynomial (4), we can construct a cubic spline S(x) in terms of the second derivatives. The following is another representation of S(x), for x 2 [xi1, xi]; i = 1(1)k, in terms of Mi1, Mi, Si1 and Si: 1 SðxÞ ¼ fðx xi1 Þ3 M i ðx xi Þ3 M i1 g 6h 1 h2 h2 ðM i1 M i Þ þ ðS i S i1 Þ ðx xi Þ þ S i M i . ð14Þ þ h 6 6 It is easy to see that S(x) and S00 (x) are continuous at the interior points xi; i = 1(1)k 1. In order to have a continuous S 0 (x) at xi, we must have the next consistency relations which is in terms of Mi1, Mi, Mi+1, and Si1, Si and Si+1. 6 ð15Þ M i1 þ 4M i þ M iþ1 ¼ 2 ðS i1 2S i þ S iþ1 Þ; i ¼ 1ð1Þk 1. h By using the integrals of (14) in (3), we obtain the following relations between MÕs, SÕs and IÕs: h3 h ðM i1 þ M i Þ þ ðS i1 þ S i Þ ¼ I i ; i ¼ 1ð1Þk; 2 24 h3 h ðM i þ M iþ1 Þ þ ðS i þ S iþ1 Þ ¼ I iþ1 ; i ¼ 0ð1Þk 1. 2 24
ð16Þ ð17Þ
12
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
By adding and subtracting (16) and (17) we obtain h3 h ðM i1 þ 2M i þ M iþ1 Þ þ ðS i1 þ 2S i þ S iþ1 Þ ¼ I iþ1 þ I i ; ð18Þ 2 24 h3 h ð19Þ ðM iþ1 M i1 Þ þ ðS iþ1 S i1 Þ ¼ I iþ1 I i . 2 24 Unlike the first derivative representation, here, we cannot eliminate SÕs between (15)–(19) to obtain a relation similar to (11) without SÕs. So, the advantage of the Hermite representation or using the relations (11) is that, to find the parameters, we need only to solve a simple (k 1) by (k 1) tridiagonal equations, not a system with a full matrix of order (2k + 2) by (2k + 2). Even using the natural cubic spline end conditions M0 = Mk = 0, or knowing the second derivative values at x = a and x = b and using them as two additional end conditions M0 = y00 (a) and Mk = y00 (b), do not lead us to solving a system of (k 1) by (k 1) tridiagonal equations. So, to construct S(x) by (14) with any end conditions, we have to solve a system of linear equations with a full matrix of order (2k + 2) by (2k + 2).
4. Other end conditions In practice we do not expect to know the values of the first derivative of the function at two end points. But these two end conditions m0 = a and mk = b reduce the (2k + 2) by (2k + 2) system to an easy system of (k 1) by (k 1) tridiagonal equations. Furthermore, if y 2 C4[a, b], then the cubic spline with these two end conditions has the highest order of convergence O(h4). In the literature, we have so many other end conditions for this purpose. Even with the same accuracy and with the same order of convergence O(h4). See for example, the end conditions that are mentioned in Behforooz [1] and the references therein. Most of those end conditions are in terms of the derivatives of the spline functions and the function values. But here, the problem is that we do not have the function values at the knots and when we use those types of end conditions, we have to solve a system of equations with full matrix of order (2k + 2) by (2k + 2). 5. Minimum property The traditional cubic spline S(x) with interpolation conditions S(xi) = y(xi); i = 0(1)k with either end conditions M0 = Mk = 0 (for natural cubic spline) or the first derivative end conditions m0 = y 0 (a), mk = y 0 (b), minimizes the integral Z b ½f 00 ðxÞ2 dx ð20Þ a
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
13
over all the functions y = f(x) 2 C2[a, b]. To investigate this, we start as follows: Z b 2 ½f 00 ðxÞ S 00 ðxÞ dx a Z b Z b Z b 2 2 ½f 00 ðxÞ dx þ ½S 00 ðxÞ dx 2 f 00 ðxÞS 00 ðxÞ dx ¼ a a a Z b Z b Z b ½f 00 ðxÞ2 dx ½S 00 ðxÞ2 dx 2 ½f 00 ðxÞ S 00 ðxÞS 00 ðxÞ dx. ¼ a
a
a
Applying the integration by parts on the last integral we obtain Z b ½f 00 ðxÞ S 00 ðxÞS 00 ðxÞ dx a Z b b 0 00 0 ½f 0 ðxÞ S 0 ðxÞS 000 ðxÞ dx ¼ ½f ðxÞ S ðxÞS ðxÞja a b
¼ ½f 0 ðxÞ S 0 ðxÞS 00 ðxÞja
k X i¼1
0
0
¼ ½f ðxÞ S ðxÞS
00
b ðxÞja
k X i¼1
Z
xi
Ci
½f 0 ðxÞ S 0 ðxÞ dx
xi1
xi C i ½f ðxÞ SðxÞ ;
ð21Þ
xi1
where, Ci = S000 (x) is a constant in each interval [xi1, xi]. So, for the traditional cubic spline S(x) with interpolation conditions S(xi) = y(xi); i = 0(1)k with either end conditions {M0 = Mk = 0} or {m0 = y 0 (a), mk = y 0 (b)} we have Z b ½f 00 ðxÞ S 00 ðxÞS 00 ðxÞ dx ¼ 0 a
and Z
b 2
½S 00 ðxÞ dx ¼ a
Z
b 2
½f 00 ðxÞ dx
a
Therefore for all f 2 C2[a, b] Z b Z b ½S 00 ðxÞ2 dx 6 ½f 00 ðxÞ2 dx; b
Z
b 2
½f 00 ðxÞ S 00 ðxÞ dx 6 a
Z
b 2
½f 00 ðxÞ dx. a
ð22Þ
a
which is called the minimum property of S(x) with either one of those two sets of end conditions. In other words, S(x) is the smoothest function which has the minimum total curvature. In our integro cubic spline with end conditions m0 = y 0 (a) and mk = y 0 (b), the first bracket in (21) is equal to zero. But since our interpolating and fitting conditions are Z xi Z xi yðxÞ dx ¼ SðxÞ dx ¼ I i ; i ¼ 1ð1Þk xi1
xi1
14
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
not S(xi) = y(xi), then we do not make the second bracket equal to zero. But the cubic spline with end conditions m0 = y 0 (a) and mk = y 0 (b) has the highest order of convergence and for f 2 C4[a, b], jf ðxÞ SðxÞj ¼ Oðh4 Þ. The numerical examples below, show that the order of convergence is high enough to see that jf(x) S(x)j is almost zero at xi; i = 0(1)k and this implies that the second bracket in Eq. (21) is almost zero and our integro cubic spline has minimum total curvature.
6. Numerical examples We examined this method on several functions and the accuracy of the approximation is very high. I compared the results of this integro cubic spline with my other ordinary higher degree splines and the results are reasonably comparable, see the numerical results in Behforooz [3,4]. Here we show the results of one of those functions. Let S be the integro cubic spline to interpolate the function y(x) = exp(x) at the knots xi = 0.05i, i = 0(1)20, and satisfies the end conditions m0 = exp(0) = 1 and m20 = exp(1). In Table 1 we have listed the values of the error E(x) = jS(x) exp(x)j computed at some points in the interval [0, 1]. For comparison, we also listed the errors which are computed by using the ordinary natural cubic spline (NCS) and the famous not-a-knot cubic spline. This example shows that our cubic spline S is more accurate than the other two cubic splines (specially at the two end subintervals). It is worth mentioning that, to find the inverse of the coefficient matrix of (11), one may use an algorithm similar to the algorithm of Moon [6], which reduces the number of multiplication from O(n2) to O(n). Also, Delhez in [5], has introduced the mc-spline of degree four, but to construct this spline we have to solve a system Table 1 Comparison of jexp(x) s(x)j for three cubic splines x
N.C.S.
Not-a-knot
E(x)
0.01 0.02 0.07 0.09 0.22 0.36 0.62 0.93 0.96 0.98 0.99
9.95E5 1.23E4 3.36E5 1.34E5 6.17E7 2.04E8 5.43E8 8.92E5 1.37E4 3.32E4 2.74E4
1.54E6 1.84E6 7.02E7 2.82E7 1.83E7 9.95E8 2.84E7 1.74E6 1.76E6 4.54E6 3.87E7
9.00E9 1.62E8 1.21E9 1.04E8 2.65E9 1.07E8 5.61E9 2.30E8 3.61E8 2.35E8 4.32E8
H. Behforooz / Appl. Math. Comput. 175 (2006) 8–15
15
of (2k + 2) by (2k + 2) blocked matrix linear equations which is more difficult than our (k + 1) by (k + 1) diagonal equations (11).
References [1] H. Behforooz, End conditions for cubic spline interpolation derived from integration, Appl. Math. Comput. 29 (1989) 231–246. [2] H. Behforooz, A new approach to spline functions, Appl. Numer. Math. 13 (1993) 271–276. [3] H. Behforooz, The not-a-knot interpolatory cubic polynomial, Appl. Math. Comput. 52 (1992) 29–35. [4] H. Behforooz, A comparison of the E(3) and not-a-knot cubic spline, Appl. Math. Comput. 72 (1992) 219–223. [5] E.J.M. Delhez, A spline interpolation technique that preserve mass budget, Appl. Math. Lett. 16 (2003) 17–26. [6] B.S. Moon, An explicit solution for the cubic spline interpolation for functions of single variable, J. Appl. Math. Comput. 117 (2001) 251–255.