Splitting-based schemes for numerical solution of ... - Semantic Scholar
Applied Mathematics and Computation 219 (2013) 8467–8485
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Splitting-based schemes for numerical solution of nonlinear diffusion equations on a sphere Yuri N. Skiba a, Denis M. Filatov b,⇑ a Centro de Ciencias de la Atmósfera (CCA), Universidad Nacional Autónoma de México (UNAM), Av. Universidad #3000, Cd. Universitaria, C.P. 04510 México D.F., Mexico b Centro de Investigación en Computación (CIC), Instituto Politécnico Nacional (IPN), Av. Juan de Dios Bátiz s/n, C.P. 07738 México D.F., Mexico
a r t i c l e
i n f o
Keywords: Nonlinear diffusion equation Dissipative balanced finite difference schemes Operator splitting on closed nonperiodic manifolds
a b s t r a c t We provide an advanced study of our recently developed method for the numerical solution of nonlinear diffusion equations on a sphere. In particular, we analyse the method in detail when applied to solving diverse diffusion phenomena, with specific conditions on the smoothness of the solution, the degree of nonlinearity, and the initial data and sources. The main idea of the method consists in splitting the original differential operator by coordinates and subsequent constructing finite difference schemes for the split one-dimensional problems using different coordinate maps for the sphere at the two split time intervals. The essential advantage of this technique is that each split 1D equation can be equipped with a periodic boundary condition, despite the sphere being not a doubly periodic domain. Therefore, unlike the existing methods, this one does not require applying special numerical procedures for careful computing the solution near the poles, which is always a challenge. Each split 1D equation is approximated by a second- or a fourth-order finite difference scheme that keeps all the substantial properties of the differential problem: it is balanced and dissipative, since the spatial finite difference operator is negative definite. The developed algorithm is cheap-to-implement from the computational point of view. The theoretical results are confirmed numerically by simulating various nonlinear diffusion processes. A numerical example, in particular, shows that the competition of the three basic mechanisms—the nonlinear interaction, forcing and dissipation—can generate wave solutions, whose spatial structures on the sphere are subjected to the alternating influence of the processes of self-organisation and self-destruction. The accuracy of the method is evaluated by comparing the numerical solutions versus the analytical ones obtained with specially chosen forcings. The convergence of the numerical solution to the analytics is verified by refining spatial grids. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction A large number of important natural phenomena, e.g., heat transfer in ionised gases, gas percolation through porous media, concentration waves in distributed chemical reactors, combustion, viscous processes in diverse media, as well as many others are described by (nonlinear) diffusion equations [1–10]. Possibly, the most obvious and well-studied applications of the diffusion equation lie in atmospheric and hydrodynamic sciences. Hundreds of papers providing researches on some or other issues of the environment are published every year. ⇑ Corresponding author. E-mail addresses: [email protected] (Y.N. Skiba), denisfi[email protected] (D.M. Filatov). 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.02.066
8468
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Among the applications there are air and water pollution, viscous flows in river channels—just to mention a couple of them [9,10]. The governing equations are specific for each application, but in general they contain a hyperbolic part responsible for the advection process (velocity field) and a parabolic term aimed to simulate fluctuation from the mainstream, or diffusion. Another important application of the diffusion equation, especially involving nonlinearity, is blow-up. In [8], among other issues, a nonlinear diffusion Cauchy problem in Rn is studied and numerical experiments are performed, simulating blow-up in a bounded domain. A not less interesting, though a little less obvious problem described by the diffusion equation is the simulation of time evolution of the orientational distribution function of non-spherical particles due to interactions with each other and external fields [11]. In many papers on diffusion equations boundary value problems are studied [12–15], while much less papers are dedicated to the numerical simulation of nonlinear diffusion processes on closed manifolds. The latters differ from bounded domains as having no boundaries, which makes it nontrivial to develop efficient numerical methods for solving diffusion equations in such domains. Although for some particular cases these problems can be reduced to boundary value problems with boundary conditions of special types (e.g., a torus can be represented as a doubly periodic domain, so that periodic boundary conditions in both directions can be employed), this cannot be done for a more or less complicated manifold. A simple, but important example of such a manifold is a two-dimensional sphere. The case of the sphere can be seen as a step before more general manifolds, which can be interesting, e.g., in nanoelectronics [16]. The point is that the sphere is not a doubly periodic domain, since it is periodic in the longitude, but is not in the latitude due to the presence of two poles. Therefore, a numerical procedure designed the solve the original 2D problem will be computationally cumbersome, because the matrix of the resulting linear system will be of a general type, hence not permitting to apply some or other fast linear solvers. For instance, in [11] a Voronoi tessellation of the sphere is used and then a finite volume method is applied. As for making a preliminary splitting of the 2D equation [17], this will bring to the necessity of constructing mathematically and physically correct boundary conditions or involving special numerical procedures near the poles (e.g., matrix bordering) when computing in the latitudinal direction. These are always a problem, as the poles represent an artificial boundary appearing exclusively due to the latitudinal-longitudinal coordinate system, and the construction of proper artificial boundary conditions is a serious independent question [18]. The nonlinear diffusion equation on a two-dimensional sphere S ¼ fðk; uÞ : k 2 ½0; 2pÞ; u 2 p2 ; p2 g has the form
@T 1 @ lT a @T @ lT a cos u @T þ þ f: ¼ AT þ f @t R cos u @k R cos u @k @u @u R
ð1Þ
It is equipped with a suitable initial condition. Here A is the diffusion operator, T ¼ Tðk; u; tÞ P 0 is the unknown function, lT a is the diffusion coefficient, l ¼ lðk; u; tÞ P 0 is the amplification factor, f ¼ f ðk; u; tÞ is the source function, R is the radius of the sphere with k as the longitude (positive eastward) and u as the latitude (positive northward). The parameter a is usually a positive integer number that determines the degree of nonlinearity of the diffusion process; if a ¼ 0 then we are dealing with linear diffusion. Introducing on S the scalar product as
hf; gi ¼
Z
fgdS;
ð2Þ
S
multiplying (1) by g and integrating over the domain S, we find
@T ; g ¼ hAT; gi þ hf ; gi: @t
ð3Þ
In the case of g ¼ 1 the first summand on the right-hand side of (3) vanishes [19] and the equality turns into the balance equation
d dt
Z
TdS ¼ S
Z
fdS;
ð4Þ
S
which, if f ¼ 0, provides the mass conservation law
d dt
Z
TdS ¼ 0:
ð5Þ
S
In the case of g ¼ T due to the negative definiteness of A it holds
hAT; gi 6 0;
ð6Þ
and so
1 d 2 dt
Z S
T 2 dS 6
Z
fTdS; S
from where under f ¼ 0 we obtain the solution’s dissipation in the L2 ðSÞ-norm
ð7Þ
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
1 d 2 dt
Z
T 2 dS 6 0:
8469
ð8Þ
S
Our aim is to develop an accurate and efficient numerical method for solving Eq. (1). The accuracy implies that the method should provide physically correct numerical solutions which would possess properties (4) and (8). The efficiency means that the method should be computationally inexpensive. First results on this issue were presented by the authors in [20], so in this paper we study the developed method in detail, analysing how it works when applied to solving diverse diffusion phenomena, with specific conditions on the smoothness of the solution, the degree of nonlinearity, and the initial data and sources. Besides, we explicitly demonstrate the constructed finite difference operators inherit important properties of the original diffusion operator. The paper is organised as follows. In Section 2 we first linearise the original nonlinear differential equation and then apply the method of splitting [21] to the resulting equation on the sphere. Employing the procedure of alternation of two coordinate maps for representing the same sphere at each split time step, we become able to develop a numerical algorithm as if we dealt with a periodic domain in both spatial directions (although, repeat again, the sphere is not a doubly periodic manifold). This leads to using simple periodic boundary conditions while computing the solutions in k and in u. Next, involving the Crank–Nicolson approximation, we construct two second-order finite difference schemes for the split 1D diffusion equations. The schemes are shown to be dissipative and balanced, according to the properties of the original differential problem, as well as easy-to-implement from the computational standpoint. Then, because of the periodicity of the boundary conditions we increase the approximation’s accuracy to the fourth order in space. In Section 3 we provide results of numerical experiments performed to test the developed method. In several examples, from the simplest linear case to a complex nonlinear one, we demonstrate the schemes are accurately simulating diffusion phenomena. We prove the schemes inherit the essential properties of the differential operators and hence yield physically adequate numerical solutions. We comprehensively study the developed method, providing a comparison of the numerical solutions with the analytics. We also give a numerical example (see Experiment 4) demonstrating the synergistic action of the external forcing, nonlinearity and dissipation. These three mechanisms produce solutions, whose spatial structures assert the presence of the processes of selforganisation and self-destruction, which phenomenon is of interest in nonlinear science. In Section 4 we give a conclusion. 2. Mathematical foundation of the method 2.1. Operator splitting on the sphere In the standard manner we define the grid spacing in time as s ¼ tnþ1 t n and in every time interval ðtn ; tnþ1 Þ we linearise (1) and then split the linearised equation by coordinates [22,23]
@T f 1 @ D @T f þ ; ¼ Ak T þ @t 2 R cos u @k R cos u @k 2 @T f 1 @ D cos u @T f þ ; ¼ Au T þ @t 2 R cos u @ u @u 2 R
ð9Þ ð10Þ
where D ¼ lðT n Þa ; T n ¼ Tðk; u; t n Þ. Hereafter the split equations will be considered in time successively: the solution to Eq. (9), solved in k, will be used as the initial condition for Eq. (10); the latter, solved in u, will supply the initial condition for (9) in the next time interval ðtnþ1 ; tnþ2 Þ, and so on. Because the metric term R cos u vanishes at u ¼ p=2, both equations are meaningless at the poles. Therefore, for covering the entire sphere we assume Dk ¼ kkþ1 kk and Du ¼ ulþ1 ul and define on S a grid, making a half step shift in u [24], i.e. ð1Þ
SDk;Du ¼
Dk Dk p Du p Du ; ul 2 þ : ; 2p þ ðkk ; ul Þ : kk 2 ; 2 2 2 2 2 2
ð11Þ
Thereby we exclude the pole singularities (Fig. 1), so that the subsequent finite difference equations will have sense everyð1Þ where on SDk;Du . For solving (9) one has, evidently, to use the periodic boundary condition in k, as the sphere S is a periodic domain in the longitude. (Note that a half step grid shifting was also applied in [14,25], and, compared to the classical approach [12], it provided benefits in the solution.) As for solving Eq. (10), one faces the difficulty since the sphere is not a periodic domain in u. Several approaches can be used to resolve it, e.g., one could involve a matrix bordering procedure or try to paste the solution from the opposite meridians at the poles. A disadvantage of these techniques is that in the case of an implicit temporal approximation on the righthand side of (10) the resulting matrix will be of a general type, and so no fast algorithms of linear algebra can be used for computing the solution. Alternatively, an explicit temporal discretisation can be used, but this will impose serious restrictions on the timestep, especially when a is large. However, this is exactly the method of splitting that allows to avoid these undesired procedures which, besides, if performed inaccurately, may easily result in introducing nonphysical modes into the solution. Specifically, due to the splitting, for computing the solution in u we change the coordinate map from (11) to
8470
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 1. Shift of the original grid (dotted lines) a half step in u (solid lines) allows excluding the pole singularities, thereby keeping the equations to have sense on the entire sphere. The selected black parallel on the sphere corresponds to that on the plane.
Dk Dk p Du 3p Du ;p ; ul 2 þ ðkk ; ul Þ : kk 2 ; þ : 2 2 2 2 2 2
ð2Þ
SDk;Du ¼
ð12Þ
Obviously, grid (12) contains the same nodes as (11) (Fig. 2). The use of the two coordinate maps allows employing the same numerical algorithm with the periodic boundary conditions in both directions, k and u [26,27,20]. The only change we have to make in (10) if using (12) is to replace cos u with j cos uj, as well. 2.2. Second-order finite difference schemes According to the splitting, in each time interval ðtn ; tnþ1 Þ we discretise the temporal derivatives in (9) and (10) as
nþ1 T kl 2 T nkl @T ; @t t¼tn s
ð13Þ
nþ1 T nþ1 T kl 2 @T kl ; @t t¼tnþ1=2 s
ð14Þ
and
respectively. In its turn, for the spatial derivatives we take the second-order central finite difference stencils—
Ak Tjk¼kk ADk Tjk¼kk
R2 cos2 u
Au T u¼u ADu T u¼u l
l
Dkþ1=2
1
1 R2 j cos ul j
T kþ1 T k Dk
k1 Dk1=2 T k T Dk ; Dk
ðDj cos ujÞlþ1=2
T lþ1 T l Du
l1 ðDj cos ujÞl1=2 T l T Du
Du
ð15Þ
;
ð16Þ
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
8471
Fig. 2. The map swap, possible due to the splitting, allows using the periodic boundary condition in u, avoiding the question of constructing boundary conditions at the poles. The selected black meridian on the sphere corresponds to that on the plane.
where
Dkþp=2 :¼
Dkþðpþ1Þ=2 þ Dkþðp1Þ=2 2
ðDj cos ujÞlþp=2 :¼
for p ¼ 1;
ðDj cos ujÞlþðpþ1Þ=2 þ ðDj cos ujÞlþðp1Þ=2 2
ð17Þ
for p ¼ 1:
ð18Þ
(In order not to overload the formulas, the nonvarying index l in (15) and the nonvarying index k in (16) are omitted.) The 1 external forcing f is computed at the intermediate time moment tn þt2nþ1 , i.e. f nþ2 ¼ f k; u; tn þt2nþ1 . The finite difference operators in (15) and (16) are correct in the sense of keeping the property of negative definiteness (cf. (6)). Indeed, for (15), multiplying the right-hand side by T k (recall, l is fixed while computing in k) and summing all over the k’s, we obtain
! k1 X X Dk1=2 T k T 1 Dk ¼ D ðT T ÞT D ðT T ÞT T k kþ1=2 kþ1 k k k1=2 k k1 k 2 2 Dk R2 cos2 ul Dk2 k R cos ul k k ! X X X 1 1 Dkþ1=2 ðT kþ1 T k ÞT k Dk0 þ1=2 ðT k0 þ1 T k0 ÞT k0 þ1 ¼ 2 Dkþ1=2 ðT kþ1 T k ÞðT k T kþ1 Þ ¼ 2 2 u Dk2 R cos2 ul Dk2 R cos 0 l k k k X 1 Dkþ1=2 ðT kþ1 T k Þ2 6 0: ð19Þ ¼ 2 2 2 R cos ul Dk k
X
Dkþ1=2
1
0
T kþ1 T k Dk
ð1Þ
Here we denoted k ¼ k 1 and used the periodicity of the solution in k on SDk;Du to return to the original index k. Analoð2Þ gously, for (16), due to the periodicity of the solution in u on SDk;Du we have
8472
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
1
l
R2 j cos ul j ¼
T lþ1 T l Du
ðDj cos ujÞlþ1=2
X
l1 ðDj cos ujÞl1=2 T l T Du
Du
T l j cos ul j ¼ . . .
X ðDj cos ujÞlþ1=2 ðT lþ1 T l Þ2 6 0:
1 2
R Du 2
ð20Þ
l
Temporal discretisation of the function T kl in (15) and (16) can be done via the Crank–Nicolson approximation nþp
nþp1 2
T 2 þ T kl T kl :¼ kl 2
for p ¼ 1; 2:
ð21Þ
Therefore, each of the constructed split one-dimensional second-order finite difference schemes is dissipative (cf. (7) and (8))—
1 2
nþ2
T k
ð1Þ
L2 SDk;u
L2
l
nþ1 2
T l
ð2Þ k ;D u
L2 Sk
* nþ1 + 2 þ 2s f k ; T k ð1Þ 2 SDk;u
6
T n 2 k
* nþ1 + 2 þ 2s fl ; T l ð2Þ Sk ;Du 2
1 2 nþ 6
T l 2 L2
l
ð2Þ k ;D u
L2 Sk
k
ð22Þ
;
ð1Þ
L2 SDk;u
l
ð23Þ
:
Gathering all together, eventually for the second-order finite difference scheme in k we obtain nþ1
nþ12
T kþ12 mkþ1 þ T k
1
nþ1 f 2 1 nþ1 þ mk T k12 mk1 ¼ T nkþ1 mkþ1 þ T nk mk þ T nk1 mk1 þ k ; s s 2
ð24Þ
where
mk ¼
Dkþ1=2 þ Dk1=2 2R
2
cos2
u l Dk
2
;
mkþp ¼
Dkþp=2 2
2R cos2 ul Dk2
for p ¼ 1;
ð25Þ
while in u we have
T nþ1 lþ1 mlþ1
þ
T nþ1 l
1
s
þ ml
T nþ1 l1 ml1
¼
nþ1 T lþ12 mlþ1
þ
nþ1 Tl 2
1
s
ml
þ
nþ1 T l12 ml1
nþ12
þ
fl
2
;
ð26Þ
where
ml ¼ mlþp
ðDj cos ujÞlþ1=2 þ ðDj cos ujÞl1=2 2R2 j cos ul jDu2 ðDj cos ujÞlþp=2 ¼ 2 for p ¼ 1: 2R j cos ul jDu2
; ð27Þ
Supplementing them with the bicyclic splitting (or Strang splitting) nþ1
nþ2
T kl 4 T nkl f 4 ¼ ADk T kl þ k ; s=2 2 nþ2
nþ1
nþ2
T kl 4 T kl 4 f 4 ¼ ADu T kl þ l ; s=2 2 nþ3
ð28Þ
nþ2
nþ2
nþ34
nþ24
T kl 4 T kl 4 f 4 ¼ ADu T kl þ l ; s=2 2 T nþ1 T kl f kl ¼ ADk T kl þ k ; s=2 2
ð29Þ ð30Þ ð31Þ
where nþ
p
nþ
T 4 þ T kl T kl :¼ kl 2
p1 4
for p ¼ 1; 4;
ð32Þ
we achieve the second order accuracy in time in the interval ðtn ; tnþ1 Þ [23,28]. Since both the approximation and stability are fulfilled, the discrete solutions to (28)–(31) converge to the solutions to the corresponding 1D differential problems. Hence, due to the periodic boundary conditions in k and in u, provided by the map swap, the eventual 2D solution to (28)–(31) converges to the solution to the unsplit linearised differential problem. The issue
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
8473
of convergence of the numerical solution to the original nonlinear problem (1) will be addressed in the next section (see Experiment 3). Aside from being dissipative, the constructed schemes are also balanced, according to (4). Specifically, multiplying the left-hand side of (24) by sR2 DkDu cos ul , we have
nþ1 nþ1 1 nþ1 T kþ12 mkþ1 þ T k 2 þ mk T k12 mk1 sR2 DkDu cos ul
s
nþ1 nþ1 ¼ T kþ12 Dkþ1=2 þ T k12 Dk1=2 sR2 DkDu cos ul
1 2R2 cos2 ul Dk2
þ
nþ1 T k 2 R 2 Dk D
1 2R
2
cos2
nþ1
2
ul Dk
þ T k 2 ðDkþ1=2 þ Dk1=2 ÞsR2 DkDu cos ul
u cos ul :
ð33Þ
Summing all over the k; l’s and using the periodicity of the solution in k, we find the terms with Dkþ1=2 and Dk1=2 vanish, and P P nþ1 there only remains R2 DkDu l cos ul k T kl 2 . Doing the same with the right-hand side of (24), we shall eventually obtain
R2 DkDu
X
cos ul
l
X
nþ1 T kl 2 T nkl ¼
k
s 2
R2 DkDu
X
cos ul
X nþ1 fkl 2 ;
l
ð34Þ
k
which is the discrete analogue of the balance Eq. (4). Calculations for (26) and (27) are similar. Under f ¼ 0 we get the mass conservation law (cf. (5))
R2 DkDu
X
cos ul
l
X nþ1 X X T kl 2 ¼ R2 DkDu cos ul T nkl : k
l
ð35Þ
k
It is to emphasise that all the constructed finite difference schemes are systems of linear algebraic equations with tridiagonal positive definite matrices. Hence, the solution can easily be computed by the Sherman–Morrison formula [29]. 2.3. Fourth-order finite difference schemes Due to the map swap (11) and (12) and periodic boundary conditions in k and in u we can involve, generally speaking, arbitrary order finite difference approximations when discretising the differential operators Ak ; Au in (9) and (10). In particular, applying the standard fourth-order stencils [30] to
@ @T @D @T @2T ¼ D þD 2; @k @k @k @k @k 1 @ @T @D @T @2T ¼ D cos u D tan u þD ; cos u @ u @u @u @u @ u2
ð36Þ ð37Þ
we shall eventually obtain the fourth-order finite difference scheme in k nþ1
nþ1
nþ12
T kþ22 mkþ2 T kþ12 mkþ1 þ T k
1
s
þ mk
¼ T nkþ2 mkþ2 þ T nkþ1 mkþ1 þ T nk
1
s
nþ1
nþ1
T k12 mk1 þ T k22 mk2
mk
nþ1
þ T nk1 mk1 T nk2 mk2 þ
fk 2 ; 2
where
Fig. 3. Experiment 1: graph of the L2 -norm of the second-order solution.
ð38Þ
8474
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 4. Experiment 1: second-order numerical solution at several time moments.
8475
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
mk ¼ mkþp mkþp Mk ¼
30Dk
; 24R cos2 ul Dk2 1 16Dk 8M k for p ¼ 1; ¼ 2 þ sgnðpÞ 24Dk R cos2 ul 24Dk2 1 Dk Mk ¼ 2 þ sgnðpÞ for p ¼ 2; 24Dk R cos2 ul 24Dk2 2
ð39Þ
Dkþ2 þ 8Dkþ1 8Dk1 þ Dk2 ; 12Dk
and in u nþ1 nþ1 T nþ1 lþ2 mlþ2 T lþ1 mlþ1 þ T l
¼
nþ1 T lþ22 mlþ2
þ
nþ1 T lþ12 mlþ1
1
þ
s
þ ml
nþ1 Tl 2
1
s
nþ1 T nþ1 l1 ml1 þ T l2 ml2
ml
þ
nþ1 T l12 ml1
nþ1 T l22 ml2
nþ12
þ
fl
2
ð40Þ
;
where
30Dl ; 24R2 Du2 1 16Dl 8M l for p ¼ 1; mlþp ¼ 2 þ sgnðpÞ 24Du R 24Du2 1 Dl Ml for p ¼ 2; þ sgnðpÞ mlþp ¼ 2 24Du R 24Du2 Dlþ2 þ 8Dlþ1 8Dl1 þ Dl2 Ml ¼ Dl tan ul : 12Du
ml ¼
ð41Þ
Using the bicyclic splitting (28)–(31), we also achieve the second approximation order in time. Analysis of the properties of the fourth-order schemes in space is more cumbersome than that done above for the second order. 3. Numerical experiments
To verify the developed schemes we have simulated a few diffusion phenomena. We used grids 2 2 ; 4 4 and 6 6 ; below, if not explicitly mentioned, the latter grid is assumed. Throughout the experiments we employed the Strang splitting that provides second-order accuracy in time.
3.1. Experiment 1—The simplest diffusion problem First of all we considered the simplest case—linear diffusion ða ¼ 0Þ with a constant factor l in the absence of sources (f ¼ 0). As the initial condition we took a smooth function in the shape of a spot located at k ¼ p6 ; u ¼ p4 . The diameter of
Fig. 5. Experiment 2: profile of the diffusion amplification factor (42) (North pole view).
8476
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 6. Experiment 2: fourth-order numerical solution at several time moments (North pole view).
the spot was chosen such that the latter did not cover the North pole. The aim of this experiment was to demonstrate the correctness of the method in the sense of accurate simulation of the diffusion process near the pole: since the initial condition is symmetric and no external forcing is present, it is evident to expect a symmetric spread of the spot over the sphere.
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
8477
This would prove that the presence of the pole singularities caused by the convergence of meridians in the lat-log coordinate system does not affect the diffusion process near the poles and that the splitting on the sphere and the map swap (11) and (12) are mathematically correct. ð1Þ In Fig. 3 we plot the temporal graph of the L2 SDk;Du -norm of the second-order numerical solution, while in Fig. 4 we show the numerical solution at a few time moments (the exact solution was expanded into a series of spherical harmonics [30]; it looked quite similar to the numerical solution, so we omit it). As it is seen, the solution’s L2 -norm is monotonically decaying, according to (22) and (23), whereas the spot is propagating uniformly all over the sphere and no any nonphysical effects (e.g., the mass accumulation near the pole, etc.) are observed. 3.2. Experiment 2—Variable diffusion coefficient Consider a more complicated problem. Let
lðk; uÞ sin4 2k sin2 u;
ð42Þ
while the initial condition be the same spot placed exactly on the North pole. Since all the directions from the North pole are equivalent, the spot is expected to spread according to the diffusion amplification factor’s profile (Fig. 5). In [20] a similar linear problem was concerned briefly, while in the current paper the process is complicated by introducing a nonlinearity with a ¼ 1. The fourth-order (in space) numerical solution is shown in Fig. 6 (we used a nonuniform colour map for better visualisation). One can see that the initial spot is propagating as it should be doing, consistent with the diffusion coefficient’s profile (42). The temporal graph of the solution’s L2 -norm shown in Fig. 7 demonstrates the property of the solution’s dissipation. 3.3. Experiment 3—Large gradients near the poles This experiment was performed to test the schemes on a function with large gradients in order to make sure that the developed method allows accurate modelling of phenomena which imply nonsmooth solutions. The schemes were tested on a series of grids 6 6 ; 4 4 and 2 2 . The numerical solution was compared with the analytical one using the relative error
kT num T exact k d dðt n Þ ¼
ð1Þ
L2 SDk;Du
n
kT exact k L2
ð1Þ SDk;Du
:
ð43Þ
As the exact solution we chose the function
Tðk; u; tÞ ¼ c1 sin n cos u cos2 t þ c2 ;
ð44Þ
n xðk #1 tan j1 uÞ þ #2 tan j2 u sin ct:
ð45Þ
where
The source function resulted to be
" ! ! 2 !# 2 @T lT a1 1 @T @2T cos u @2T @T @T þa ; f ðk; u; tÞ ¼ þT 2 þ tan u T a @t R cos u R cos u @k @ u2 @u @u R @k
ð46Þ
where
@T ¼ c1 cos u #2 c cos n tan j2 u cos ct cos2 t sin n sin 2t ; @t @T ¼ c1 x cos u cos2 t cos n; @k @2T ¼ c1 x2 cos u cos2 t sin n; @k2 @T @n ¼ c1 cos2 t cos u cos n sin n sin u ; @u @u " 2 ! 2 @ T @n @n 2 ¼ c cos t 2 sin u cos n sin n cos u 1 þ 1 @ u2 @u @u sin j1 u sin j2 u 2 ; þ # j sin c t þ 2 cos n cos u x#1 j21 2 2 cos3 j1 u cos3 j2 u @n #1 j1 #2 j2 ¼ x þ sin ct : @u cos2 j1 u cos2 j2 u
ð47Þ
8478
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 7. Experiment 2: graph of the L2 -norm of the fourth-order solution.
Fig. 8. Experiment 3: graphs of the function sðuÞ ¼ sinð#2 tan j2 uÞ at #2 ¼ 15; j2 ¼ 0:1; 0:3; 0:5; 0:7; 0:9.
Fig. 9. Experiment 3: graph of the relative error (grid 2 2 ).
Solution (44) was chosen because its spatial structures convey spiral spherical waves interesting for a variety of applications [31–33]. A smooth solution of a similar form was studied in [20], while now we substantially complicate the problem introducing in (45) the term sin n with n #2 tan j2 u. That term has the intention to simulate rather large gradients in the
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
8479
Fig. 10. Experiment 3: numerical solution at several time moments.
solution at high latitudes. Examples of such a behaviour at #2 ¼ 15 and j2 ¼ 0:1; 0:3; 0:5; 0:7; 0:9 are given in Fig. 8. As it is seen, from j2 ¼ 0:1 to 0.9 the solution is getting more oscillating, with sharp gradients. In our experiment we assumed c1 ¼ 10; c2 ¼ 100; x ¼ 7; #1 ¼ 1; j1 ¼ 0:5; #2 ¼ 15; j2 ¼ 0:7; c ¼ 4; l ¼ 1 107 ; a ¼ 2.
8480
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 11. Experiment 4: graphs of the ‘‘mass’’ (left) and the ‘‘integrated sources’’ (right) (a ¼ 3).
In Fig. 9 we plot the graph of dn in time on the grid 2 2 . According to the periodical character of the analytical solution (44) in time, the relative error is concentrated within a bounded range, no trends towards unlimited growth are observed. From (46) and (47) it follows that the forcing generates a periodic wavy solution in k and in u due to the terms sin n; cos n and various combinations thereof; yet, the wave’s phase is being changed in time due to the terms sin ct; cos ct. This is ex actly what is observed in Fig. 10 where the numerical solution on the grid 6 6 at several time moments is shown. One can see the solution has a spiral structure whose direction of rotation is being varied in time with the period 2cp (cf., e.g., t ¼ 0:4 vs. t ¼ 2:0; t ¼ 0:7 vs. t ¼ 2:3; etc.). Grid refinement from 6 6 to 4 4 and 2 2 provides dnmax ¼ 4:599 102 ; 2:005 102 3 and 6:458 10 , respectively, which indicates that the solution to the linearised split 2D finite difference problem converges to the solution to the original unsplit nonlinear differential Eq. (1). 3.4. Experiment 4—The balance property and schemes’ accuracy at a high nonlinearity This experiment had the aim to test the schemes to be balanced in the sense of (34) when dealing with a highly nonlinear diffusion process. Besides, along with the second-order schemes we tested the fourth-order ones which are expected to provide more accurate results. We considered a series of diffusion problems with different a’s and l’s. The exact (analytical) solution had the form
Tðk; u; tÞ ¼ c1 sin n cos u cos2 t þ c2 ;
ð48Þ
where
n xk þ # cos ju sin t:
ð49Þ 32a
We took c1 ¼ 2:5; c2 ¼ 100; x ¼ 9; #ðaÞ ¼ 5a; j ¼ 3; a ¼ 0; 1; 2; 3; lðaÞ ¼ 1 10 . Unlike (45), the solution is now smooth due to n cos ju instead of n tan j2 u. However, the dependence # ¼ #ðaÞ results in a forcing f that produces two totally different behaviours of the solution T: the case of a ¼ 0 corresponds to the linear diffusion equation and a wave solution determined by the steady-state term sin xk, while the case a – 0 corresponds to various nonlinear diffusion problems with nonstationary wave solutions of the form sinðxk þ 5a cos ju sin tÞ. P P nþ1 ð1Þ nþ1 ð1Þ In Fig. 11 there are the graphs of the ‘‘mass’’ T nR ¼ k;l T nkl DSDk;Du and the ‘‘integrated sources’’ fR 2 ¼ k;l fkl 2 DSDk;Du in time for a ¼ 3. The ‘‘mass’’ looks as a decaying function of time almost everywhere, except, roughly, t 2 ½1:3; 1:7 and t 2 ½4:5; 4:9 , where it is nearly constant. This perfectly matches the behaviour of the ‘‘integrated sources’’ (cf. (34)). In Fig. 12 we show the fourth-order numerical solution at a ¼ 3 (note that unlike the previous experiment the periodicity is now equal to 2p, not p2 , due to sin t in (49)). The solution is behaving in accordance with the term sin n: as the time grows, the sinusoid begins travelling to the West (clockwise) at the low latitudes, while at the high ones it is going to the East (anticlockwise); later, when the time comes to t ¼ p2 , the direction of rotation is changed to the opposite one, etc. Numerical solutions related to different a’s are shown in Fig. 13. As it is seen, the higher the degree of nonlinearity, the more complicated structure of the solution, especially when the phase # cos ju sin t in (49) achieves the maximum at t p2 and 32p. In fact, in the strongly nonlinear cases (a ¼ 2; 3) the interaction between the three basic mechanisms—the nonlinearity, forcing and dissipation—generates solutions with periodic breakups and renewals of complex structures via selforganisation. In Table 1 we summarise the maximum relative errors dnmax ’s, corresponding to different values of a. As expected, the fourth-order schemes yield more precise solutions.
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
8481
Fig. 12. Experiment 4: fourth-order numerical solution at several time moments (a ¼ 3).
3.5. Experiment 5—Modelling of combustion The phenomenon of combustion was studied in [8] on the real line as a Cauchy problem described by the nonlinear diffusion equation with an appropriate parameter a and a source function f. We tested the developed method having simulated this phenomenon within a bounded area on a sphere. This problem has an interesting application in metal flaming and burning under the laser action [34]. For the initial condition we took a spot with the centre at k ¼ u ¼ p4 . We also had l ¼ 1 103 ; a ¼ 1 and f ¼ cT b , where c ¼ 4:1. The parameter b > 0 introduces an additional, aside from a, nonlinearity, now with respect to the sources, and determines two different regimes of combustion: the case b < a þ 1 corresponds to the HS-regime—the area of combustion is getting expanded while burning; the case b > a þ 1 corresponds to the LS-regime—the area of combustion is getting narrow [8].
8482
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 13. Experiment 4: fourth-order numerical solutions at several time moments: a ¼ 0 (top), a ¼ 1 (middle), a ¼ 2 (bottom).
Table 1 Experiment 4: maximum relative errors dnmax ’s at different a’s.
a
dnmax 2nd order
4th order
0
4:650 103
6:587 104
1
3
3:639 10
5:397 104
2
1:314 102
7:939 103
3
2
2:591 102
3:496 10
The point is that in both cases the given source leads to an infinite increase of the temperature T while the time grows, that is a blow-up occurs [8]. In Fig. 14 we present the numerical solutions related to b ¼ 1 (HS-regime) and b ¼ 4 (LS-regime). In Fig. 15 we also plot the L2 -norms of the second-order numerical solutions. The experiment confirms the theory: in the HS-regime we observe the increase of the area of combustion, while in the LS-regime the area is reducing to a singularity. The graphs of the energy prove the solutions tend to infinity, that is a blow-up happens. Further analysis allows to conclude that in the HS-regime the entire character of the blow-up asymptotics is slower than that in the LS-regime—the curve of the solution’s L2 -norm at b ¼ 1 is growing bit by bit all the time from t ¼ 0 to t ¼ 1, while at b ¼ 4 it is going up slowly only up to t 0:1, and
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
8483
Fig. 14. Experiment 5: second-order numerical solutions at several time moments for HS-regime (first two rows) and LS-regime (last two rows).
it rapidly changes the growth rate after that. In its turn, the growth of the combustion area observed in Fig. 14 is smooth at b ¼ 1 all the time, whereas it is much sharper at b ¼ 4 after t 0:1 (cf., e.g., t ¼ 0:09 and t ¼ 0:11).
8484
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 15. Experiment 5: graphs of the L2 -norm of the second-order solutions in time for HS-regime (left) and LS-regime (right).
4. Conclusions A new method for the numerical solution of nonlinear diffusion equations on a sphere has been developed. Due to the splitting of the differential operator by coordinates and the use of different coordinate maps for the sphere at two split time intervals, both one-dimensional split problems obey periodic boundary conditions. Thus, although the sphere is not a doubly periodic manifold, our method allows to construct a numerical algorithm for the sphere in the same way as for a torus. Hence, the method does not require to impose any artificial boundary conditions at the poles (which, in addition, must be physically and mathematically correct)—and it must be so indeed, since the sphere has no boundaries and preferred points, and so all directions on the sphere are equivalent. This resulted in constructing second- and fourth-order finite difference schemes in space without using computationally cumbersome numerical procedures for computing the numerical solution at the poles. The schemes keep the properties required by the original differential problem: they appear balanced and dissipative, thus providing physically adequate numerical solutions. The schemes are also computationally inexpensive and can be solved by direct band linear solvers. The numerical experiments proved the approach to allow accurate simulating diverse diffusion phenomena with constant and variable diffusion coefficients on the entire sphere, including those that imply high nonlinearity, nonsmoothness and infinite growth of the solutions. The accuracy of the method was evaluated by comparing the numerical results with the analytical solutions obtained under forcings of special types. The convergence of the numerical solution to the analytics was verified on a series of spatial grids 6 6 ; 4 4 and 2 2 . As an interesting nonlinear effect, it was also constructed an example showing that the competition of the three basic mechanisms of the model—the nonlinear interaction, external forcing and dissipation—generates a wave solution, whose spatial structure varies periodically due to the alternating influence of the processes of self-organisation and self-destruction. Acknowledgements This research was partially supported by grants No. 14539 and 26073 of the National System of Researchers of Mexico (SNI), and is part of the project PAPIIT-UNAM IN104811, Mexico. The authors express gratitude to the anonymous reviewer for the remarks resulted in improvement of the paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
J. Bear, Dynamics of Fluids in Porous Media, Courier Dover Publications, 1988. F. Catté et al, Image selective smoothing and edge detection by nonlinear diffusion, SIAM J. Numer. Anal. 29 (1992) 182–193. M.E. Glicksman, Diffusion in Solids: Field Theory, Solid-State Principles and Applications, John Wiley & Sons, 2000. J.R. King, ‘Instantaneous source’ solutions to a singular nonlinear diffusion equation, J. Eng. Math. 27 (1993) 31–72. A.A. Lacey, J.R. Ockerdon, A.B. Tayler, ‘Waiting-time’ solutions of a nonlinear diffusion equation, SIAM J. Appl. Math. 42 (1982) 1252–1264. G.A. Rudykh, E.I. Semenov, Non-self-similar solutions of multidimensional nonlinear diffusion equations, Math. Notes 67 (2000) 200–206. A.Kh. Vorob’yov, Diffusion Problems in Chemical Kinetics, Moscow University Press, 2003 (in Russian). A.A. Samarskii et al, Blow-up in Quasilinear Parabolic Equations, Walter de Gruyter, 1995. Yu.N. Skiba, Dual oil concentration estimates in ecologically sensitive zones, Environ. Monit. Assess. 43 (1996) 139–151. V.I. Agoshkov, F. Saleri, Recent developments in the numerical simulation of shallow water equations Part. III: boundary conditions and finite element approximations in the river flow calculations, Math. Model. 8 (1996) 3–24. X. Zheng, P. Palffy-Muhoray, Nonlinear diffusion on a sphere, Bull. Am. Phys. Soc. 57 (2012). W.-Z. Huang, D.M. Sloan, Pole condition for singular problems: the pseudospectral approximation, J. Comput. Phys. 107 (1993) 254–261. X.-Y. Chen et al, Stabilization of vortices in the Ginzburg–Landau equation with a variable diffusion coefficient, SIAM J. Math. Anal. 29 (1998) 903–912. M.-C. Lai, A note on finite difference discretizations for Poisson equation on a disk, Numer. Methods Partial Diff. Eq. 17 (2001) 199–203. M.-C. Lai, Y.-H. Tseng, A fast iterative solver for the variable coefficient diffusion equation on a disk, J. Comput. Phys. 208 (2005) 196–205.
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485 [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]
8485
P. Gérard, F. Méhats, The Schrödinger–Poisson system on the sphere, SIAM J. Math. Anal. 43 (2011) 1232–1268. W. Hundsdorfer, J.G. Verwer, Numerical Solution of Time-Dependent Advection–Diffusion–Reaction Equations, Springer, 2003. S.V. Tsynkov, Numerical solution of problems on unbounded domains. A review, Appl. Numer. Math. 27 (1998) 465–532. Z. Wu et al, Nonlinear Diffusion Equations, World Scientific Publishing, Singapore, 2001. Yu.N. Skiba, D.M. Filatov, Simulation of nonlinear diffusion on a sphere, in: Proceedings of the 1st International Conference on Simulation and Modelling Methodologies, Technologies and Applications (SIMULTECH 2011), Noordwijkerhout, the Netherlands, 2011, pp. 24–30. D.W. Peaceman, H.H. Rachford Jr., The numerical solution of parabolic and elliptic differential equations, J. Soc. Ind. Appl. Math. 3 (1955) 28–41. N.N. Yanenko, The Method of Fractional Steps: Solution of Problems of Mathematical Physics in Several Variables, Springer, 1971 (translated from Russian: Nauka, Novosibirsk, 1967). G.I. Marchuk, Methods of Computational Mathematics, Springer, 1982 (translated from Russian: Nauka, Moscow, 1977). D.L. Williamson, Difference approximations of flow equations on the sphere, in: Numerical Methods Used in Atmospherical Models, vol. II. GARP Publ. Series, No 17, 1979. K. Mohseni, T. Colonius, Numerical treatment of polar coordinate singularities, J. Comput. Phys. 157 (2000) 787–795. Yu.N. Skiba, D.M. Filatov, On splitting-based mass and total energy conserving arbitrary order shallow-water schemes, Numer. Methods Partial Diff. Eq. 23 (2007) 534–552. Yu.N. Skiba, D.M. Filatov, Conservative arbitrary order finite difference schemes for shallow-water flows, J. Comput. Appl. Math. 218 (2008) 579–591. G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 506–517. W.H. Press et al, Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, Cambridge, 2007. G.A. Korn, T.M. Korn, Mathematical Handbook for Scientists and Engineers, Courier Dover Publications, 2000. V.S. Zykov, S.C. Müller, Spiral waves on circular and spherical domains of excitable medium, Physica D 97 (1996) 322–332. H. Yagisita, M. Mimura, M. Yamada, Spiral wave behaviors in an excitable reaction-diffusion system on a sphere, Physica D 124 (1998) 126–136. K. Rohlf, L. Glass, R. Kapral, Spiral wave dynamics in excitable media with spherical geometries, Chaos 16 (2006) 037115 (10 pp.). A.A. Samarskii, Nonlinear effects of blow-up and localization processes in burning problems, in: C.-M. Brauner, C. Schmidt-Lainé (Eds.), Mathematical Modeling in Combustion and Related Topics, Martinus Nijhoff Publishers, 1988, pp. 217–231.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Splitting-based schemes for numerical solution of nonlinear diffusion equations on a sphere Yuri N. Skiba a, Denis M. Filatov b,⇑ a Centro de Ciencias de la Atmósfera (CCA), Universidad Nacional Autónoma de México (UNAM), Av. Universidad #3000, Cd. Universitaria, C.P. 04510 México D.F., Mexico b Centro de Investigación en Computación (CIC), Instituto Politécnico Nacional (IPN), Av. Juan de Dios Bátiz s/n, C.P. 07738 México D.F., Mexico
a r t i c l e
i n f o
Keywords: Nonlinear diffusion equation Dissipative balanced finite difference schemes Operator splitting on closed nonperiodic manifolds
a b s t r a c t We provide an advanced study of our recently developed method for the numerical solution of nonlinear diffusion equations on a sphere. In particular, we analyse the method in detail when applied to solving diverse diffusion phenomena, with specific conditions on the smoothness of the solution, the degree of nonlinearity, and the initial data and sources. The main idea of the method consists in splitting the original differential operator by coordinates and subsequent constructing finite difference schemes for the split one-dimensional problems using different coordinate maps for the sphere at the two split time intervals. The essential advantage of this technique is that each split 1D equation can be equipped with a periodic boundary condition, despite the sphere being not a doubly periodic domain. Therefore, unlike the existing methods, this one does not require applying special numerical procedures for careful computing the solution near the poles, which is always a challenge. Each split 1D equation is approximated by a second- or a fourth-order finite difference scheme that keeps all the substantial properties of the differential problem: it is balanced and dissipative, since the spatial finite difference operator is negative definite. The developed algorithm is cheap-to-implement from the computational point of view. The theoretical results are confirmed numerically by simulating various nonlinear diffusion processes. A numerical example, in particular, shows that the competition of the three basic mechanisms—the nonlinear interaction, forcing and dissipation—can generate wave solutions, whose spatial structures on the sphere are subjected to the alternating influence of the processes of self-organisation and self-destruction. The accuracy of the method is evaluated by comparing the numerical solutions versus the analytical ones obtained with specially chosen forcings. The convergence of the numerical solution to the analytics is verified by refining spatial grids. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction A large number of important natural phenomena, e.g., heat transfer in ionised gases, gas percolation through porous media, concentration waves in distributed chemical reactors, combustion, viscous processes in diverse media, as well as many others are described by (nonlinear) diffusion equations [1–10]. Possibly, the most obvious and well-studied applications of the diffusion equation lie in atmospheric and hydrodynamic sciences. Hundreds of papers providing researches on some or other issues of the environment are published every year. ⇑ Corresponding author. E-mail addresses: [email protected] (Y.N. Skiba), denisfi[email protected] (D.M. Filatov). 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.02.066
8468
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Among the applications there are air and water pollution, viscous flows in river channels—just to mention a couple of them [9,10]. The governing equations are specific for each application, but in general they contain a hyperbolic part responsible for the advection process (velocity field) and a parabolic term aimed to simulate fluctuation from the mainstream, or diffusion. Another important application of the diffusion equation, especially involving nonlinearity, is blow-up. In [8], among other issues, a nonlinear diffusion Cauchy problem in Rn is studied and numerical experiments are performed, simulating blow-up in a bounded domain. A not less interesting, though a little less obvious problem described by the diffusion equation is the simulation of time evolution of the orientational distribution function of non-spherical particles due to interactions with each other and external fields [11]. In many papers on diffusion equations boundary value problems are studied [12–15], while much less papers are dedicated to the numerical simulation of nonlinear diffusion processes on closed manifolds. The latters differ from bounded domains as having no boundaries, which makes it nontrivial to develop efficient numerical methods for solving diffusion equations in such domains. Although for some particular cases these problems can be reduced to boundary value problems with boundary conditions of special types (e.g., a torus can be represented as a doubly periodic domain, so that periodic boundary conditions in both directions can be employed), this cannot be done for a more or less complicated manifold. A simple, but important example of such a manifold is a two-dimensional sphere. The case of the sphere can be seen as a step before more general manifolds, which can be interesting, e.g., in nanoelectronics [16]. The point is that the sphere is not a doubly periodic domain, since it is periodic in the longitude, but is not in the latitude due to the presence of two poles. Therefore, a numerical procedure designed the solve the original 2D problem will be computationally cumbersome, because the matrix of the resulting linear system will be of a general type, hence not permitting to apply some or other fast linear solvers. For instance, in [11] a Voronoi tessellation of the sphere is used and then a finite volume method is applied. As for making a preliminary splitting of the 2D equation [17], this will bring to the necessity of constructing mathematically and physically correct boundary conditions or involving special numerical procedures near the poles (e.g., matrix bordering) when computing in the latitudinal direction. These are always a problem, as the poles represent an artificial boundary appearing exclusively due to the latitudinal-longitudinal coordinate system, and the construction of proper artificial boundary conditions is a serious independent question [18]. The nonlinear diffusion equation on a two-dimensional sphere S ¼ fðk; uÞ : k 2 ½0; 2pÞ; u 2 p2 ; p2 g has the form
@T 1 @ lT a @T @ lT a cos u @T þ þ f: ¼ AT þ f @t R cos u @k R cos u @k @u @u R
ð1Þ
It is equipped with a suitable initial condition. Here A is the diffusion operator, T ¼ Tðk; u; tÞ P 0 is the unknown function, lT a is the diffusion coefficient, l ¼ lðk; u; tÞ P 0 is the amplification factor, f ¼ f ðk; u; tÞ is the source function, R is the radius of the sphere with k as the longitude (positive eastward) and u as the latitude (positive northward). The parameter a is usually a positive integer number that determines the degree of nonlinearity of the diffusion process; if a ¼ 0 then we are dealing with linear diffusion. Introducing on S the scalar product as
hf; gi ¼
Z
fgdS;
ð2Þ
S
multiplying (1) by g and integrating over the domain S, we find
@T ; g ¼ hAT; gi þ hf ; gi: @t
ð3Þ
In the case of g ¼ 1 the first summand on the right-hand side of (3) vanishes [19] and the equality turns into the balance equation
d dt
Z
TdS ¼ S
Z
fdS;
ð4Þ
S
which, if f ¼ 0, provides the mass conservation law
d dt
Z
TdS ¼ 0:
ð5Þ
S
In the case of g ¼ T due to the negative definiteness of A it holds
hAT; gi 6 0;
ð6Þ
and so
1 d 2 dt
Z S
T 2 dS 6
Z
fTdS; S
from where under f ¼ 0 we obtain the solution’s dissipation in the L2 ðSÞ-norm
ð7Þ
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
1 d 2 dt
Z
T 2 dS 6 0:
8469
ð8Þ
S
Our aim is to develop an accurate and efficient numerical method for solving Eq. (1). The accuracy implies that the method should provide physically correct numerical solutions which would possess properties (4) and (8). The efficiency means that the method should be computationally inexpensive. First results on this issue were presented by the authors in [20], so in this paper we study the developed method in detail, analysing how it works when applied to solving diverse diffusion phenomena, with specific conditions on the smoothness of the solution, the degree of nonlinearity, and the initial data and sources. Besides, we explicitly demonstrate the constructed finite difference operators inherit important properties of the original diffusion operator. The paper is organised as follows. In Section 2 we first linearise the original nonlinear differential equation and then apply the method of splitting [21] to the resulting equation on the sphere. Employing the procedure of alternation of two coordinate maps for representing the same sphere at each split time step, we become able to develop a numerical algorithm as if we dealt with a periodic domain in both spatial directions (although, repeat again, the sphere is not a doubly periodic manifold). This leads to using simple periodic boundary conditions while computing the solutions in k and in u. Next, involving the Crank–Nicolson approximation, we construct two second-order finite difference schemes for the split 1D diffusion equations. The schemes are shown to be dissipative and balanced, according to the properties of the original differential problem, as well as easy-to-implement from the computational standpoint. Then, because of the periodicity of the boundary conditions we increase the approximation’s accuracy to the fourth order in space. In Section 3 we provide results of numerical experiments performed to test the developed method. In several examples, from the simplest linear case to a complex nonlinear one, we demonstrate the schemes are accurately simulating diffusion phenomena. We prove the schemes inherit the essential properties of the differential operators and hence yield physically adequate numerical solutions. We comprehensively study the developed method, providing a comparison of the numerical solutions with the analytics. We also give a numerical example (see Experiment 4) demonstrating the synergistic action of the external forcing, nonlinearity and dissipation. These three mechanisms produce solutions, whose spatial structures assert the presence of the processes of selforganisation and self-destruction, which phenomenon is of interest in nonlinear science. In Section 4 we give a conclusion. 2. Mathematical foundation of the method 2.1. Operator splitting on the sphere In the standard manner we define the grid spacing in time as s ¼ tnþ1 t n and in every time interval ðtn ; tnþ1 Þ we linearise (1) and then split the linearised equation by coordinates [22,23]
@T f 1 @ D @T f þ ; ¼ Ak T þ @t 2 R cos u @k R cos u @k 2 @T f 1 @ D cos u @T f þ ; ¼ Au T þ @t 2 R cos u @ u @u 2 R
ð9Þ ð10Þ
where D ¼ lðT n Þa ; T n ¼ Tðk; u; t n Þ. Hereafter the split equations will be considered in time successively: the solution to Eq. (9), solved in k, will be used as the initial condition for Eq. (10); the latter, solved in u, will supply the initial condition for (9) in the next time interval ðtnþ1 ; tnþ2 Þ, and so on. Because the metric term R cos u vanishes at u ¼ p=2, both equations are meaningless at the poles. Therefore, for covering the entire sphere we assume Dk ¼ kkþ1 kk and Du ¼ ulþ1 ul and define on S a grid, making a half step shift in u [24], i.e. ð1Þ
SDk;Du ¼
Dk Dk p Du p Du ; ul 2 þ : ; 2p þ ðkk ; ul Þ : kk 2 ; 2 2 2 2 2 2
ð11Þ
Thereby we exclude the pole singularities (Fig. 1), so that the subsequent finite difference equations will have sense everyð1Þ where on SDk;Du . For solving (9) one has, evidently, to use the periodic boundary condition in k, as the sphere S is a periodic domain in the longitude. (Note that a half step grid shifting was also applied in [14,25], and, compared to the classical approach [12], it provided benefits in the solution.) As for solving Eq. (10), one faces the difficulty since the sphere is not a periodic domain in u. Several approaches can be used to resolve it, e.g., one could involve a matrix bordering procedure or try to paste the solution from the opposite meridians at the poles. A disadvantage of these techniques is that in the case of an implicit temporal approximation on the righthand side of (10) the resulting matrix will be of a general type, and so no fast algorithms of linear algebra can be used for computing the solution. Alternatively, an explicit temporal discretisation can be used, but this will impose serious restrictions on the timestep, especially when a is large. However, this is exactly the method of splitting that allows to avoid these undesired procedures which, besides, if performed inaccurately, may easily result in introducing nonphysical modes into the solution. Specifically, due to the splitting, for computing the solution in u we change the coordinate map from (11) to
8470
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 1. Shift of the original grid (dotted lines) a half step in u (solid lines) allows excluding the pole singularities, thereby keeping the equations to have sense on the entire sphere. The selected black parallel on the sphere corresponds to that on the plane.
Dk Dk p Du 3p Du ;p ; ul 2 þ ðkk ; ul Þ : kk 2 ; þ : 2 2 2 2 2 2
ð2Þ
SDk;Du ¼
ð12Þ
Obviously, grid (12) contains the same nodes as (11) (Fig. 2). The use of the two coordinate maps allows employing the same numerical algorithm with the periodic boundary conditions in both directions, k and u [26,27,20]. The only change we have to make in (10) if using (12) is to replace cos u with j cos uj, as well. 2.2. Second-order finite difference schemes According to the splitting, in each time interval ðtn ; tnþ1 Þ we discretise the temporal derivatives in (9) and (10) as
nþ1 T kl 2 T nkl @T ; @t t¼tn s
ð13Þ
nþ1 T nþ1 T kl 2 @T kl ; @t t¼tnþ1=2 s
ð14Þ
and
respectively. In its turn, for the spatial derivatives we take the second-order central finite difference stencils—
Ak Tjk¼kk ADk Tjk¼kk
R2 cos2 u
Au T u¼u ADu T u¼u l
l
Dkþ1=2
1
1 R2 j cos ul j
T kþ1 T k Dk
k1 Dk1=2 T k T Dk ; Dk
ðDj cos ujÞlþ1=2
T lþ1 T l Du
l1 ðDj cos ujÞl1=2 T l T Du
Du
ð15Þ
;
ð16Þ
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
8471
Fig. 2. The map swap, possible due to the splitting, allows using the periodic boundary condition in u, avoiding the question of constructing boundary conditions at the poles. The selected black meridian on the sphere corresponds to that on the plane.
where
Dkþp=2 :¼
Dkþðpþ1Þ=2 þ Dkþðp1Þ=2 2
ðDj cos ujÞlþp=2 :¼
for p ¼ 1;
ðDj cos ujÞlþðpþ1Þ=2 þ ðDj cos ujÞlþðp1Þ=2 2
ð17Þ
for p ¼ 1:
ð18Þ
(In order not to overload the formulas, the nonvarying index l in (15) and the nonvarying index k in (16) are omitted.) The 1 external forcing f is computed at the intermediate time moment tn þt2nþ1 , i.e. f nþ2 ¼ f k; u; tn þt2nþ1 . The finite difference operators in (15) and (16) are correct in the sense of keeping the property of negative definiteness (cf. (6)). Indeed, for (15), multiplying the right-hand side by T k (recall, l is fixed while computing in k) and summing all over the k’s, we obtain
! k1 X X Dk1=2 T k T 1 Dk ¼ D ðT T ÞT D ðT T ÞT T k kþ1=2 kþ1 k k k1=2 k k1 k 2 2 Dk R2 cos2 ul Dk2 k R cos ul k k ! X X X 1 1 Dkþ1=2 ðT kþ1 T k ÞT k Dk0 þ1=2 ðT k0 þ1 T k0 ÞT k0 þ1 ¼ 2 Dkþ1=2 ðT kþ1 T k ÞðT k T kþ1 Þ ¼ 2 2 u Dk2 R cos2 ul Dk2 R cos 0 l k k k X 1 Dkþ1=2 ðT kþ1 T k Þ2 6 0: ð19Þ ¼ 2 2 2 R cos ul Dk k
X
Dkþ1=2
1
0
T kþ1 T k Dk
ð1Þ
Here we denoted k ¼ k 1 and used the periodicity of the solution in k on SDk;Du to return to the original index k. Analoð2Þ gously, for (16), due to the periodicity of the solution in u on SDk;Du we have
8472
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
1
l
R2 j cos ul j ¼
T lþ1 T l Du
ðDj cos ujÞlþ1=2
X
l1 ðDj cos ujÞl1=2 T l T Du
Du
T l j cos ul j ¼ . . .
X ðDj cos ujÞlþ1=2 ðT lþ1 T l Þ2 6 0:
1 2
R Du 2
ð20Þ
l
Temporal discretisation of the function T kl in (15) and (16) can be done via the Crank–Nicolson approximation nþp
nþp1 2
T 2 þ T kl T kl :¼ kl 2
for p ¼ 1; 2:
ð21Þ
Therefore, each of the constructed split one-dimensional second-order finite difference schemes is dissipative (cf. (7) and (8))—
1 2
nþ2
T k
ð1Þ
L2 SDk;u
L2
l
nþ1 2
T l
ð2Þ k ;D u
L2 Sk
* nþ1 + 2 þ 2s f k ; T k ð1Þ 2 SDk;u
6
T n 2 k
* nþ1 + 2 þ 2s fl ; T l ð2Þ Sk ;Du 2
1 2 nþ 6
T l 2 L2
l
ð2Þ k ;D u
L2 Sk
k
ð22Þ
;
ð1Þ
L2 SDk;u
l
ð23Þ
:
Gathering all together, eventually for the second-order finite difference scheme in k we obtain nþ1
nþ12
T kþ12 mkþ1 þ T k
1
nþ1 f 2 1 nþ1 þ mk T k12 mk1 ¼ T nkþ1 mkþ1 þ T nk mk þ T nk1 mk1 þ k ; s s 2
ð24Þ
where
mk ¼
Dkþ1=2 þ Dk1=2 2R
2
cos2
u l Dk
2
;
mkþp ¼
Dkþp=2 2
2R cos2 ul Dk2
for p ¼ 1;
ð25Þ
while in u we have
T nþ1 lþ1 mlþ1
þ
T nþ1 l
1
s
þ ml
T nþ1 l1 ml1
¼
nþ1 T lþ12 mlþ1
þ
nþ1 Tl 2
1
s
ml
þ
nþ1 T l12 ml1
nþ12
þ
fl
2
;
ð26Þ
where
ml ¼ mlþp
ðDj cos ujÞlþ1=2 þ ðDj cos ujÞl1=2 2R2 j cos ul jDu2 ðDj cos ujÞlþp=2 ¼ 2 for p ¼ 1: 2R j cos ul jDu2
; ð27Þ
Supplementing them with the bicyclic splitting (or Strang splitting) nþ1
nþ2
T kl 4 T nkl f 4 ¼ ADk T kl þ k ; s=2 2 nþ2
nþ1
nþ2
T kl 4 T kl 4 f 4 ¼ ADu T kl þ l ; s=2 2 nþ3
ð28Þ
nþ2
nþ2
nþ34
nþ24
T kl 4 T kl 4 f 4 ¼ ADu T kl þ l ; s=2 2 T nþ1 T kl f kl ¼ ADk T kl þ k ; s=2 2
ð29Þ ð30Þ ð31Þ
where nþ
p
nþ
T 4 þ T kl T kl :¼ kl 2
p1 4
for p ¼ 1; 4;
ð32Þ
we achieve the second order accuracy in time in the interval ðtn ; tnþ1 Þ [23,28]. Since both the approximation and stability are fulfilled, the discrete solutions to (28)–(31) converge to the solutions to the corresponding 1D differential problems. Hence, due to the periodic boundary conditions in k and in u, provided by the map swap, the eventual 2D solution to (28)–(31) converges to the solution to the unsplit linearised differential problem. The issue
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
8473
of convergence of the numerical solution to the original nonlinear problem (1) will be addressed in the next section (see Experiment 3). Aside from being dissipative, the constructed schemes are also balanced, according to (4). Specifically, multiplying the left-hand side of (24) by sR2 DkDu cos ul , we have
nþ1 nþ1 1 nþ1 T kþ12 mkþ1 þ T k 2 þ mk T k12 mk1 sR2 DkDu cos ul
s
nþ1 nþ1 ¼ T kþ12 Dkþ1=2 þ T k12 Dk1=2 sR2 DkDu cos ul
1 2R2 cos2 ul Dk2
þ
nþ1 T k 2 R 2 Dk D
1 2R
2
cos2
nþ1
2
ul Dk
þ T k 2 ðDkþ1=2 þ Dk1=2 ÞsR2 DkDu cos ul
u cos ul :
ð33Þ
Summing all over the k; l’s and using the periodicity of the solution in k, we find the terms with Dkþ1=2 and Dk1=2 vanish, and P P nþ1 there only remains R2 DkDu l cos ul k T kl 2 . Doing the same with the right-hand side of (24), we shall eventually obtain
R2 DkDu
X
cos ul
l
X
nþ1 T kl 2 T nkl ¼
k
s 2
R2 DkDu
X
cos ul
X nþ1 fkl 2 ;
l
ð34Þ
k
which is the discrete analogue of the balance Eq. (4). Calculations for (26) and (27) are similar. Under f ¼ 0 we get the mass conservation law (cf. (5))
R2 DkDu
X
cos ul
l
X nþ1 X X T kl 2 ¼ R2 DkDu cos ul T nkl : k
l
ð35Þ
k
It is to emphasise that all the constructed finite difference schemes are systems of linear algebraic equations with tridiagonal positive definite matrices. Hence, the solution can easily be computed by the Sherman–Morrison formula [29]. 2.3. Fourth-order finite difference schemes Due to the map swap (11) and (12) and periodic boundary conditions in k and in u we can involve, generally speaking, arbitrary order finite difference approximations when discretising the differential operators Ak ; Au in (9) and (10). In particular, applying the standard fourth-order stencils [30] to
@ @T @D @T @2T ¼ D þD 2; @k @k @k @k @k 1 @ @T @D @T @2T ¼ D cos u D tan u þD ; cos u @ u @u @u @u @ u2
ð36Þ ð37Þ
we shall eventually obtain the fourth-order finite difference scheme in k nþ1
nþ1
nþ12
T kþ22 mkþ2 T kþ12 mkþ1 þ T k
1
s
þ mk
¼ T nkþ2 mkþ2 þ T nkþ1 mkþ1 þ T nk
1
s
nþ1
nþ1
T k12 mk1 þ T k22 mk2
mk
nþ1
þ T nk1 mk1 T nk2 mk2 þ
fk 2 ; 2
where
Fig. 3. Experiment 1: graph of the L2 -norm of the second-order solution.
ð38Þ
8474
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 4. Experiment 1: second-order numerical solution at several time moments.
8475
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
mk ¼ mkþp mkþp Mk ¼
30Dk
; 24R cos2 ul Dk2 1 16Dk 8M k for p ¼ 1; ¼ 2 þ sgnðpÞ 24Dk R cos2 ul 24Dk2 1 Dk Mk ¼ 2 þ sgnðpÞ for p ¼ 2; 24Dk R cos2 ul 24Dk2 2
ð39Þ
Dkþ2 þ 8Dkþ1 8Dk1 þ Dk2 ; 12Dk
and in u nþ1 nþ1 T nþ1 lþ2 mlþ2 T lþ1 mlþ1 þ T l
¼
nþ1 T lþ22 mlþ2
þ
nþ1 T lþ12 mlþ1
1
þ
s
þ ml
nþ1 Tl 2
1
s
nþ1 T nþ1 l1 ml1 þ T l2 ml2
ml
þ
nþ1 T l12 ml1
nþ1 T l22 ml2
nþ12
þ
fl
2
ð40Þ
;
where
30Dl ; 24R2 Du2 1 16Dl 8M l for p ¼ 1; mlþp ¼ 2 þ sgnðpÞ 24Du R 24Du2 1 Dl Ml for p ¼ 2; þ sgnðpÞ mlþp ¼ 2 24Du R 24Du2 Dlþ2 þ 8Dlþ1 8Dl1 þ Dl2 Ml ¼ Dl tan ul : 12Du
ml ¼
ð41Þ
Using the bicyclic splitting (28)–(31), we also achieve the second approximation order in time. Analysis of the properties of the fourth-order schemes in space is more cumbersome than that done above for the second order. 3. Numerical experiments
To verify the developed schemes we have simulated a few diffusion phenomena. We used grids 2 2 ; 4 4 and 6 6 ; below, if not explicitly mentioned, the latter grid is assumed. Throughout the experiments we employed the Strang splitting that provides second-order accuracy in time.
3.1. Experiment 1—The simplest diffusion problem First of all we considered the simplest case—linear diffusion ða ¼ 0Þ with a constant factor l in the absence of sources (f ¼ 0). As the initial condition we took a smooth function in the shape of a spot located at k ¼ p6 ; u ¼ p4 . The diameter of
Fig. 5. Experiment 2: profile of the diffusion amplification factor (42) (North pole view).
8476
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 6. Experiment 2: fourth-order numerical solution at several time moments (North pole view).
the spot was chosen such that the latter did not cover the North pole. The aim of this experiment was to demonstrate the correctness of the method in the sense of accurate simulation of the diffusion process near the pole: since the initial condition is symmetric and no external forcing is present, it is evident to expect a symmetric spread of the spot over the sphere.
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
8477
This would prove that the presence of the pole singularities caused by the convergence of meridians in the lat-log coordinate system does not affect the diffusion process near the poles and that the splitting on the sphere and the map swap (11) and (12) are mathematically correct. ð1Þ In Fig. 3 we plot the temporal graph of the L2 SDk;Du -norm of the second-order numerical solution, while in Fig. 4 we show the numerical solution at a few time moments (the exact solution was expanded into a series of spherical harmonics [30]; it looked quite similar to the numerical solution, so we omit it). As it is seen, the solution’s L2 -norm is monotonically decaying, according to (22) and (23), whereas the spot is propagating uniformly all over the sphere and no any nonphysical effects (e.g., the mass accumulation near the pole, etc.) are observed. 3.2. Experiment 2—Variable diffusion coefficient Consider a more complicated problem. Let
lðk; uÞ sin4 2k sin2 u;
ð42Þ
while the initial condition be the same spot placed exactly on the North pole. Since all the directions from the North pole are equivalent, the spot is expected to spread according to the diffusion amplification factor’s profile (Fig. 5). In [20] a similar linear problem was concerned briefly, while in the current paper the process is complicated by introducing a nonlinearity with a ¼ 1. The fourth-order (in space) numerical solution is shown in Fig. 6 (we used a nonuniform colour map for better visualisation). One can see that the initial spot is propagating as it should be doing, consistent with the diffusion coefficient’s profile (42). The temporal graph of the solution’s L2 -norm shown in Fig. 7 demonstrates the property of the solution’s dissipation. 3.3. Experiment 3—Large gradients near the poles This experiment was performed to test the schemes on a function with large gradients in order to make sure that the developed method allows accurate modelling of phenomena which imply nonsmooth solutions. The schemes were tested on a series of grids 6 6 ; 4 4 and 2 2 . The numerical solution was compared with the analytical one using the relative error
kT num T exact k d dðt n Þ ¼
ð1Þ
L2 SDk;Du
n
kT exact k L2
ð1Þ SDk;Du
:
ð43Þ
As the exact solution we chose the function
Tðk; u; tÞ ¼ c1 sin n cos u cos2 t þ c2 ;
ð44Þ
n xðk #1 tan j1 uÞ þ #2 tan j2 u sin ct:
ð45Þ
where
The source function resulted to be
" ! ! 2 !# 2 @T lT a1 1 @T @2T cos u @2T @T @T þa ; f ðk; u; tÞ ¼ þT 2 þ tan u T a @t R cos u R cos u @k @ u2 @u @u R @k
ð46Þ
where
@T ¼ c1 cos u #2 c cos n tan j2 u cos ct cos2 t sin n sin 2t ; @t @T ¼ c1 x cos u cos2 t cos n; @k @2T ¼ c1 x2 cos u cos2 t sin n; @k2 @T @n ¼ c1 cos2 t cos u cos n sin n sin u ; @u @u " 2 ! 2 @ T @n @n 2 ¼ c cos t 2 sin u cos n sin n cos u 1 þ 1 @ u2 @u @u sin j1 u sin j2 u 2 ; þ # j sin c t þ 2 cos n cos u x#1 j21 2 2 cos3 j1 u cos3 j2 u @n #1 j1 #2 j2 ¼ x þ sin ct : @u cos2 j1 u cos2 j2 u
ð47Þ
8478
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 7. Experiment 2: graph of the L2 -norm of the fourth-order solution.
Fig. 8. Experiment 3: graphs of the function sðuÞ ¼ sinð#2 tan j2 uÞ at #2 ¼ 15; j2 ¼ 0:1; 0:3; 0:5; 0:7; 0:9.
Fig. 9. Experiment 3: graph of the relative error (grid 2 2 ).
Solution (44) was chosen because its spatial structures convey spiral spherical waves interesting for a variety of applications [31–33]. A smooth solution of a similar form was studied in [20], while now we substantially complicate the problem introducing in (45) the term sin n with n #2 tan j2 u. That term has the intention to simulate rather large gradients in the
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
8479
Fig. 10. Experiment 3: numerical solution at several time moments.
solution at high latitudes. Examples of such a behaviour at #2 ¼ 15 and j2 ¼ 0:1; 0:3; 0:5; 0:7; 0:9 are given in Fig. 8. As it is seen, from j2 ¼ 0:1 to 0.9 the solution is getting more oscillating, with sharp gradients. In our experiment we assumed c1 ¼ 10; c2 ¼ 100; x ¼ 7; #1 ¼ 1; j1 ¼ 0:5; #2 ¼ 15; j2 ¼ 0:7; c ¼ 4; l ¼ 1 107 ; a ¼ 2.
8480
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 11. Experiment 4: graphs of the ‘‘mass’’ (left) and the ‘‘integrated sources’’ (right) (a ¼ 3).
In Fig. 9 we plot the graph of dn in time on the grid 2 2 . According to the periodical character of the analytical solution (44) in time, the relative error is concentrated within a bounded range, no trends towards unlimited growth are observed. From (46) and (47) it follows that the forcing generates a periodic wavy solution in k and in u due to the terms sin n; cos n and various combinations thereof; yet, the wave’s phase is being changed in time due to the terms sin ct; cos ct. This is ex actly what is observed in Fig. 10 where the numerical solution on the grid 6 6 at several time moments is shown. One can see the solution has a spiral structure whose direction of rotation is being varied in time with the period 2cp (cf., e.g., t ¼ 0:4 vs. t ¼ 2:0; t ¼ 0:7 vs. t ¼ 2:3; etc.). Grid refinement from 6 6 to 4 4 and 2 2 provides dnmax ¼ 4:599 102 ; 2:005 102 3 and 6:458 10 , respectively, which indicates that the solution to the linearised split 2D finite difference problem converges to the solution to the original unsplit nonlinear differential Eq. (1). 3.4. Experiment 4—The balance property and schemes’ accuracy at a high nonlinearity This experiment had the aim to test the schemes to be balanced in the sense of (34) when dealing with a highly nonlinear diffusion process. Besides, along with the second-order schemes we tested the fourth-order ones which are expected to provide more accurate results. We considered a series of diffusion problems with different a’s and l’s. The exact (analytical) solution had the form
Tðk; u; tÞ ¼ c1 sin n cos u cos2 t þ c2 ;
ð48Þ
where
n xk þ # cos ju sin t:
ð49Þ 32a
We took c1 ¼ 2:5; c2 ¼ 100; x ¼ 9; #ðaÞ ¼ 5a; j ¼ 3; a ¼ 0; 1; 2; 3; lðaÞ ¼ 1 10 . Unlike (45), the solution is now smooth due to n cos ju instead of n tan j2 u. However, the dependence # ¼ #ðaÞ results in a forcing f that produces two totally different behaviours of the solution T: the case of a ¼ 0 corresponds to the linear diffusion equation and a wave solution determined by the steady-state term sin xk, while the case a – 0 corresponds to various nonlinear diffusion problems with nonstationary wave solutions of the form sinðxk þ 5a cos ju sin tÞ. P P nþ1 ð1Þ nþ1 ð1Þ In Fig. 11 there are the graphs of the ‘‘mass’’ T nR ¼ k;l T nkl DSDk;Du and the ‘‘integrated sources’’ fR 2 ¼ k;l fkl 2 DSDk;Du in time for a ¼ 3. The ‘‘mass’’ looks as a decaying function of time almost everywhere, except, roughly, t 2 ½1:3; 1:7 and t 2 ½4:5; 4:9 , where it is nearly constant. This perfectly matches the behaviour of the ‘‘integrated sources’’ (cf. (34)). In Fig. 12 we show the fourth-order numerical solution at a ¼ 3 (note that unlike the previous experiment the periodicity is now equal to 2p, not p2 , due to sin t in (49)). The solution is behaving in accordance with the term sin n: as the time grows, the sinusoid begins travelling to the West (clockwise) at the low latitudes, while at the high ones it is going to the East (anticlockwise); later, when the time comes to t ¼ p2 , the direction of rotation is changed to the opposite one, etc. Numerical solutions related to different a’s are shown in Fig. 13. As it is seen, the higher the degree of nonlinearity, the more complicated structure of the solution, especially when the phase # cos ju sin t in (49) achieves the maximum at t p2 and 32p. In fact, in the strongly nonlinear cases (a ¼ 2; 3) the interaction between the three basic mechanisms—the nonlinearity, forcing and dissipation—generates solutions with periodic breakups and renewals of complex structures via selforganisation. In Table 1 we summarise the maximum relative errors dnmax ’s, corresponding to different values of a. As expected, the fourth-order schemes yield more precise solutions.
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
8481
Fig. 12. Experiment 4: fourth-order numerical solution at several time moments (a ¼ 3).
3.5. Experiment 5—Modelling of combustion The phenomenon of combustion was studied in [8] on the real line as a Cauchy problem described by the nonlinear diffusion equation with an appropriate parameter a and a source function f. We tested the developed method having simulated this phenomenon within a bounded area on a sphere. This problem has an interesting application in metal flaming and burning under the laser action [34]. For the initial condition we took a spot with the centre at k ¼ u ¼ p4 . We also had l ¼ 1 103 ; a ¼ 1 and f ¼ cT b , where c ¼ 4:1. The parameter b > 0 introduces an additional, aside from a, nonlinearity, now with respect to the sources, and determines two different regimes of combustion: the case b < a þ 1 corresponds to the HS-regime—the area of combustion is getting expanded while burning; the case b > a þ 1 corresponds to the LS-regime—the area of combustion is getting narrow [8].
8482
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 13. Experiment 4: fourth-order numerical solutions at several time moments: a ¼ 0 (top), a ¼ 1 (middle), a ¼ 2 (bottom).
Table 1 Experiment 4: maximum relative errors dnmax ’s at different a’s.
a
dnmax 2nd order
4th order
0
4:650 103
6:587 104
1
3
3:639 10
5:397 104
2
1:314 102
7:939 103
3
2
2:591 102
3:496 10
The point is that in both cases the given source leads to an infinite increase of the temperature T while the time grows, that is a blow-up occurs [8]. In Fig. 14 we present the numerical solutions related to b ¼ 1 (HS-regime) and b ¼ 4 (LS-regime). In Fig. 15 we also plot the L2 -norms of the second-order numerical solutions. The experiment confirms the theory: in the HS-regime we observe the increase of the area of combustion, while in the LS-regime the area is reducing to a singularity. The graphs of the energy prove the solutions tend to infinity, that is a blow-up happens. Further analysis allows to conclude that in the HS-regime the entire character of the blow-up asymptotics is slower than that in the LS-regime—the curve of the solution’s L2 -norm at b ¼ 1 is growing bit by bit all the time from t ¼ 0 to t ¼ 1, while at b ¼ 4 it is going up slowly only up to t 0:1, and
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
8483
Fig. 14. Experiment 5: second-order numerical solutions at several time moments for HS-regime (first two rows) and LS-regime (last two rows).
it rapidly changes the growth rate after that. In its turn, the growth of the combustion area observed in Fig. 14 is smooth at b ¼ 1 all the time, whereas it is much sharper at b ¼ 4 after t 0:1 (cf., e.g., t ¼ 0:09 and t ¼ 0:11).
8484
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485
Fig. 15. Experiment 5: graphs of the L2 -norm of the second-order solutions in time for HS-regime (left) and LS-regime (right).
4. Conclusions A new method for the numerical solution of nonlinear diffusion equations on a sphere has been developed. Due to the splitting of the differential operator by coordinates and the use of different coordinate maps for the sphere at two split time intervals, both one-dimensional split problems obey periodic boundary conditions. Thus, although the sphere is not a doubly periodic manifold, our method allows to construct a numerical algorithm for the sphere in the same way as for a torus. Hence, the method does not require to impose any artificial boundary conditions at the poles (which, in addition, must be physically and mathematically correct)—and it must be so indeed, since the sphere has no boundaries and preferred points, and so all directions on the sphere are equivalent. This resulted in constructing second- and fourth-order finite difference schemes in space without using computationally cumbersome numerical procedures for computing the numerical solution at the poles. The schemes keep the properties required by the original differential problem: they appear balanced and dissipative, thus providing physically adequate numerical solutions. The schemes are also computationally inexpensive and can be solved by direct band linear solvers. The numerical experiments proved the approach to allow accurate simulating diverse diffusion phenomena with constant and variable diffusion coefficients on the entire sphere, including those that imply high nonlinearity, nonsmoothness and infinite growth of the solutions. The accuracy of the method was evaluated by comparing the numerical results with the analytical solutions obtained under forcings of special types. The convergence of the numerical solution to the analytics was verified on a series of spatial grids 6 6 ; 4 4 and 2 2 . As an interesting nonlinear effect, it was also constructed an example showing that the competition of the three basic mechanisms of the model—the nonlinear interaction, external forcing and dissipation—generates a wave solution, whose spatial structure varies periodically due to the alternating influence of the processes of self-organisation and self-destruction. Acknowledgements This research was partially supported by grants No. 14539 and 26073 of the National System of Researchers of Mexico (SNI), and is part of the project PAPIIT-UNAM IN104811, Mexico. The authors express gratitude to the anonymous reviewer for the remarks resulted in improvement of the paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
J. Bear, Dynamics of Fluids in Porous Media, Courier Dover Publications, 1988. F. Catté et al, Image selective smoothing and edge detection by nonlinear diffusion, SIAM J. Numer. Anal. 29 (1992) 182–193. M.E. Glicksman, Diffusion in Solids: Field Theory, Solid-State Principles and Applications, John Wiley & Sons, 2000. J.R. King, ‘Instantaneous source’ solutions to a singular nonlinear diffusion equation, J. Eng. Math. 27 (1993) 31–72. A.A. Lacey, J.R. Ockerdon, A.B. Tayler, ‘Waiting-time’ solutions of a nonlinear diffusion equation, SIAM J. Appl. Math. 42 (1982) 1252–1264. G.A. Rudykh, E.I. Semenov, Non-self-similar solutions of multidimensional nonlinear diffusion equations, Math. Notes 67 (2000) 200–206. A.Kh. Vorob’yov, Diffusion Problems in Chemical Kinetics, Moscow University Press, 2003 (in Russian). A.A. Samarskii et al, Blow-up in Quasilinear Parabolic Equations, Walter de Gruyter, 1995. Yu.N. Skiba, Dual oil concentration estimates in ecologically sensitive zones, Environ. Monit. Assess. 43 (1996) 139–151. V.I. Agoshkov, F. Saleri, Recent developments in the numerical simulation of shallow water equations Part. III: boundary conditions and finite element approximations in the river flow calculations, Math. Model. 8 (1996) 3–24. X. Zheng, P. Palffy-Muhoray, Nonlinear diffusion on a sphere, Bull. Am. Phys. Soc. 57 (2012). W.-Z. Huang, D.M. Sloan, Pole condition for singular problems: the pseudospectral approximation, J. Comput. Phys. 107 (1993) 254–261. X.-Y. Chen et al, Stabilization of vortices in the Ginzburg–Landau equation with a variable diffusion coefficient, SIAM J. Math. Anal. 29 (1998) 903–912. M.-C. Lai, A note on finite difference discretizations for Poisson equation on a disk, Numer. Methods Partial Diff. Eq. 17 (2001) 199–203. M.-C. Lai, Y.-H. Tseng, A fast iterative solver for the variable coefficient diffusion equation on a disk, J. Comput. Phys. 208 (2005) 196–205.
Y.N. Skiba, D.M. Filatov / Applied Mathematics and Computation 219 (2013) 8467–8485 [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]
8485
P. Gérard, F. Méhats, The Schrödinger–Poisson system on the sphere, SIAM J. Math. Anal. 43 (2011) 1232–1268. W. Hundsdorfer, J.G. Verwer, Numerical Solution of Time-Dependent Advection–Diffusion–Reaction Equations, Springer, 2003. S.V. Tsynkov, Numerical solution of problems on unbounded domains. A review, Appl. Numer. Math. 27 (1998) 465–532. Z. Wu et al, Nonlinear Diffusion Equations, World Scientific Publishing, Singapore, 2001. Yu.N. Skiba, D.M. Filatov, Simulation of nonlinear diffusion on a sphere, in: Proceedings of the 1st International Conference on Simulation and Modelling Methodologies, Technologies and Applications (SIMULTECH 2011), Noordwijkerhout, the Netherlands, 2011, pp. 24–30. D.W. Peaceman, H.H. Rachford Jr., The numerical solution of parabolic and elliptic differential equations, J. Soc. Ind. Appl. Math. 3 (1955) 28–41. N.N. Yanenko, The Method of Fractional Steps: Solution of Problems of Mathematical Physics in Several Variables, Springer, 1971 (translated from Russian: Nauka, Novosibirsk, 1967). G.I. Marchuk, Methods of Computational Mathematics, Springer, 1982 (translated from Russian: Nauka, Moscow, 1977). D.L. Williamson, Difference approximations of flow equations on the sphere, in: Numerical Methods Used in Atmospherical Models, vol. II. GARP Publ. Series, No 17, 1979. K. Mohseni, T. Colonius, Numerical treatment of polar coordinate singularities, J. Comput. Phys. 157 (2000) 787–795. Yu.N. Skiba, D.M. Filatov, On splitting-based mass and total energy conserving arbitrary order shallow-water schemes, Numer. Methods Partial Diff. Eq. 23 (2007) 534–552. Yu.N. Skiba, D.M. Filatov, Conservative arbitrary order finite difference schemes for shallow-water flows, J. Comput. Appl. Math. 218 (2008) 579–591. G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 506–517. W.H. Press et al, Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, Cambridge, 2007. G.A. Korn, T.M. Korn, Mathematical Handbook for Scientists and Engineers, Courier Dover Publications, 2000. V.S. Zykov, S.C. Müller, Spiral waves on circular and spherical domains of excitable medium, Physica D 97 (1996) 322–332. H. Yagisita, M. Mimura, M. Yamada, Spiral wave behaviors in an excitable reaction-diffusion system on a sphere, Physica D 124 (1998) 126–136. K. Rohlf, L. Glass, R. Kapral, Spiral wave dynamics in excitable media with spherical geometries, Chaos 16 (2006) 037115 (10 pp.). A.A. Samarskii, Nonlinear effects of blow-up and localization processes in burning problems, in: C.-M. Brauner, C. Schmidt-Lainé (Eds.), Mathematical Modeling in Combustion and Related Topics, Martinus Nijhoff Publishers, 1988, pp. 217–231.
