Computational Engineering and Technical Physics Technomathematics
Ibrahim Khalifa
COMPARING NUMERICAL METHODS FOR SOLVING THE FISHER EQUATION
Master’s Thesis
Examiners: Professor Heikki Haario Professor Lassi Roininen Supervisor: Professor Heikki Haario
Lappeenranta-Lahti University of Technology LUT School of Engineering Science
Computational Engineering and Technical Physics Technomathematics
Ibrahim Khalifa
Comparing Numerical Methods for Solving the Fisher Equation
Master’s Thesis 2020
49 pages.
Examiners: Professor Heikki Haario Professor Lassi Roininen
Keywords: Fisher’s equation, Finite Difference Methods, The Method of Lines
This thesis provides a comparison of results from Finite Difference Methods, pdepe solver and the Method of Lines when solving the Fisher equation and the heat equation. We study the Fisher equation, which is an example of the parabolic equation and the heat equation as a special case of the Fisher equation. We compute the numerical solutions of the two parabolic equations using the FDMs, the pdepe MATLAB solver and the MOLs and compare their results with the analytical solutions at different times. Finally the CPU time and the accuracy were computed and presented when solving these parabolic equations.
I would like to thank my supervisor Professor Heikki Haario for his continuous support and encouragement during the course of this project. His valuable contribution to the success of this study cannot be quantified. Many thanks to Nosiba Elfadel for her support and cooperation during this work. I also want to thank my parents and everyone who helped me and encouraged me to finish this project.
Lappeenranta, June 25, 2020
Ibrahim Khalifa
CONTENTS
1 Introduction 6
1.1 Related work . . . 7
1.2 Objectives . . . 8
1.3 Scope . . . 8
1.4 Structure of the thesis . . . 8
2 Analytical Solutions of the Parabolic Equations 9 2.1 Analytical solution of the Fisher equation . . . 9
2.1.1 Dimensional analysis . . . 9
2.2 Analytical solution of the heat equation . . . 15
3 Finite Difference Methods 18 3.1 Taylor’s theorem . . . 18
3.2 Finite difference approximations . . . 18
3.3 The heat equation . . . 20
3.3.1 Explicit method . . . 21
3.3.2 Implicit method . . . 22
3.4 Fisher’s equation . . . 23
3.4.1 Explicit FD method . . . 23
3.4.2 A Semi implicit FD method . . . 24
3.5 The Method of Lines . . . 26
3.5.1 Five point approximations . . . 26
4 MATLAB’s pdepe Function 27 4.1 The MATLAB pdepe solution of the Fisher equation . . . 28
4.2 The MATLAB pdepe solution of the heat equation . . . 30
5 Numerical Experiments 33 5.1 Fisher’s equation results . . . 33
5.1.1 Explicit FDM results . . . 33
5.1.2 Semi Implicit FDM results . . . 35
5.1.3 The Method of Lines results . . . 36
5.1.4 pdepe solver results . . . 38
5.1.5 Comparing the CPU time and the accuracy . . . 40
5.2 The heat equation results . . . 40
6 Discussion and Conclusions 47
LIST OF ABBREVIATIONS
ASD Accurate Space Derivative BCs Boundary Conditions DEs Differential Equations
EFDM Explicit Finite Difference Method FD Finite Difference
FDM Finite Difference Method FDMs Finite Difference Methods FEM Finite Element Method IC Initial Condition
IFDM Implicit Finite Difference Method IVP Initial Value Problem
MOL Method of Lines
ODE Ordinary Differential Equation ODEs Ordinary Differential Equations PDE Partial Differential Equation PDEs Partial Differential Equations RDE Reaction Diffusion Equation
1 Introduction
In explaining a wide variety of natural phenomena undergoing transition, differential equations (DEs) play a major role. Partial and ordinary differential equations occur in different applications like fluid flow, mechanical processes, image processing, earth sci- ences, physics, and mathematical finance [1]. Mathematicians have studied the essence of those equations for many years, and there are many well-developed methods for comput- ing the solutions of the DEs. For a given DE, the analytical solution is usually preferable but many of them can not be resolved exactly by using analytical methods. It is important to establish the numerical methods in order to solve differential equations in this case.
However, there’s a growing demand today for simpler and more effective approaches due to the increasing demand of applications [1]. In science and engineering, a signifi- cant number of basic physical phenomena also are modelled using nonlinear PDEs. Large number of researchers have applied analytical and computational methods to get solutions to the nonlinear differential equations.
Finite Difference Methods, The Method of Line (MOL) and Finite Element Method (FEM), are among the main numerical techniques for solving PDEs. There are other approaches also, such as the Crank-Nicolson method, Monte Carlo Methods, and so on.
In this thesis we employ some numerical methods for solving PDEs in particular, the Fisher equation and the heat equation. Fisher’s equation is PDE with nonlinear term.
Fisher’s equation has been studied first by R. A. Fisher in propagation of a gene within a population [2]. It’s also known as Fisher-Kolmogorov equation, and Fisher-KPP equa- tion (FKPP). Fisher’s equation is special case of reaction diffusion equation (RDE) and it’s useful in applications like biological phenomena [4], physical, chemical reaction pro- cesses, heat and mass transfer [80], optics, and combustion [3]. Many RDEs have travel- ling wave fronts that play a key role in understanding physical, chemical and biological phenomena [10], [11], [12]. Consider the following reaction-diffusion equation:
ut = (D(x, u)ux)x+g(x)f(u) (1) whereD(x, u)andf(u)are the diffusivity and source term, respectively. Note thatDis dependent onxanduandfonu. Furthermore, the coefficient off is spatially dependent.
Assume that the diffusion coefficient is constant this suggestsD(x, u) = k wherek is a constant, andg(x) = 1. Hence,
ut= (kux)x+f(u) (2)
Givenf(u) = αu(1−u)equation(2) become Fisher’s equation which is heat equation with additional non-linear term. The following is FKPP equation
∂u
∂t =k∂2u
∂x2 +αu(1−u) (3)
wherekis constant andαu(1−u)is that the nonlinear term its represent logistic growth.
Equation (3) has been studied first by R. A. Fisher in propagation of a gene within a population [2].
1.1 Related work
In [9] Ablowitz,M.J.and Zeppetella found that the Fisher KPP equation gives explicit solution for special speed c = ±5/√
6. In [19] Al-Khaled applied the sinc collocation method to get a solution for such model equation. Abdullaev [14] has described the stability of the travelling wave solutions. In [16] Jenö Gazdag and José Canosa have used the ASD method for PDEs to the Fisher equation solution. Logan has studied the FKPP by employing a perturbation technique and found an approximate answer by representing the answer as a power series [15]. Wang [10] provides the concepts of analytical and explicit solutions which are related to the general form, also Dag [20] found solution of Fisher’s equation, he found the answer by B-spline Galerkin technique, also Sloan and Qio [21] discover approximate solutions of Fisher KPP equation by the technigue of moving the mesh, in the mean time using element free Galerkin method Ting [21] studied generalized Fisher’s equation. Also modified cubic B-spline collocation approach is hired through Mittal and Jain [21] to take a look at numerical answers of nonlinear Fisher’s equation. In [25] Joseph G and charles R have studied the stochastically perturbed Fisher’s equation, in particular its wave velocity properties. Chandraker found also the solution of the Fisher equation numerically and Tomasiello discussed the stability in numerical method of differential quadrature solutions of wave problems [21].
In [26] M. F. K. Abur-Robb looked at the explicit solution of the Fisher quation with 3 zeros. This paper [21] implemented polynomial differential quadrature technique to look at numerical answers of non-linear equations. In [22] Kaysar carried out FD based methods to look for numerical answers of the Burger equation and the Fisher equation.
1.2 Objectives
Our objectives are as follows:
1. To introduce the FDMs, pdepe MATLAB solver and the Method of Lines (MOL).
2. To apply the methods in 1 to solve the Fisher equation and the heat equation.
3. To compare the CPU time and the accuracy of the methods in 1.
1.3 Scope
Our work is focusing on comparing the CPU time and the accuracy of the methods in 1.
The comparison is made between the FDMs, pdepe solver and the Method of Lines.
1.4 Structure of the thesis
This thesis contains six sections. In this introductory section we point to the importance of the differential equations and we define the Fisher equation. This section also contains a literature review of the related work of the Fisher equation, the objectives of the study and the scope of this work. In section 2 we looked at the analytical solution of the FKPP equation and the heat equation. The homotopy perturbation method is applied to find the analytical solution of FKPP equation. The method of separation of variables is used to find analytic solution for the heat equation. In section 3 we introduced the FDMs and its approximations and we gave the solution of the heat equation by using the EFDM and IFDM and we solved the Fisher equation using the explicit, semi implicit FDMs and in the end of the this section the Method of Lines is introduced. In section 4 we introduced the MATLAB pdepe built in function and used it to solve the Fisher equation and the heat equation. Section 5 discusses the results of the exact solution and the numerical solution of Fisher’s equation and the heat equation. Finally, Section 6 summarizes the results obtained from the methods and describes the possible future work.
2 Analytical Solutions of the Parabolic Equations
In this section we present in details the analytical solution of the Fisher equation using the homotopy perturbation method. The separation of variables technique is employed to find out an analytical solution for the heat equation
2.1 Analytical solution of the Fisher equation
The reaction-diffusion equation in the commonly known form can be written as shown in the following formula
∂u
∂t =k∂2u
∂x2 + +f(u) (4)
where f is a nonlinear function of u. Having the reaction term αu(1− u), equation(4) gives Fisher’s equation
∂u
∂t =k∂2u
∂x2 +αu(1−u) (5)
k is the diffusion coefficient andαis the growth rate. u =u(x, t)is the state variable in locationxand in timet.
2.1.1 Dimensional analysis
To nondimensionalize equation (5) we substitute equation(6) in equation (5). Let tˆ=αt xˆ=x α
k 1/2
(6) Then
x= k α
1/2
ˆ
x (7)
t = 1
αˆt (8)
By substituting (7) and (8) into (5), we get
∂u
∂t =α∂u
∂ˆt,
∂u
∂x = α k
1/2∂u
∂xˆ
∂2u
∂x2 = ∂
∂x
∂u
∂x = ∂
∂x α
k
1/2∂u
∂xˆ
(9)
= α k
1/2 ∂
∂xˆ α
k
1/2∂u
∂xˆ
(10)
=α k
∂2u
∂xˆ2 (11)
Hence equation (5) becomes α∂u
∂ˆt =k α k
∂2u
∂xˆ2
+αu(1−u) (12) Or,
∂u
∂ˆt = ∂2u
∂xˆ2 +u(1−u) (13)
Now drop the superscript hat notation, we get
∂u
∂t = ∂2u
∂x2 +u(1−u) (14)
The above Fisher equation with the IC
u(0, x) = λ (15)
can be solved using the homotopy perturbation method [24].
Now by constructing the homotopy we get
∂u
∂t − ∂u0
∂t =p ∂2u
∂x2 −u(1−u)−∂u0
∂t
(16)
Suppose that the solution of Fisher’s equation can be written as
u=u0+pu1+p2u2+... (17) By substituting (15), (17) into (16) and equating the terms of the same power ofp,we get
∂
∂t(u0+pu1+p2u2+...)− ∂u0
∂t =p ∂2
∂x2(u0+pu1+p2u2+...)−(u0+pu1+p2u2+...) + (u0+pu1+p2u2 +...)2−∂u0
∂t
Equate the terms of thep0on both sides we get
∂u0
∂t −∂u0
∂t = 0, u0(0, x) =λ,
Equate the terms of thep1on both sides we get
∂u1
∂t = ∂2u0
∂x2 +u0 −u20 , u1(0, x) = 0, Similarly we compare the terms ofp2, p3, ...on both sides we have
∂u2
∂t = ∂2u1
∂x2 −2u0u1+u1, u2(0, x) = 0,
∂u3
∂t = ∂2u2
∂x2 −2u0u2−u21+u2, u3(0, x) = 0,
∂u4
∂t = ∂2u3
∂x2 −2u0u3−2u1u2+u3, u4(0, x) = 0, ...
By solving the above system of equations, we obtain
∂u0
∂t −∂u0
∂t = 0andu(0, x) = λ=⇒u0(t, x) =λ
∂u1
∂t = ∂2u0
∂x2 +u0−u20
= ∂2λ
∂x2 +λ−λ2 By integrating the both sides we obtain
u1(t, x) =λ(1−λ)t,
∂u2
∂t = ∂2u1
∂x2 −u0u1+u1(1−u0)
= ∂2u1
∂x2 −u0u1−u0u1+u1
= ∂2
∂x2(λ(1−λ)t)
| {z }
=0
−2λ(λ(1−λ)t) +λ(1−λ)t
= −2λ2(1−λ)t+λ(1−λ)t
Therefore, integrating the both sides we obtain u2(t, x) = −2λ2(1−λ)t2
2 +λ(1−λ)t2 2
= (1−λ)(λ−2λ2)t2 2!
= λ(1−λ)(1−2λ)t2 2!
∂u3
∂t = ∂2u2
∂x2 −2u0u2−u21+u2
= ∂2
∂x2
λ(1−λ)(1−2λ)t2 2!
| {z }
=0
−2λ
λ(1−λ)(1−2λ)t2 2!
+
λ(1−λ)(1−2λ)t2 2!
2
+λ(1−λ)(1−2λ)t2 2!
=λ(1−λ)−2λ(1−2λ)
2 −λ(1−λ) + 1−2λ 2
t2
=λ(1−λ)
1−6λ+ 6λ2 t2
2!
By integrating both sides we get
u3(t, x) = λ(1−λ)(1−6λ+ 6λ2)t3 3!
By doing the same we have
u4(t, x) = λ(1−λ)(1−2λ)(1−12λ+ 12λ2)t4 4!, ...
(18)
Now by substituting p = 1 in (17) we can express the solution of equation (14) in the form
u=u0+u1+u2+u3+· · · (19) Hence, the solutionu(t, x)can be as follows
u(t, x) =λ+λ(1−λ)t+λ(1−λ)(1−2λ)t2
2!+λ(1−λ)(1−6λ+ 6λ2)t3 3!
+λ(1−λ)(1−2λ)(1−12λ+ 12λ2)t4 4!+· · ·
Therefor, the equation (14) has the following analytic close form solution u(t, x) = λet
1−λ+λet (20)
In the following we study the analytic close form solution of equation (5) when k = 1 andα = 6Therefore, we have
∂u
∂t = ∂2u
∂x2 + 6u(1−u) (21)
with the initial condition
u(0, x) = 1
(1 +ex)2 (22)
By constructing the homotopy we have
∂u
∂t −∂u0
∂t =p ∂2u
∂x2 −6u(1−u)− ∂u0
∂t
(23)
Substituting (22), (17) into (23) and equating the terms of the same power ofp,we get
∂
∂t(u0+pu1+p2u2+...)− ∂u0
∂t =p ∂2
∂x2(u0+pu1+p2u2+...)−6((u0+pu1+p2u2+...) + (u0+pu1+p2u2 +...)2)−∂u0
∂t
Equating the terms of thep0on both sides we get
∂u0
∂t − ∂u0
∂t = 0, u0(0, x) = 1 (1 +ex)2, Equating the terms of thep1on both sides we get
∂u1
∂t = ∂2u0
∂x2 + 6(u0−u20)− ∂u0
∂t , u1(0, x) = 0, Similarly
∂u2
∂t = ∂2u1
∂x2 −12u0u1+ 6u1, u2(0, x) = 0,
∂u3
∂t = ∂2u2
∂x2 −12u0u2−6u21+ 6u2, u3(0, x) = 0, ...
By solving the above system of equations, we get
∂u0
∂t − ∂u0
∂t = 0 =⇒u0(t, x) = 1 (1 +ex)2
∂u1
∂t = ∂2u0
∂x2 + 6(u0−u20)−∂u0
∂t
= ∂2
∂x2
1 (1 +ex)2
+ 6 1
(1 +ex)2 − 1 (1 +ex)4
− ∂
∂t
1 (1 +ex)2
| {z }
=0
= 4e2x−2ex
(1 +ex)2 +12ex+ 6e2x (1 +ex)4
= 10ex (1 +ex)4
By integrating both sides we obtain
u1(t, x) = 10ex(1 +ex)−3t
∂u2
∂t = ∂2u1
∂x2 −12u0u1+ 6u1
= ∂2
∂x2
10ex (1 +ex)3t
−12 1 (1 +ex)2
10ex (1 +ex)3t
+ 6 10ex (1 +ex)3t
= 10ex−70e2x+ 40e3x
(1 +ex)5 − 120ex
(1 +ex)5 + 60ex (1 +ex)3
t
= 50ex(−1 + 2ex) (1 +ex)4 t By integrating both sides we get
u2(t, x) = 25ex(−1 +e2x)(1 +ex)−4t2 Doing the same we have
u3(t, x) = 125
3 (−1 + 7ex−4e2x)(1 +ex)−5t3, ...
(24)
Substitutingp= 1in (17) we obtain
u=u0+u1+u2+u3+· · · (25)
Therefore, the solution can be written of the form u(t, x) = 1
(1 +ex)2 + 10ex
(1 +ex)3t+ 25ex(−1 +e2x) (1 +ex)4 t2+ +125(−1 + 7ex−4e2x)
3(1 +ex)5 t3+· · ·
Hence the equation (21) has the following the analytic close form solution u(t, x) = 1
(1 +ex−5t)2 (26)
Thus, applying the homotopy perturbation technique to the Fisher equation
∂u
∂t = ∂2u
∂x2 +αu(1−u) (27)
gives the following exact close form solution
u(t, x) = 1
[1 + exp(pα
6x− 5α6 t)]2 x∈[−∞,∞], and t ∈[0, τ] (28) whereαis a constant andc= 5p
α/6.
2.2 Analytical solution of the heat equation
Consider an IVP for the heat equation :
ut=kuxx, 0≤x≤1, t >0 (29) with
BCs: u(t,0) = u(t,1) = 0 ∀t >0 (30) IC: u(0, x) = f(x) ∀x∈[0,1] (31) Let us attempt to find the solution of equation (29) that satisfies the BCs by using the separation of variables method. We guess that we have solutions on the form
u(x, t) =X(t)T(t) (32)
X(t)andT(t)depend only onxandt. Substituting (32) back into (29), we get X00(x)
X(x) = 1 k
T0(t)
T(t) (33)
However, the above equation is equal to some constant value−λ.
X00(x) +λX(x) = 0, (34)
T0(t) +kλT(t) = 0. (35)
The BCs imply thatu(t,0) = 0 =X(0)T(t)andu(t,1) = 0 =X(1)T(t).
Then,X(x)must satisfy
X(0) = 0andX(1) = 1
withT(t)not equal to zero fort >0. For three distinct cases ofλwe get
I. λ <0 :
X(x) =C1e
√−λx+C2e
√−λx
=⇒C1 =C2 = 0which only provides a trivial solution.
II. λ= 0 :
X(x) =C1x+C2
Again, from the boundary conditionsC1 = C2 = 0,one gets only trivial solution of the problem.
III. λ >0 :
X(x) = C1cos(√
λx) +C2sin(√ λx) By substituting the BCs ,we obtain
X(0) =C1 = 0 X(1) =C2sin(√
λ) = 0⇒sin(√
λ) = 0 ⇒λn= (nπ)2, n= 1,2, ...
Therefore,
X(x) =Cnsin(nπx)
Now, the second equation forT(t)takes the form:
T0(t) +k(πn)2·T(t) = 0⇒T(t) =Bnexp(−k(nπ)2t), whereBnis constant.
Hence, the solution equation (29) can be written as u(t, x) =
∞
X
n=1
Ansin(nπx) exp(−k(nπ)2t), Anis const. (36)
Ancan be determined by the IC
u(x,0) =f(x), f(x)can be written as a Fourier series, as follows
f(x) =
∞
X
n=1
Fnsin(nπx) =
∞
X
n=1
Ansin(nπx)
An=Fn = 2
1
Z
0
f(η) sin(nπη)dη
Therefore, the general solution is
u(t, x) =
∞
X
n=1
2
1
Z
0
f(η) sin(nπη)dη sin
nπx exp
−k(nπ)2t
. (37)
3 Finite Difference Methods
The aim of this method is to create a numerical scheme. It’s based on differential operator approximations. In this section the FDMs are discussed and used to solve two parabolic equations. The MOL method is also introduced at the end of the section.
3.1 Taylor’s theorem
Letu(x)which is defined over the interval(a, b)havencontinuous derivatives. Then for a < x0, x0+h < b
u(x0+h) = u(x0) +hux(x0) +h2uxx(x0)
2! +...+hn−1u(n−1)(x0)
(n−1)! +O(hn) (38) where,
ux = du
dx, uxx = d2u
dx2, ..., u(n−1) = dn−1u dxn−1.
ux(x0)is the derivative ofuw.r.txandO(hn)its the order ofhton.
This theorem states that if the value of uand its derivatives values at point x0 is known then equation(38)can be written at the pointx0+h.If we cancelO(h)in (38) we obtain an approximation ofu(x0+h).The termO(hn)is related to the error order.
3.2 Finite difference approximations
Now, we are constructing common FD approximations for common partial derivatives.
let us takeu = u(x, t) and approximate the partial derivatives ofuwith reference tox.
Consideringtis constant, ubecome a function ofxand so in equation (38)and replace the step size h by∆xhence(38)becomes,
u(t, x0+∆x) = u(t, x0)+(∆x)ux(t, x0)+(∆x)2
2! uxx(t, x0)+...+(∆x)n−1
(n−1)!un−1(t, x0)+O(∆x)n (39) Truncating equation(39)toO(∆x)2 gives,
u(t, x0+ ∆x) =u(t, x0) + (∆x)ux(t, x0) +O(∆x)2 (40)
Rearranging(40)gives,
ux(t, x0) = u(t, x0+ ∆x)−u(t, x0)
∆x −O(∆x)2
∆x Hence,
ux(t, x0) = u(t, x0+ ∆x)−u(t, x0)
∆x −O(∆x) (41)
When using numerical schemes to solve PDEs we’re restrained to a grid of discrete x values,x1, x2, x3, ..., xN and discretetlevels,tn, n= 0,1,2, .... Assume a constant grid space, so that∆x=xi+1−xi.Evaluating equation(41)at the point(tn, xi)gives,
ux(tn, xi) = u(tn, xi+1+ ∆x)−u(tn, xi)
∆x −O(∆x) (42)
For simplicity we will use the notation Uin for evaluating u at the point (tn, xi). Now using this new notation and dropping the error termO(∆x), equation(42)becomes,
ux(tn, xi) = Ui+1n −Uin
∆x (43)
Equation(43)is forward difference approximation toux(tn, xi)of order one. Substitute
−∆xinstead of∆xin(40)we get,
u(t, x0−∆x) =u(t, x0)−(∆x)ux(t, x0) +O(∆x)2 (44) Rearranging(44)as previously we get
ux(tn, xi) = Uin−Ui−1n
∆x (45)
Equation(45)is called backward difference approximation toux(tn, xi)of order one. Our equations(43)and(45)are of the first order but we can increase the order by taking more terms in the Taylor expansion. Many of the partial differential equations contains higher order derivatives terms so we have to derive their approximations. We shall restrict our focus to the PDEs of second order i.e. uxx.
Approximating(39)toO(∆x)4 gives,
u(t, x0+∆x) =u(t, x0)+(∆x)ux(t, x0)+(∆x)2
2! uxx(t, x0)+(∆x)3
3! uxxx(t, x0)+O(∆x)4 (46) Substitute−∆xinstead of∆xin(46)we get
u(t, x0−∆x) =u(t, x0)−(∆x)ux(t, x0)+(∆x)2
2! uxx(t, x0)−(∆x)3
3! uxxx(t, x0)+O(∆x)4 (47) Adding(46)and(47)gives,
u(t, x0+ ∆x) +u(t, x0−∆x) = 2u(t, x0) + ∆x2uxx(t, x0) +O(∆x)4 (48) Using our discrete notation, evaluating equation(48)at the point(tn, xi)gives,
Ui+1n +Ui−1n = 2Uin+ ∆x2uxx(t, x0) +O(∆x4) (49) drop the error termO(∆x4)and rewrite equation49, we obtain
uxx(t, x0) = Ui+1n −2Uin+Ui−1n
∆x2 (50)
Equation (50)is symmetric FD approximation of the second order. We can also derive FD approximations with respect totin the same way as we did onx.
3.3 The heat equation
The diffusion equation (or heat equation) is important for scientific fields and engineering problems. Its considered as special case of the Fisher equation when takingf(u) = 0.
Consider
ut =kuxx+f(u)
| {z }
=0
0< x < 1, 0< t < τ (51)
whereu=u(x, t)andkis constant. Having the IC
u(0, x) = f(x) (52)
and the following BCs
u(t,0) =g1(t) u(t,1) =g2(t)
(53)
The solution of(51),(52)and(53)is to findu(t, x)satisfying the boundary conditions [6].
We want to solve the heat equation numerically, so let us first partition the spatial interval [0,1]and[0, τ]as follows:
xi =i∆x i= 0,1, ..., N where 1
N = ∆x (54)
tn=n∆t n= 0,1, ..., M where τ
M = ∆t (55)
In the following we recall the Finite Difference Approximations in section 3.2 and use them to solve equation (51) by using the explicit and implicit methods.
3.3.1 Explicit method
The explicit method of finite difference based on the approximation of first derivative forward difference.
ut(tn, xi) = Uin+1−Uin
∆t , (56)
and the central difference to second order derivative.
uxx(tn, xi) = Ui+1n −2Uin+Ui−1n
∆x2 , (57)
Substituting(56)and(57)in equation(51)results Uin+1−Uin
∆t =kUi+1n −2Uin+Ui−1n
(∆x)2 (58)
The temperature at timen+ 1is determined by the temperature at timen, in the explicit finite difference method. However, solving forUin+1in equation(58)we get
Uin+1−Uin= k∆t
(∆x)2 Ui+1n −2Uin+Ui−1n
(59)
wherer = k∆t (∆x)2
Uin+1 =Uin+r Ui+1n −2Uin+Ui−1n
(60) Therefore,
Uin+1 =rUi+1n + (1−2r)Uin+rUi−1n (61)
Equation (61) is known as the explicit FD approximation to the heat equation given in equation(51)[7]. Also we can rewrite equation (61) in the matrix form as
U1n+1 U2n+1
... ... UN−1n+1
=
1−2r r 0 . . . 0
r 1−2r r 0
... . .. . .. ... ...
0 r
0 0 . . . r 1−2r
U1n+rU0n U2n
... ... UNn−1+rUNn
(62)
3.3.2 Implicit method
To derive the implicit form, we substitute (56) and (50) at timen+ 1in equation (51), we get
Uin+1−Uin
∆t =kUi+1n+1−2Uin+1+Ui−1n+1
(∆x)2 (63)
Equation (63) can be written as follows
Uin =−rUi−1n+1+ (1 + 2r)Uin+1−rUi+1n+1 (64) wherer = k∆t
(∆x)2 andiis from1toN −1.
Equation (64) is called the Implicit Method, this method can not be reorganized to achieve a simple algebraic formula similar to the explicit method. Although this is an implicit method drawback, this method has the benefit of being unconditionally consistent [8].
Equation ( 61) can be in the matrix form as follows
1 + 2r −r 0 . . . 0
−r 1 + 2r −r 0
... . .. . .. ... ...
0 −r
0 0 . . . −r 1 + 2r
U1n+1 U2n+1
... ... UN−1n+1
=
U1n+rU0n U2n
... ... UN−1n +rUNn
(65)
3.4 Fisher’s equation
In the following we want to study Fisher’s equation using a standard explicit and semi implicit finite difference schemes. The equation is given by
∂u
∂t = ∂2u
∂x2 +αu(1−u) (66)
with the BCs and IC
u(t,0) = 1
(1 +e−5t)2, 0≤t≤τ, (67)
u(t,1) = 1
(1 +e1−5t)2, 0≤t≤τ, (68) u(0, x) = 1
(1 +ex)2, 0≤x≤1. (69)
3.4.1 Explicit FD method
Partition the spatial interval[0,1]and[0, τ]as in subsubsection 3.3. Recall (56) and (57) and plug them into (66) we get
Uin+1−Uin
∆t = Ui+1n −2Uin+Ui−1n
(∆x)2 +αUin(1−Uin) (70) The above equation can be written as
Uin+1−Uin = ∆t
(∆x)2 Ui+1n −2Uin+Ui−1n
+α∆tUin(1−Uin) (71) Or,
Uin+1 =rUi−1n + (1−2r)Uin+rUi+1n +α∆tUin(1−Uin) (72) Uin+1 =rUi−1n + (1−2r+α∆t−α∆tUin)Uin+rUi+1n (73) wherer = (∆x)∆t2, n≥0, i= 1,2, ..., N −1.
We initialize the scheme by
UN0 = 1
(1 +ex)2, 0≤x≤1
The BCs atx= 0, x= 1are
U0M = 1
(1 +e−5t)2, 0≤t≤τ, U1M = 1
(1 +e1−5t)2, 0≤t ≤τ,
Equation (72) is called explicit FD approximation to the Fisher equation [7]. We can rewrite equation (72) in the matrix form as follows
U1n+1 U2n+1
... ... UN−1n+1
=
Q1 r 0 . . . 0
r Q2 r 0
... . .. ... ... ...
0 r
0 0 . . . r QN−1
U1n+rU0n U2n
... ... UN−1n +rUNn
(74)
where,
Qi = (1−2r+α∆t−α∆tUin), i= 1,2, ...N−1
3.4.2 A Semi implicit FD method
The domain is also discreetized on a fixed mesh in this method. The[0,1]space interval and the[0, τ]time interval are divided as follows
xi =i∆x i= 0,1, ..., N where 1
N = ∆x (75)
tn=n∆t n= 0,1, ..., M where τ
M = ∆t (76)
Here we calculate the derivatives at the time level ofn+ 1, but in the explicit FD method the derivatives are determined at the time level ofn. In this scheme, we use the forward difference in time and the central difference in space. In this scheme, the nonlinear term is linearized using the method of lagging. The method of lagging idea is that the solution is determined at different time levels, for example one is calculated at time stepn+ 1, and the other is calculated at time stepn.
Consider the Fisher equation
ut =uxx+αu(1−u) (77)
By substituting (56) and central FD approximation at n+ 1 and use method of lagging into the above Fisher equation we get,
Uin+1−Uin
∆t = Ui+1n+1−2Uin+1+Ui−1n+1
(∆x)2 +αUin+1(1−Uin) (78) Uin+1 = ∆t
(∆x)2 Ui+1n+1−2Uin+1+Ui−1n+1
+α∆tUin+1(1−Uin) +Uin (79) Or,
Uin+1 =r Ui+1n+1−2Uin+1+Ui−1n+1
+ηUin+1(1−Uin) +Uin (80) where,
r= ∆t
(∆x)2 and α∆t =η (81)
Equation (80) may be expressed as below
AiUi−1n+1+BiUin+1+SiUi+1n+1 =Di (82) such that Ai ,Si are constants. Ai = Si = −r and Bi = 1 + 2r−η(1−Ui). When the whole scheme is formed and the BCs are applied, the matrix for Ax = b where A is representing a tridiagonal matrix ofM ×N. Thexvector contains the values we are looking for. Alsox, bare vectors , andbcontains the all known values.
A=
B1 −r . . . .
−r B2 −r . . . ... ... ... . .. ... ... ... ... ... . . . .. ... −r ... ... . . . −r BN−2 −r . . . . . . −r BN−1 −r
x=
U1n+1 U2n+1
... ... UNn+1−1 UNn+1
, b=
U1n U2n ... ... UN−1n
UNn
(83)
In section 5, we provide the results of the semi-implicit scheme for the solution of the Fisher equation to confirm the theoretical results of the scheme. The MATLAB software is used to program the Fisher equation and to generate numerical results and graphs.
3.5 The Method of Lines
This Method is a numerical method for solving PDEs. The simple concept is to discretize in space only and use an ODE package to solve the resulting ODEs in time. Discretization can be done by using FD approximations, but more research is currently going on to use new techniques, such as Galerkin or collocation methods. In this method we discretize the space using the five point FD approximations.
3.5.1 Five point approximations
Consider a square grid in 1-dimension, the five point stencil of a pointxin the grid is (x−2∆x, x−∆x, x, x+ ∆x, x+ 2∆x)
The first and second approximations of a functionuat the pointxcan be determined using the following formulas
ux = −u(x+ 2∆x) + 8u(x+ ∆x)−8u(x−∆x) +u(x−2∆x) 12∆x
uxx = −u(x+ 2∆x) + 16u(x+ ∆x)−30u(x) + 16u(x−∆x)−u(x−2∆x) 12(∆x)2
it is called a five-point approximations, because it is derived from the Lagrange polyno- mials for the five points(x−2∆x, x−∆x, x, x+ ∆x, x+ 2∆x).In section 5 we give the solution of the Fisher equation and the heat equation by using the MOL method and compare the results with the exact solution.
4 MATLAB’s pdepe Function
MATLAB has a built in function called pdepe that can solve PDEs with initial and bound- ary values. It solves parabolic and elliptic PDEs in spacexand timet. Its built on FDMs and the FEM is used in the PDE Toolbox. The pdepe function solves PDEs of the follow- ing form
c(x, t, u,∂u
∂x)∂u
∂t =x−m ∂
∂x xmf(x, t, u,∂u
∂x)
+s(x, t, u,∂u
∂x)
where m = 0,1 or 2; f and s are given functions of x, t, u and ∂u∂x; c(x, t, u,∂u∂x) is a diagonal matrix with identically zero or positive coefficients.
with the IC and the BCs
u(t= 0, x) = u0(x) p(x, t, u) +q(x, t)f(x, t, u,∂u
∂x) = 0
wherepis a function depending onx, t, uandqis a function depending onx, t.
In order to solve any parabolic or elliptic PDEs system using pdepe solver we need to specify the following steps refer to Matlab help:
1. Write the PDEs as follows c(x, t, u,∂u
∂x)∂u
∂t =x−m ∂
∂x xmf(x, t, u,∂u
∂x)
+s(x, t, u,∂u
∂x) 2. Write the code for the PDE
function[c, f, s] =pdefunction(t, x, u, DuDx)
where x is the position, t represent the time, u is dependent variable and DuDx is uxwhich is the first derivative of u with respect to x.
3. Write the code for the IC
functionu0 = pdeinitial(x) which specifies the functionu(0, x) =f(x).
4. Write the code for the BCs
[pl, ql, pr, qr] =pdeboundary(t, xl, ul, xr, ur)
In the function pdeboundary,plandqlrepresent the left BCs, andprandqrrepre- sent the right BCs.
5. Select mesh points and time steps. For instance x=linspace(0,1,11) t=linspace(0,0.1,101) 6. Apply the PDE solver
m= 0,1or2; sol=pdepe(m,@pdefunction,@pdeinitial,@pdeboundary, x, t,options)
In the following the steps to solve the Fisher equation and heat equation when using the MATLAB pdepe function are given in details. Fisher’s equation solution is given in Fig.1 and the heat equation solution which is special case of the Fisher equation is shown when s(x, t, u,∂u∂x) = 0in Fig.2.
4.1 The MATLAB pdepe solution of the Fisher equation
This part gives in details the solution steps of the Fisher equation (66) using the built in function pdepe.
Firstly, write Fisher’s equation
∂u
∂t = ∂2u
∂x2 +αu(1−u)
= ∂
∂x(∂u
∂x) +αu(1−u) in the form of PDEs in step 1, this gives
c∂u
∂t =x−m ∂
∂x(xmf) +s Comparing the above equation with the Fisher’s equation
(1)∂u
∂t = ∂
∂x(∂u
∂x) +αu(1−u)
this means
m= 0, c(t, x) = 1, f(t, x, u,∂u
∂x) = ∂u
∂x ands(t, x, u,∂u
∂x) = αu(1−u) Write the initial condition
u0 = 1
(1 + exp(x))2 and the boundary conditions
u(t,0) = 1
(1 +e−5t)2
u(t,1) = 1
(1 +e1−5t)2 takes the following forms
pl(t, xl, u) +ql(t, xl)f(t, xl, u,∂u
∂x) = 0 pr(t, xr, u) +qr(t, xr)f(t, xr, u,∂u
∂x) = 0 By comparing the boundary conditions with the above forms we get
pl(t, xl, u) = ul− 1 (1 +e−5t)2 ql = 0
pr(t, xr, u) = ur− 1 (1 +e1−5t)2 qr = 0
Note thatulindicates theuat the left boundary andurindicates theuin the right boundary Step five take the mesh points
For instance
x=linspace(0,1,30)taking 30 meshes for (xmesh) t=linspace(0,2,10)taking 30 meshes for (tmesh) In the last solve using the function pdepe
m= 0; sol=pdepe(m,@pdefunction,@pdeinitial,@pdeboundary, x, t,options) where@pdefunction, @pdeinitial and@pdeboundary are the M-files for our PDE, the IC and the BCs respectively.
The following is the mesh plot solution to the Fisher Equation
Figure 1. pdepe function solution of the Fisher equation (66) computed with α = 6,∆x = 0.1, and∆t= 0.001.
4.2 The MATLAB pdepe solution of the heat equation
The heat equation is a specific case of Fisher’s equation, given the term αu(1−u) = 0.
In the following we give the solution steps of the heat equation using the built in function pdepe. Write the heat equation
∂u
∂t = ∂2u
∂x2
= ∂
∂x(∂u
∂x) in the form of PDEs this gives
c∂u
∂t =x−m ∂
∂x(xmf) +s comparing the above equation with the heat equation
(1)∂u
∂t = ∂
∂x(∂u
∂x)
this means
m= 0, c(t, x) = 1, f(t, x, u,∂u
∂x) = ∂u
∂x, s(t, x, u,∂u
∂x) = 0.
The only difference from the Fisher equation is thats(t, x, u,∂u
∂x) = 0 Now write the IC
u0 = 1
(1 + exp(x))2 and the BCs
u(t,0) = 1
(1 +e−5t)2
u(t,1) = 1
(1 +e1−5t)2 takes the following forms
pl(t, xl, u) +ql(t, xl)f(t, xl, u,∂u
∂x) = 0 pr(t, xr, u) +qr(t, xr)f(t, xr, u,∂u
∂x) = 0 By comparing the boundary conditions with the above forms we get
pl(t, xl, u) = ul− 1 (1 +e−5t)2 ql = 0
pr(t, xr, u) = ur− 1 (1 +e1−5t)2 qr = 0
Note thatul, ur indicates respectively the left and the right boundary conditions.
Take the mesh points For example take
x=linspace(0,1,30)taking 30 meshes for (xmesh) t=linspace(0,2,10)taking 30 meshes for (tmesh) In the last solve using the function pdepe
m= 0; sol=pdepe(m,@pdefunction,@pdeinitial,@pdeboundary, x, t,options)
where@pdefunction, @pdeinitial and@pdeboundary are the M-files for our PDE, the IC and the BCs respectively.
Figure 2. pdepe function Solution of the equation (66) computed with f(u) = s = 0,∆x = 0.1, and∆t= 0.001.
5 Numerical Experiments
This section compare the numerical results and the analytical solutions of the Fisher’s equation and the heat equation. We present the exact and the numerical solutions of the Fisher equation and the heat equation using the explicit and implicit FDMs, the method of lines and the pdepe MATLAB solver. Also we compare the CPU time and the accuracy of these methods when solving the Fisher equation and the heat equation. In the all figures we are taking the distancex∈[0,1]and the timet∈[0,0.1].
5.1 Fisher’s equation results
Recall the Fisher equation
∂u
∂t = ∂2u
∂x2 +αu(1−u) with boundary and initial conditions
u(t,0) = 1
(1 +e−5t)2, 0≤t≤τ,
u(t,1) = 1
(1 +e1−5t)2, 0≤t≤τ, u(0, x) = 1
(1 +ex)2, 0≤x≤1.
5.1.1 Explicit FDM results
Table 1 gives the EFDM solution and the exact solution of the Fisher equation and the ab- solute error between the two. The graphs represents the comparison between the explicit FDM solution and the exact solution of Fisher’s equation in different times.
Table 1.The Explicit FDM solution and the exact solution of the Fisher equation whenx∈[0,1]andt∈[0,0.1]at the final time t=0.1.
Distance x Approx. Solution Exact Solution Absolute Error
0 0.3875 0.3875 0
0.1 0.3584 0.3584 5.3785e-05
0.2 0.3299 0.3300 9.9794e-05
0.3 0.3022 0.3023 1.3615e-04
0.4 0.2754 0.2756 1.6145e-04
0.5 0.2498 0.2500 1.7444e-04
0.6 0.2255 0.2256 1.7375e-04
0.7 0.2025 0.2026 1.5779e-04
0.8 0.1810 0.1811 1.2470e-04
0.9 0.1610 0.1611 7.2620e-05
1 0.1425 0.1425 0
(a) (b)
(c) (d)
Figure 3. Comparison between the results obtained by the explicit FDM Scheme and exact solu- tion at different times : (a)t=0.005; (b)t=0.0490; (c)t=0.08; (d) at the finalt=0.1.
5.1.2 Semi Implicit FDM results
Table 2 gives the semi implicit solution, the exact solution of the Fisher equation and the absolute error between two. The graphs represents the comparison between the implicit FDM solution and the exact solution Fisher’s equation in different times.
Table 2.The Implicit FDM solution and the exact solution of the Fisher equation whenx∈[0,1]andt∈[0,0.1]at the final time t=0.1.
Distance x Approx. Solution Exact Solution Absolute Error
0 0.3875 0.3875 0
0.1 0.3585 0.3584 1.0684e-04
0.2 0.3302 0.3300 1.8115e-04
0.3 0.3025 0.3023 2.2825e-04
0.4 0.2759 0.2756 2.5246e-04
0.5 0.2503 0.2500 2.5690e-04
0.6 0.2259 0.2256 2.4332e-04
0.7 0.2029 0.2026 2.1213e-04
0.8 0.1813 0.1811 1.6247e-04
0.9 0.1611 0.1611 9.2527e-05
1 0.1425 0.1425 0
(a) (b)
(c) (d)
Figure 4. Comparison between the results obtained by the semi implicit FDM Scheme and exact solution at different times: (a)t=0.005; (b)t=0.0490; (c)t=0.08; (d) at final timet= 0.1
5.1.3 The Method of Lines results
Table 3 gives the MOL solution, the exact solution of the Fisher equation and the absolute error between two. The graphs represents the comparison between the MOL solution and the exact solution of the Fisher equation in different times.
Table 3.The MOL solution and the exact solution of the Fisher equation whenx∈[0,1]andt∈[0,0.1]at the final time t=0.1.
Distance x Approx. Solution Exact Solution Absolute Error
0 0.3875 0.3875 0
0.1 0.3591 0.3584 7.1559e-04
0.2 0.3307 0.3300 6.9152e-04
0.3 0.3029 0.3023 5.8317e-04
0.4 0.2761 0.2756 4.9226e-04
0.5 0.2504 0.2500 4.3369e-04
0.6 0.2261 0.2256 4.1043e-04
0.7 0.2031 0.2026 4.2087e-04
0.8 0.1816 0.1811 4.5476e-04
0.9 0.1615 0.1611 4.4794e-04
1 0.1425 0.1425 0
(a) (b)
(c) (d)
Figure 5. Comparison between the results obtained by the explicit FDM Scheme and exact solu- tion at different times: (a)t=0.005; (b)t=0.0490; (c)t=0.08; (d) at the finalt=0.1.
5.1.4 pdepe solver results
Table 4.The pdepe solution and the exact solution of the Fisher equation whenx∈[0,1]andt∈[0,0.1]at the final time t=0.1.
Distance x Approx. Solution Exact Solution Absolute Error
0 0.3875 0.3875 0
0.1 0.3584 0.3584 8.2098e-07
0.2 0.3300 0.3300 7.4960e-06
0.3 0.3023 0.3023 2.1840e-05
0.4 0.2756 0.2756 3.8433e-05
0.5 0.2501 0.2500 5.3360e-05
0.6 0.2257 0.2256 6.2959e-05
0.7 0.2027 0.2026 6.4103e-05
0.8 0.1812 0.1811 5.4493e-05
0.9 0.1611 0.1611 3.2995e-05
1 0.1425 0.1425 0
(a) (b)
(c) (d)
Figure 6. Comparison between the results obtained by the pdepe MATLAB solver and exact solution at different time where ∆x= 0.1, ∆t= 0.001and (a)t= 0.005; (b)t= 0.0490; (c) t= 0.08; (d) at final timet= 0.1.
5.1.5 Comparing the CPU time and the accuracy
In the following table we give the CPU time and the accuracy of the methods when solving the Fisher equation. In this computation x ∈ [0,1], t ∈ [0,0.1],∆x = 0.1 and∆t = 0.001.The accuracy computed by comparing the numerical solution and the exact solution using theL2-norm of the error. TheL2-norm is defined as
L2 =||uexact−ucomputed||2 = v u u t
n
X
i=1
|uexacti −ucomputedi |2 N
whereuexact ,ucomputed represents the exact and the numerical solution andN is the mesh nodes.
Explicit FDM Semi Implicit FDM MOL Pdepe Solver
CPU time in seconds 0.0023 0.0111 0.1385 0.2107
L2−norm 0.0026 0.0037 0.0113 0.0015
The above table shows that the all methods gives accurate results in short CPU time. We see that the pdepe solver gives the smallest error compare to the other methods. The explicit FDM gives good result in smallest time compare to the other methods.
5.2 The heat equation results
In this part we solve the heat equation :
ut=kuxx, 0≤x≤1, t >0 (84) with the boundary and initial condition
u(t,0) = u(t,1) = 0 ∀t >0 (85)
u(0, x) = f(x) ∀x∈[0,1] (86)
for which the exact solution is given as
u(t, x) = exp(−4π2t)·sin(2πx) (87) The solutions of the equation (84) have been obtained with the above initial and boundary conditions by using the explicit, implicit FDM, the Method of Lines and the pdepe solver.
The following figures are the comparison between the exact and the numerical results of these methods. We mainly compute the error for these methods when ∆t = 0.001 and∆t = 0.01. In the semi implicit method we found that considering∆x = 0.1with
∆t = 0.001 and 0.01 gives large error compare to the explicit method. Therefore we studied further when∆x= 0.01with∆t = 0.001and∆t= 0.00001and this gives better error up toO(10−3)andO(10−4)respectively.
(a) (b)
(c) (d)
Figure 7. Comparison between the results obtained by the explicit FDM Scheme and exact so- lution when x ∈ [0,1], t ∈ [0,0.1]with ∆x = 0.1 and ∆t = 0.001in the upper figures and
∆x= 0.1and ∆t= 0.01in the lower figures
(a) (b)
(c) (d)
Figure 8. Comparison between the results obtained by the implicit FDM Scheme and exact so- lution when x ∈ [0,1], t ∈ [0,0.1]with ∆x = 0.1 and ∆t = 0.001in the upper figures and
∆x= 0.1and ∆t= 0.01in the lower figures
(a) (b)
(c) (d)
Figure 9. Comparison between the results obtained by the implicit FDM Scheme and exact so- lution whenx ∈ [0,1], t ∈ [0,0.1]with ∆x = 0.01and ∆t = 0.001in the upper figures and
∆x= 0.01and ∆t= 0.01in the lower figures
(a) (b)
(c) (d)
Figure 10. Comparison between the results obtained by the implicit FDM Scheme and exact solution whenx ∈[0,1], t ∈ [0,0.1]with∆x = 0.01and ∆t= 0.001in the upper figures and
∆x= 0.01and ∆t= 0.00001in the lower figures
(a) (b)
(c) (d)
Figure 11. Comparison between the results obtained by the pdepe solver and exact solution when x ∈[0,1], t∈ [0,0.1]withM = 11, N = 101in the upper figures andM = 100, N = 100in the lower figures.
