- Email: [email protected]
Chebyshev spectral methods for solving two-point boundary value problems arising in heat transfer
COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 70 (1988) 103-121 NORTH-HOLLAND
CHEBYSHEV SPECTRAL METHODS FOR SOLVING TWO-POINT BOUNDARY VALUE PROBLEMS ARISING IN HEAT TRANSFER Andreas KARAGEORGHIS Department of Mathematics, The University College of Wales, Aberystwyth, Dyfed SY23 3BZ, United Kingdom Received 18 January 1988 The performance of Chebyshev spectral methods for solving nonlinear two-point boundary value problems is tested on several such problems arising in heat transfer. The collocation and the Galerkin spectral formulations are examined in detail and their various computational aspects compared. The accuracy of these methods is compared with that of standard software packages on a number of problems which appear in the literature.
1. Introduction Problems which arise in heat transfer can often be reformulated as nonlinear two-point boundary value problems for ordinary differential equations. These can be solved numerically via traditional numerical methods which have been implemented in general-purpose software packages. Three such packages were discussed in detail in a paper by Fairweather and Vedha-Nayagam [1]. In it, the authors analyse the performance of the packages DVCPR[2], COLSYS [3], and D ~ [2] for solving two-point boundary value problems arising in heat transfer. The aim of this paper is to investigate the use of Chebyshev spectral methods, in both a collocation and a Galerkin framework, for solving nonlinear two-point boundary value problems arising in heat transfer. In Chebyshev spectral methods, the solution of a differential equation is represented by a truncated series of Chebyshev polynomials. These polynomials are the eigenfunctions of a singular Sturm-Liouville problem and for problems with smooth solutions, such expansions converge faster than any fmite power of 1/n as n--, oo. Chebyshev collocation and Chebyshev-Galerkin methods have been used successfully in two recent studies [4,5] for the solution of a highly nonlinear two-point boundary value problem in modelling viscoelastic flows. In the present work, Chebyshev collocation methods are compared to Chebyshev-Galerkin methods on the set of problems solved in [1] using the packages DVCPR,COLSYS,and D~re (which are based on the trapezoidal rule witb deferred corrections, spline collocation, and multiple shooting, respectively). 2. Chebyshev S l ~ r a i methods Con~ider the two-point boundary value problem
F(x, u(x), u
on the finite interval [0, L],
0045-7825188/$3.50 ~ 1988, Elsevier Science Publishers B.V. (North-Holland)
(2.1)
104
A. Karageorghis, Chebyshev spectral methods
subject to the boundary conditions (2.2)
li[u(x)][x= o = 0 ,
i= 1,2,...,m
I,
rj[u(x)lL= L = O,
j = 1, 2 , . . . , m 2 ( f m - m , ) ,
(2.3)
where F(t~, h , - - . , tin+2) is a nonlinear function of the arguments {ti}i=~.2..... m+2, and {li}i--~,2, ml and {r~},=~, ,m are linear differential operators. The al~P'roaches used in2the" pr2esent paper are formulated in the spirit of the work of Gottlieb and Orszag [6] and are based on the idea of approximating the solution of a two-point boundary value problem by a truncated series expansion of the form N
ti(x) = ~ a , T ~ ( x )
(2.4)
nffi0
The trial functions T ~ ( x ) , n = O, 1 , . . . ,
N , are the shifted Chebyshev polynomials
(2.5) (2.6)
(y) = cos(n cos- y).
In the approximation (2.4) the coefficients a, are unknown and are chosen so that (i) ti(x) satisfies the ordinary differential equation in either the collocation sense (see Section 2.1) or the Galerkin sense (see Section 2.2), and (ii) ti(x) satisfies the associated boundary conditions. 2.1. Collocation f o r m u l a t i o n
In the collocation framework, the approximation (2.4) is forced to satisfy the differential equation at a finite set of points. In accordance with the recommendations of Gottlieb et al. [7] these points are chosen to be a subset of the points in the interval [0, L], where T ~ ( x ) reaches its extrema, i.e. the points Xk, k ffi O, 1 , . . . , N , where
L(
Xk= ~
1+COS'a" 1 - - ~
,
(2.7)
k=0,1,...,N.
The reason for this choice of collocation points is that, in principle, one can make use of FFTs in constructing the approximation ([7]). The ordinary differential equation (2.1) is satisfied in the collocation sense if F(x k, ~(Xk), f~t~)(Xk), . . . , f~tm)(Xk) ) = O,
k = m~, . . . , N -
m 2.
(2.8)
The satisfaction of the two sets of boundary conditions (2.2) and (2.3) provides m ffi m~ + m 2 extra linear equations, the total number of equations therefore being N + 1, which is sufficient to determine the N + 1 unknown coefficients in (2.4). By omitting the m I points in (2.7) nearest 0 (including 0) and the m 2 points in (2.7) nearest L (including L) we wish the boundary conditions to dominate over the differential equation near the end points (see also
[41).
A. Karageorghis, Chebyshev spectral methods
105
2.2. Galerkin formulation
In the Galerkin framework (see [8]) the approximation (2.4) is required to satisfy (1) (F(x, ~ ( x ) , . . . , ~(")(x)), r~(x)) = 0,
n = 0, 1 , . . . , N - m, - m2,
(2.9)
where
1 fo~(x-~ ( 1 - -~)) x~ -''2 v ( x ) w ( x ) d x ,
(V(X), W(X)) = -~
(2.10)
and the boundary conditions (2.2), (2.3) l,[~(x)]l~:o =
o,
i
=
1, 2 , . . . ,
m, ,
and rj[~(x)]lx=,=O,
1= 1,2,
....
,m
2 .
In general, equations (2.9) give rise to complicated expressions involving integrals such as (2.11)
( d(P)(x)~(q)(x), TL, (x)) .
It is often possible, however, to avoid performing such integrations, as was demonstrated in [5]. By using the orthogonality properties of Chebyshev polynomials and the property (see
[9]):
r~(x)r~;(x) ffi ½[r~÷,(x)
+ T~_jl(x)],
(2.12)
and in expression (2.11) assuming that N-q
d(q)(x)= ~
a (q)''L" ~jtx), "
qfO,...,m,
(2.13)
j=O
the integrals involved in (2.11) can be expressed as combinations of integrals of the form
(r,~(~)T~(x), r~(,)).
(2.14)
Using property (2.12), expression (2.14) can be rewritten as a combination of terms of the form
( r~ (,), r~ (x)).
(2.15)
Rather than evaluating integrals of the form (2.15) one may use the orthogonality properties of Chebyshev polynomials which, in the case of (2.15), are equivalent to equating the coefficients of Chebyshev polynomials of order n, n = O, 1 , . . . , N - m~ - m 2 . It is therefore possible to express equations (2.9) in terms of expressions of the form (2.15) and then equate the coefficients of the Chebyshev polynomials of the same order. This process
106
A. Karageorghis, Chebyshev spectral methods
gives rise to the set of equations G,(a, a ~'), a¢'-), . . . , a (m)) = O,
i-O,...
(2.16)
,N,
where a Cq)~ - (a¢J)~ a~q), a~q), . . . , U-
=
4q ~ (j+q-2)! (q-1)!j=, (j-l)!
where
(n+j+q-2)! (n+j-1)!
(n + 2j +
q - 2)a,,+2j+q_ 2 ,
1 for n > O , c,,= 0 for n < 0 .
c0 = 2 ,
(2.17)
For the purposes of this study we extend this for a general (finite) L. COROLLARY
2.2.
1
nun
,~ (j+q-2)! ( q - 1)! ~.= ( j - 1)l
(n+j+q-2)! (n + j - 1)!
(n + 2j + q - 2)a,+2j+q_ 2 .
(2.18)
A proof of the above corollan., is given in Appendix A. Application of the corollary greatly facilitates the solution of the system (2.16) since each of the a
Both the collocation and Galerkin formulations described in Sections 2.1 and 2.2 lead to systen~s of equations which are mixed, i.e. some of the equations are linear and the rest are nonlinear, in the unknown coefficients of the approximation (2.4). In the present study such systems were solved using the NAG subroutine C0SPBF[11] which is based on a modification of the hybrid Powell method [12]. The subroutine c0sPBr requires the user to provide the exact Jacobian of the system. This is not a major difficulty in the case of collocation, but can lead to some tedious programming in the case of the Galerkin formulation due to the often complicated nature of (2.16) and (2.18). An alternative approach is to use the N A G subroutine C05NBFwhich is the same as c05PBr except that the user is not required to provide the exact Jacobian, which is calculated internally. This, however, could lead to a substantial increase in the cost of the process.
A. Karageorghis, Chebyshev spectral methods
107
The numerical experiments in this work were carried out in double precision on a Honeywell 6080 computer at The University College of Wales, Aberystwyth, and on a CDC 7600 computer located at the University of Manchester. The subroutine C05pBF stops when the relative increment in the nonlinear iteration is less than a user specified tolerance, which in this case was taken to be equal to 10 -6
3. Numerical examples
3.1. Example 1 3.1.1. The problem This example describes free convective flow past a vertical plate embedded in a saturated porous medium. The governing nonlinear ordinary differential equation is
f(3)(X) + (2 +
~m)f(x)f(2)(x)-( 1 + 2m)(f(t)(x))2-0
f(0) - 0,
f(2)(0) --- - 1 ,
on
0 ~ x
(3.1)
subject to /(1)(oo)--0.
(3.2)
This problem was solved by Na and Pop [13] by replacing (3.1) by a first-order system and then applying a finite difference technique.
3.1.2. The numerical methods In this problem, infinity is replaced by a finite value L, say, and the problem solved on [0, L]. Chebyshev collocation and Chebyshev-Galerkin methods are applied as described in Section 2. In order to illustrate the Galerkin formulation for this specific problem, application of (2.9)-(2.15) leads to a system of the form (2.16)
a(3) + ½(32_+ lm)[ m~__uiu ....(2) i + n
i
n
a i u~ +
i+j=n
s
E
li-jl
=n
E
a(2)]j
aii
..(l)_(l~i u, uj i - - 0 ,
n=0,1,...,N-3,
(3.3)
ii_/l= v
with
f(O) = ~ ( - 1)'a,, = O, n=O
N-2 J~(2)(0)---- -~ E (--1)nn2( n 2 - 1 ) a n - - 1 ,
(3.4)
n=O N-1
](t)(L)= 2 E n2an --'0. n=O
In order to deal with the system (3.3) one requires C~rollary 2.2 to express u,,-(3),an(2), and a(~) in terms of the a n Direct application leads to n
108
A.
Karageorghis, Chebyshevspectral methods
4 n+2~I~N (n + 2j-1)a~+zj_l ,
- -£
a(2) = 16 Cn n a(3) _-- 32
Cn
n
.3
n+2j~N
~
jffil
j(n + ])(n + 2j)a,,+zj,
(3.5)
n+2j+l~N
]=1
j ( j + 1)(n + j)(n + j + 1)(n + 2j + 1)a,,+2j+~ •
3.1.3. Numerical results The quantity of interest in this problem is the dimensionless excess surface temperature fcl)(0) and its values were calculated for m = 0, 1, 2. An often crucial point in problems where infinity is replaced by a finite length L is the choice of L. Following [1] the initial choice was L = 20. In Table 1 one can see the convergence of the approximations to f ° ) ( 0 ) obtained using a collocation and a Galerkin formulation for the case m = 0. From this table it is evident that the Galerkin formulation converges more rapidly. Further, it was observed that although continuation (in N) was required for both formulations, Galerkin proved more robust in that continuation was started for a larger value of N. In order to confirm that L = 20 is sufficiently large, both the collocation method and the Galerkin formulation (with N = 38 and 28, respectively) were applied for the values of L ffi 22, 24, 26, 28, and 30. This had negligible effect on the value of f'(0) thus confirming the initial choice of L. Results for m = 1, 2 were also computed (using continuation in m) and are compared to the results of [13] and [1] in Table 2. The Chebyshev approximations are in extremely good agreement with the results of [1]. The only result which is comparable in [13] is for m = 0. For m - 1 and m = 2 the results obtained in [13] correspond to different problems (see [1]). The results presented in [1] were obtained using the packages DVCPR and COLSYS. The package COLSYS required 40 degrees of freedom for m = 0, 1 and 64 degrees of freedom for m = 2, Table 1 Dimensionless excess surface temperature yCl)(0) for m ffi0, L = 20.00 N 12 14 16 18 20 22 24 26 28 30 32 34 36 38
Chebyshev collocation 0.12416632(01) 0.12760410(01) 0.12966599(01) 0.13000872(01) 0.12979998(01) 0.12965069(01) 0.12961028(01) 0.12961056(01) 0.12961534(01) 0.12961736(01) 0.12961771(01) 0.12961765(01) 0.12961760(01) 0.12961758(01)
Cbebyshev-Galerkin 0.12968270(01) 0.12961553(01) 0.12961279(01) 0.12961618(01) 0.12961741(01) 0.12961761(01) 0.12961760(01) 0.12961758(01) 0.12961758(01)
A. Karageorghis, Chebyshev spectral methods
109
Table 2 Comparison of results for f(l)(0) for m ffi 0,1, and 2 m
Na and Pop [13]
Fairweather and Vedha-Nayayam [1]
Chebyshev collocation, N - 38
Chebyshev-Galerkin, N ffi 28
0 1 2
0.12955(01) ---
0.129618(01) 0.100000(01) 0.864977
0.12961758(01) 0.10000000(01) 0.86497615
0.12961758(01) 0.10000000(01) 0.86497615
whereas DvcPB required 63 degrees of freedom for convergence. The package vTvre failed for this problem.
3.2. Example 2 3.2.1. The problem The second example considered is
u(2)(r) + s_u(,(r) + # e.(, ) = o. r
s=1,2,
0
(3.6)
subject to
u(')(o) = o ,
u(1) = o .
(3.7)
The two-point boundary value problem (3.6), (3.7) describes the steady-state temperature distribution in the radial direction in the interior of a cylinder of unit radius (when s = 1) and a sphere of unit radius (when s = 2). This problem was examined by Na and Tang [14] and solved numerically as an initial value problem with a Runge-Kutta technique. The same problem was also examined in [1] for a variety of values of the parameter/3 using each of the packages mentioned in Section 1.
3.2.2. The numerical method The presence of the exponential term in (3.6) makes the use of a Galerkin formulation for this problem extremely complicated and therefore impractical (see [15]). In applying a Chebyshev collocation method the following useful modification is applied, often used in spectral methods (see e.g. [16]). Rather than considering an approximation of the form (2.4), one considers N
fi(r) = ~ anP*(r) ,
(3.8)
n=2
where the polynomials P* (r), n - 2 , . . . , N, are modified Chebyshev polynomials of the form e*(r) = r * ( 0 + ~.r~(r) + ~ . r ; ( 0 , with T*(r) = TL(r) for L = '..
(3.9)
A. Karageorghis, Chebyshev spectral methods
110
The coefficients a. and/3, are chosen so that the polynomials P*(r) satisfy the boundary conditions (3.7); that is (3.10)
P*<'(o) = P * O ) = o ,
which leads to an = ( - 1 ) " n
2,
fin = - 1 - 1 ( - 1 ) " n 2 .
(3.11)
The modified Chebyshev approximation (3.8) has the advantage of satisfying the boundary conditions, which therefore need not be added to the system, which reduces the number of unknowns by two (m ! + m 2 in general).
3.2.3. Numerical results Case s = 1. In this case there is no solution for values of/3 > 1.7 (see [14]). As in [1], numerical experiments were performed for several values of fl in (0, 1.7]. The quantity of interest u(0) (the dimensionless centre temperature), is tabulated in Table 3, from which it can be seen that relatively very few degrees of freedom produce results which are in excellent agreement with the results of [1]. In fact, a maximum of 11 degrees of freedom was required in Chebyshev collocation whereas in [1] coLsYs required a maximum of 60 degrees of freedom.
Case s = 2. In this case no solution exists for/3 > 3.0 ([14]) and the values of u(0) for various degrees of freedom and values of/3 are tabulated in Table 4. Although convergence was slower for values of/3 approaching 3.0, very few degrees of freedom were required and the results are almost indistinguishable from the results reported in [1]. The maximum number of degrees of freedom required was 11 (for/3 ffi 3.0) whereas COLSYSrequired 60 degrees of freedom for/3 = 2.8. In neither case was continuation required, that is, for any value of/3 and N, the nonlinear solver, with zero starting values for all the coefficients a,, converged without difficulty.
3.3. Example 3 3.3.1. The problem Takhar and Soundalgekar [17] show that the boundary-layer equations governing a two-dimensional steady flow of a micropolar fluid past a continuously moving semi-infinite porous plate can be modelled by the two-point boundary value problem f(3)(x) + f(x)f{2)(x) + kg(1)(x) = O,
(3.12)
GgC2)(x) - 2(2g(x) + fc2)(x)) = 0,
(3.13)
v¢2)(x) + Prf(x)v (1) (x) + PrEc(f(2)(x))2 subject
to the boundary conditions
---- 0
on 0 ~
(3.14)
0.25482216(-1) 0.25482214(-1) 0.25482214(-1)
0.254822(-1)
0.252(-1)
5 6 7 8 9 10 11 12
[1]
[14]
0.519(-1)
0.519865(-1)
0.51986548(-1) 0.51986536(-1) 0.51986536(-1)
0.2
0.796(-1)
0.796091(-1)
0.79609185(-1) 0.79609144(-1) 0.79609144(-1)
0.3
0.16864486(--1) 0.16864485(--1) 0.16864485(--1)
0.168645(--1)
5 6 7 8 9 10 11 12
[1]
[14] 0.168(--1)
0.1
N
0.343(--1)
0.341387(--1)
0.34138674(--1) 0.34138671(--1) 0.34138671(--1)
0.2
0.522(--1)
0.518450(--1)
0.51845019(--1) 0.51845008(--1) 0.51845008(--1)
0.3
0.706(--1)
0.700079(--1)
0.70007967(--1) 0.70007941(--1) 0.70007941(--1)
0.4
0.5
0.5
B
0.138
0.894(--1) .
0.886541(--1)
0.6
0.107
0.107813
0.6
0.7
0.127
0.127516
0.12751624 0.12751610 0.12751610
0.169
0.170397
0.17039728 0.17039693 0.17039692 0.17039692
0.10781293 0.10781284 0.10781284
0.138673
0.13867313 0.13867293 0.13867292 0.13867292
0.88654200(--1) 0.88654149(--1) 0.88654149(--1)
0.109
0.108461
0.10846142 0.10846132 0.10846132
0.4
Table 4 Dimensionless centre temperature u(0); sphere, s - 2
0.1
N
Table 3 Dimensionless centre temperature u(0); cylinder, s = 1 0.7
0.147
0.147799
0.14779947 0.14779928 0.14779928
0.8
0.2035
0.203815
0.20381603 0.20381546 0.20381544 0.20381544
0.8
0.259448
0.25944881 0.25944814 0.25944816 0.25944816
1.3
0.238
0.239148
1.0
1.8
0.65
0.393966
2.3
1.5
1.1
0.565667
0.812248
1.7'
3.0
--
0.95840045 0.95841487 0.95842562 0.95842217 0.95842234 0.95842239 0.95842238 0.95842238
2.0
0.731577
0.73158548 0.73158144 0.73|57851 0.73157788 0.73157792 0.73157794 0.73157794 0.73157794
0.81222808 0.81224669 0.81225074 0.81224943 0.81224946 0.81224949 0.81224948 0.81224948
2.8
0.575364
0.57537024 0.57536537 0.57536426 0.57536413 0.57536414 0.57536415 0.57536414 0.57536414
0.56566360 0.56566657 0.56566718 0.56566703 0.56566703
0.316694
0.31669619 0.31669447 0.31669437 0.31669437
0.39396636 0.39396560 0.39396572 0.39396570 0.39396570
0.23914892 0.23914806 0.23914802 0.23914802
b-t
A. Karageorghis, Chebyshev spectral methods
112
f(O) =
= 1,
f<'(®) - 0 ,
o(®) =0,
v(O) = 1 ,
g(o) =o, (3.15)
8(00) =0,
where k is a coupling parameter, G the microrotation parameter, Pr and Ec the Prandtl and Eckert numbers, respectively, and f~ represents the suction or injection at the plate. Problem (3.12)-(3.15) was solved numerically by Takhar and Soundalgekar in [17], and in [18] for the case f ~ - 0 using a Runge-Kutta method, and by Fairweather and VedhaNayagam [1] using the package DVCPR(the other two packages were also used but proved costly).
3.3.2. The numerical methods The Chebyshev spectral methods of Section 2 can be readily applied to boundary value problems of the type (3.12)-(3.15). After replacing infinity by L, say, one simply determines the approximations N
f ( x ) - ~ a.r~(x),
(3.16)
nffi0 N
~(x)= ~ b.TL(x),
(3.17)
n--0 N
(3.18) n--0
In this problem the collocation and Galerkin spectral methods lead to mixed systems of order 3(N + 1) in the 3(N + 1) unknown coefficients of (3.!6)-(3.18).
3.3.3. Numerical results Following the considerations of [1], L was taken to be equal to 5. The quantity of interest in this problem is the dimensionless temperature gradient at the plate -v~l)(0). The problem with fw ffi 0 describes the flow of micropolar fluid past an impermeable, continuously moving, semi-infinite plate and was studied in [18]. For this problem the performance of the Chebyshev collocation and the Chebyshev-Galerkin formulations with respect to the number of degrees of freedom used is examined in Table 5. It is clear that the Chebyshev-Galerkin formulation converges more rapidly than the Chebyshev collocation method, but the difference is not as pronounced as in the first example. A comparison of the results obtained in [1,18] and the present study (for N = 24) is given in Table 6 and indicates extremely good agreement of our results with results of [1]. The problem with various values of fw is examined in Table 7 where the results obtained using the spectral methods are compared to those obtained in [1, 17] for different values of k and G. There is excellent agreement with the results of [1] whereas the results reported in [17] appear to be quite inaccurate. In [1], as many as 195 degrees of freedom were required for convergence using the package Dvcea, whereas as few as 72 were required using the spectral methods. Continuation in N had to be performed in both the collocation method and the Galerkin formulation. As was observed in the previous example, however, the continuation requirements of the Galerkin formulation were not as demanding as for collocation.
A. Karageorghis, Chebyshev spectral methods Table 5 Values of - v ° ) ( 0 ) obtained using the Chebyshev collocation and Chebyshev-Galerkin methods for fw -- 0, P r = 7 , Ec = 0.02 N 10 12 14 16 16 20 22 24
Collocation
Table 6 Comparison of values of -v(l)(0) obtained using the present formulations and the ones of [1, 17] (fw = 0, Pr = 7, Ec = 0.02)
Galerkin
0.190~.~A.8(01) 0.19218366(01) 0.19321768(01) 0.19316883(01) 0.19307409(01) 0.19306859(01) 0.19307468(01) 0.19307542(01)
0.19371829(01) 0.19329357(01) 0.19306532(01) 0.19306231(01) 0.19307460(01) 0.19307567(01) 0.19307514(01) 0.19307506(01)
113
k
G
[18]
[1]
0.1 0.5 0.1 0.5
2 2 4 4
1.944 1.801 1.946 --
1.93075 1.93017 1.93093 1.93129
Collocation
Galerkin
N=24
N=24
1.9307542 1.9301748 1.9309349 1.9312888
1.9307506 1.9301709 1.9309312 1.9312850
3.4. Example 4 3.4.1. The problem The buoyancy-induced flow in a variable-porosity medium adjacent to a horizontal heated plate can be modelled [19] by the two-point boundary value problem f ( 2 ) ( ' r l ) - - - ( 1 + de-'a)[Av(@)+ -~(A-2)7/v°)('¢/)]- de-'~(1 + de-'~)-lf°)('q),
(3.19)
D(2)(~) ~-[80(1 "l" d * e - ' ~ ) ( l - b)+ b]-~[Av('rl)f(1)(71)- ~(A + l)f(@)v(~)O/) + eod*e-'~(1 - b)v(~)(v/)] on 0~
(3.20)
subject to the boundary conditions
f(O) =0,
v(O)=1,
f(~)(oo)=O,
v(oo)=O,
(3.21)
wl~ere b is the ratio of the thermal conductivity of the solid to that of the fluid, A describes the temperature variation of the heated ~!ate, and 80, d, and d* are constants chosen to be 0.4, 3, and 1.5 respectively (see [19]).
3.4.2. The numerical method The presence of the exponential terms in (3.19) and (3.20) makes the use of a Galerkin formulation for this problem impractical. Instead, a collocation formulation is used after replacing infinity by a finite value L.
3.4.3. Numerical results One of the quantities of interest in this problem is - v ° ) ( 0 ) , and is examined in [1, 19]. In [1, 19] the problem is solved for b = 2 and b = 100 for A = 0.5, 1.0, 1.5, and 2.0. For b = 2.0, A = 0.5, the performance of the Chebyshev collocation methods is examined in Table 8 for L = 18.0. As few as 30 terms in the expansions for the approximations produced very satisfactory results. This is confirmed in Table 9 where the results of the present study are in
(a)
0.5742 0.4544 0.1946 0.3767 0.1344
f~
0.7 0.5 0.0 0.5 0.7
¢c)
5.75418 5.7541845 4.55283 4.5528343 1.93093 1.9309349 0.375206 0.3752085f __0"133095 0.13309052
(b)
k =0.1, G = 4
5.7541771 4.5528311 1.9309312 0.37520622 0.13308492
(d) 0.5733 0.4538 -0.3776 0.1353
(a) 5.75311 4.55193 1.93067 0.375546 0.133411
(b) 5.7531148 4.5519334 1.9306748 0.37554813 0.13341693
(c)
k = 0.02, G = 2 (a)
(b)
5.7531072 -5.75387 4.5519301 0.4538 4.55255 1.9306711 -1.93075 0.37554596 0.3769 0.375167 0.13341091 0.1348 0.133084
(d)
5.7538732 4.5525522 1.9307542 0.37516932 0.13308898
(c)
k =0.1, G = 2
5.7538657 4.5525489 1.9307506 0.37516701 0.13308331
(d)
Table 7 Values of - v (1)(0): (a) from [17], (b) from [ 1], (c) Chebyshev collocation, and (d) Chebyshev-Galerkin method, N = 24, for Pr = 7 and Ec = 0.02
l
4~
A. Karageorghis, Chebyshev spectral methods Table 8 Convergence of v(l)(0) for L 18.0, b = 2.0, A = 0.5
N
v(')(O)
14 16 18
-0.13936116(01) -0.13988288(01) -0.14007927(01)
20
-0.14014816(01)
22 24
-0.14017115(01) -0.14017847(01)
115
Table 9 Values of -v(l)(0): (a) from [19], (b) from [1], and (c) computed in present study
b=2
26
-0.1401s068(01)
28 29
-0.14018131(01) -0.14018140(01)
A
(a)
b=lO0
(b)
(c)
(a)
(b)
(c)
0.5
1.4013 1.40183 1.4018140-2.959 8.03640 7.6667223
1.0 1.5 2.0
1.890 2.307 2.682
1.89080 1.8907711-4.89 2.30743 2.3072697-6.726 2.68250 2.6820309-8.406
9.77060 9.3153150 11.0519 10.531564 12.0883 11.513973
excellent agreement with the results presented in [1, 19]. Graphs of f against T/for b = 2, A =0.5,1.0, 1.5, and 2.0 obtained in this study are in very good agreement with the corresponding graphs produced in [1] (Fig. 1). For b = 100, L must be taken to be as large as 60 and the problem becomes more difficult to solve. Because of the lack of structure in the spectral method Jacobian a large number of degrees of freedom requires a large amount of storage. As a result we could only reach N = 114, that is 115 degrees of freedom in each Chebyshev polynomial expansion. Consequently the results obtained are not in good agreement with the ones obtained in [1] using the packages DVCPRand co~Ys, each of which, however, required over 200 degrees of freedom to determine each unknown function. Plots of
3.0
~:O.S 2-S
X:I 2.0
},=2
f 1.5
1.0
O.S
0
0
/ I
I
I
I
I
I
I
2
3
4
5
6
7
Fig. 1. Graph of f against 7/for b = 2, A = 0.5, 1.0, 1.5, and 2.0.
116
A. Karageorghis, Chebyshev spectral methods
f against ~? for b : 100, A = 0.5, 1.0, 1.5, and 2.0 are, however, in good agreement with the corresponding graphs of [1] (Fig. 2). While the results presented in [19] for the case b = 2 are in good agreement with those presented in [1] and in the present study, the results for b - 100 are rather questionable. 3.5. Example 5 3.5.1. The problem Heat conduction with temperature-dependent thermal properties may be modelled by the two-point boundary value problem [20] (1 + Hv(~'l))v(2)(vl)+ H(v(1)(~?)) 2 +
2~v(1)(vt) = 0 on 0<~vt
(3.22)
subject to the boundary conditions v(0} = 1 and
v(ao)= 0.
(3.23)
This problem was solved by Aziz and Na [21] with a Runge-Kutta technique and in [1] with the package COLSYS. 3.5.2. The numerical method For this problem we use a different technique. Instead of replacing infinity by a finite value L, we map the semi-infinite interval [0, o0) onto the interval [0, 1]. This is done in the spirit of the work of Grosch and Orszag [22] and the ideas of Boyd [23] (see also [15]). 25
2O X:I
f
~L:1.5 - -
15
I0
5
O0
i
l
I
I
i
I
I
2
6
6
O
I0
12
16
16
Fig. 2. Graph of f against ~ for b - 100, A - 0.5, 1.0, 1.5, and 2.0.
117
A. Karageorghis, Chebyshev spectral methods
The algebraic-type mapping z =
'/
L>O
,/+L'
(3.24)
maps the semi-infinite interval onto [0, 1]. Since
= ( l - z ) 2 dr(z) L
dz
1 [ v(2)(~) = ~ ( 1 - z) 3 ( 1 - z)
'
d2o(z) dz 2
dr(z) ] 2 dz
' (3.25)
by setting
dr(z) 0(z)=
dz
O(z) =
'
d2v(z) dz 2
equation (3.22) becomes (1 - z ) 4
L2
//(z) + 2(1 -
2H ~-~ ( 1 -
z)[z
z)3v(z)O(z) +
(1 - z) 2 H '~'~ ]t~(z) + ~ (1 H(1 - z) 4
L2
z)4v(z)iJ(z) (3.26)
(O(z)y =o,
and the boundary conditions (3.23) become v ( 0 ) - 1,
v(1)-0.
(3.27)
The transformed problem (3.26), (3.27) was subsequently solved using Chebyshev collocation on the interval [0, 1] with a Chebyshev approximation of the form N
~(z) = 1 -
z + ~, a.P*(z),
(3.28)
n----2
where the modified Chebyshev polynomials were chosen as in (3.9).
3.5.3. Numerical results The degree of difficulty of this problem increased considerably as the value of H approached - 1 . For other values of H the results obtained were independent of the value of the mapping parameter L. As H approached - 1 , however, a reduction in L proved necessary as well as an increase in the number of degrees of freedom. This is due to the fact that as H approaches - 1 , (3.22), (3.23) become a singular perturbation problem and this results in a thickening of the boundary layer at the origin. The effect of reducing L is to concentrate more collocation points near the origin ~hich is necessary as H approaches - 1 . Table 10 shows the convergence in the values of the wall temperature gradient vtl)(0) as N increases and L decreases, for H - -0.96. Table 11 shows a comparison of the results obtained in the present
A. Karageorghis, Chebyshev spectral methods
118
Table 10 Wall temperature gradient for H--0.96 N 44 44 44 44 49 59 61 63 65 67
L
v~l~(O)
0.30 0.17211038(02) 0.29 0.17214930(02) 0.28 0.17218631(02) 0 . 2 4 5 0.17230011(02) 0.2 0.172~6355(02) 0.16 0.17250591(02) 0.16 0.17250646(02) 0.16 0.17250686(02) 0.16 0.17250727(02) 0.16 0.17250765(02)
Table 11 Wall temperature gradient -v°~(0) for values of H approaching -1, (a) computed in the present study, (b) presented in [1]
x
N
L
(a)
(b)
-0.95 -0.96 -0.97 -0.98 -0.99 -0.995
59 67 73 99 122 122
0.16 0.16 0.155 0.15 0.11 0.10
0.13926499(02) 0.17250765(02) 0.22788320(02) 0.33859194(02) 0.67005346(02) 0.12988865(03)
0.139265(02) 0.172508(02) 0.227890(02) 0.338615(02) 0.670711(02) 0.133483(03)
study and the results of [1] computed with COLSVS.The sharp increase in N and decrease in L show the degree of difficulty of the problem as H tends to -1. Continuation in both N and L and to be employed and continuation in L became more difficult as L decreased, i.e. smaller and smaller steps in the value of L had to be taken as L was decreased. The use of an algebraic mapping for such problems has the advantage that via L one can regulate the concentration of collocation points near the origin. This results in considerable savings in the number of degrees of freedom required (see [1]).
4. Comments and conclusions
The numerical experiments performed in this study indicate that the Chebyshev collocation formulation offers considerable advantages over the Chebyshev-Galerkin formulation. From the implementational aspect, the Chebyshev collocation approach is extremely simple in contrast to the Galerkin approach which can be very tedious, especially if one uses a nonlinear equations solver which requires the provision of the 3acobian. Failure to use such a solver, in favour of a solver not requiring the user to provide the Jacobian, such as C05NBF,results in loss of efficiency. On the other hand, it was observed that the Galerkin formulation led to systems of equations, the numerical solution of which required less effort when continuation was needed. This advantage was, however, rarely exploited, because the Galerkin formulation (even when the exact 3acobian was provided) proved to be computationally expensive. As is to be expected, even in cases when continuation was not required it was nonetheless employed in order to reduce the computational cost of the process. In terms of accuracy it was noted that there is little to choose between the two formulations although, in general, the Galerkin approach produced more accurate results with relatively fewer degrees of freedom than the collocation approach. The Chebyshev collocation approach proved to be a more general numerical tool than the Galerkin approach as it could be applied with equal ease to all the problems examined herein, regardless of the complexity of the differential equations, whereas application of the Galerkin formulation to problems involving complicated differential equations proved very difficult. From these observations one feels that the advantages of the
A. Karageorghis, Chebyshev spectral methods
119
Chebyshev collocation formulation make it much more attractive than the Galerkin formulation. Most of these points are confirmed in the conclusions reported in [4, 5]. In all examples considered in this study, it was observed that spectral Chebyshev methods required less degrees of freedom than the ones required by the packages of [1] to obtain comparable accuracy. On the negative side, it was noticed that large amounts of storage were required by the spectral Chebyshev methods when relatively many degrees of freedom were used. This is a direct result of the fact that in spectral methods the Jacobians are dense. A way of alleviating this problem would be to use a spectral element method, i.e. two or more elements to describe the region rather than one. Such a process would produce block-diagonal Jacobians, the structure of which could be exploited to reduce the storage requirements. Appendix A
In order to prove the Corollary 2.2 one requires the following elementary results, under the assumptions of Section 2.
RESULT I L{
c,
TL'(x)= 4
...LtI~, ,
(rT1) l,+tq, x)
dr-2 L(I)} (r---i) T,_, ,
(A.1)
where
Co=2,
c,=l
for r > O ,
c,=O
for r < O ,
and
d,=l
d,=O for r < O .
for r ~ O ,
Result I can be easily proved by considering the relevant formula on [0, 1], i.e. L - 1 [9] and then mapping from [0, 1] to [0, L].
RESULT H Assuming that
(n~_0
oo
~.(q-l)
L
\(1)
~---
_ ,.
n(q)
-
then 4n < q - l ) L an
--
Vn-l'~n-I
n ~'~ u-(q)T~(x) , n----O
dn
,,(q)
-l'*n+l,
n ~>1.
(A.2)
(A.3)
Result II can be proved by using (A.1) in the right-hand side of (A.2) and equating the coefficients of Y~Cl~(x) (see [24]).
RESULT III (q 1) _ _(q) _ _4 ~ (n + 2i - 1)a.~+2~_ 1. UnUn L i=o
Formula (A.4) is the solution of the recurrence relation (A.3) (see [25]).
(A.4)
A. Karageorghis, Chebyshev spectral methods
120
The following lemma, which is proved in [i0], is also required.
L E M M A IV (from [10]). i
(n+2i-1)(J+q-2)l ( j - 1)l
,,/=,
(n + 2i + ] + q - 3)l (n + 2i + j - 2)t
i+j=i+l
=-
1 (/+q-l)! q ( l - 1)l
(n+l+q-1) t " Vn~O, (n + l - 1)!
(A.5)
Vq~l.
PROOF OF C O R O L L A R Y 2.2. Following the same steps as in [10], for q = 1 equation (A.4) with q = 1 yields the required formula. Proceeding by induction, assuming that (2.15) is valid for q, one needs to show that (n+j+q-1)l (n + j - 1)!
c..(q+l) (4~q+1 1 ~ , ( j + q - 1 ) l """ = L / • q'-~f"~.= ( j - 1)!
(n + 2j + q - 1)a,,+2j+q-1.
(A.6)
Application of (A.4) with q + 1 instead of q and assumption of the corollary for q, yields
c"aq+l = (4)q+l 1--q'.i=~ (n +
-1){~.=
(j+q--2)[ (j -- 1)l
(n + 2i + 2j + q
(n+2i+j+q--3)l (n + 2i + j - 2)1 -
3)an+2~+2i+q_3} .
(A.7)
Replacement of i + j - 1 by I and subsequent application of the lemma leads to (A.6).
R E M A R K V. It can be easily shown that Corollary 2•2 is also valid for any finite interval [a, b] with L ffi b - a. The proof is identical and based on a simple translation of the variable•
Acknowledgment The a,thor is grateful to the U.K. Science and Engineering Research Council for its financial support during the course of this work and to his colleagues, Dr. T.N. Phillips, Dr. B. Rao (UCW Aberystwyth) and Dr. M. Vedha-Nayagam (University of Kentucky) for their useful comments. I am particularly indebted to Professor G. Fairweather of the University of Kentucky for his constructive criticisms and helpful suggestions. References [1] G. Fairweather and M. Vedha-Nayagam, ~,n assessment of numerical software for solving two-point boundary-value problems arising in heat transfer, Numer. Heat Transfer 11 (1987) 281-293.
A. Karageorghis, Chebyshev spectral methods
121
[21 IMSL Library Edition 9, IMSL Inc., Houston, TX, 1982. V. Ascher, J. Cristiansen and R.D. Russell, Collocation software for boundary-value ODEs, ACM Trans. Math. Software 7 (1981) 209-229. [41 A. Karageorghis, T.N. Phillips and A.R. Davies, Spectral collocation methods for the primary two-point boundary-value problem in modelling viscoelastic flows, Internat. J. Numer. Meths. Engrg. 26 (1988) 805-813. [51 A.R. Davies, A. Karageorghis and T.N. Phillips, Spectral Galerkin methods for the primary two-point boundary-value problem in modelling viscoelastic flows, Internat. J. Numer. Meths. Engrg. 26 (1988) 647-662. [6] D. Gottlieb and S.A. Orszag, Numerical Analysis of Spectral Methods; Theory and Applications, CBMS-NSF Regional Conference Series in Applied Mathematics 26 (SIAM, Philadelphia, PA, 1977). [7l D. Gottlieb, M.Y. Hussaini and S.A. Orszag, Introduction: theory and applications of spectral methods, in: R.G. Voight, D. Gottlieb and M.Y. Hussaini, eds., Spectral Methods for Partial Differential Equations (SIAM, Philadelphia, PA, 1984) 1-54. [81 C.A.J. Fletcher, Computational Galerkin Methods (Springer, Berlin, 1984). [9] L. Fox and I.B. Parker, Chebyshev Polynomials in Numerical Analysis (Oxford University Press, London, 1968). [lOl A. Ka:ageorghis, A note on the Chebyshev coefficients of the general order derivative of an infinitely differentiable function, J. Comput. Appl. Math. 21 (1988) 129-132. [11] NAG Library Mark 11, NAG (UK) Ltd., Oxford. [121 M.J.D. Powell, A hybrid method for nonfinear equations, in: E Rabinowitz, ed., Numerical Methods for Nonlinear Algebraic Equations (Gordon and Breach, London, 1970) 87-114. [13] T.Y. Na and I. Pop, Free convection flow past a vertical flat plate embedded in a saturated porous medium, Internat. J. Engrg. Sci. 41 (1983) 517-526. [14l T.Y. Na and S.C. Tang, A method for the solution of conduction heat transfer with non-finear heat generation, Z. Angew. Math. Mech. 49 (1969) 45-52. [15] S.A. Orszag, Spectral methods for problems in complex geometries, in: S.V. Parter, ed., Numerical Methods for Partial Differential Equations (Academic Press, New York, 1979) 273-305. [16] A. Karageorghis and T.N. Phillips, Spectral collocation methods for Stokes flow in contraction geometries and unbounded domains, J. Comput. Phys. (to appear). [17] H.S. Takhar and V.M. Soundalgekar, Flow and heat transfer of a microplanar fluid past a continuously moving porous plate, Intemat. J. Engrg. Sci. 23 (1985) 201-205. [18] V.M. Soundalgekar and H.S. Takhar, Flow of micropolar fluid past a continuously moving plate, Intemat. J. Engrg. Sci. 21 (1983)961-965. [19] B.C. Chandrasekhara, EM.S. Namboodiri and A.R. Hanumanthappa, Similarity conditions for buoyancy induced flows in a saturated porous medium adjacent to impermeable horizontal surfaces, W~irme- und Stoffubertragung 18 (1984) 17-23. [2o1 S.H. Cho and J.E. Sunderland, Phase change problems with temperature-dependent thermal conductivity, Trans. ASME J. Heat Transfer, 96 (1974) 214-217. [211 A. Aziz and T.Y. Na, A noniterative numerical solution for step-heated semi-infinite solid with temperature dependent thermal properties, Comput. Chem. Engrg. 5 (1981) 115-117. [221 C.E. Grosch and S.A. Orszag, Numerical solutions of problems in unbounded regions: coordinate transforms, J. Comput. Phys. 25 (1977) 273-296. ~" [231 J.P. Boyd, The optimization of convergence for Chebyshev polynomial methods in an unbounded domain, J. Comput. Phys. 45 (1982) 43-79. [241 S.A. Orszag, Accurate solution of the Orr-sommerfeld stability equation, J. Fluid Mech. 50(4) (1971) 689-703. [251 T.N. Phillips, T.A. Zang and M.Y. Hussaini, Preconditioners for the spectral multigrid method, IMA J. Numer. Anal. 6 (1986) 273-293.
CHEBYSHEV SPECTRAL METHODS FOR SOLVING TWO-POINT BOUNDARY VALUE PROBLEMS ARISING IN HEAT TRANSFER Andreas KARAGEORGHIS Department of Mathematics, The University College of Wales, Aberystwyth, Dyfed SY23 3BZ, United Kingdom Received 18 January 1988 The performance of Chebyshev spectral methods for solving nonlinear two-point boundary value problems is tested on several such problems arising in heat transfer. The collocation and the Galerkin spectral formulations are examined in detail and their various computational aspects compared. The accuracy of these methods is compared with that of standard software packages on a number of problems which appear in the literature.
1. Introduction Problems which arise in heat transfer can often be reformulated as nonlinear two-point boundary value problems for ordinary differential equations. These can be solved numerically via traditional numerical methods which have been implemented in general-purpose software packages. Three such packages were discussed in detail in a paper by Fairweather and Vedha-Nayagam [1]. In it, the authors analyse the performance of the packages DVCPR[2], COLSYS [3], and D ~ [2] for solving two-point boundary value problems arising in heat transfer. The aim of this paper is to investigate the use of Chebyshev spectral methods, in both a collocation and a Galerkin framework, for solving nonlinear two-point boundary value problems arising in heat transfer. In Chebyshev spectral methods, the solution of a differential equation is represented by a truncated series of Chebyshev polynomials. These polynomials are the eigenfunctions of a singular Sturm-Liouville problem and for problems with smooth solutions, such expansions converge faster than any fmite power of 1/n as n--, oo. Chebyshev collocation and Chebyshev-Galerkin methods have been used successfully in two recent studies [4,5] for the solution of a highly nonlinear two-point boundary value problem in modelling viscoelastic flows. In the present work, Chebyshev collocation methods are compared to Chebyshev-Galerkin methods on the set of problems solved in [1] using the packages DVCPR,COLSYS,and D~re (which are based on the trapezoidal rule witb deferred corrections, spline collocation, and multiple shooting, respectively). 2. Chebyshev S l ~ r a i methods Con~ider the two-point boundary value problem
F(x, u(x), u
on the finite interval [0, L],
0045-7825188/$3.50 ~ 1988, Elsevier Science Publishers B.V. (North-Holland)
(2.1)
104
A. Karageorghis, Chebyshev spectral methods
subject to the boundary conditions (2.2)
li[u(x)][x= o = 0 ,
i= 1,2,...,m
I,
rj[u(x)lL= L = O,
j = 1, 2 , . . . , m 2 ( f m - m , ) ,
(2.3)
where F(t~, h , - - . , tin+2) is a nonlinear function of the arguments {ti}i=~.2..... m+2, and {li}i--~,2, ml and {r~},=~, ,m are linear differential operators. The al~P'roaches used in2the" pr2esent paper are formulated in the spirit of the work of Gottlieb and Orszag [6] and are based on the idea of approximating the solution of a two-point boundary value problem by a truncated series expansion of the form N
ti(x) = ~ a , T ~ ( x )
(2.4)
nffi0
The trial functions T ~ ( x ) , n = O, 1 , . . . ,
N , are the shifted Chebyshev polynomials
(2.5) (2.6)
(y) = cos(n cos- y).
In the approximation (2.4) the coefficients a, are unknown and are chosen so that (i) ti(x) satisfies the ordinary differential equation in either the collocation sense (see Section 2.1) or the Galerkin sense (see Section 2.2), and (ii) ti(x) satisfies the associated boundary conditions. 2.1. Collocation f o r m u l a t i o n
In the collocation framework, the approximation (2.4) is forced to satisfy the differential equation at a finite set of points. In accordance with the recommendations of Gottlieb et al. [7] these points are chosen to be a subset of the points in the interval [0, L], where T ~ ( x ) reaches its extrema, i.e. the points Xk, k ffi O, 1 , . . . , N , where
L(
Xk= ~
1+COS'a" 1 - - ~
,
(2.7)
k=0,1,...,N.
The reason for this choice of collocation points is that, in principle, one can make use of FFTs in constructing the approximation ([7]). The ordinary differential equation (2.1) is satisfied in the collocation sense if F(x k, ~(Xk), f~t~)(Xk), . . . , f~tm)(Xk) ) = O,
k = m~, . . . , N -
m 2.
(2.8)
The satisfaction of the two sets of boundary conditions (2.2) and (2.3) provides m ffi m~ + m 2 extra linear equations, the total number of equations therefore being N + 1, which is sufficient to determine the N + 1 unknown coefficients in (2.4). By omitting the m I points in (2.7) nearest 0 (including 0) and the m 2 points in (2.7) nearest L (including L) we wish the boundary conditions to dominate over the differential equation near the end points (see also
[41).
A. Karageorghis, Chebyshev spectral methods
105
2.2. Galerkin formulation
In the Galerkin framework (see [8]) the approximation (2.4) is required to satisfy (1) (F(x, ~ ( x ) , . . . , ~(")(x)), r~(x)) = 0,
n = 0, 1 , . . . , N - m, - m2,
(2.9)
where
1 fo~(x-~ ( 1 - -~)) x~ -''2 v ( x ) w ( x ) d x ,
(V(X), W(X)) = -~
(2.10)
and the boundary conditions (2.2), (2.3) l,[~(x)]l~:o =
o,
i
=
1, 2 , . . . ,
m, ,
and rj[~(x)]lx=,=O,
1= 1,2,
....
,m
2 .
In general, equations (2.9) give rise to complicated expressions involving integrals such as (2.11)
( d(P)(x)~(q)(x), TL, (x)) .
It is often possible, however, to avoid performing such integrations, as was demonstrated in [5]. By using the orthogonality properties of Chebyshev polynomials and the property (see
[9]):
r~(x)r~;(x) ffi ½[r~÷,(x)
+ T~_jl(x)],
(2.12)
and in expression (2.11) assuming that N-q
d(q)(x)= ~
a (q)''L" ~jtx), "
qfO,...,m,
(2.13)
j=O
the integrals involved in (2.11) can be expressed as combinations of integrals of the form
(r,~(~)T~(x), r~(,)).
(2.14)
Using property (2.12), expression (2.14) can be rewritten as a combination of terms of the form
( r~ (,), r~ (x)).
(2.15)
Rather than evaluating integrals of the form (2.15) one may use the orthogonality properties of Chebyshev polynomials which, in the case of (2.15), are equivalent to equating the coefficients of Chebyshev polynomials of order n, n = O, 1 , . . . , N - m~ - m 2 . It is therefore possible to express equations (2.9) in terms of expressions of the form (2.15) and then equate the coefficients of the Chebyshev polynomials of the same order. This process
106
A. Karageorghis, Chebyshev spectral methods
gives rise to the set of equations G,(a, a ~'), a¢'-), . . . , a (m)) = O,
i-O,...
(2.16)
,N,
where a Cq)~ - (a¢J)~ a~q), a~q), . . . , U-
=
4q ~ (j+q-2)! (q-1)!j=, (j-l)!
where
(n+j+q-2)! (n+j-1)!
(n + 2j +
q - 2)a,,+2j+q_ 2 ,
1 for n > O , c,,= 0 for n < 0 .
c0 = 2 ,
(2.17)
For the purposes of this study we extend this for a general (finite) L. COROLLARY
2.2.
1
nun
,~ (j+q-2)! ( q - 1)! ~.= ( j - 1)l
(n+j+q-2)! (n + j - 1)!
(n + 2j + q - 2)a,+2j+q_ 2 .
(2.18)
A proof of the above corollan., is given in Appendix A. Application of the corollary greatly facilitates the solution of the system (2.16) since each of the a
Both the collocation and Galerkin formulations described in Sections 2.1 and 2.2 lead to systen~s of equations which are mixed, i.e. some of the equations are linear and the rest are nonlinear, in the unknown coefficients of the approximation (2.4). In the present study such systems were solved using the NAG subroutine C0SPBF[11] which is based on a modification of the hybrid Powell method [12]. The subroutine c0sPBr requires the user to provide the exact Jacobian of the system. This is not a major difficulty in the case of collocation, but can lead to some tedious programming in the case of the Galerkin formulation due to the often complicated nature of (2.16) and (2.18). An alternative approach is to use the N A G subroutine C05NBFwhich is the same as c05PBr except that the user is not required to provide the exact Jacobian, which is calculated internally. This, however, could lead to a substantial increase in the cost of the process.
A. Karageorghis, Chebyshev spectral methods
107
The numerical experiments in this work were carried out in double precision on a Honeywell 6080 computer at The University College of Wales, Aberystwyth, and on a CDC 7600 computer located at the University of Manchester. The subroutine C05pBF stops when the relative increment in the nonlinear iteration is less than a user specified tolerance, which in this case was taken to be equal to 10 -6
3. Numerical examples
3.1. Example 1 3.1.1. The problem This example describes free convective flow past a vertical plate embedded in a saturated porous medium. The governing nonlinear ordinary differential equation is
f(3)(X) + (2 +
~m)f(x)f(2)(x)-( 1 + 2m)(f(t)(x))2-0
f(0) - 0,
f(2)(0) --- - 1 ,
on
0 ~ x
(3.1)
subject to /(1)(oo)--0.
(3.2)
This problem was solved by Na and Pop [13] by replacing (3.1) by a first-order system and then applying a finite difference technique.
3.1.2. The numerical methods In this problem, infinity is replaced by a finite value L, say, and the problem solved on [0, L]. Chebyshev collocation and Chebyshev-Galerkin methods are applied as described in Section 2. In order to illustrate the Galerkin formulation for this specific problem, application of (2.9)-(2.15) leads to a system of the form (2.16)
a(3) + ½(32_+ lm)[ m~__uiu ....(2) i + n
i
n
a i u~ +
i+j=n
s
E
li-jl
=n
E
a(2)]j
aii
..(l)_(l~i u, uj i - - 0 ,
n=0,1,...,N-3,
(3.3)
ii_/l= v
with
f(O) = ~ ( - 1)'a,, = O, n=O
N-2 J~(2)(0)---- -~ E (--1)nn2( n 2 - 1 ) a n - - 1 ,
(3.4)
n=O N-1
](t)(L)= 2 E n2an --'0. n=O
In order to deal with the system (3.3) one requires C~rollary 2.2 to express u,,-(3),an(2), and a(~) in terms of the a n Direct application leads to n
108
A.
Karageorghis, Chebyshevspectral methods
4 n+2~I~N (n + 2j-1)a~+zj_l ,
- -£
a(2) = 16 Cn n a(3) _-- 32
Cn
n
.3
n+2j~N
~
jffil
j(n + ])(n + 2j)a,,+zj,
(3.5)
n+2j+l~N
]=1
j ( j + 1)(n + j)(n + j + 1)(n + 2j + 1)a,,+2j+~ •
3.1.3. Numerical results The quantity of interest in this problem is the dimensionless excess surface temperature fcl)(0) and its values were calculated for m = 0, 1, 2. An often crucial point in problems where infinity is replaced by a finite length L is the choice of L. Following [1] the initial choice was L = 20. In Table 1 one can see the convergence of the approximations to f ° ) ( 0 ) obtained using a collocation and a Galerkin formulation for the case m = 0. From this table it is evident that the Galerkin formulation converges more rapidly. Further, it was observed that although continuation (in N) was required for both formulations, Galerkin proved more robust in that continuation was started for a larger value of N. In order to confirm that L = 20 is sufficiently large, both the collocation method and the Galerkin formulation (with N = 38 and 28, respectively) were applied for the values of L ffi 22, 24, 26, 28, and 30. This had negligible effect on the value of f'(0) thus confirming the initial choice of L. Results for m = 1, 2 were also computed (using continuation in m) and are compared to the results of [13] and [1] in Table 2. The Chebyshev approximations are in extremely good agreement with the results of [1]. The only result which is comparable in [13] is for m = 0. For m - 1 and m = 2 the results obtained in [13] correspond to different problems (see [1]). The results presented in [1] were obtained using the packages DVCPR and COLSYS. The package COLSYS required 40 degrees of freedom for m = 0, 1 and 64 degrees of freedom for m = 2, Table 1 Dimensionless excess surface temperature yCl)(0) for m ffi0, L = 20.00 N 12 14 16 18 20 22 24 26 28 30 32 34 36 38
Chebyshev collocation 0.12416632(01) 0.12760410(01) 0.12966599(01) 0.13000872(01) 0.12979998(01) 0.12965069(01) 0.12961028(01) 0.12961056(01) 0.12961534(01) 0.12961736(01) 0.12961771(01) 0.12961765(01) 0.12961760(01) 0.12961758(01)
Cbebyshev-Galerkin 0.12968270(01) 0.12961553(01) 0.12961279(01) 0.12961618(01) 0.12961741(01) 0.12961761(01) 0.12961760(01) 0.12961758(01) 0.12961758(01)
A. Karageorghis, Chebyshev spectral methods
109
Table 2 Comparison of results for f(l)(0) for m ffi 0,1, and 2 m
Na and Pop [13]
Fairweather and Vedha-Nayayam [1]
Chebyshev collocation, N - 38
Chebyshev-Galerkin, N ffi 28
0 1 2
0.12955(01) ---
0.129618(01) 0.100000(01) 0.864977
0.12961758(01) 0.10000000(01) 0.86497615
0.12961758(01) 0.10000000(01) 0.86497615
whereas DvcPB required 63 degrees of freedom for convergence. The package vTvre failed for this problem.
3.2. Example 2 3.2.1. The problem The second example considered is
u(2)(r) + s_u(,(r) + # e.(, ) = o. r
s=1,2,
0
(3.6)
subject to
u(')(o) = o ,
u(1) = o .
(3.7)
The two-point boundary value problem (3.6), (3.7) describes the steady-state temperature distribution in the radial direction in the interior of a cylinder of unit radius (when s = 1) and a sphere of unit radius (when s = 2). This problem was examined by Na and Tang [14] and solved numerically as an initial value problem with a Runge-Kutta technique. The same problem was also examined in [1] for a variety of values of the parameter/3 using each of the packages mentioned in Section 1.
3.2.2. The numerical method The presence of the exponential term in (3.6) makes the use of a Galerkin formulation for this problem extremely complicated and therefore impractical (see [15]). In applying a Chebyshev collocation method the following useful modification is applied, often used in spectral methods (see e.g. [16]). Rather than considering an approximation of the form (2.4), one considers N
fi(r) = ~ anP*(r) ,
(3.8)
n=2
where the polynomials P* (r), n - 2 , . . . , N, are modified Chebyshev polynomials of the form e*(r) = r * ( 0 + ~.r~(r) + ~ . r ; ( 0 , with T*(r) = TL(r) for L = '..
(3.9)
A. Karageorghis, Chebyshev spectral methods
110
The coefficients a. and/3, are chosen so that the polynomials P*(r) satisfy the boundary conditions (3.7); that is (3.10)
P*<'(o) = P * O ) = o ,
which leads to an = ( - 1 ) " n
2,
fin = - 1 - 1 ( - 1 ) " n 2 .
(3.11)
The modified Chebyshev approximation (3.8) has the advantage of satisfying the boundary conditions, which therefore need not be added to the system, which reduces the number of unknowns by two (m ! + m 2 in general).
3.2.3. Numerical results Case s = 1. In this case there is no solution for values of/3 > 1.7 (see [14]). As in [1], numerical experiments were performed for several values of fl in (0, 1.7]. The quantity of interest u(0) (the dimensionless centre temperature), is tabulated in Table 3, from which it can be seen that relatively very few degrees of freedom produce results which are in excellent agreement with the results of [1]. In fact, a maximum of 11 degrees of freedom was required in Chebyshev collocation whereas in [1] coLsYs required a maximum of 60 degrees of freedom.
Case s = 2. In this case no solution exists for/3 > 3.0 ([14]) and the values of u(0) for various degrees of freedom and values of/3 are tabulated in Table 4. Although convergence was slower for values of/3 approaching 3.0, very few degrees of freedom were required and the results are almost indistinguishable from the results reported in [1]. The maximum number of degrees of freedom required was 11 (for/3 ffi 3.0) whereas COLSYSrequired 60 degrees of freedom for/3 = 2.8. In neither case was continuation required, that is, for any value of/3 and N, the nonlinear solver, with zero starting values for all the coefficients a,, converged without difficulty.
3.3. Example 3 3.3.1. The problem Takhar and Soundalgekar [17] show that the boundary-layer equations governing a two-dimensional steady flow of a micropolar fluid past a continuously moving semi-infinite porous plate can be modelled by the two-point boundary value problem f(3)(x) + f(x)f{2)(x) + kg(1)(x) = O,
(3.12)
GgC2)(x) - 2(2g(x) + fc2)(x)) = 0,
(3.13)
v¢2)(x) + Prf(x)v (1) (x) + PrEc(f(2)(x))2 subject
to the boundary conditions
---- 0
on 0 ~
(3.14)
0.25482216(-1) 0.25482214(-1) 0.25482214(-1)
0.254822(-1)
0.252(-1)
5 6 7 8 9 10 11 12
[1]
[14]
0.519(-1)
0.519865(-1)
0.51986548(-1) 0.51986536(-1) 0.51986536(-1)
0.2
0.796(-1)
0.796091(-1)
0.79609185(-1) 0.79609144(-1) 0.79609144(-1)
0.3
0.16864486(--1) 0.16864485(--1) 0.16864485(--1)
0.168645(--1)
5 6 7 8 9 10 11 12
[1]
[14] 0.168(--1)
0.1
N
0.343(--1)
0.341387(--1)
0.34138674(--1) 0.34138671(--1) 0.34138671(--1)
0.2
0.522(--1)
0.518450(--1)
0.51845019(--1) 0.51845008(--1) 0.51845008(--1)
0.3
0.706(--1)
0.700079(--1)
0.70007967(--1) 0.70007941(--1) 0.70007941(--1)
0.4
0.5
0.5
B
0.138
0.894(--1) .
0.886541(--1)
0.6
0.107
0.107813
0.6
0.7
0.127
0.127516
0.12751624 0.12751610 0.12751610
0.169
0.170397
0.17039728 0.17039693 0.17039692 0.17039692
0.10781293 0.10781284 0.10781284
0.138673
0.13867313 0.13867293 0.13867292 0.13867292
0.88654200(--1) 0.88654149(--1) 0.88654149(--1)
0.109
0.108461
0.10846142 0.10846132 0.10846132
0.4
Table 4 Dimensionless centre temperature u(0); sphere, s - 2
0.1
N
Table 3 Dimensionless centre temperature u(0); cylinder, s = 1 0.7
0.147
0.147799
0.14779947 0.14779928 0.14779928
0.8
0.2035
0.203815
0.20381603 0.20381546 0.20381544 0.20381544
0.8
0.259448
0.25944881 0.25944814 0.25944816 0.25944816
1.3
0.238
0.239148
1.0
1.8
0.65
0.393966
2.3
1.5
1.1
0.565667
0.812248
1.7'
3.0
--
0.95840045 0.95841487 0.95842562 0.95842217 0.95842234 0.95842239 0.95842238 0.95842238
2.0
0.731577
0.73158548 0.73158144 0.73|57851 0.73157788 0.73157792 0.73157794 0.73157794 0.73157794
0.81222808 0.81224669 0.81225074 0.81224943 0.81224946 0.81224949 0.81224948 0.81224948
2.8
0.575364
0.57537024 0.57536537 0.57536426 0.57536413 0.57536414 0.57536415 0.57536414 0.57536414
0.56566360 0.56566657 0.56566718 0.56566703 0.56566703
0.316694
0.31669619 0.31669447 0.31669437 0.31669437
0.39396636 0.39396560 0.39396572 0.39396570 0.39396570
0.23914892 0.23914806 0.23914802 0.23914802
b-t
A. Karageorghis, Chebyshev spectral methods
112
f(O) =
= 1,
f<'(®) - 0 ,
o(®) =0,
v(O) = 1 ,
g(o) =o, (3.15)
8(00) =0,
where k is a coupling parameter, G the microrotation parameter, Pr and Ec the Prandtl and Eckert numbers, respectively, and f~ represents the suction or injection at the plate. Problem (3.12)-(3.15) was solved numerically by Takhar and Soundalgekar in [17], and in [18] for the case f ~ - 0 using a Runge-Kutta method, and by Fairweather and VedhaNayagam [1] using the package DVCPR(the other two packages were also used but proved costly).
3.3.2. The numerical methods The Chebyshev spectral methods of Section 2 can be readily applied to boundary value problems of the type (3.12)-(3.15). After replacing infinity by L, say, one simply determines the approximations N
f ( x ) - ~ a.r~(x),
(3.16)
nffi0 N
~(x)= ~ b.TL(x),
(3.17)
n--0 N
(3.18) n--0
In this problem the collocation and Galerkin spectral methods lead to mixed systems of order 3(N + 1) in the 3(N + 1) unknown coefficients of (3.!6)-(3.18).
3.3.3. Numerical results Following the considerations of [1], L was taken to be equal to 5. The quantity of interest in this problem is the dimensionless temperature gradient at the plate -v~l)(0). The problem with fw ffi 0 describes the flow of micropolar fluid past an impermeable, continuously moving, semi-infinite plate and was studied in [18]. For this problem the performance of the Chebyshev collocation and the Chebyshev-Galerkin formulations with respect to the number of degrees of freedom used is examined in Table 5. It is clear that the Chebyshev-Galerkin formulation converges more rapidly than the Chebyshev collocation method, but the difference is not as pronounced as in the first example. A comparison of the results obtained in [1,18] and the present study (for N = 24) is given in Table 6 and indicates extremely good agreement of our results with results of [1]. The problem with various values of fw is examined in Table 7 where the results obtained using the spectral methods are compared to those obtained in [1, 17] for different values of k and G. There is excellent agreement with the results of [1] whereas the results reported in [17] appear to be quite inaccurate. In [1], as many as 195 degrees of freedom were required for convergence using the package Dvcea, whereas as few as 72 were required using the spectral methods. Continuation in N had to be performed in both the collocation method and the Galerkin formulation. As was observed in the previous example, however, the continuation requirements of the Galerkin formulation were not as demanding as for collocation.
A. Karageorghis, Chebyshev spectral methods Table 5 Values of - v ° ) ( 0 ) obtained using the Chebyshev collocation and Chebyshev-Galerkin methods for fw -- 0, P r = 7 , Ec = 0.02 N 10 12 14 16 16 20 22 24
Collocation
Table 6 Comparison of values of -v(l)(0) obtained using the present formulations and the ones of [1, 17] (fw = 0, Pr = 7, Ec = 0.02)
Galerkin
0.190~.~A.8(01) 0.19218366(01) 0.19321768(01) 0.19316883(01) 0.19307409(01) 0.19306859(01) 0.19307468(01) 0.19307542(01)
0.19371829(01) 0.19329357(01) 0.19306532(01) 0.19306231(01) 0.19307460(01) 0.19307567(01) 0.19307514(01) 0.19307506(01)
113
k
G
[18]
[1]
0.1 0.5 0.1 0.5
2 2 4 4
1.944 1.801 1.946 --
1.93075 1.93017 1.93093 1.93129
Collocation
Galerkin
N=24
N=24
1.9307542 1.9301748 1.9309349 1.9312888
1.9307506 1.9301709 1.9309312 1.9312850
3.4. Example 4 3.4.1. The problem The buoyancy-induced flow in a variable-porosity medium adjacent to a horizontal heated plate can be modelled [19] by the two-point boundary value problem f ( 2 ) ( ' r l ) - - - ( 1 + de-'a)[Av(@)+ -~(A-2)7/v°)('¢/)]- de-'~(1 + de-'~)-lf°)('q),
(3.19)
D(2)(~) ~-[80(1 "l" d * e - ' ~ ) ( l - b)+ b]-~[Av('rl)f(1)(71)- ~(A + l)f(@)v(~)O/) + eod*e-'~(1 - b)v(~)(v/)] on 0~
(3.20)
subject to the boundary conditions
f(O) =0,
v(O)=1,
f(~)(oo)=O,
v(oo)=O,
(3.21)
wl~ere b is the ratio of the thermal conductivity of the solid to that of the fluid, A describes the temperature variation of the heated ~!ate, and 80, d, and d* are constants chosen to be 0.4, 3, and 1.5 respectively (see [19]).
3.4.2. The numerical method The presence of the exponential terms in (3.19) and (3.20) makes the use of a Galerkin formulation for this problem impractical. Instead, a collocation formulation is used after replacing infinity by a finite value L.
3.4.3. Numerical results One of the quantities of interest in this problem is - v ° ) ( 0 ) , and is examined in [1, 19]. In [1, 19] the problem is solved for b = 2 and b = 100 for A = 0.5, 1.0, 1.5, and 2.0. For b = 2.0, A = 0.5, the performance of the Chebyshev collocation methods is examined in Table 8 for L = 18.0. As few as 30 terms in the expansions for the approximations produced very satisfactory results. This is confirmed in Table 9 where the results of the present study are in
(a)
0.5742 0.4544 0.1946 0.3767 0.1344
f~
0.7 0.5 0.0 0.5 0.7
¢c)
5.75418 5.7541845 4.55283 4.5528343 1.93093 1.9309349 0.375206 0.3752085f __0"133095 0.13309052
(b)
k =0.1, G = 4
5.7541771 4.5528311 1.9309312 0.37520622 0.13308492
(d) 0.5733 0.4538 -0.3776 0.1353
(a) 5.75311 4.55193 1.93067 0.375546 0.133411
(b) 5.7531148 4.5519334 1.9306748 0.37554813 0.13341693
(c)
k = 0.02, G = 2 (a)
(b)
5.7531072 -5.75387 4.5519301 0.4538 4.55255 1.9306711 -1.93075 0.37554596 0.3769 0.375167 0.13341091 0.1348 0.133084
(d)
5.7538732 4.5525522 1.9307542 0.37516932 0.13308898
(c)
k =0.1, G = 2
5.7538657 4.5525489 1.9307506 0.37516701 0.13308331
(d)
Table 7 Values of - v (1)(0): (a) from [17], (b) from [ 1], (c) Chebyshev collocation, and (d) Chebyshev-Galerkin method, N = 24, for Pr = 7 and Ec = 0.02
l
4~
A. Karageorghis, Chebyshev spectral methods Table 8 Convergence of v(l)(0) for L 18.0, b = 2.0, A = 0.5
N
v(')(O)
14 16 18
-0.13936116(01) -0.13988288(01) -0.14007927(01)
20
-0.14014816(01)
22 24
-0.14017115(01) -0.14017847(01)
115
Table 9 Values of -v(l)(0): (a) from [19], (b) from [1], and (c) computed in present study
b=2
26
-0.1401s068(01)
28 29
-0.14018131(01) -0.14018140(01)
A
(a)
b=lO0
(b)
(c)
(a)
(b)
(c)
0.5
1.4013 1.40183 1.4018140-2.959 8.03640 7.6667223
1.0 1.5 2.0
1.890 2.307 2.682
1.89080 1.8907711-4.89 2.30743 2.3072697-6.726 2.68250 2.6820309-8.406
9.77060 9.3153150 11.0519 10.531564 12.0883 11.513973
excellent agreement with the results presented in [1, 19]. Graphs of f against T/for b = 2, A =0.5,1.0, 1.5, and 2.0 obtained in this study are in very good agreement with the corresponding graphs produced in [1] (Fig. 1). For b = 100, L must be taken to be as large as 60 and the problem becomes more difficult to solve. Because of the lack of structure in the spectral method Jacobian a large number of degrees of freedom requires a large amount of storage. As a result we could only reach N = 114, that is 115 degrees of freedom in each Chebyshev polynomial expansion. Consequently the results obtained are not in good agreement with the ones obtained in [1] using the packages DVCPRand co~Ys, each of which, however, required over 200 degrees of freedom to determine each unknown function. Plots of
3.0
~:O.S 2-S
X:I 2.0
},=2
f 1.5
1.0
O.S
0
0
/ I
I
I
I
I
I
I
2
3
4
5
6
7
Fig. 1. Graph of f against 7/for b = 2, A = 0.5, 1.0, 1.5, and 2.0.
116
A. Karageorghis, Chebyshev spectral methods
f against ~? for b : 100, A = 0.5, 1.0, 1.5, and 2.0 are, however, in good agreement with the corresponding graphs of [1] (Fig. 2). While the results presented in [19] for the case b = 2 are in good agreement with those presented in [1] and in the present study, the results for b - 100 are rather questionable. 3.5. Example 5 3.5.1. The problem Heat conduction with temperature-dependent thermal properties may be modelled by the two-point boundary value problem [20] (1 + Hv(~'l))v(2)(vl)+ H(v(1)(~?)) 2 +
2~v(1)(vt) = 0 on 0<~vt
(3.22)
subject to the boundary conditions v(0} = 1 and
v(ao)= 0.
(3.23)
This problem was solved by Aziz and Na [21] with a Runge-Kutta technique and in [1] with the package COLSYS. 3.5.2. The numerical method For this problem we use a different technique. Instead of replacing infinity by a finite value L, we map the semi-infinite interval [0, o0) onto the interval [0, 1]. This is done in the spirit of the work of Grosch and Orszag [22] and the ideas of Boyd [23] (see also [15]). 25
2O X:I
f
~L:1.5 - -
15
I0
5
O0
i
l
I
I
i
I
I
2
6
6
O
I0
12
16
16
Fig. 2. Graph of f against ~ for b - 100, A - 0.5, 1.0, 1.5, and 2.0.
117
A. Karageorghis, Chebyshev spectral methods
The algebraic-type mapping z =
'/
L>O
,/+L'
(3.24)
maps the semi-infinite interval onto [0, 1]. Since
= ( l - z ) 2 dr(z) L
dz
1 [ v(2)(~) = ~ ( 1 - z) 3 ( 1 - z)
'
d2o(z) dz 2
dr(z) ] 2 dz
' (3.25)
by setting
dr(z) 0(z)=
dz
O(z) =
'
d2v(z) dz 2
equation (3.22) becomes (1 - z ) 4
L2
//(z) + 2(1 -
2H ~-~ ( 1 -
z)[z
z)3v(z)O(z) +
(1 - z) 2 H '~'~ ]t~(z) + ~ (1 H(1 - z) 4
L2
z)4v(z)iJ(z) (3.26)
(O(z)y =o,
and the boundary conditions (3.23) become v ( 0 ) - 1,
v(1)-0.
(3.27)
The transformed problem (3.26), (3.27) was subsequently solved using Chebyshev collocation on the interval [0, 1] with a Chebyshev approximation of the form N
~(z) = 1 -
z + ~, a.P*(z),
(3.28)
n----2
where the modified Chebyshev polynomials were chosen as in (3.9).
3.5.3. Numerical results The degree of difficulty of this problem increased considerably as the value of H approached - 1 . For other values of H the results obtained were independent of the value of the mapping parameter L. As H approached - 1 , however, a reduction in L proved necessary as well as an increase in the number of degrees of freedom. This is due to the fact that as H approaches - 1 , (3.22), (3.23) become a singular perturbation problem and this results in a thickening of the boundary layer at the origin. The effect of reducing L is to concentrate more collocation points near the origin ~hich is necessary as H approaches - 1 . Table 10 shows the convergence in the values of the wall temperature gradient vtl)(0) as N increases and L decreases, for H - -0.96. Table 11 shows a comparison of the results obtained in the present
A. Karageorghis, Chebyshev spectral methods
118
Table 10 Wall temperature gradient for H--0.96 N 44 44 44 44 49 59 61 63 65 67
L
v~l~(O)
0.30 0.17211038(02) 0.29 0.17214930(02) 0.28 0.17218631(02) 0 . 2 4 5 0.17230011(02) 0.2 0.172~6355(02) 0.16 0.17250591(02) 0.16 0.17250646(02) 0.16 0.17250686(02) 0.16 0.17250727(02) 0.16 0.17250765(02)
Table 11 Wall temperature gradient -v°~(0) for values of H approaching -1, (a) computed in the present study, (b) presented in [1]
x
N
L
(a)
(b)
-0.95 -0.96 -0.97 -0.98 -0.99 -0.995
59 67 73 99 122 122
0.16 0.16 0.155 0.15 0.11 0.10
0.13926499(02) 0.17250765(02) 0.22788320(02) 0.33859194(02) 0.67005346(02) 0.12988865(03)
0.139265(02) 0.172508(02) 0.227890(02) 0.338615(02) 0.670711(02) 0.133483(03)
study and the results of [1] computed with COLSVS.The sharp increase in N and decrease in L show the degree of difficulty of the problem as H tends to -1. Continuation in both N and L and to be employed and continuation in L became more difficult as L decreased, i.e. smaller and smaller steps in the value of L had to be taken as L was decreased. The use of an algebraic mapping for such problems has the advantage that via L one can regulate the concentration of collocation points near the origin. This results in considerable savings in the number of degrees of freedom required (see [1]).
4. Comments and conclusions
The numerical experiments performed in this study indicate that the Chebyshev collocation formulation offers considerable advantages over the Chebyshev-Galerkin formulation. From the implementational aspect, the Chebyshev collocation approach is extremely simple in contrast to the Galerkin approach which can be very tedious, especially if one uses a nonlinear equations solver which requires the provision of the 3acobian. Failure to use such a solver, in favour of a solver not requiring the user to provide the Jacobian, such as C05NBF,results in loss of efficiency. On the other hand, it was observed that the Galerkin formulation led to systems of equations, the numerical solution of which required less effort when continuation was needed. This advantage was, however, rarely exploited, because the Galerkin formulation (even when the exact 3acobian was provided) proved to be computationally expensive. As is to be expected, even in cases when continuation was not required it was nonetheless employed in order to reduce the computational cost of the process. In terms of accuracy it was noted that there is little to choose between the two formulations although, in general, the Galerkin approach produced more accurate results with relatively fewer degrees of freedom than the collocation approach. The Chebyshev collocation approach proved to be a more general numerical tool than the Galerkin approach as it could be applied with equal ease to all the problems examined herein, regardless of the complexity of the differential equations, whereas application of the Galerkin formulation to problems involving complicated differential equations proved very difficult. From these observations one feels that the advantages of the
A. Karageorghis, Chebyshev spectral methods
119
Chebyshev collocation formulation make it much more attractive than the Galerkin formulation. Most of these points are confirmed in the conclusions reported in [4, 5]. In all examples considered in this study, it was observed that spectral Chebyshev methods required less degrees of freedom than the ones required by the packages of [1] to obtain comparable accuracy. On the negative side, it was noticed that large amounts of storage were required by the spectral Chebyshev methods when relatively many degrees of freedom were used. This is a direct result of the fact that in spectral methods the Jacobians are dense. A way of alleviating this problem would be to use a spectral element method, i.e. two or more elements to describe the region rather than one. Such a process would produce block-diagonal Jacobians, the structure of which could be exploited to reduce the storage requirements. Appendix A
In order to prove the Corollary 2.2 one requires the following elementary results, under the assumptions of Section 2.
RESULT I L{
c,
TL'(x)= 4
...LtI~, ,
(rT1) l,+tq, x)
dr-2 L(I)} (r---i) T,_, ,
(A.1)
where
Co=2,
c,=l
for r > O ,
c,=O
for r < O ,
and
d,=l
d,=O for r < O .
for r ~ O ,
Result I can be easily proved by considering the relevant formula on [0, 1], i.e. L - 1 [9] and then mapping from [0, 1] to [0, L].
RESULT H Assuming that
(n~_0
oo
~.(q-l)
L
\(1)
~---
_ ,.
n(q)
-
then 4n < q - l ) L an
--
Vn-l'~n-I
n ~'~ u-(q)T~(x) , n----O
dn
,,(q)
-l'*n+l,
n ~>1.
(A.2)
(A.3)
Result II can be proved by using (A.1) in the right-hand side of (A.2) and equating the coefficients of Y~Cl~(x) (see [24]).
RESULT III (q 1) _ _(q) _ _4 ~ (n + 2i - 1)a.~+2~_ 1. UnUn L i=o
Formula (A.4) is the solution of the recurrence relation (A.3) (see [25]).
(A.4)
A. Karageorghis, Chebyshev spectral methods
120
The following lemma, which is proved in [i0], is also required.
L E M M A IV (from [10]). i
(n+2i-1)(J+q-2)l ( j - 1)l
,,/=,
(n + 2i + ] + q - 3)l (n + 2i + j - 2)t
i+j=i+l
=-
1 (/+q-l)! q ( l - 1)l
(n+l+q-1) t " Vn~O, (n + l - 1)!
(A.5)
Vq~l.
PROOF OF C O R O L L A R Y 2.2. Following the same steps as in [10], for q = 1 equation (A.4) with q = 1 yields the required formula. Proceeding by induction, assuming that (2.15) is valid for q, one needs to show that (n+j+q-1)l (n + j - 1)!
c..(q+l) (4~q+1 1 ~ , ( j + q - 1 ) l """ = L / • q'-~f"~.= ( j - 1)!
(n + 2j + q - 1)a,,+2j+q-1.
(A.6)
Application of (A.4) with q + 1 instead of q and assumption of the corollary for q, yields
c"aq+l = (4)q+l 1--q'.i=~ (n +
-1){~.=
(j+q--2)[ (j -- 1)l
(n + 2i + 2j + q
(n+2i+j+q--3)l (n + 2i + j - 2)1 -
3)an+2~+2i+q_3} .
(A.7)
Replacement of i + j - 1 by I and subsequent application of the lemma leads to (A.6).
R E M A R K V. It can be easily shown that Corollary 2•2 is also valid for any finite interval [a, b] with L ffi b - a. The proof is identical and based on a simple translation of the variable•
Acknowledgment The a,thor is grateful to the U.K. Science and Engineering Research Council for its financial support during the course of this work and to his colleagues, Dr. T.N. Phillips, Dr. B. Rao (UCW Aberystwyth) and Dr. M. Vedha-Nayagam (University of Kentucky) for their useful comments. I am particularly indebted to Professor G. Fairweather of the University of Kentucky for his constructive criticisms and helpful suggestions. References [1] G. Fairweather and M. Vedha-Nayagam, ~,n assessment of numerical software for solving two-point boundary-value problems arising in heat transfer, Numer. Heat Transfer 11 (1987) 281-293.
A. Karageorghis, Chebyshev spectral methods
121
[21 IMSL Library Edition 9, IMSL Inc., Houston, TX, 1982. V. Ascher, J. Cristiansen and R.D. Russell, Collocation software for boundary-value ODEs, ACM Trans. Math. Software 7 (1981) 209-229. [41 A. Karageorghis, T.N. Phillips and A.R. Davies, Spectral collocation methods for the primary two-point boundary-value problem in modelling viscoelastic flows, Internat. J. Numer. Meths. Engrg. 26 (1988) 805-813. [51 A.R. Davies, A. Karageorghis and T.N. Phillips, Spectral Galerkin methods for the primary two-point boundary-value problem in modelling viscoelastic flows, Internat. J. Numer. Meths. Engrg. 26 (1988) 647-662. [6] D. Gottlieb and S.A. Orszag, Numerical Analysis of Spectral Methods; Theory and Applications, CBMS-NSF Regional Conference Series in Applied Mathematics 26 (SIAM, Philadelphia, PA, 1977). [7l D. Gottlieb, M.Y. Hussaini and S.A. Orszag, Introduction: theory and applications of spectral methods, in: R.G. Voight, D. Gottlieb and M.Y. Hussaini, eds., Spectral Methods for Partial Differential Equations (SIAM, Philadelphia, PA, 1984) 1-54. [81 C.A.J. Fletcher, Computational Galerkin Methods (Springer, Berlin, 1984). [9] L. Fox and I.B. Parker, Chebyshev Polynomials in Numerical Analysis (Oxford University Press, London, 1968). [lOl A. Ka:ageorghis, A note on the Chebyshev coefficients of the general order derivative of an infinitely differentiable function, J. Comput. Appl. Math. 21 (1988) 129-132. [11] NAG Library Mark 11, NAG (UK) Ltd., Oxford. [121 M.J.D. Powell, A hybrid method for nonfinear equations, in: E Rabinowitz, ed., Numerical Methods for Nonlinear Algebraic Equations (Gordon and Breach, London, 1970) 87-114. [13] T.Y. Na and I. Pop, Free convection flow past a vertical flat plate embedded in a saturated porous medium, Internat. J. Engrg. Sci. 41 (1983) 517-526. [14l T.Y. Na and S.C. Tang, A method for the solution of conduction heat transfer with non-finear heat generation, Z. Angew. Math. Mech. 49 (1969) 45-52. [15] S.A. Orszag, Spectral methods for problems in complex geometries, in: S.V. Parter, ed., Numerical Methods for Partial Differential Equations (Academic Press, New York, 1979) 273-305. [16] A. Karageorghis and T.N. Phillips, Spectral collocation methods for Stokes flow in contraction geometries and unbounded domains, J. Comput. Phys. (to appear). [17] H.S. Takhar and V.M. Soundalgekar, Flow and heat transfer of a microplanar fluid past a continuously moving porous plate, Intemat. J. Engrg. Sci. 23 (1985) 201-205. [18] V.M. Soundalgekar and H.S. Takhar, Flow of micropolar fluid past a continuously moving plate, Intemat. J. Engrg. Sci. 21 (1983)961-965. [19] B.C. Chandrasekhara, EM.S. Namboodiri and A.R. Hanumanthappa, Similarity conditions for buoyancy induced flows in a saturated porous medium adjacent to impermeable horizontal surfaces, W~irme- und Stoffubertragung 18 (1984) 17-23. [2o1 S.H. Cho and J.E. Sunderland, Phase change problems with temperature-dependent thermal conductivity, Trans. ASME J. Heat Transfer, 96 (1974) 214-217. [211 A. Aziz and T.Y. Na, A noniterative numerical solution for step-heated semi-infinite solid with temperature dependent thermal properties, Comput. Chem. Engrg. 5 (1981) 115-117. [221 C.E. Grosch and S.A. Orszag, Numerical solutions of problems in unbounded regions: coordinate transforms, J. Comput. Phys. 25 (1977) 273-296. ~" [231 J.P. Boyd, The optimization of convergence for Chebyshev polynomial methods in an unbounded domain, J. Comput. Phys. 45 (1982) 43-79. [241 S.A. Orszag, Accurate solution of the Orr-sommerfeld stability equation, J. Fluid Mech. 50(4) (1971) 689-703. [251 T.N. Phillips, T.A. Zang and M.Y. Hussaini, Preconditioners for the spectral multigrid method, IMA J. Numer. Anal. 6 (1986) 273-293.