- Email: [email protected]
A modification of pseudo-spectral method for solving a linear ODEs with singularity
Applied Mathematics and Computation 188 (2007) 1260–1266 www.elsevier.com/locate/amc
A modification of pseudo-spectral method for solving a linear ODEs with singularity E. Babolian a, M. Bromilow b, R. England b, M. Saravi a
b,c,*
Department of Mathematics, Teacher Training University, 599 Taleghani Avenue, Tehran 15618, Iran b Department of Mathematics, OU Walton Hall, Milton Keynes, England, UK c Islamic Azad University-Noor Branch, Noor, Iran
Abstract In this paper, first we introduce, briefly, pseudo-spectral method for numerical solution of ODE’s and focus on those problems in which some of coefficient functions or solution function are singular. Then, by expressing weak and strong aspects of spectral methods to solve these kind of problems, a modified pseudo-spectral method which is more efficient than other spectral methods is suggested. We compare the methods with some numerical examples. 2006 Elsevier Inc. All rights reserved. Keywords: Spectral method; ODE; Singularity; Pseudo-spectral method
1. Introduction The spectral methods arise from the fundamental problem of approximation of a function by interpolation on an interval, and are very much successful for the numerical solution of ordinary or partial differential equations [1]. Since the time of Fourier (1882), spectral representations in analytic study of differential equations are used and their applications for numerical solution of ordinary differential equations refers, at least, to the time of Lanczos [2]. A survey of some applications is given in [3]. Spectral methods may be viewed as an extreme development of the class of discretization scheme for differential equations known generally as the method of weighted residuals (MWR) [5]. The key elements of the MWR are the trial functions (also called approximating functions) which are used as basis functions for a truncated series expansion of the solution, and the test functions (also known as weight functions) which are used to ensure that the differential equation is satisfied as closely as possible by the truncated series expansion. The choice of such functions distinguishes between the three most commonly used spectral schemes, namely, Galerkin, collocation (also called pseudo-spectral) and Tau version. The Tau approach is a modification of Galerkin method that is applicable to problems with non-periodic boundary conditions. In broad *
Corresponding author. Address: Islamic Azad University-Noor Branch, Noor, Iran. E-mail addresses: [email protected] (E. Babolian), [email protected] (M. Saravi).
0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.10.079
E. Babolian et al. / Applied Mathematics and Computation 188 (2007) 1260–1266
1261
terms, Galerkin and Tau methods are implemented in terms of the expansion coefficients [4], whereas collocation methods are implemented in terms of physical space values of unknown function. The basic idea of spectral methods to solve differential equations is to expand the solution function as a finite series of very smooth basis functions, as given y N ðxÞ ¼
N X
ak /k ðxÞ;
ð1:1Þ
k¼0
in which, the best choice of /k, are the eigenfunctions of a singular Sturm–Liouville problem. If the function y belongs to C1[a, b], the produced error of approximation (1.1), when N tends to infinity, approaches zero with exponential rate [1]. This phenomenon is usually referred to as ‘‘spectral accuracy’’. [3]. The accuracy of derivatives obtained by direct, term by term, differentiation of such truncated expansion naturally deteriorates [1], but for low order derivatives and sufficiently high- order truncations this deterioration is negligible. So, if solution function and coefficient functions are analytic on [a, b], spectral methods will be very efficient and suitable. 2. Pseudo-spectral method In this section, we briefly introduce pseudo-spectral method. For this reason, first we consider the following differential equation: m X Ly ¼ fmi ðxÞDi y ¼ f ðxÞ; x 2 ½1; 1; ð2:1Þ i¼0
Ty ¼ C;
ð2:2Þ i
where fi, i = 0, 1, . . . , m, f, are known real functions of x, D denotes ith order of differentiation with respect to x, T is a linear functional of rank N and C 2 Rm. Here (2.2) can be initial, boundary or mixed conditions. The basic of spectral methods to solve this class of equations is to expand the solution function, y, in (2.1) and (2.2), as a finite series of very smooth basis functions, as given below: y N ðxÞ ¼
N X
ak T k ðxÞ;
ð2:3Þ
k¼0
where, fT k ðxÞgk0 is the sequence of Chebyshev polynomials of the first kind. With replacing y by yN in (2.2), we define residual term by rN(x) as follows: rN ðxÞ ¼ Ly N f :
ð2:4Þ N
In spectral methods, the main target is to minimize r (x) as much as possible with regard to (2.2). Implementation of these methods lead to a system of linear equations with N + 1 equations and N + 1 unknowns a0, a 1, . . . , a N . In the rest of this section, we discuss only one of the three spectral methods, namely, collocation (also known as pseudo-spectral) method, also we use Tau method for numerical solution of second order linear differential equations to compare the results with pseudo-spectral method. It is to be noted that this discussion can be extended to the general form Eqs. (2.1) and (2.2). Consider the following differential equation: P ðxÞy 00 þ QðxÞy 0 þ RðxÞy ¼ SðxÞ; yð1Þ ¼ a;
yð1Þ ¼ b:
x 2 ð1; 1Þ;
ð2:5Þ
First, for an arbitrary natural number N, we suppose that the approximate solution of Eq. (2.5) is given by k (2.3), where a = (a0, a1, . . . , aN)t 2 RN+1 is the coefficients vector and fT k ðxÞg0 is the sequence of Chebyshev polynomials of the first kind. Now if we put V ðxÞ ¼
N X k¼0
ak T k ðxÞ;
ð2:6Þ
1262
E. Babolian et al. / Applied Mathematics and Computation 188 (2007) 1260–1266
then corresponding to functions V, V 0 and V00 , we can define matrices A(0), A(1) and A(2) as follows [7]: V ffi Að0Þ ;
Að0Þ ¼ I ðN þ1ÞðN þ1Þ ; ( 1 2i; for j > i; i þ j ¼ odd; ð1Þ ð1Þ 0 ci V ffi A ; Aij ¼ 0 otherwise: 3 2 0 1 0 3 0 ... 6 0 4 0 8 ...7 7 6 7 6 ð1Þ 7 6 0 6 0 . . . A ¼6 7 6 0 8 ...7 5 4 .. . ( 1 ðj iÞjðj þ iÞ; for j > i; i þ j ¼ even; ci V 0 ffi Að2Þ ; ðAð2Þ Þij ¼ 0; otherwise: 3 2 0 0 4 0 32 . . . 6 0 0 24 0 ...7 7 6 7 6 0 0 48 . . . 7; Að2Þ ¼ 6 7 6 6 0 0 ...7 5 4 .. . 2; i ¼ 0; where N P j, i P 0 and ci ¼ 1; i > 0: Now using differential equation (2.5), we define matrix AA(x) as follows: AAðxÞ ¼ P ðxÞAð2Þ þ QðxÞAð1Þ þ RðxÞA; that is,
2
RðxÞ QðxÞ 4P ðxÞ 6 RðxÞ 4QðxÞ 6 6 6 RðxÞ 6 AAðxÞ ¼ 6 6 6 6 O 4
3 3QðxÞ 32P ðxÞ . . . 24P ðxÞ 8QðxÞ . . . 7 7 7 6QðxÞ 48P ðxÞ . . . 7 7 RðxÞ 8QðxÞ . . . 7 7 7 RðxÞ ...7 5 .. .
Hence, differential equation (2.5) converts to N X ai /i ðxÞ SðxÞ;
ð2:7Þ
:
ðN þ1ÞðN þ1Þ
ð2:8Þ
i¼0
in which, /i ðxÞ ¼
i X ðAAÞki T k ðxÞ; k¼0
that is, /0 ðxÞ ¼ RðxÞT 0 ðxÞ; /1 ðxÞ ¼ QðxÞT 0 ðxÞ þ RðxÞT 1 ðxÞ; /2 ðxÞ ¼ 4P ðxÞT 0 ðxÞ þ 4QðxÞT 1 ðxÞ þ RðxÞT 2 ðxÞ; .. .
ð2:9Þ
E. Babolian et al. / Applied Mathematics and Computation 188 (2007) 1260–1266
1263
It must be noted that, if A(k+2); k P 1 be corresponding matrix of (k + 2)th order differentiation of V(x), it follows that [7]: ðkþ2Þ
Aij
ðkÞ
¼ Aij
ðj i kÞðj þ i þ kÞðj þ i kÞ ; 4kðk þ 1Þ
0 6 i; j 6 N :
Now, if we impose boundary condition from (2.5) on V(x) we will have; V ð1Þ ¼ a )
N X
N X
ak T k ð1Þ ¼
k¼0
V ð1Þ ¼ b )
N X
k
ak ð1Þ ¼ a;
k¼0
ak T k ð1Þ ¼
k¼0
N X
ak ¼ b:
k¼0
So "
1 1
1 1 1
1
2 a0 #6 a1 ð1ÞN 6 6 . 6 1 4 ..
... ...
3 7 7 a 7¼ : 7 b 5
ð2:10Þ
aN Relation (2.10) forms a system with two equations and N + 1 unknowns. To construct the remaining N 1 equations we substitute points xj = cos(pj/N), j = 1, 2, . . . , N 1, in (2.8) and put: N X
ai /i ðxj Þ ¼ Sðxj Þ;
j ¼ 1; 2; . . . ; N 1;
ð2:11Þ
i¼0
to obtain N 1 equations. 3. Numerical experiments Now, we consider some ordinary differential equations with Tau method (as representative of Tau and Galerkian methods) and discuss the results. Problem 1. Let us consider y 00 ðxÞ þ xy 0 ðxÞ þ y ¼ x cosðxÞ;
yð1Þ ¼ sinð1Þ;
yð1Þ ¼ sinð1Þ;
with the exact solution y(x) = sin(x). We solved it by Runge–Kutta with different orders and all had nearly a maximum error of 3.0 · 101. We also solved it by Tau method with N = 5, 8, 16 and the maximum error produced using this method is given in Table 1. This table shows the power of spectral method. Problem 2. Consider 1 y 00 ðxÞ þ y 0 ðxÞ ¼ x
8 8 x2
2 ;
x 2 ð0; 1Þ; yð1Þ ¼ 0; y 0 ð0Þ ¼ 0;
7 with the exact solution yðxÞ ¼ 2 lnð8x 2 Þ. This problem was chosen from [6]. It was solved by extrapolation method with maximum error of 108. Here we solved it by Tau method for different values of N, and the results are given in Table 2.
Table 1 N
kyðxÞ ~y s ðxÞk1
5 8 16
2.11 · 105 5.71 · 108 1.11 · 1016
1264
E. Babolian et al. / Applied Mathematics and Computation 188 (2007) 1260–1266
Table 2 N
Maximum error
5 15 16 30 95
5 · 105 2 · 106 8 · 107 5 · 107 8 · 108
As we can see, when N increases, the rate of improvement of accuracy is very low. This is because of smoothless of coefficient function. But, when we solved it by pseudo-spectral method, since coefficient functions are not expanded, the produced error will be better. In the above example by Tau method with N = 95, maximum error was about 8 · 108, but by pseudo-spectral method with N = 18 we come to a maximum error about 4 · 1019. Therefore, in this case, this method is more successful than Tau method. We may have different cases. One in which, at least, one of coefficient functions or solution function is not analytic, or they are analytic. We consider two other cases. Problem 3 (Case 1). Consider 0 1 s s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 3 3 1 1 A expðxÞ; yðxÞ ¼ @1 þ jxj þ x2 y 00 ðxÞ þ jxjy 0 ðxÞ þ x2 4 4 yð1Þ ¼ expð1Þ;
x 2 ½1; 1;
yð1Þ ¼ expð1Þ
with the exact solution y(x) = exp(x). This problem was chosen from [8]. This problem was solved by Tau method, the error produced from expansion of, at least, one of the coefficient functions follows from Table 3 [7]. So, in making system of equations (comparing with approximate error solution function) the error is considerable and causes undesirable effect on final solution. Now, if this problem is solved by pseudo-spectral, since solution function, y, is analytic and there is no need to expand coefficient functions, the produced error from making system is less than system for Tau method. Therefore, in this case, pseudo-spectral method is more successful than the Tau method. Since, for Problem 3 coefficient functions are not analytic it has been solved by pseudo-spectral method. For N = 8, 11, 16 the results are given in Table 4. Now we are going to consider an ODE with solution function and, at least, one of the coefficient functions non-analytic. This problem, also, was chosen from [8].
Table 3 Function
Order of infinity norm of the error Close to irregular points
Far from irregular points
f, discontinuous
1
f 0 , discontinuous
1 N 1 N2
1 N 1 N2 1 N3
f 00 , discontinuous .. .
.. .
.. .
Table 4 N 8 11 16
ky(x) yt(x)k1 6
3.13 · 10 6.40 · 108 3.92 · 108
ky(x) yps(x)k1 3.24 · 108 2.52 · 1012 3.50 · 1018
E. Babolian et al. / Applied Mathematics and Computation 188 (2007) 1260–1266
1265
Problem 4 (Case 2). Consider y 00 ðxÞ þ jxjy 0 ðxÞ þ yðxÞ ¼ j6xj þ jx3 j þ 3x3 ;
x 2 ½1; 1; yð1Þ ¼ yð1Þ ¼ 1;
with the exact solution y(x) = jx3j. In this case, since the order of infinity norm of the producing error from approximation of solution function and, at least, one of the coefficient functions follows from Table 3, the errors are not accurate enough. Due to the fact that the representing matrices A, A(1), A(2), and AA, because of existing a non-analytic solution function, accompany the error which is indispensable. A modified spectral method is suggested in [8], in such a way that it considerably decreases the error in setting up the system of equations. Before we go further, let’s consider another problem. Problem 5. Consider 1
1
y 00 ðxÞ þ ex y 0 ðxÞ þ yðxÞ ¼ 6x þ x3 þ 3x2 ex ; with the exact solution y(x) = x3. Here, we have essential singularity. If we choose N as an even number Tau method, pseudo- spectral method and even modified method introduced in [8], do not work. Here we are going to introduce another modification which produces better results and works very well whenever those methods do not work. The idea comes from the point that we checked coefficients matrix we found that the difficulty arises from N/2-row. This is because of initial condition which produces (2.10) and elements in this row, which are infinity. If we take the infinity, as a factor, out of determinant of this matrix we get a matrix with two identical rows. Hence determinant becomes zero. To avoid this difficulty we choose one of these elements close to 1, for example, 0.99999 and continue the process. We examined this choice and solved this problem with this modification and obtained the results given in Table 5. Now, we go back to Problem 4. Here, again we applied this modification to this example and obtained the results given in Table 6. Here, this method gives slightly better results than given in [8], where was notified by ms1. The notations ms1 and ms2 mean modified methods in [8] and here, respectively. Main question is, why we will have better results? The answer is: Because of condition number of the coefficient matrix. When we considered condition number of the coefficient matrix before and after using this change, we found substaintial difference. For example, in Problem 5 after using this method condition number changes from infinity to 4.6058 · 109. We also checked other problems, and in all cases condition numbers reduces, at least, by a factor of 107. Let’s consider another example. Problem 6. Consider 1 1 y 00 ðxÞ y 0 ðxÞ þ yðxÞ ¼ jxj; x x with the exact solution y(x) = xjxj. We tested it with different values of N, and the results are given in Table 7. Table 5 N
kyðxÞ ~y ðxÞkmps
5 8 12
2.2204 · 1015 1.6653 · 1015 1.3843 · 1015
Table 6 N 8 15 20
kyðxÞ ~y t ðxÞk1 2
8.98 · 10 1.54 · 102 1.68 · 102
kyðxÞ ~y ps ðxÞk1 1
1.21 · 10 1.76 · 102 1.92 · 102
kyðxÞ ~y ms1 ðxÞk1 3
3.82 · 10 5.61 · 104 1.81 · 104
kyðxÞ ~y ms2 ðxÞk1 2.41 · 103 1.85 · 104 1.28 · 104
1266
E. Babolian et al. / Applied Mathematics and Computation 188 (2007) 1260–1266
Table 7 N
kyðxÞ ~y t ðxÞk1
kyðxÞ ~y ps ðxÞk1
kyðxÞ ~y ms2 ðxÞk1
5 8 9 17
8.31 · 102 8.75 · 101 1.54 · 102 1.12 · 102
7.64 · 102 8.86 · 101 3.97 · 102 2.05 · 102
1.92 · 102 3.23 · 102 1.00 · 103 5.47 · 104
N
kyðxÞ ~y t ðxÞk1
kyðxÞ ~y ps ðxÞk1
kyðxÞ ~y ms2 ðxÞk1
5 8 16
3.26 · 101 2.42 · 101 2.05 · 101
5.99 · 101 9.59 · 101 5.27 · 101
3.52 · 101 1.52 · 101 7.43 · 102
Table 8
As we can see the results obtained by the proposed modification are better. In particular, when N is even. In this example again we get a lower condition number. Introduction of this method will be ended by representing another example with non-analytical solution. Problem 7. Consider y 00 ðxÞ þ
1 0 1 y ðxÞ þ yðxÞ ¼ þ jxj; jxj x
with the exact solution y(x) = jxj. We used these methods for different N, and the results are given in Table 8. The results are not very accurate, but still the manner we described above governs. Here condition number changes from 3.3466 · 1017 to 4.1590 · 109. Acknowledgements The authors would like to thank Professor R. Saadati for giving useful comments and suggestions for the improvement of this paper. References [1] [2] [3] [4] [5] [6]
C. Canuto, M. Hussaini, A. Quarteroni, T. Zang, Spectral Methods in Fluid Dynamics, Springer, Berlin, 1988. C. Lanczos, Trigonometric interpolation of empirical and analytical functions, J. Math. Phys. 17 (1938) 123–129. D. Gottlieb, S. Orszag, Numerical Analysis of Spectral Methods, Theory and Applications, SIAM, Philadelphia, PA, 1977. L.M. Delves, J.L. Mohamed, Computational Methods for Integral Equations, Cambridge University Press, 1985. A. Finlason, L.E. Scrinven, The method of weighted residuals, Appl. Mech. Rev. 19 (1996) 735–748. U.M. Ascher, R.M. Mahheij, R.D. Russell, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, Perentice-Hall Inc., 1988. [7] B. Fornberg, A Practical Guide to Pseudo-Spectral Methods, Cambridge University Press, Cambridge, 1996. [8] E. Babolian, M.M. Hosseini, A modified spectral method for numerical solution of ordinary differential equations with non-analytical solution, Appl. Math. Comput. 132 (2002) 341–351.
A modification of pseudo-spectral method for solving a linear ODEs with singularity E. Babolian a, M. Bromilow b, R. England b, M. Saravi a
b,c,*
Department of Mathematics, Teacher Training University, 599 Taleghani Avenue, Tehran 15618, Iran b Department of Mathematics, OU Walton Hall, Milton Keynes, England, UK c Islamic Azad University-Noor Branch, Noor, Iran
Abstract In this paper, first we introduce, briefly, pseudo-spectral method for numerical solution of ODE’s and focus on those problems in which some of coefficient functions or solution function are singular. Then, by expressing weak and strong aspects of spectral methods to solve these kind of problems, a modified pseudo-spectral method which is more efficient than other spectral methods is suggested. We compare the methods with some numerical examples. 2006 Elsevier Inc. All rights reserved. Keywords: Spectral method; ODE; Singularity; Pseudo-spectral method
1. Introduction The spectral methods arise from the fundamental problem of approximation of a function by interpolation on an interval, and are very much successful for the numerical solution of ordinary or partial differential equations [1]. Since the time of Fourier (1882), spectral representations in analytic study of differential equations are used and their applications for numerical solution of ordinary differential equations refers, at least, to the time of Lanczos [2]. A survey of some applications is given in [3]. Spectral methods may be viewed as an extreme development of the class of discretization scheme for differential equations known generally as the method of weighted residuals (MWR) [5]. The key elements of the MWR are the trial functions (also called approximating functions) which are used as basis functions for a truncated series expansion of the solution, and the test functions (also known as weight functions) which are used to ensure that the differential equation is satisfied as closely as possible by the truncated series expansion. The choice of such functions distinguishes between the three most commonly used spectral schemes, namely, Galerkin, collocation (also called pseudo-spectral) and Tau version. The Tau approach is a modification of Galerkin method that is applicable to problems with non-periodic boundary conditions. In broad *
Corresponding author. Address: Islamic Azad University-Noor Branch, Noor, Iran. E-mail addresses: [email protected] (E. Babolian), [email protected] (M. Saravi).
0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.10.079
E. Babolian et al. / Applied Mathematics and Computation 188 (2007) 1260–1266
1261
terms, Galerkin and Tau methods are implemented in terms of the expansion coefficients [4], whereas collocation methods are implemented in terms of physical space values of unknown function. The basic idea of spectral methods to solve differential equations is to expand the solution function as a finite series of very smooth basis functions, as given y N ðxÞ ¼
N X
ak /k ðxÞ;
ð1:1Þ
k¼0
in which, the best choice of /k, are the eigenfunctions of a singular Sturm–Liouville problem. If the function y belongs to C1[a, b], the produced error of approximation (1.1), when N tends to infinity, approaches zero with exponential rate [1]. This phenomenon is usually referred to as ‘‘spectral accuracy’’. [3]. The accuracy of derivatives obtained by direct, term by term, differentiation of such truncated expansion naturally deteriorates [1], but for low order derivatives and sufficiently high- order truncations this deterioration is negligible. So, if solution function and coefficient functions are analytic on [a, b], spectral methods will be very efficient and suitable. 2. Pseudo-spectral method In this section, we briefly introduce pseudo-spectral method. For this reason, first we consider the following differential equation: m X Ly ¼ fmi ðxÞDi y ¼ f ðxÞ; x 2 ½1; 1; ð2:1Þ i¼0
Ty ¼ C;
ð2:2Þ i
where fi, i = 0, 1, . . . , m, f, are known real functions of x, D denotes ith order of differentiation with respect to x, T is a linear functional of rank N and C 2 Rm. Here (2.2) can be initial, boundary or mixed conditions. The basic of spectral methods to solve this class of equations is to expand the solution function, y, in (2.1) and (2.2), as a finite series of very smooth basis functions, as given below: y N ðxÞ ¼
N X
ak T k ðxÞ;
ð2:3Þ
k¼0
where, fT k ðxÞgk0 is the sequence of Chebyshev polynomials of the first kind. With replacing y by yN in (2.2), we define residual term by rN(x) as follows: rN ðxÞ ¼ Ly N f :
ð2:4Þ N
In spectral methods, the main target is to minimize r (x) as much as possible with regard to (2.2). Implementation of these methods lead to a system of linear equations with N + 1 equations and N + 1 unknowns a0, a 1, . . . , a N . In the rest of this section, we discuss only one of the three spectral methods, namely, collocation (also known as pseudo-spectral) method, also we use Tau method for numerical solution of second order linear differential equations to compare the results with pseudo-spectral method. It is to be noted that this discussion can be extended to the general form Eqs. (2.1) and (2.2). Consider the following differential equation: P ðxÞy 00 þ QðxÞy 0 þ RðxÞy ¼ SðxÞ; yð1Þ ¼ a;
yð1Þ ¼ b:
x 2 ð1; 1Þ;
ð2:5Þ
First, for an arbitrary natural number N, we suppose that the approximate solution of Eq. (2.5) is given by k (2.3), where a = (a0, a1, . . . , aN)t 2 RN+1 is the coefficients vector and fT k ðxÞg0 is the sequence of Chebyshev polynomials of the first kind. Now if we put V ðxÞ ¼
N X k¼0
ak T k ðxÞ;
ð2:6Þ
1262
E. Babolian et al. / Applied Mathematics and Computation 188 (2007) 1260–1266
then corresponding to functions V, V 0 and V00 , we can define matrices A(0), A(1) and A(2) as follows [7]: V ffi Að0Þ ;
Að0Þ ¼ I ðN þ1ÞðN þ1Þ ; ( 1 2i; for j > i; i þ j ¼ odd; ð1Þ ð1Þ 0 ci V ffi A ; Aij ¼ 0 otherwise: 3 2 0 1 0 3 0 ... 6 0 4 0 8 ...7 7 6 7 6 ð1Þ 7 6 0 6 0 . . . A ¼6 7 6 0 8 ...7 5 4 .. . ( 1 ðj iÞjðj þ iÞ; for j > i; i þ j ¼ even; ci V 0 ffi Að2Þ ; ðAð2Þ Þij ¼ 0; otherwise: 3 2 0 0 4 0 32 . . . 6 0 0 24 0 ...7 7 6 7 6 0 0 48 . . . 7; Að2Þ ¼ 6 7 6 6 0 0 ...7 5 4 .. . 2; i ¼ 0; where N P j, i P 0 and ci ¼ 1; i > 0: Now using differential equation (2.5), we define matrix AA(x) as follows: AAðxÞ ¼ P ðxÞAð2Þ þ QðxÞAð1Þ þ RðxÞA; that is,
2
RðxÞ QðxÞ 4P ðxÞ 6 RðxÞ 4QðxÞ 6 6 6 RðxÞ 6 AAðxÞ ¼ 6 6 6 6 O 4
3 3QðxÞ 32P ðxÞ . . . 24P ðxÞ 8QðxÞ . . . 7 7 7 6QðxÞ 48P ðxÞ . . . 7 7 RðxÞ 8QðxÞ . . . 7 7 7 RðxÞ ...7 5 .. .
Hence, differential equation (2.5) converts to N X ai /i ðxÞ SðxÞ;
ð2:7Þ
:
ðN þ1ÞðN þ1Þ
ð2:8Þ
i¼0
in which, /i ðxÞ ¼
i X ðAAÞki T k ðxÞ; k¼0
that is, /0 ðxÞ ¼ RðxÞT 0 ðxÞ; /1 ðxÞ ¼ QðxÞT 0 ðxÞ þ RðxÞT 1 ðxÞ; /2 ðxÞ ¼ 4P ðxÞT 0 ðxÞ þ 4QðxÞT 1 ðxÞ þ RðxÞT 2 ðxÞ; .. .
ð2:9Þ
E. Babolian et al. / Applied Mathematics and Computation 188 (2007) 1260–1266
1263
It must be noted that, if A(k+2); k P 1 be corresponding matrix of (k + 2)th order differentiation of V(x), it follows that [7]: ðkþ2Þ
Aij
ðkÞ
¼ Aij
ðj i kÞðj þ i þ kÞðj þ i kÞ ; 4kðk þ 1Þ
0 6 i; j 6 N :
Now, if we impose boundary condition from (2.5) on V(x) we will have; V ð1Þ ¼ a )
N X
N X
ak T k ð1Þ ¼
k¼0
V ð1Þ ¼ b )
N X
k
ak ð1Þ ¼ a;
k¼0
ak T k ð1Þ ¼
k¼0
N X
ak ¼ b:
k¼0
So "
1 1
1 1 1
1
2 a0 #6 a1 ð1ÞN 6 6 . 6 1 4 ..
... ...
3 7 7 a 7¼ : 7 b 5
ð2:10Þ
aN Relation (2.10) forms a system with two equations and N + 1 unknowns. To construct the remaining N 1 equations we substitute points xj = cos(pj/N), j = 1, 2, . . . , N 1, in (2.8) and put: N X
ai /i ðxj Þ ¼ Sðxj Þ;
j ¼ 1; 2; . . . ; N 1;
ð2:11Þ
i¼0
to obtain N 1 equations. 3. Numerical experiments Now, we consider some ordinary differential equations with Tau method (as representative of Tau and Galerkian methods) and discuss the results. Problem 1. Let us consider y 00 ðxÞ þ xy 0 ðxÞ þ y ¼ x cosðxÞ;
yð1Þ ¼ sinð1Þ;
yð1Þ ¼ sinð1Þ;
with the exact solution y(x) = sin(x). We solved it by Runge–Kutta with different orders and all had nearly a maximum error of 3.0 · 101. We also solved it by Tau method with N = 5, 8, 16 and the maximum error produced using this method is given in Table 1. This table shows the power of spectral method. Problem 2. Consider 1 y 00 ðxÞ þ y 0 ðxÞ ¼ x
8 8 x2
2 ;
x 2 ð0; 1Þ; yð1Þ ¼ 0; y 0 ð0Þ ¼ 0;
7 with the exact solution yðxÞ ¼ 2 lnð8x 2 Þ. This problem was chosen from [6]. It was solved by extrapolation method with maximum error of 108. Here we solved it by Tau method for different values of N, and the results are given in Table 2.
Table 1 N
kyðxÞ ~y s ðxÞk1
5 8 16
2.11 · 105 5.71 · 108 1.11 · 1016
1264
E. Babolian et al. / Applied Mathematics and Computation 188 (2007) 1260–1266
Table 2 N
Maximum error
5 15 16 30 95
5 · 105 2 · 106 8 · 107 5 · 107 8 · 108
As we can see, when N increases, the rate of improvement of accuracy is very low. This is because of smoothless of coefficient function. But, when we solved it by pseudo-spectral method, since coefficient functions are not expanded, the produced error will be better. In the above example by Tau method with N = 95, maximum error was about 8 · 108, but by pseudo-spectral method with N = 18 we come to a maximum error about 4 · 1019. Therefore, in this case, this method is more successful than Tau method. We may have different cases. One in which, at least, one of coefficient functions or solution function is not analytic, or they are analytic. We consider two other cases. Problem 3 (Case 1). Consider 0 1 s s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 3 3 1 1 A expðxÞ; yðxÞ ¼ @1 þ jxj þ x2 y 00 ðxÞ þ jxjy 0 ðxÞ þ x2 4 4 yð1Þ ¼ expð1Þ;
x 2 ½1; 1;
yð1Þ ¼ expð1Þ
with the exact solution y(x) = exp(x). This problem was chosen from [8]. This problem was solved by Tau method, the error produced from expansion of, at least, one of the coefficient functions follows from Table 3 [7]. So, in making system of equations (comparing with approximate error solution function) the error is considerable and causes undesirable effect on final solution. Now, if this problem is solved by pseudo-spectral, since solution function, y, is analytic and there is no need to expand coefficient functions, the produced error from making system is less than system for Tau method. Therefore, in this case, pseudo-spectral method is more successful than the Tau method. Since, for Problem 3 coefficient functions are not analytic it has been solved by pseudo-spectral method. For N = 8, 11, 16 the results are given in Table 4. Now we are going to consider an ODE with solution function and, at least, one of the coefficient functions non-analytic. This problem, also, was chosen from [8].
Table 3 Function
Order of infinity norm of the error Close to irregular points
Far from irregular points
f, discontinuous
1
f 0 , discontinuous
1 N 1 N2
1 N 1 N2 1 N3
f 00 , discontinuous .. .
.. .
.. .
Table 4 N 8 11 16
ky(x) yt(x)k1 6
3.13 · 10 6.40 · 108 3.92 · 108
ky(x) yps(x)k1 3.24 · 108 2.52 · 1012 3.50 · 1018
E. Babolian et al. / Applied Mathematics and Computation 188 (2007) 1260–1266
1265
Problem 4 (Case 2). Consider y 00 ðxÞ þ jxjy 0 ðxÞ þ yðxÞ ¼ j6xj þ jx3 j þ 3x3 ;
x 2 ½1; 1; yð1Þ ¼ yð1Þ ¼ 1;
with the exact solution y(x) = jx3j. In this case, since the order of infinity norm of the producing error from approximation of solution function and, at least, one of the coefficient functions follows from Table 3, the errors are not accurate enough. Due to the fact that the representing matrices A, A(1), A(2), and AA, because of existing a non-analytic solution function, accompany the error which is indispensable. A modified spectral method is suggested in [8], in such a way that it considerably decreases the error in setting up the system of equations. Before we go further, let’s consider another problem. Problem 5. Consider 1
1
y 00 ðxÞ þ ex y 0 ðxÞ þ yðxÞ ¼ 6x þ x3 þ 3x2 ex ; with the exact solution y(x) = x3. Here, we have essential singularity. If we choose N as an even number Tau method, pseudo- spectral method and even modified method introduced in [8], do not work. Here we are going to introduce another modification which produces better results and works very well whenever those methods do not work. The idea comes from the point that we checked coefficients matrix we found that the difficulty arises from N/2-row. This is because of initial condition which produces (2.10) and elements in this row, which are infinity. If we take the infinity, as a factor, out of determinant of this matrix we get a matrix with two identical rows. Hence determinant becomes zero. To avoid this difficulty we choose one of these elements close to 1, for example, 0.99999 and continue the process. We examined this choice and solved this problem with this modification and obtained the results given in Table 5. Now, we go back to Problem 4. Here, again we applied this modification to this example and obtained the results given in Table 6. Here, this method gives slightly better results than given in [8], where was notified by ms1. The notations ms1 and ms2 mean modified methods in [8] and here, respectively. Main question is, why we will have better results? The answer is: Because of condition number of the coefficient matrix. When we considered condition number of the coefficient matrix before and after using this change, we found substaintial difference. For example, in Problem 5 after using this method condition number changes from infinity to 4.6058 · 109. We also checked other problems, and in all cases condition numbers reduces, at least, by a factor of 107. Let’s consider another example. Problem 6. Consider 1 1 y 00 ðxÞ y 0 ðxÞ þ yðxÞ ¼ jxj; x x with the exact solution y(x) = xjxj. We tested it with different values of N, and the results are given in Table 7. Table 5 N
kyðxÞ ~y ðxÞkmps
5 8 12
2.2204 · 1015 1.6653 · 1015 1.3843 · 1015
Table 6 N 8 15 20
kyðxÞ ~y t ðxÞk1 2
8.98 · 10 1.54 · 102 1.68 · 102
kyðxÞ ~y ps ðxÞk1 1
1.21 · 10 1.76 · 102 1.92 · 102
kyðxÞ ~y ms1 ðxÞk1 3
3.82 · 10 5.61 · 104 1.81 · 104
kyðxÞ ~y ms2 ðxÞk1 2.41 · 103 1.85 · 104 1.28 · 104
1266
E. Babolian et al. / Applied Mathematics and Computation 188 (2007) 1260–1266
Table 7 N
kyðxÞ ~y t ðxÞk1
kyðxÞ ~y ps ðxÞk1
kyðxÞ ~y ms2 ðxÞk1
5 8 9 17
8.31 · 102 8.75 · 101 1.54 · 102 1.12 · 102
7.64 · 102 8.86 · 101 3.97 · 102 2.05 · 102
1.92 · 102 3.23 · 102 1.00 · 103 5.47 · 104
N
kyðxÞ ~y t ðxÞk1
kyðxÞ ~y ps ðxÞk1
kyðxÞ ~y ms2 ðxÞk1
5 8 16
3.26 · 101 2.42 · 101 2.05 · 101
5.99 · 101 9.59 · 101 5.27 · 101
3.52 · 101 1.52 · 101 7.43 · 102
Table 8
As we can see the results obtained by the proposed modification are better. In particular, when N is even. In this example again we get a lower condition number. Introduction of this method will be ended by representing another example with non-analytical solution. Problem 7. Consider y 00 ðxÞ þ
1 0 1 y ðxÞ þ yðxÞ ¼ þ jxj; jxj x
with the exact solution y(x) = jxj. We used these methods for different N, and the results are given in Table 8. The results are not very accurate, but still the manner we described above governs. Here condition number changes from 3.3466 · 1017 to 4.1590 · 109. Acknowledgements The authors would like to thank Professor R. Saadati for giving useful comments and suggestions for the improvement of this paper. References [1] [2] [3] [4] [5] [6]
C. Canuto, M. Hussaini, A. Quarteroni, T. Zang, Spectral Methods in Fluid Dynamics, Springer, Berlin, 1988. C. Lanczos, Trigonometric interpolation of empirical and analytical functions, J. Math. Phys. 17 (1938) 123–129. D. Gottlieb, S. Orszag, Numerical Analysis of Spectral Methods, Theory and Applications, SIAM, Philadelphia, PA, 1977. L.M. Delves, J.L. Mohamed, Computational Methods for Integral Equations, Cambridge University Press, 1985. A. Finlason, L.E. Scrinven, The method of weighted residuals, Appl. Mech. Rev. 19 (1996) 735–748. U.M. Ascher, R.M. Mahheij, R.D. Russell, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, Perentice-Hall Inc., 1988. [7] B. Fornberg, A Practical Guide to Pseudo-Spectral Methods, Cambridge University Press, Cambridge, 1996. [8] E. Babolian, M.M. Hosseini, A modified spectral method for numerical solution of ordinary differential equations with non-analytical solution, Appl. Math. Comput. 132 (2002) 341–351.