Spectral methods were mostly developed in the

1970s. They have played an important role in re-

cent studies of numerical solutions of diﬀerential equa-

tions in regular domains. They are considered to be

eﬃcient due to their high accuracy [1]. According to

Vilhena et al. [2], spectral methods encompass repre-

senting the solution to the problem as truncated series

of known functions of independent variables. Spectral

methods that make use of collocation methods are usu-

ally called pseudo-spectral methods [3], and these meth-

ods have been widely used [4]. The collocation points

are normally the zeros of the polynomial chosen for ap-

proximation [5]. In the case of Chebyshev collocation

method, the collocation points chosen are the Gauss-

Lobatto points [6], [7], [8], [9]. The Chebyshev Gauss-

Lobatto nodes have also been used in Rashad [10]. More

recently, Legendre and Chebyshev spectral approxima-

tions have been used in PDEs in bounded domains [1].

Some progress have been made in solving problems in

unbounded domains. In particular, Tatari and Haghighi

[1] investigated Laguerre and Hermite polynomials and

used appropriate choices for semi-inﬁnite and inﬁnite do-

mains. Speciﬁcally, an accurate Christov-Galerkin spec-

tral technique for the solutions of interacting localized

wave solutions of fourth and sixth order generalized wave

equations was achieved by application of spectral meth-

ods on an inﬁnite domain [11]. The use of spectral meth-

ods has been motivated by their accuracy and robustness

in solving incompressible Navier-Stokes equations [12].

Some approximations like the Galerkin approximations

or collocation schemes have been described by [13]. The

main advantage of the method is that there is no need for

numerical integration [14]. Thus, spectral methods pro-

vide more accurate solution approximations with a small

number of unknowns, and so play important roles in op-

timizing engineering designs and other scientiﬁc compu-

tations [15].

Spectral methods are applied in numerical solutions

for neutron transport [16], Darcy and coupled Stokes

equations [17], modelling of bridges [18]. Pozrikidis

[19] applied spectral collocation method with triangular

boundary elements for solving integral equations arising

from boundary integral formulations over surfaces dis-

cretized into ﬂat or curved triangular elements. Dlamini

et al. [20] compared the pseudo-spectral method and

the compact ﬁnite diﬀerence (CFD) method for solv-

ing boundary layer problems and found that the spectral

method were better than the CFD in terms of computa-

tional speed.

Other applications of spectral methods are in

the spectral local linearization method (SLLM), spec-

tral relaxation method (SRM), the spectral quasi-

linearization method (SQLM), the spectral perturbation

method (SPM), the bivariate quasi-linearization method

(BQLM/BSQLM). Most of these methods are based on

linearization methods. The spectral relaxation method

(SRM) requires converting the equations into a system

of ﬁrst order diﬀerential equations [21] or arranging the

equations in a particular order, placing the equations

with least number of unknowns at the top of the equa-

tion list [22]. The resulting system is then decoupled

using ideas from the Gauss-Seidel method, which is nor-

Numerical analysis of fluid flow problems using spectral relaxation

method (SRM)

GILBERT MAKANDA

Central University of Technology, Free State,

Private Bag X20539, Bloemfontein, 9300

SOUTH AFRICA

Abstract: The paper presents four examples arising from mathematical models in fluid flow. Examples

1 and 2 illustrate the implementation of the spectral relaxation method (SRM) on problems involving

ordinary differential equations. Examples 3 and 4 illustrate the application of the SRM on partial

differential equations. The SRM is accurate and robust if it is used together with the successive over-

relaxation (SOR) technique. The method is easy to implement and requires less computational time

than similar methods that can be used to solve similar problems. The method converges after a few

iterations and is stable. The method can be used as an alternative method to solve problems arising in

fluid flow.

Keywords: Successive relaxation, spectral relaxation method, ordinary and partial differential

equations

Received: July 18, 2021. Revised: August 5, 2022. Accepted: September 5, 2022. Published: October 11, 2022.

1. Introduction

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

199

Volume 21, 2022

mally used to solve linear algebraic systems of equations

[9]. The decoupled system is numerically integrated us-

ing the Chebyshev pseudo spectral method [6]. Unlike

other iterative schemes for solving nonlinear diﬀerential

equations, the SRM does not require any evaluation of

derivatives and perturbation [9]. Developments of the

SRM are noted in [20] where a multistage spectral re-

laxation method is used for solving problems of chaos

control and synchronization. Thus, the SRM is an eﬃ-

cient, reliable, convergent, numerically stable and very

easy method to implement that has a great potential as

very useful tool for solving boundary layer ﬂow problems

arising from ﬂuid dynamics applications [9]. In this study

we implement the spectral relaxation method (SRM).

From the literature mentioned above, it is clear that this

method is highly accurate as well as easy to implement.

In this section we give a detailed description of the

spectral relaxation method as described in Motsa et

al. [32]. We consider a system of nnonlinear ordi-

nary diﬀerential equations in nunknown functions fi(η),

i= 1,2,...n where η∈[a, b] is the dependent variable.

We deﬁne a vector Fito be a vector of derivatives of the

variable fiwith respect to η

Fi(η) = hf(0)

i, f(1)

i, . . . , f(mi)

ii.(1)

Where f(0)

i=fi,f(p)

i, is the pth derivative of fi

with respect to ηand mi(i= 1,2,3...,n) is the highest

derivative order of the variable fiwhich is in the system

of equations. The system can be written in terms of Fi

as the sum of linear (Li) and nonlinear components (Ni)

Li[F1, F2, . . . , Fn] + Ni[F1, F2, . . . , Fn]

=Gi(η), i = 1, . . . n. (2)

Where G(η) is a known function of η.

Eq. (2) is solved subject to two point boundary con-

ditions which are expressed as

j=1

mj−1

p=0

α(p)

ν,j f(p)

j(a) = Ka,ν , ν = 1,2, . . . , na,(3)

j=1

mj−1

p=0

γ(p)

ν,j f(p)

j(b) = Kb,σ, σ = 1,2, . . . , nb,(4)

where α(p)

ν,j ,γ(p)

σ,j are the constant coeﬃcients of f(p)

in the boundary conditions, and naand nbare the to-

tal number of prescribed boundary conditions at η=a

and η=brespectively. Starting from the initial ap-

proximation F1,0,F2,0,...,Fn,0, the iterative method

is obtained as

L1[F1,r+1, F2,r , . . . , Fn,r ] = G1

+N1[F1,r, F2,r , . . . , Fn,r ],

L2[F1,r+1, F2,r+1, . . . , Fn,r] = G2

+N2[F1,r+1, F2,r , . . . , Fn,r ],

Ln−1[F1,r+1, F2,r+1, . . . , Fn−1,r+1, Fn,r] = Gn−1

+Nn−1[F1,r+1, . . . , Fn−2,r+1, Fn−1,r , Fn,r],

Ln[F1,r+1, F2,r+1, . . . , Fn−1,r+1, Fn,r+1] = Gn

+Nn[F1,r, F2,r , . . . , Fn,r ].(5)

Where Fi,r+1 and Fi,r are the approximation of Fi

at the current and the previous iterations respectively.

We state that Eqs. (5) form a system of nlinear de-

coupled equations which can be solved iteratively for

r= 1,2. . . . We start from a an appropriate initial ap-

proximation Fi,0which satisfy boundary conditions. The

iterations are repeated until convergence is reached. The

decoupling error can be used to assess the desired conver-

gence. The decoupling error Erat the (r+1)th iteration

is deﬁned by

The idea incorporated in this method is the Gauss-

Seidel relaxation method which is normally used for solv-

ing large systems of algebraic equations. To implement

the spectral collocation method, we deﬁne the diﬀeren-

tiation matrix

dfi(ηl)

dη =

k=0

Dl,kfi(τk) = D Fi, l = 0, . . . , N.(6)

Where N+ 1 is the number collocation points, D=

2D/(b−a) and

F= [f(τ0), f(τ1), . . . , f(τN)]Tis the vector function of

the collocation points and higher order derivatives are

obtained in powers of Dgiven by

f(p)

j=DpFj.(7)

We then apply the Chebyshevpseudo spectral method to

the iteration scheme shown in Eqs. (5)-(5). This then

gives

j=1

p=0

β[p]

i,j f(p)

j+Ni[F1,F2,...,Fn] = Gi,(8)

where β(p)

i,j are constants coeﬃcients of f(p)

j, the

derivative of fj(j= 1,2, . . . , n) that is in the ith equa-

tion for i= 1,2, . . . , n. The iteration scheme used in

Eqs.(5-5) can be expressed as

j=1

p=0

β[p]

i,j f(p)

j,r+1 =Gi−

j=1+1

p=0

β[p]

i,j f(p)

j,r+1

−Ni[F1,r+1, . . . , Fi−1,r+1, Fi,r , . . . , Fm,r] (9)

1.1 Spectral relaxation method (SRM)

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

200

Volume 21, 2022

for i= 1,2, . . . , m. Using the Eq. (8) on Eq. (9) and

the boundary conditions we otain the spectral Gauss-

Seidel relaxation method iteration scheme given by

j=1

p=0

β[p]

i,j D(p)Fj,r+1 =Gi−

j=1+1

p=0

β[p]

i,jD(p)Fj,r

− Ni[F1,r+1,...,Fi−1,r+1,Fi,r,...,Fn,r],(10)

subject to

j=1

mj−1

p=0

α(p)

ν,j

k=0

N,kfj,r+1(τk) = Ka,ν ,

ν= 1,2, . . . , na,(11)

j=1

mj−1

p=0

γ[p]

σ,j

k=0

N,kfj,r+1(τk) = Kb,σ ,

σ= 1,2, . . . , nb,(12)

The substitution of previously known functions de-

couples the system of equations and an eﬃcient iteration

scheme is created giving accurate results. The spectral

relaxation method (SRM) will be implemented in the

next section.

Numerical examples

In this section we demonstrate the implementation

of the spectral relaxation method, we give two examples

of applications in ordinary diﬀerential equations and two

examples of applications in partial diﬀerential equations.

We begin by considering the problem of Ece [24]; free

convection ﬂow about a vertical spinning cone under a

magnetic ﬁeld. The problem is fully described in Ece

[24]. The governing equations are written as

∂(ru)

∂x +∂(rv)

∂y = 0 (13)

u∂u

∂x +v∂u

∂y −Re2

r0w2

r=∂2u

∂y2+ Θ −MΛ2u,(14)

u∂w

∂x +v∂w

∂y −uw r0

r=∂2w

∂y2−MΛ2u, (15)

u∂Θ

∂x +v∂Θ

∂y =1

P r

∂2Θ

∂y2,(16)

where the rotational Reynolds number Re, the magnetic

ﬁeld function Λ and the magnetic parameter Mare given

Re =ΩL2

ν,Λ = b(x)

r√1−r2, M =σB2

Uρ ,(17)

where σis the electrical conductivity and ρis the den-

sity of the ﬂuid. Boundary conditions considered are as

follows;

u(x, 0) = 0, v(x, 0) = 0, w(x, 0) = r,

Θ(x, 0) = a(x),(18)

u(x, y)→0, w(x, y)→0,

Θ(x, y)→0 as y→ ∞,(19)

where

a(x)=[A(x)−T0]/(Tr−T0),(20)

c(x) = Gr−1

4C(x)/k(Tr−T0) (21)

Using the stream function and boundary layer vari-

ables deﬁned as

ru =∂ψ

∂y , rv =−∂ψ

∂x ,(22)

ψ(x, y) = xrF (y), w =rG(y),Θ = xH(y).(23)

where r=xsin φ.

we obtain the following ordinary diﬀerential equa-

tions as in Ece[24],

F000 + 2F F 00 −F02−MF 0+G2+H= 0,(24)

G00 + 2F G0−2F0G−MG = 0,(25)

P r H00 + 2F H0−F0H= 0,(26)

where = (Re sin φ)2/Gr is the spin parameter, the

boundary conditions considered in this example are given

as;

F(0) = 0, F 0(0) = 0, G(0) = 1, H(0) = 1,(27)

F0(y)→0, G(y)→0, H(y)→0 as y→ ∞.(28)

Applying the SRM as described in section A. in equa-

tions (24) - (26) and ﬁrst reducing it to a system of sec-

ond order equations we obtain,

F0=K, (29)

K00 + 2fK0−K2−M K +G2+H= 0,(30)

G00 + 2F G0−2F0G−MG = 0,(31)

P r H00 + 2F H0−F0H= 0,(32)

with boundary conditions

F(0) = 0, K(0) = 0, G(0) = 1, H(0) = 1,(33)

K(y)→0, G(y)→0, H(y)→0 as y→ ∞.(34)

The iterative scheme for the system (29) -(32) be-

comes

r+1 =Kr, Fr+1(0) = 0 (35)

K00

r+1 + 2fr+1K0

r+1 −MKr+1

=K2

r−G2

r−Hr,(36)

G00

r+1 + 2Fr+1G0

r+1 −2Kr+1Gr+1

−MGr+1 = 0,(37)

P r H00

r+1 + 2Fr+1H0

r+1 −Kr+1Hr+1 = 0,(38)

2. Example 1: Boundary layer free

convection flow from a spinning cone

under magnetic field

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

201

Volume 21, 2022

with boundary conditions

Fr+1(0) = 0, Kr+1(0) = 0, Gr+1(0) = 1,

Hr+1(0) = 1,(39)

Kr+1(y)→0, Gr+1(y)→0,

Hr+1(y)→0 as y→ ∞.(40)

Applying the Chebyshev pseudospectral method on

system (35) -(38) we obtain

A1Fr+1 =B1, Fr+1(τN) = 0,(41)

A2Kr+1 =B2, Kr+1(τN) = 0, Kr+1(τ0) = 0,(42)

A3Gr+1 =B3, Gr+1(τN)=1, Gr+1(τ0)=0,(43)

A4Hr+1 =B4, Hr+1(τN)=0, Hr+1(τ0)=0,(44)

where

A1=D,B1=Kr,(45)

A2=D2+ 2diagFr+1D−M I, (46)

B2=K2

r−G2

r−Hr,(47)

A3=D2+ 2diagFr+1D−2diagKr+1,(48)

B3= 0,(49)

A4=1

P r D2+ 2diagFr+1D−diagKr+1,(50)

B4= 0,(51)

We use the successive over-relaxation (SOR) method

to accelerated the convergence rates of the spectral re-

laxation method. This method is normally used in nu-

merical linear algebra to control the convergence of the

Gauss-Seidel method for linear systems of equations. In

this section we propose a similar method to improve the

convergence of the spectral relaxation method. If the

SRM scheme for solving a function xat the (r+ 1)th

iteration is

axr+1 =b,(52)

We deﬁne the modiﬁed version of the SRM as

axr+1 = (1 −ω)axr+ωb,(53)

where a,bare matrices and ωis the convergence con-

trolling relaxation parameter. The case ω= 1 reduce the

system (53) to (52). when ω < 1 convergence speeds up

and slows down when ω > 1. The eﬀect of changing ω

will be discussed in the results and discussion section.

We considering the problem of Awad et al. [25] which

investigated free convection ﬂow from an inverted cone in

porous medium with cross diﬀusion, the governing equa-

tions were presented as;

∂(ru)

∂x +∂(rv)

∂y = 0 (54)

u∂u

∂x +v∂u

∂y =ν∂2u

∂y2−ν

−gβ cos Ω(T−T∞) + gβ∗cos Ω(C−C∞),(55)

u∂T

∂x +v∂T

∂y =α∂2T

∂y2+Dk

CsCp

∂2C

∂y2,(56)

u∂C

∂x +v∂C

∂y =D∂2C

∂y2+Dk

CsCp

∂2T

∂y2,(57)

where r=xsin Ω, uand vare velocity components in

the xand ydirections respectively. ρis the ﬂuid density,

gis the acceleration due to gravity, Kis the permeabil-

ity, νis the kinematic viscosity of the ﬂuid. βand β∗

are the thermal expansion and concentration expansion

coeﬃcients respectively. αand Dare thermal and mass

diﬀusivities, kis the thermal diﬀusion ratio. cpis the

speciﬁc heat capacity at constant pressure and csis the

concentration susceptibility. The boundary conditions

are given by

u= 0, v = 0, T =Tw=T∞+Axλ,

C=Cw=C∞+Bxλ,on y= 0, x > 0,

u= 0, T =T∞, C =C∞,as y→ ∞ (58)

where A, B > 0 are constants and λis the power-law in-

dex. By using the similarity transformations and stream

functions described in Awad et.al [25], we obtain the fol-

lowing system of ordinary diﬀerential equations.

f000 +λ+ 7

4ff00 −λ+ 1

2f02−Λf0

+θ+N1φ= 0,(59)

θ00 +Dfφ00 +P r λ+ 7

4fθ0−P rλf 0θ= 0,(60)

φ00 +Srθ00 +Sc λ+ 7

4fφ0−Scλf0φ= 0,(61)

with boundary conditions given as

f= 0, f0= 0, θ =φ= 1,on η= 0,

f0= 0, θ = 0, φ = 0,as η→ ∞ (62)

where Dfis the Dufour number, Sr is the Soret num-

ber, N1is the buoyancy parameter, P r is the Prandtl

number, Sc is the Schmidt number. To solve the system

(59) - (62), we ﬁrst substitute for θ00 in equation (61) and

for φ00 in equation (60) reduces to the following system,

3. Improving the convergence of the

spectral relaxation method

([DPSOH)UHHFRQYHFWLRQIURP

DQLQYHUWHGFRQHLQSRURXVPHGLXP

ZLWKFURVVGLIIXVLRQHIIHFWV

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

202

Volume 21, 2022

f000 +λ+ 7

4ff00 −λ+ 1

2f02−Λf0

+θ+N1φ= 0,(63)

(1 −DfSr)θ00 −ScDfλ+ 7

4fφ0+ScDfλf0φ

+P r λ+ 7

4fθ0−P rλf 0θ= 0,(64)

(1 −DfSr)φ00 −SrP r λ+ 7

4fθ0+SrP rλf 0θ

+Sc λ+ 7

4fφ0−Scλf0φ= 0,(65)

subject to boundary conditions

f= 0, f0= 0, θ =φ= 1,on η= 0,

f0= 0, θ = 0, φ = 0,as η→ ∞ (66)

We now apply the spectral relaxation technique as

follows, the system is converted yo a system of second

order ordinary diﬀerential equations. The iterative sys-

tem as written as;

r+1 =gr, fr+1(0) = 0,(67)

g00

r+1 +λ+ 7

4fr+1g0

r+1 −Λgr+1

=λ+ 1

2g2

r−θr−N1φr,(68)

(1 −DfSr)θ00

r+1 −P r λ+ 7

4fr+1θ0−P rλgr+1θr+1

=ScDfλ+ 7

4fr+1φ0

r−λgr+1φr,(69)

(1 −DfSr)φ00

r+1 −Sc λ+ 7

4fr+1φ0−Scλgr+1φr+1

=SrP r λ+ 7

4fr+1θ0

r+1 −λgr+1θr+1,(70)

subject to

gr+1(0) = 0, gr+1(∞) = 0, θr+1(0) = 1,

φr+1(0) = 1, φr+1(∞)=0.(71)

Applying the Chebyshev pseudospectral method on

system (67) -(71) we obtain

A1fr+1 =B1, fr+1(τN)=0,(72)

A2gr+1 =B2, gr+1(τN)=0, gr+1(τ0) = 0,(73)

A3θr+1 =B3, θr+1(τN)=1, θr+1(τ0)=0,(74)

A4φr+1 =B4, φr+1(τN)=0, φr+1(τ0) = 0,(75)

where

A1=D,B1=gr, fr+1(τ¯

N)=0,(76)

A2=D2+λ+ 7

4diagfr+1D−ΛI,(77)

B2=g2

r−θ2

r−N1φr,(78)

gr+1(τ¯

N) = 0, gr+1(τ0)=0.(79)

A3= (1 −DfSr)D2+P r λ+ 7

4diagfr+1D

−P rλdiaggr+1,(80)

B3=ScDfλ+ 7

4fr+1φ0

r−λgr+1φr,(81)

θr+1(τ¯

N) = 0, θr+1(τ0)=0.(82)

A4= (1 −DfSr)D2+Sc λ+ 7

4diagfr+1D

−Scλdiaggr+1,(83)

B4=SrP r λ+ 7

4fr+1θ0

r+1 −λgr+1θr+1,(84)

φr+1(τ¯

N)=0, φr+1(τ0) = 0.(85)

A semi-inﬁnite vertical plate placed in a saturated

porous medium with uniform ambient temperature T∞,

at time t= 0 the ﬂuid is impulsively moved with velocity

Ueand the surface temperature is suddenly raised.

The governing equations in this ﬂuid ﬂow are given

∂u

∂x +∂v

∂y = 0,(86)

∂u

∂t +u∂u

∂x +v∂u

∂y =Ue

∂Ue

∂x +1

∂

∂y µ∂u

∂y 

+gβT(T−T∞) + 

ρK (µ∞U∞−µu)

+F12

(U2

e−u2) (87)

∂T

∂t +u∂T

∂x +v∂T

∂y =α∂2T

∂y2−1

∂qr

∂y .(88)

with initial conditions

u(x, y) = v(x, y)=0, T (x, y) = T∞, t < 0,(89)

The boundary conditions for t≥0

u(x, 0) = v(x, 0) = 0, u(x, ∞) = Ue=ax, a > 0,

T x, ∞=T∞, T (x, 0) = Tw(x)bxn, b > 0, n ≥0,(90)

The corresponding transformations used are

. Example 3: Non-Darcy unsteady mixed

convection flow near the stagnation

point on a heated vertical surface

embedded in porous medium with thermal

radiation and variable viscosity

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

203

Volume 21, 2022

η=aρ

µ∞

yξ−1

2, ξ = 1 −e−t∗,(91)

t∗=at, a > 0, U(x, y, t) = axf0,(92)

v(x, y, t) = −aµ∞

ρ

ξ1

2f(η, ξ),(93)

T(x, y, t) = T∞+ (Tw−T∞)θ(η, ξ),(94)

P r =µ∞

ρα , λ =Grx

Re2

,(95)

Grx=gβ(Tw−T∞)x3ρ2

µ2

∞

,(96)

Rex

ax2ρ

µ∞

,(97)

Using equation (91)-(97) in equations (86)-(88) together

with boundary conditions equations (90) we get

(1 + σθ)f000 +σθ0f00 +1

2η(1 −ξ)f00

+ξf f00 +ξ(1 −(f0)2) + (98)

+λθξ +γξ(1 −(1 + σθ)f0)+∆ξ(1 −(f0)2)

=ξ(1 −ξ)∂f0

∂ξ ,(99)

P r θ00 −Rξ(θ−1) + 1

2η(1 −ξ)θ0+ξ(fθ0−nf0θ)

=ξ(1 −ξ)∂θ

∂ξ ,(100)

where γ=µ∞/aKρ is the ﬁrst order resistance param-

eter, ∆ = F14Rexµinf ty/aKρ

2, is the second order

parameter,R= 4αI/ka is the radiation parameter, λ > 0

for buoyancy assisting ﬂow and λ < 0 for buoyancy op-

posing ﬂow. The boundary conditions reduce to

f(0, ξ) = f0(0, ξ)=0, θ(0,∞)=1,

f0(∞, ξ)=1, θ(∞, ξ) = 0.(101)

For the case ξ= 0 and σ= 0 the system admits to the

solution

f=ηerfc(η

2)−2

√π1−exp(−η2

4),

θ=erfc(P r 1

2η

2),(102)

The spectral relaxation method is applied to the sys-

tem (99) and (101) and becomes;

a1,rU00

r+1 +a2,rU0

r+1 +a3,rUr+1

=ξ(1 −ξ)∂Ur+1

∂ξ ,(103)

r+1 =Ur+1,(104)

P r θ00

r+1 +b1,rθ0

r+1 +b2,rθr+1

=ξ(1 −ξ)∂θr+1

∂ξ ,(105)

with boundary conditions;

η= 0, fr+1 =fw,

Ur+1 =1 + 1

βSfU0

r+1, θr+1 = 1 + STθ0

r+1,

η→ ∞, Ur+1 →0, θr+1 →0.(106)

where

a1,r = 1 + σθr, a2,r =σθ0

r+1

2η(1 −ξ) + ξfr,

a3,r =−γξ(1 + σθr),(107)

a4,r =ξ(γ+λθr) + (1 −U2

r)(ξ+ ∆ξ),(108)

b1,r =1

2η(1 −ξ) + ξfr,

b2,r =ξ(R−nUr), b3,r =Rξ. (109)

The initial approximations for solving the equations

(99) - (101) are obtained by setting ξ= 0.

f0= 1 −ηe(−η)−e(−η), θ0=e−η,(110)

The iterative schemes (99) - (101) can be solved it-

eratively for ur+1 when r= 0,1,2. . . , the solution for

ur+1 is used to solve for fr+1. To solve the equations

we use the Chebyshev spectral method in the η−di-

rection and use the implicit ﬁnite diﬀerence method in

the ξ−direction. The ﬁnite diﬀerence scheme is applied

with centering about the midpoint halfway between ξn+1

and ξn, which is deﬁned as ξn+1

2= (ξn+1 +ξn)/2. The

derivatives with respect to ηare deﬁned in terms of the

Chebyshev diﬀerentiation matrices. Using the centering

about ξn+1

2to any function, for example f(η, ξ) and its

corresponding derivative we get,

f(ηj, ξn+1

2) = fn+1

j=fn+1

j+fn

∂f

∂ξ n+1

=fn+1

j−fn

∆ξ.(111)

Writing the system (99) - (101) in terms of the diﬀeren-

tiation matrix we obtain,

a1,rD2+a2,r D+a3,r Ur+1 +a4,r

=ξ(1 −ξ)∂Ur+1

∂ξ ,(112)

Dfr+1 =Ur+1, fr+1(ηNx, ξ) = fw(113)

1

P r D2+b1,rD+b2,r θr+1 +b3,r

=ξ(1 −ξ)∂θr+1

∂ξ ,(114)

with boundary conditions;

Ur+1(ηNx, ξ) = 1 + 1

βSfU0

r+1(ηNx, ξ),

θr+1(ηNx, ξ) = 1 + STθ0

r+1(ηNx, ξ),

Ur+1(η0, ξ)→0, θr+1(η0, ξ)→0.(115)

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

204

Volume 21, 2022

We now apply the ﬁnite diﬀerence scheme on (112) -

(115) in the ξ−direction with centering about the mid-

point ξn+1

2to get,

A1Un+1

r+1 =B1Un

r+1 +K1,(116)

subject to

Ur+1(η0, ξn) = (1 + 1

β)SfU0

r+1(η0, ξn),

Ur+1(ηNx, ξn)→0, n = 0,1,2, . . . , Nt(117)

Uηj,0= erfc ηj

2, j = 0,1,2, . . . , Nx.(118)

where

A1=1

2(a1,r)n+1

2D2+ (a2,r)n+1

2D+ (a3,r)n+1

2

−ξn+1

2(1 −ξn+1

∆ξ,(119)

B1=−1

2(a1,r)n+1

2D2+ (a2,r)n+1

2D+ (a3,r)n+1

2

−ξn+1

2(1 −ξn+1

∆ξ,(120)

K1=−an+1

4,r ,(121)

AFn+1

r+1 =Un

r+1,(122)

subject to

Fr+1(ηNx, ξ) = fw.(123)

A2Θn+1

r+1 =B2Θn

r+1 +K2,(124)

subject to

Θr+1(ηNx, ξ) = 1 + STΘ0

r+1(ηNx, ξ),

Θr+1(η0, ξ)→0,(125)

where

A2=1

21

P r D2+ (b1,r)n+1

2D+b2,r)n+1

2

−ξn+1

2(1 −ξn+1

∆ξ,(126)

B2=−1

21

P r D2+ (b1,r)n+1

2D+b2,r)n+1

2

−ξn+1

2(1 −ξn+1

∆ξ,(127)

K2=bn+1

3,r ,(128)

The boundary conditions in (118) are imposed on the

ﬁrst and last rows of equation (116), those in equation

(123) are imposed on the last row of equation (122) and

those in equation (125) are imposed in the ﬁrst and last

rows of equation (124).

In this section we make use of the example from Patil

et al. [26] mixed convection from a vertical stretching

surface and viscous dissipation. The governing equation

is given by;

∂u

∂x +∂v

∂y = 0,(129)

u∂u

∂x +v∂u

∂y =Ue

dUe

dx +ν∂2u

∂y2

+gβ(T−T∞) + gβ∗(C−C∞),(130)

u∂T

∂x +v∂T

∂y =ν

P r

∂2T

∂y2+ν

Cp∂u

∂y 2

+QT−T∞

ρCp

,(131)

u∂C

∂x +v∂C

∂y =ν

∂2C

∂y2,(132)

In the equations, Bis the magnetic ﬁeld strength, β

and β∗are the thermal and concentration coeﬃcients

of volumetric expansion respectively, σis the electrical

conductivity. The boundary conditions are given as;

y=0:, u =Uw(x), v =vw,

T=Tw,=T∞+ (Tw0−T∞)exp(2x

L),

C=Cw,=C∞+ (Cw0−C∞)exp(2x

L),

y→ ∞ :, u →Ue(x), T →T∞, C →C∞.(133)

Using the following dimensionless variables,

ξ=x

L, η =U0

νxexp(x

2L)y, (134)

ψ(x, y) = (νU0x)1

2exp(x

2L)f(ξ, η),(135)

T−T∞= (Tw−T∞)G(ξ, η),(136)

T−T∞= (Tw−T∞)e(x

2L),(137)

C−C∞= (Cw−C∞)H(ξ, η),(138)

C−C∞= (Cw−C∞)e(x

2L),(139)

u=∂ψ

∂y , v =−∂ψ

∂x , u =U0e(x

L)f0,(140)

v=−(νU0

x)1

2e(x

2L)(1 + ξ)f

2+ξ∂f

∂ξ +η

∂f

∂η ,(141)

The following system of partial diﬀerential equations

is obtained,

f000 +1

2(1 + ξ)ff00 −ξf02+ξRi(G+NH) + ξ2

= 4ξf0∂f0

∂ξ −f00 ∂f

∂ξ ,(142)

. Example 4: Double diffusive mixed

convection flow from a vertical

exponentially stretching surface in the

presence of viscous dissipation

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

205

Volume 21, 2022

G00 +P r

2(1 + ξ)fG0−2P rξf0G+P rEc(f00 )2

+ReΓP rξG =P rξ f0∂θ

∂ξ −θ0∂f

∂ξ ,(143)

H00 +Sc

2(1 + ξ)fH0−2Scξf0H

=Scξ f0∂H

∂ξ −H0∂f

∂ξ ,(144)

with boundary conditions;

f0(ξ, 0) = 0, f(ξ, 0) = 0, θ(ξ, 0) = 1, H(ξ, 0) = 1,

f0(ξ, ∞) = 1, θ(ξ, ∞)=0, H(ξ, ∞) = 0.(145)

In the above equations, for the case ξ= 0, the equa-

tions do not reduce to those of Magyari and Keller [33],

Bidin and Nazar [29] and Mukhopadhyay [30] as men-

tioned in Patil et al. [26], but rather to those of Watan-

abe and Pop [27]. The spectral relaxation method is

applied to the system (142) and (145) and becomes;

U00

r+1 +a1,rU0

r+1 +a2,r =a3,r

∂Ur+1

∂ξ ,(146)

Ur+1(ξ, 0) = 1, Ur+1(ξ, ∞)=0,(147)

r+1 =Ur+1, fr+1(ξ, 0) = 0,(148)

G00

r+1 +b1,rG0

r+1 +b2,rGr+1 +b3,r

=b4,r

∂Gr+1

∂ξ ,(149)

Gr+1(ξ, 0) = 1, Gr+1(ξ, ∞)=0,(150)

H00

r+1 +c1,rH0

r+1 +c2,rHr+1 =c3,r

∂Hr+1

∂ξ ,(151)

Hr+1(ξ, 0) = 1, Hr+1(ξ, ∞) = 0,(152)

where

a1,r =1

2(1 + ξ)fr+ξ∂fr

∂ξ ,(153)

a2,r =−ξU 2

r+ξRi(Gr+NHr) + ξ2,(154)

a3,r =ξUr,(155)

b1,r =P r

2(1 + ξ)fr+P rξ ∂fr

∂ξ ,(156)

b2,r =−ξP r(2Ur−ReΓ),(157)

b3,r =P rEc(U0

r)2, b4,r =ξP rUr,(158)

c1,r =Sc

2(1 + ξ)fr+Scξ ∂fr

∂ξ ,(159)

c2,r =−2ScξUr,(160)

c3,r =ξScUr,(161)

We now apply the ﬁnite diﬀerence scheme on (146) -

(152) in the ξ−direction with centering about the mid-

point ξn+1

2to obtain,

A1Un+1

r+1 =B1Un

r+1 +K1,(162)

subject to

Ur+1(ηNx, ξ)=1, Ur+1(η0, ξ)→0,(163)

where

A1=1

2D2+ (a1,r)n+1

2D

−(a3,r)n+1

∆ξ,(164)

B1=−1

2D2+ (a1,r)n+1

2D

−(a3,r)n+1

∆ξ,(165)

K1=−an+1

2,r ,(166)

AFn+1

r+1 =Un

r+1,(167)

subject to

Fr+1(ηNx, ξ)=0.(168)

A2θn+1

r+1 =B2θn

r+1 +K2,(169)

subject to

θr+1(ηNx, ξ)=1, θr+1(η0, ξ)→0,(170)

where

A2=1

2D2+ (b1,r)n+1

2D+b2,r)n+1

2

−(b4,r)n+1

∆ξ,(171)

B2=−1

2D2+ (b1,r)n+1

2D+b2,r)n+1

2

−(b4,r)n+1

∆ξ,(172)

K1=−bn+1

3,r ,(173)

A3Hn+1

r+1 =B3Hn

r+1 +K3,(174)

subject to

Hr+1(ηNx, ξ)=1, Hr+1(η0, ξ)→0,(175)

where

A3=1

2D2+ (c1,r)n+1

2D+ (c2,r)n+1

2

−(c3,r)n+1

∆ξ,(176)

B3=−1

2D2+ (c1,r)n+1

2D+ (c2,r)n+1

2

−(c3,r)n+1

∆ξ,(177)

K3=0.(178)

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

206

Volume 21, 2022

The spectral relaxation method (SRM) for solving

boundary value problems is presented in this paper, four

examples of boundary value problems were presented in

the previous sections. In this section we discuss the re-

sults for all numerical examples considered. The method

depends on the length of the governing domain (b−a)

and the number of collocation points Nsometimes re-

ferred to as grid points. The domain of the given nu-

merical examples are deﬁned on [0,∞), to implement

the SRM scheme we need to ﬁnd the appropriate ﬁnite

value η∞which is large enough to approximate inﬁnity

for the given example. We start with an initial guess and

solve the SRM scheme over [0, η∞] so we obtain solutions

for the ﬂuid ﬂow functions f(η), g(η), h(η) for example

1, f(η), θ(η), φ(η) for example 2, f(η), θ(η) for example

3 and f(η), θ(η), φ(η) for example 4. The derivatives of

these functions at η= 0 are also computed to obtain the

important values f00 (0) (skin friction coeﬃcient), −θ0(0)

(heat transfer coeﬃcient) and −φ0(0) (mass transfer co-

eﬃcient).

The numerical scheme is implemented by increasing

the value of η∞by one and the solution process repeated

until the diﬀerence in values of f00 (0),−θ0(0) and −φ0(0)

between current and previous solution is less that a spec-

iﬁed tolerance level. If the optimal value of η∞is iden-

tiﬁed for a particular set governing constants, the min-

imum number of grid points is obtained using a similar

process. Starting with a small value of N, say N= 40,

the value of Nis increased by 10 and the solution is

solved until the results do not change within the spec-

iﬁed tolerance level. The accuracy of the SRM results

presented in this Chapter validated by comparison with

those previously published in the literature. The results

were compared against the MATLAB bvp4c solver. The

comparisons show a good agreement. It is also observed

that the error decrease at every iteration showing the

convergence and stability of the solutions produced by

the SRM.

Table 1: Comparison of the values of f00 (0) obtained by

SRM against those of Ece [24] and bvp4c for example 1

for =M= 0, 90 iterations and η∞= 10

Pr f00 (0)Ece [24] f00 (0) bvp4c f00 (0) SRM

1 0.68150212 0.68147996 0.68147995

10 0.43327726 0.43327802 0.43327802

In Table 1, the SRM results for skin friction coeﬃcient

agrees with the results of Ece [24] and MATLAB bvp4c

to six decimal places, this shows the accuracy of the SRM

for example 1.

Table 2: Comparison of the values of f00 (0) obtained by

SRM against those of Awad [25] and bvp4c for example

2 with Sc = 0.2, N1= 0.5, Df= 0.1, Sr= 0.3, P r =

0.71, λ = 1.90 iterations and η∞= 12

Λf00 (0) Awad [25] f00 (0) bvp4c f00 (0) SRM

0.3 1.0160906 1.0160906 1.0160906

0.5 0.9748755 0.9748755 0.9748755

In Table 2, the SRM results for skin friction coeﬃcient

agrees with the results of Awad [25] and MATLAB bvp4c

to six decimal places, this shows the accuracy of the SRM

for example 2.

Table 3: Comparison of the values of −θ0(0) obtained by

SRM against those of Hassanien and Al-Arabi. [31] and

QLM for example 3 with γ= Λ = σ=R=ξ= 0, n =

0, P r = 0.7 90 iterations and η∞= 10

ξ−θ0(0) [31] −θ0(0) QLM −θ0(0) SRM

0 0.47204 0.472035 0.472035

0.5 0.46347 0.463471 0.463471

In Table 3, the SRM results for the heat transfer coef-

ﬁcient agrees with the results of Hassanien and Al-Arabi.

[31] and QLM to six decimal places, this shows the ac-

curacy of the SRM for example 3.

Table 4: Comparison of the values of −θ0(0) obtained

by SRM against those of Watanabe [27] and bvp4c for

example 4 with M=N= 0, Sc = 1, P r = 1.

ξ−θ0(0) [27] −θ0(0) QLM −θ0(0) SRM

0 0.33026 0.33025734 0.33025734

0.5 0.40280 0.40279671 0.40279671

In Table 4, the SRM results for the heat transfer coef-

ﬁcient agrees with the results of Watanabe [27] and QLM

to eight decimal places, this shows the accuracy of the

SRM for example 4.

. Results and discussion

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

207

Volume 21, 2022

0 10 20 30 40 50 60 70 80

−18

−16

−14

−12

−10

−8

−6

−4

−2

Iterations

ln(Ed)

(a)

0 10 20 30 40 50 60 70 80 90

−25

−20

−15

−10

−5

Iterations

ln(Ed)

(b)

Fig. 1: Logarithm of the SRM decoupling error for (a)

Example 1 and (b) Example 2.

In Figure 1, the logarithm of the SRM decoupling er-

ror decreases with increasing iteration showing the rapid

convergence and stability of the SRM for (a) Example

1 and (b) Example 2. This is the standard SRM with-

out the successive over-relaxation (SOR). Convergence

is also aﬀected by the parameter values. These aﬀect in

many diﬀerent ways, they can slow convergence or cause

complete divergence. When certain parameters assume

higher or even lower values, the method might require

more iterations to converge. In the graph, the larger the

negative values on the vertical axis, the smaller the er-

ror. The SRM requires the use of the SOR to control its

convergence. In these two examples, the SRM requires

more than 70 iterations. Example 1 required more iter-

ations to achieve a larger error than that of example 2

that required less number of iterations.

0 10 20 30 40 50 60 70 80

−25

−20

−15

−10

−5

Iterations

ln(Ed)

ω = 0.9

ω = 1

ω = 1.1

(a)

0 10 20 30 40 50 60 70 80 90 100

−25

−20

−15

−10

−5

Iterations

ln(Ed)

ω = 0.9

ω = 1

ω = 1.1

(b)

Fig. 2: Logarithm of the SRM with SOR decoupling

error for (a) Example 1 and (b) Example 2.

The convergence of the SRM can be controlled by us-

ing the successive over-relaxation (SOR) method. This is

achieved by varying the relaxation parameter ω. In Fig-

ure 2, when ω= 0.9 the SRM converges faster and when

ω= 1.1 slows down convergence. The value ω= 1 cor-

responds to the standard SRM (without SOR). The ap-

plication of the successive over-relaxation was consistent

in both examples 1 and 2. The eﬀect of over-relaxation

was is more enhanced in example 2 than in example 1.

This is caused by the parameter values in each of these

examples. The SRM will require a smaller value of ωto

achieve a faster convergence. For values of ω > 1 (purple

line) more iterations are required. For ω= 1 (Red line),

similar to Figure 1. ω < 1 (black line), the number of

iterations is reduced signiﬁcantly.

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

208

Volume 21, 2022

0123456789

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

H(η)

ε=0,M=0

ε=0,M=1

ε=0,M=2

ε=0,M=5

ε=2,M=0

ε=2,M=1

ε=2,M=2

ε=2,M=5

(a)

0123456789

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

H(η)

Pr=1,M=0

Pr=1,M=1

Pr=1,M=2

Pr=1,M=5

Pr=10,M=0

Pr=10,M=1

Pr=10,M=2

Pr=10,M=5

(b)

Fig. 3: Variation of (a) spin parameter and (b) magnetic

parameter on temperature proﬁles for Example 1

Figure 3 shows the variation the spin parameter 

and Prandtl number P r in the presence of the magnetic

ﬁeld on temperature proﬁles. It is observed that increas-

ing both the spin parameter and Prandtl number reduce

temperature proﬁles. Increasing the spin parameter has

an eﬀect of reducing the velocity proﬁles, this is because

the direction of rotation is at right angles with the direc-

tion of ﬂuid ﬂow. The reduction in ﬂuid velocity in the

direction of ﬂow cause reduction in heat transfer, thereby

reducing temperature proﬁles. Increasing the magnetic

parameter will have the same eﬀect of reducing velocity

proﬁles and also reducing heat transfer. Increasing the

Prandtl number has an eﬀect of reducing thermal diﬀu-

sivity as shown in this ﬁgure. These results are consistent

with those reported by Ece [24]. This conﬁrms that the

SRM is accurate and can also be used as an alternative

method for solving

0123456789

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

ε=0,M=0

ε=0,M=1

ε=0,M=2

ε=0,M=5

ε=2,M=0

ε=2,M=1

ε=2,M=2

ε=2,M=5

(a)

0123456789

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

Pr=1,M=0

Pr=1,M=1

Pr=1,M=2

Pr=1,M=5

Pr=10,M=0

Pr=10,M=1

Pr=10,M=2

Pr=10,M=5

(b)

Fig. 4: Variation of (a) spin parameter and (b) magnetic

parameter on velocity proﬁles for Example 1

Figure 4 shows the variation the spin parameter 

and Prandtl number P r in the presence of the magnetic

ﬁeld on velocity proﬁles. It is observed that increasing

the spin parameter result in the decrease in velocity pro-

ﬁles and increasing the Prandtl number reduce velocity

proﬁles. The increase in spin parameter result in the de-

crease in velocity proﬁles, this caused by the spin velocity

which is tangential to the velocity of the ﬂuid ﬂow. This

has a dragging eﬀect on the ﬂuid ﬂow. The Increase in

the magnetic parameter result in the retardation of ﬂuid

ﬂow as it is directed perpendicular to the ﬂuid ﬂow. In-

creasing the Prandtl number result in the direct increase

of momentum diﬀusivity. This agrees with the result re-

ported in Ece [24]. This shows that the SRM is accurate

and can be used to solve boundary value problems in

ﬂuid ﬂow.

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

209

Volume 21, 2022

012345678

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

θ(η)

N1=0.1

N1=0.3

N1=0.7

N1=1.3

(a)

012345678

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

θ(η)

λ=0

λ=0.3

λ=0.6

λ=1

(b)

Fig. 5: Variation of (a) concentration buoyancy param-

eter and (b) power-law index on tempearture proﬁles for

Example 2

Figure 5 shows the variation of concentration buoy-

ancy parameter N1and power-law index λon temper-

ature proﬁles. Increasing both the concentration buoy-

ancy parameter and power-law index result in the de-

crease in temperature proﬁles. Increasing the concentra-

tion buoyancy parameter result in the decrease in tem-

perature proﬁles, the high buoyancy aids ﬂuid motion.

The presence of the solute cause a decrease in temper-

ature proﬁles as shown in Figure 5 (a). Also increasing

the power-law index result in the decrease in tempera-

ture proﬁles as shown in Figure 5 (b). The results are

consistent with those of Awad et. al [25].

012345678

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

f′(η)

N1=0.1

N1=0.3

N1=0.7

N1=1.3

(a)

012345678

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

f′(η)

λ=0

λ=0.3

λ=0.6

λ=1

(b)

Fig. 6: Variation of (a) concentration buoyancy param-

eter and (b) power-law index on velocity proﬁles for Ex-

ample 2

Figure 6 shows the eﬀect of increasing the buoy-

ancy parameter N1and power-law index λ.Increasing the

buoyancy parameter result in the increase in velocity pro-

ﬁles, concentration buoyancy tend to aid ﬂuid motion,

the solute is pushed by the buoyancy force and increase

ﬂuid motion as shown in 6 (a). Increase in power-law

index result in the decrease in velocity proﬁles as shown

in Figure 6 (b). These results agree with those of Awad

et. al [25].

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

210

Volume 21, 2022

0 5 10 15 20 25 30 35 40

10−12

10−10

10−8

10−6

10−4

10−2

100

Iterations

Ef, E θ

Eθ

(a)

(b)

Fig. 7: (a) Convergence graph for f(η, ξ), θ(η, ξ), (b)

Residual graph for f(η, ξ) for Example 3

Figure 7 indicates the error graphs for example 3.

Figure 7 (a) shows the rapid convergence of the SRM

after ten iterations, thereafter it is observed that the

SRM is stable. The stability of the SRM was observed for

both momentum and thermal equations. The errors are

as small as the order of 10−12. In Figure 7 (b) we observe

the increase in the error as ξincrease. The error decrease

with increasing iterations. The most important result is

the eﬀect of the time dimension ξon error propagation.

0 5 10 15 20 25 30 35 40

10−10

10−8

10−6

10−4

10−2

100

102

Iterations

Ef, E G, E H

(a)

0 0.2 0.4 0.6 0.8 1

10−7

10−6

10−5

10−4

10−3

10−2

10−1

kR e s(f)k∞

Ite r.

(b)

Fig. 8: (a) Convergence graph for

f(η, ξ), G(η, ξ), H(η, ξ), (b) Residual graph for f(η, ξ)

for Example 4

The accuracy of the SRM for example 4 is shown in

Figure 8 (a) and Figure 8 (b), showing the reduction

of the error with increasing number of iterations. It is

observed that there is rapid convergence of the SRM in

about 10 iterations, the error being of the order of 10−11.

The stability of the SRM is shown after 10 iterations,

the variation of the error does not change much showing

the accuracy of the method. Figure 8 (b) shows the

plot of residual error of the function f(η, ξ) against ξ

with increasing iteration. It is observed that the rate

of increasing the residual error is almost linear and is

minimum at ξ= 0. The error increases sharply and

reach a level where it is almost constant. We further

observe that the error is small after two iterations in the

entire range of the values of ξ. The small size of the

residual error indicate the accuracy of the method.

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

211

Volume 21, 2022

Table 5: Heat transfer coeﬃcient −θ0(0, ξ) for various

values of ξcomputed using SRM for example 3.

ξ/∆ξ0.004 0.001 0.0005

0.1 0.475693 0.475694 0.475694

0.3 0.483407 0.483407 0.483407

0.5 0.491688 0.491688 0.491688

0.7 0.500488 0.500488 0.500488

CPU Time 1.64 5.98 11.79

Tables 5 and 6 show the numerical values of values of

the heat transfer coeﬃcient −θ0(0) for diﬀerent values of

∆ξ, computed using the SRM and the QLM respectively

for equation (99)-(101).The total computational time to

execute the integration is also shown. The computation

was done using the same number of collocation points

Nxand η∞. Reducing the step size ∆ξimproves the

accuracy of the results until they are consistent to within

six decimal places.

Table 6: Heat transfer coeﬃcient −θ0(0, ξ) for various

values of ξcomputed using QLM for example 3.

ξ/∆ξ0.004 0.001 0.0005

0.1 0.465195 0.467389 0.467922

0.3 0.459767 0.461433 0.461856

0.5 0.455776 0.457277 0.456466

0.7 0.457627 0.456340 0.500488

CPU Time 5.16 20.1 40.5

The SRM takes less computational time than the

QLM, we observe that the SRM converge more rapidly

than the QLM when the step size is reduced. Full con-

vergence to six decimal digits is arrived when ∆ξ= 0.001

in the SRM compared to the QLM which showed con-

vergence at ∆ξ= 0.0002.

Table 7: Skin friction coeﬃcient f00 (0, ξ) for various val-

ues of ξcomputed using SRM for example 4.

ξ/∆ξ0.004 0.001 0.0005

0.1 0.287824 0.287832 0.287832

0.2 0.249535 0.249542 0.249542

0.3 0.216382 0.216387 0.216387

0.4 0.187672 0.187677 0.187677

CPU Time 4.03 15.8 35.9

Table 8: Skin friction coeﬃcient f00 (0, ξ) for various val-

ues of ξcomputed using QLM for example 4.

ξ/∆ξ0.004 0.001 0.0005

0.1 0.286511 0.287828 0.287832

0.2 0.248912 0.249538 0.249542

0.3 0.215822 0.216291 0.216387

0.4 0.186538 0.187482 0.187677

CPU Time 5.16 20.1 40.9

The results shown in Tables 7 and 8 show the ac-

curacy of the SRM, the results for the skin friction co-

eﬃcient become consistent to six decimal places when

∆ξ= 0.001 compared to the QLM which become con-

sistent at ∆ξ= 0.0005. This shows that the SRM con-

verges more rapidly than the QLM. The computation

times shown also reveal that the SRM requires less com-

putation time to execute the integration, these results

were also observed in example 3. For further reading,

readers are referred to the works of Magagula et al. [33],

Motsa et al. [34], Kamwswaran et al. [35], Agbaje and

Motsa [36], Shateyi [37] and Motsa and Makukula [38].

The objective of this paper was to describe the spec-

tral relaxation method (SRM) and consider diﬀerent ex-

amples in implementing it. Two examples involving or-

dinary diﬀerential equations and two examples involv-

ing partial diﬀerential equations were considered. The

method is easy to implement as compared to ﬁnite dif-

ference methods such as the Keller-box method. The

SRM is a robust and accurate method to solve diﬀeren-

tial equations and requires less computation time than

the quasi-linearization method. The method can be used

to solve boundary value problems. The method give re-

sults that are accurate in the speciﬁed space and time do-

mains. The residual error quickly approaches zero within

a few iterations. The method make use of well-known

Chebyshev spectral collocation method for discretiza-

tion. The method has also been improved by considering

discretization in both space and time. The challenge in

the method is the determination of the appropriate num-

ber of collocation points or grid points.

I would like to thank the Faculty Research and Inno-

vation Committee (FRIC) of the Central University of

Technology, Free State for the sponsorship of this article

registration fees, travel and accommodation.

[1] M. Tatari, M. Haghighi, A generalized Laguerre

spectral collocation method for solving initial-

boundary value problems, Applied Mathematical

Modelling,38, 1351-1364, 2014

. Conclusion

Acknowledgment

References

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

212

Volume 21, 2022

[2] M. T. Vilhena, L. B. Barichello, J. R. Zabadal, C. F.

Segatto, A. V. Cardona, R. P. Pazos, Solutions to

the multidimensional linear transport equation by

the spectral method, Progress in Nuclear Energy,

35, 275-291, 1999

[3] B. Orel, A. Perne, Chebyshev-Fourier spectral

methods for nonperiodic boundary value problems.

Journal of Applied Mathematics,2014, Article ID

572694, 10 pages, doi: 10.1155/2014/572694, 2014

[4] F. Fakhar-Izadi, M. Dehghan, Spectral methods

for parabolic Volterra integro-diﬀerential equations,

Journal of Computation and Applied Mathematics,

235, 4032-4046, 2011

[5] R. Moulla, L. Lefevre,B. Maschke, Pseudo-spectral

methods for the spatial sympletic reduction of open

systems of conservation laws, Journal of Computa-

tional Physics,231, 1272-1292, 2012

[6] G. Makanda, O. D. Makinde, P. Sibanda, Natural

convection of viscoelastic ﬂuid from a cone embed-

ded in a porous medium with viscous dissipation,

Mathematical Problems in Engineering,2013, Arti-

cle ID 934712, 11 pages, doi: 10.1155/2013/934712,

2013

[7] S. Motsa, S. Shateyi, A successive linearization

method approach to solve Lane-Emdem type of

equations, Mathematical Problems in Engineer-

ing,2012, Article ID 280702, 14 pages, doi:

10.1155/2012/280702, 2012

[8] S. Motsa, P. Sibanda, A note on the solutions of

the Van der Pol and Duﬃng equations using a lin-

earization method, Mathematical Problems in En-

gineering,2012, Article ID 693453, 10 pages, doi:

10.1155/2012/693453, 2012

[9] S. S. Motsa, V. M. Magagula, P. Sibanda, A bivarate

Chebyshev spectral collocation quasilinearization

method for nonlinear evolution parabolic equations,

The Scientiﬁc World Journal,2014, Article ID

581987, 13 pages, doi: 10.1155/2014/581987, 2014

[10] A. Rashid, M. Abbas, A.I. Ismail, A. A. Md

Majid, Numerical solution of the coupled vis-

cous Burgers equations by Chebyshev-Legendre-

Legendre pseudo-spectral method, Applied Mathe-

matics and Computation,245, 372-381, 2014

[11] M. A. Christou, N. C. Papanicolaou, Kawahara

solutions in Boussinesq equations using a robust

Christov-Galerkin spectral method, Applied Math-

ematics and Computation,243, 245-257, 2014

[12] D. Gottlieb, J. S. Hesthaven, Spectral methods for

hyperbolic problems. Journal of Computation and

Applied Mathematics,128, 83-131, 2001

[13] D. Gottlieb, S. A. Orszag, Numerical analysis of

spectral methods: Theory and applications, Capital

City Press, Montpelier, Vermont, USA, 1977

[14] R. Jovanovic, S. Kais, F. H. Alharbi, Spectral

method for solving the nonlinear Thomas-Fermi

equation based on exponential functions, Journal

of Applied Mathematics,2014, Article ID 168568, 8

pages, doi: 10.1155/2014/168568, 2014

[15] J. Zhou, The constants in a posteriori error in-

dicator for state-constrained optimal control prob-

lems with spectral methods, Abstracts and Applied

Analysis,2014, Article ID 946026, 8 pages, doi:

10.1155/2014/946026, 2014

[16] A. Kadem, Analytical solutions for the neutron

transport using the spectral methods, International

Journal of Mathematical Science,2006, Article ID

16214, 11 pages, doi:10.1155/IJMMS/2006/16214,

2006

[17] W. Wang, C. Xu, Spectral methods based on new

formulations of coupled Stokes and Darcy equations,

Journal of Computational Physics,257, 126-142,

2014

[18] I. Kozar ,N. T. Malic, Spectral method in realistic

modelling of bridges under moving vehicles, Engi-

neering Structures,50, 149-157, 2013

[19] C. Pozrikidis, A spectral collocation method with

triangular boundary elements, Engineering Analysis

with Boundary Elements,30, 315-324, 2006

[20] P. G. Dlamini, M. Khumalo, S. S. Motsa, A Note on

the multi-stage spectral relaxation method for chaos

control and synchronization, International Journal

of Nonlinear Sciences and Numerical Simulation,

15, 289-298, 2014

[21] S. Shateyi, G. T. Marewo, ”Numerical analysis of

MHD stagnation point ﬂow of Casson ﬂuid, heat

and mass transfer over a stretching sheet” In Bal-

icki, J. (Ed) Advances in Applied and Pure Math-

ematics, WSEAS,Proceedings of the 7th interna-

tional conference on ﬁnite diﬀerences, ﬁnite ele-

ments, ﬁnite volumes, boundary elements (F-and-

B’14), Gdansk, Poland, 128-132, 2014 ISBN: 978-

960-474-380-3

[22] S. S. Motsa, Z. G. Makukula, On spectral relax-

ation method approach for steady Von Karman ﬂow

of a Riner-Rivlin ﬂuid with Joule heating, viscous

dissipation and suction/injection, Central European

Journal of Physics,11, 363-374, 2013

[23] S. S. Motsa, Z. G. Makukula, S. Shateyi, Spectral

local linearization approach for natural convection

boundary layer ﬂow, Mathematical Problems in En-

gineering,2013, Article ID 765013, 7 pages, doi:

10.1155/2013/765013, 2013

[24] C. M. Ece, Free convection ﬂow about a vertical

spinning cone under a magnetic ﬁeld, Applied Math-

ematics and Computation,179, 231-242, 2006

[25] F. G. Awad, P. Sibanda, S. S. Motsa, O. D.

Makinde, Convection from an inverted cone in a

porous medium with cross-diﬀusion eﬀects. Comput-

ers and Mathematics with Applications,61, 1431-

1441, 2011

[26] P. M. Patil, D. N. Latha, S. Roy,E. Momoniat, Dou-

ble diﬀusive convection ﬂow from a vertical expo-

nentially stretching surface in the presence of vis-

cous dissipation, International Journal of Heat and

Mass Transfer,112, 758-766, 2017

[27] T. Watanabe, I. Pop, Eﬀects of viscous dissipation

and stress work on the steady two dimensional lam-

inar magnetohydrodynamic forced convection ﬂow

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

213

Volume 21, 2022

over a non-isothermal wedge, International Com-

munications in Heat and Mass Transfer,20, 871,

1993

[28] E. Magyari,B. Keller, Heat and mass transfer in the

boundary layers on exponentially stretching con-

tinuous surface, Journal of Physics D: Applied

Physics,32, 577-585, 1999

[29] B. Bidin, R. Nazar, Numerical solutions of bound-

ary layer ﬂow over an exponentially stretching sur-

face with thermal radiation, European Journal of

Scientiﬁc Research,33, 710-717, 2009

[30] S. Mukhopadhyay, Slip eﬀects on MHD boundary

layer ﬂows over an exponentially stretching sheet

with suction/blowing and thermal radiation, Ain

Shams Engineering Journal,4, 485-491, 2013

[31] I. A. Hassanien, T. H. Al-arabi, Non-Darcy un-

steady mixed convection ﬂow near the stagnation

point on a heated vertical surface embedded in a

porous medium with thermal radiation and vari-

able viscosity. Commun Nonliner Sci Numer Sim-

ulat,14, 1366-1376, 2009

[32] S. S. Motsa, A new spectral relaxation method for

similarity variable nonlinear boundary layer ﬂow

systems, Chemical Engineering Communications,

201, 241-256, 2014

[33] V. M. Magagula, S. S. Motsa, P. Sibanda, P.

G. Dlamini, On the bivariate spectral relaxation

method for unsteady magneto-hydrodynamic

ﬂow in porous media. Springer plus,doi:

10.1186/s40064-016-2053-4, 2016

[34] S. S. Motsa, P. G. Dlamini, M. Khumalo,

Spectral relaxation method and spectral quasi-

lineariation method for solving unsteady bound-

ary layer ﬂow problems. Advances in Mathematical

Physics,doi:10.1155/2014/341964, 2014

[35] P. Kameswaran, P. Sibanda, S. S. Motsa, A spectral

relaxation method for thermal dispersion and radi-

ation eﬀects in a nanoﬂuid ﬂow, Boundary Value

Problems,2013, 2013:242, 2013

[36] T. M. Agbaje, S. S. Motsa, Comparison between

spectral perturbation and spectral relaxation ap-

proach for unsteady heat and mass transfer by MHD

mixed convection ﬂow over an impulsively stretched

vertical surface with chemical reaction eﬀect. Jour-

nal of Interpolation and approximation in Scientiﬁc

Computing,2015, 48-83, 2015

[37] S. Shateyi, On spectral relaxation method for an

MHD ﬂow and heat transfer of Maxwell ﬂuid, Inter-

national Conference on Mechanics, Heat and Mass

Transfer,ISBN: 978-1-61804-220-0, 102-106, 2014

[38] S. S. Motsa, Z. G. Makukula, On spectral relax-

ation method approach for steady von-Karman ﬂow

of a Reiner-Rivlin ﬂuid with Joule heating, viscous

dissipation and suction/injection, Central European

Journal of Physics,11: 363-374, 2013

Creative Commons Attribution License 4.0

(Attribution 4.0 International, CC BY 4.0)

This article is published under the terms of the Creative

Commons Attribution License 4.0

https://creativecommons.org/licenses/by/4.0/deed.en_US

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2022.21.22

Gilbert Makanda

E-ISSN: 2224-2678

214

Volume 21, 2022