Numerical Solution of Parabolic in Partial Differential Equations (PDEs) in One and Two Space Variable ()
1. Introduction
Partial differential equations (PDEs) form the basis of very many mathematical models of physical, chemical and biological phenomena, and more recently their use has spread into economics, financial forecasting, image processing and other fields. There are other explicit numerical methods that can be applied to the 1-D heat or diffusion equation such as the Method of Lines which is used by Matlab and Mathematica [1] [2]. In calculating the quantities to a good approximation, there is a thin boundary layer near the wing surface where viscous forces are important and that outside this an inviscid flow can be assumed. Which we will assume is locally flat, we can model the flow by
(1)
where u is the flow velocity in the direction of the tangential co-ordinate x, y is the normal co-ordinate, ν is the viscosity, ρ is the density and p the pressure; we have here neglected the normal velocity. This is a typical parabolic equation for
u with
treated as a forcing term [3] [4]. Away from the wing, considered
just as a two-dimensional cross section, we can suppose the flow velocity to be inviscid and of the form
where u and v are small compared with the flow speed at infinity,
in the x-direction. One can often assume that the flow is irrational so that we have
(2)
then combining the conservation laws for mass and the x-component of momentum, and retaining only first-order quantities while assuming homentropic flow, we can deduce the simple model
(3)
where
is the Mach number at infinity,
, and
is the sound speed [5] [6].
2. Explicit Methods for 1-D Heat or Diffusion Equation
2.1. Difference Approximations for Derivative Terms in PDEs
We consider
for
,
.
Discrete time and spatial variable x:
Let
Consider Taylor series expansion for
:
(4)
If we only consider
terms in Equation (4) then we arrive at the forward difference in time approximation for
:
We can also derive a higher-order approximation for
if we consider the Taylor series expansion for
as well:
(5)
(4) – (5)
leap-frog
(or center difference) in time. This gives higher-order accuracy than forwarding difference. We can also perform similar manipulations to arrive at approximations for the second derivative
:
(4) + (5)
central difference
The finite difference method makes use of the above approximations to solve Partial Differential Equations, PDEs numerically [7] [8] [9].
2.2. Central Differences
The central differences solve Partial Differential Equations, as well:
(6-a)
(6-b)
when the central difference operator is applied twice we obtain the very useful second-order central difference
(7)
For first differences, it is often convenient to use the double interval central difference
(8)
A Taylor series expansion of the forward difference in t gives for the solution of:
for
(9)
(10)
By adding together, the Taylor series expansions in the x variable for
, and
, we see that all the odd powers of
cancel, giving
(11)
We can now define the truncation error of the scheme
(12)
We first multiply the difference equation throughout by a factor, if necessary, so that each term is an approximation to the corresponding derivative in the differential equation. Here this step is unnecessary, provided that we use the form
(13)
rather than (12) [10] [11] [12]. The truncation error is then the difference between the two sides of the equation, when the approximation
is replaced throughout by the exact solution
of the differential equation. Indeed, at any point away from the boundary we can define the truncation error
[13] [14].
2.3. Definition
The truncation error
is define by:
(14)
so that:
(15)
where these leading terms are called the principal part of the truncation error, and we have used the fact that u satisfies the differential equation [5] [15].
We have used Taylor series expansions to express the truncation error as an infinite series [16]. It is often convenient to truncate the infinite Taylor series, introducing a remainder term, for instance:
(16)
where
lies somewhere between t and
. If we do the same thing for the x expansion the truncation error becomes
(17)
where
, from which it follows that:
(18)
(19)
where
is a bound for
and
is a bound for
. It is now clear why we assumed that the initial and boundary data for u were consistent, and why it is helpful if we can also assume that the initial data are sufficiently smooth.
For then we can assume that the bounds
and
hold uniformly over the closed domain
[9] [17].
For example, suppose the boundary conditions specify that u must vanish on the boundaries
and
, and that u must take the value 1 on the initial line, where
. Then the solution
is obviously discontinuous at the corners, and in the full domain defined by
all its derivatives are unbounded, so our bound for the truncation error is useless over the full domain [18] [19].
3. Implicit Backward Euler Method for 1-D Heat Equation
Unconditionally stable (but usually slower than explicit methods). Implicit because it evaluates difference approximations to derivatives at next time step
and not current time step we are solving for
.
(20)
(21)
becomes:
(22)
where
as before. We still need to solve for
given
is known. This requires solving a tridiagonal linear system of n equations [1] [20].
Again we let
,
and
.
3.1. Numerical Implementation of the Implicit Back Ward Euler Method
We are solving the same problem:
3.2. Example
We are solving the same system again with the method of lines:
where the initial conditions are
boundary conditions are
and
.
Again we have:
we replace
where
with boundary conditions:
(23)
In matrix form for
elements:
(24)
The Dirichlet boundary conditions is:
For simplicity we consider only 4 elements in x in this example to find the matrix system we need to solve for:
We writing
as a matrix equation:
(25)
4. Heat Conservation Properties
Assume that in our model heat flow problem
we define the total heat in the system at time t by:
(26)
Then from the differential equation we have:
(27)
This is not very helpful if we have Dirichlet boundary conditions: but suppose we are given Neumann boundary conditions at each end, say,
and
. Then we have:
(28)
so thatn is given by integrating an ordinary differential equation [21] [22]. Now suppose we carry out a similar manipulation for the
-method equations:
(29)
introducing the total heat by means of a summation over the points for which (29) holds:
(30)
Then, recalling from the definitions of the finite difference notation that
(31)
we have:
(32)
So, recalling from the definitions of the finite difference notation that
We have
(33)
The rest of the analysis will depend on how the boundary condition is approximated. Consider the simplest case as in
(34)
namely we set
(35)
as an approximation to (28) this approximation may be very accurate, even though we have seen that Un may not give a good pointwise approximation to un.
In particular, if
and
are independent of t the change in H in one time step exactly equals that in
[5] [17] [23].
To make the most of this matching we should relate (30) as closely as possible to (27) [13].
If u and U were constants that would suggest we take
, rather than
as we have been assuming; and we should compare
with
(36)
And that it is centered at
and we have
. (37)
Note that this interpretation matches very closely the scheme that we were led to in
by analyzing the truncation error. It would also mean that for initial condition we should take
as defined by (36) [10] [12]. Then for time-independent boundary conditions, we have
for all n. Moreover, it is easy to see that the function:
(38)
with any constant C satisfies the differential equation, and the two boundary conditions [7] [11] [17].
5. Numerical Test Problems
This section is concerning to the numerical test problems and their visualization.
Problem. Let the coupled system of FPDEs given as
such that the external functions
and
are given as
The exact solution of the above system is
Table 1. Absolute error at various values of
for
in
and
of problem.
We evaluate the approximate solution of Problem with our proposed method. [8] [18]. The comparison between exact and approximate solution, while the absolute error corresponds to scale level
. We have also computed absolute error at various scale levels and different points of the spaces as given in Table 1.
6. Conclusion
In this paper, we have developed an efficient numerical technique. We shall be concerned with the numerical solution of parabolic equations in one space variable and the time variable t. We begin with the simplest model problem, for heat conduction in a uniform medium. For this model problem, an explicit difference method is very straightforward in use, and the analysis of its error is easily accomplished by the use of a maximum principle, or by Euler Method, the numerical solution becomes unstable unless the time step is severely restricted, so we shall go on to consider other, more elaborate, numerical methods which can avoid such a restriction. The additional complication in the numerical calculation is more than offset by the smaller number of time steps needed. We then extended the methods to problems with more general boundary conditions, then to more general linear parabolic equations. Finally, we shall discuss the more difficult problem of the solution of nonlinear equations.
Acknowledgements
I would like to thank my supervisor, Dr. Muhsin Hassan Abdallah who was a great help to me and also I thank my husband Bashir Alfadol Albdawi without whose help, I could not have written this paper.