# Wavelet Method for Numerical Solution of Parabolic Equations

Authors: A. H. Choudhury
Journal: Journal of Computational Engineering
 (2014)

---

## Abstract

We derive a highly accurate numerical method for the solution of parabolic partial differential equations in one space dimension using semidiscrete approximations. The space direction is discretized by wavelet-Galerkin method using some special types of basis functions obtained by integrating Daubechies functions which are compactly supported and differentiable. The time variable is discretized by using various classical finite difference schemes. Theoretical and numerical results are obtained for problems of diffusion, diffusion-reaction, convection-diffusion, and convection-diffusion-reaction with Dirichlet, mixed, and Neumann boundary conditions. The computed solutions are highly favourable as compared to the exact solutions.


---

## Body

## 1. Introduction

In this paper, we seek a numerical solution of linear parabolic partial differential equation (PDE) in one space dimension given by(1)
∂
u
∂
t
-
∂
∂
x
(
α
∂
u
∂
x
)
+
β
u
+
γ
∂
u
∂
x
=
f
(
x
,
t
)
,
a
<
x
<
b
,
t
>
0
,
with initial condition
(2)
u
(
x
,
0
)
=
g
(
x
)
,
a
<
x
<
b
,
and with boundary conditions of Dirichlet, mixed, or Neumann type, where α, β, and γ are functions of x and t, in general. Here we assume that α
>
0 for a
<
x
<
b and t
>
0. PDE (1) is classified as follows.(i)
Ifβ
=
0 and γ
=
0, PDE (1) is called diffusion(or heat)equation.
(ii)
Ifβ
≠
0 and γ
=
0, PDE (1) is called diffusion-reaction equation.
(iii)
Ifβ
=
0 and γ
≠
0, PDE (1) is called convection-diffusion equation.
(iv)
Ifβ
≠
0 and γ
≠
0, PDE (1) is called convection-diffusion-reaction equation.These equations have applications in a number of fields, for example, heat transfer, fluid mechanics, ground water pollutants, financial mathematics, and so on.Several methods exist for the solution of parabolic problems, for example, [1–7]. But still there is a need for modification of the solution methodology in case of (i) Neumann and mixed boundary conditions, (ii) time-dependent boundary conditions, and (iii) time-dependent source term f. The algorithm in this paper is suitable for such situations.In the present paper, we apply semidiscrete approximation for the solution of initial-boundary value problem (IBVP) governed by the PDE (1). Here we use wavelet-Galerkin (variational) method to discretize the space direction and the time variable is discretized by using various classical finite difference schemes. Wavelets in consideration here are Daubechies compactly supported wavelets [8] which are differentiable.Wavelet applications to the solution of PDE problems are relatively new. Some recent applications are [5, 6, 9–12] among many more. To discretize a PDE problem by wavelet-Galerkin method, the Galerkin bases are constructed from orthonormal bases of compactly supported wavelets. This can be done in a number of ways. If wavelets are used directly in such construction such as one in [4], the resulting basis functions do not satisfy various types of boundary conditions. The basis functions obtained by integrating wavelet functions in a suitable manner can be made to satisfy various types of boundary conditions; see Deka and Choudhury [12], which is a generalization of [13] for general domains. Also, this approach makes one free from computing connection coefficients [14] for the solution of higher dimensional problems and is convenient to use low order wavelets resulting in low computational efforts. The present paper is a generalization of [12] to the solution of parabolic problems. Basis functions satisfying various types of boundary conditions can also be obtained without integrating wavelet functions; see Choudhury and Deka [10]. For the demonstration of the algorithm presented in the paper, some numerical results are computed which show that the computed solutions compare highly favourably with the exact solutions. All the numerical tests are performed in MATLAB. The paper is organized as follows.Section2 describes wavelet integral functions and their approximation properties elaborately including a short introduction to wavelet functions. In Section 3, we formulate the problems to be solved. Sections 4 and 5 discuss the spatial and time approximations of the problems using wavelet-Galerkin method and classical finite difference schemes, respectively. In Section 6, we perform some numerical tests to demonstrate the method presented. Section 7 contains the necessary conclusion.


## 2. Wavelet Integrals and Their Approximation Properties

In this section, we discuss the approximations of various function (Sobolev) spaces in finite dimension, useful for wavelet-Galerkin solutions of PDE problems, using integrals of Daubechies scaling functions (the wavelet integrals) [12], after a short account of Daubechies scaling and wavelet functions. For details about Daubechies functions, we refer [8, 15].


### 2.1. Basic Properties of Daubechies Wavelets

For a positive integerN, consider two functions ϕ
,
ψ
∈
L
2
(
ℝ
) defined by
(3)
ϕ
(
x
)
=
∑
k
‍
a
k
ϕ
(
2
x
-
k
)
,
ψ
(
x
)
=
∑
k
‍
b
k
ϕ
(
2
x
-
k
)
,
where {
a
k
}
k
∈
ℤ and {
b
k
}
k
∈
ℤ are two specific sequences [8] such that a
k
=
b
k
=
0 for k
∉
{
0,1
,
…
,
S
}
,
S
=
2
N
-
1. The functions ϕ and ψ are called d
b
N scaling function and d
b
N wavelet function, respectively, and N is called their order. These functions are compactly supported with supp
⁡
(
ϕ
)
=
supp
⁡
(
ψ
)
=
[
0
,
S
] and they are available in wavelet toolbox of MATLAB for 1
⩽
N
⩽
45. They satisfy the following properties:
(4)
∫
ℝ
‍
ϕ
(
x
)
d
x
=
1
,
∑
k
‍
ϕ
(
x
-
k
)
=
1
,
∫
ℝ
‍
x
m
ψ
(
x
)
d
x
=
0
,
0
⩽
m
<
N
.
The integer dilates and translates of ϕ and ψ are defined as
(5)
ϕ
n
,
k
(
x
)
=
2
n
/
2
ϕ
(
2
n
x
-
k
)
,
ψ
n
,
k
(
x
)
=
2
n
/
2
ψ
(
2
n
x
-
k
)
,
for n
,
k
∈
ℤ.Now, for alln
∈
ℤ, define
(6)
V
n
=
L
2
-closure
(
span
⁡
{
ϕ
n
,
k
:
k
∈
ℤ
}
)
.
Then we have
(7)
⋃
n
∈
ℤ
‍
V
n
=
L
2
(
ℝ
)
,
in the sense that
(8)
lim
⁡
n
→
∞
P
n
f
=
f
,
where P
n is the L
2-orthogonal projection of L
2
(
ℝ
) onto V
n.


### 2.2. Approximation of Function Spaces Using Wavelet Integrals

We mention here that for an open interval(
a
,
b
) and for an integer m
⩾
1, the space
(9)
H
m
(
a
,
b
)
=
{
u
∈
H
m
-
1
(
a
,
b
)
:
u
′
∈
H
m
-
1
(
a
,
b
)
}
is called the Sobolev space of order m, which is a Hilbert space with inner product
(10)
〈
u
,
v
〉
m
=
∑
i
=
0
m
‍
∫
a
b
‍
u
(
i
)
v
(
i
)
d
x
and associated norm ∥
·
∥
m, where the derivatives are in weak (or distributional) sense. It is to be noted here that H
0
(
a
,
b
)
=
L
2
(
a
,
b
) is the space of all square (Lebesegue) integrable functions which is also a Hilbert space with inner product 〈
u
,
v
〉
=
∫
a
b
‍
u
v
d
x and associated norm ∥
·
∥.Now, we define the following subspaces ofH
m
(
a
,
b
) which will be used for the solution of PDE problems in this paper:
(11)
H
0
m
(
a
,
b
)
=
{
u
∈
H
m
(
a
,
b
)
:
u
(
a
)
=
u
(
b
)
=
0
}
,
H
*
m
(
a
,
b
)
=
{
u
∈
H
m
(
a
,
b
)
:
u
(
a
)
=
0
}
,
H
+
m
(
a
,
b
)
=
{
u
∈
H
m
(
a
,
b
)
:
u
(
b
)
=
0
}
.LetN be any positive integer and let ϕ and ψ be the d
b
N scaling function and wavelet function, respectively. Then there exists an integer m
,
0
⩽
m
<
N such that ϕ
,
ψ
∈
H
m
(
ℝ
) and the Sobolev space H
m
(
a
,
b
) can be approximated by the restrictions of translates and dilates of ϕ to (
a
,
b
) [4]. Hence H
m
+
1
(
a
,
b
) can be approximated by their integrals.Now we have, supp(
ϕ
)
=
supp
⁡
(
ψ
)
=
[
0
,
S
]. Let n
⩾
0 be any integer and let V
n be defined by ϕ as in (6). Now we define V
n
(
0
,
S
) to be the set of restrictions of all functions in V
n to (
0
,
S
). In fact, we take
(12)
V
n
(
0
,
S
)
=
span
⁡
{
ϕ
n
,
k
|
(
0
,
S
)
:
supp
⁡
(
ϕ
n
,
k
)
∩
(
0
,
S
)
≠
∅
}
.
We shift the support of ϕ from [
0
,
S
] to [
a
,
b
] by using the transformation
(13)
y
=
b
-
a
S
x
+
a
,
where 0
⩽
x
⩽
S and a
⩽
y
⩽
b.Now, consideringx
∈
[
a
,
b
] instead of x
∈
[
0
,
S
], we define the space V
n
(
a
,
b
) as
(14)
V
n
(
a
,
b
)
=
span
⁡
{
ϕ
n
,
k
|
(
a
,
b
)
:
k
∈
I
n
}
,
where
(15)
I
n
=
{
k
∈
ℤ
:
supp
⁡
(
ϕ
n
,
k
)
∩
(
a
,
b
)
≠
∅
}
.
Since (
a
,
b
) is bounded, the space V
n
(
a
,
b
) is finite dimensional and is a closed subspace of H
m
(
a
,
b
). By Proposition  4.1 in [12], dim(V
n
(
a
,
b
)
)
=
2
n
S
+
S
-
1 and a basis for V
n
(
a
,
b
) can be constructed as
(16)
{
ϕ
n
,
k
∈
V
n
(
a
,
b
)
:
1
-
S
⩽
k
⩽
2
n
S
-
1
}
.
Now, for ϕ
n
,
k
∈
V
n
(
a
,
b
), we define the spaces(i)
S
n
0
(
a
,
b
)
=
span
⁡
{
Φ
n
,
k
:
k
∈
I
n
}, where
(17)
Φ
n
,
k
(
x
)
=
∫
a
x
‍
ϕ
n
,
k
(
s
)
d
s
-
x
-
a
b
-
a
∫
a
b
‍
ϕ
n
,
k
(
x
)
d
x
,
(ii)
S
n
*
(
a
,
b
)
=
span
⁡
{
Φ
n
,
k
:
k
∈
I
n
}, where
(18)
Φ
n
,
k
(
x
)
=
∫
a
x
‍
ϕ
n
,
k
(
s
)
d
s
,and(iii)
S
n
+
(
a
,
b
)
=
span
⁡
{
Φ
n
,
k
:
k
∈
I
n
}, where
(19)
Φ
n
,
k
(
x
)
=
∫
x
b
‍
ϕ
n
,
k
(
s
)
d
s
.Now by Proposition  4.2 in [12], the spaces S
n
0
(
a
,
b
), S
n
*
(
a
,
b
), and S
n
+
(
a
,
b
) are finite dimensional subspaces of H
0
m
+
1
(
a
,
b
), H
*
m
+
1
(
a
,
b
), and H
+
m
+
1
(
a
,
b
), respectively, and the sets
(20)
{
Φ
n
,
k
∈
S
n
0
(
a
,
b
)
:
2
-
S
⩽
k
⩽
2
n
S
-
1
}
,
{
Φ
n
,
k
∈
S
n
*
(
a
,
b
)
:
1
-
S
⩽
k
⩽
2
n
S
-
1
}
,
{
Φ
n
,
k
∈
S
n
+
(
a
,
b
)
:
1
-
S
⩽
k
⩽
2
n
S
-
1
}
are bases for S
n
0
(
a
,
b
), S
n
*
(
a
,
b
), and S
n
+
(
a
,
b
), respectively.Figure1 shows the basis functions for S
n
0
(
0,1
), S
n
*
(
0,1
), and S
n
+
(
0,1
) for N
=
2 and n
=
2.Basis functions forS
2
0
(
0,1
), S
2
*
(
0,1
), and S
2
+
(
0,1
) for N
=
2.
(a)
(b)
(c)The derivatives of the above basis functions are the interpolated scaling functions themselves, with a little difference, as follows:(i)
Φ
n
,
k
′
(
x
)
=
ϕ
n
,
k
(
x
), for Φ
n
,
k
∈
S
n
*
(
a
,
b
),
(ii)
Φ
n
,
k
′
(
x
)
=
-
ϕ
n
,
k
(
x
), for Φ
n
,
k
∈
S
n
+
(
a
,
b
),
(iii)
Φ
n
,
k
′
(
x
)
=
ϕ
n
,
k
(
x
)
-
(
1
/
(
b
-
a
)
)
∫
a
b
‍
ϕ
n
,
k
(
x
)
d
x, for Φ
n
,
k
∈
S
n
0
(
a
,
b
).By Theorem  4.2 in [12] and a standard duality argument [13], the error estimates of the above approximations in H
m norm can be obtained as follows.For anyu
∈
H
0
t
+
1
(
a
,
b
) (resp., H
*
t
+
1
(
a
,
b
) or H
+
t
+
1
(
a
,
b
)) and u
n
∈
S
n
0
(
a
,
b
) (resp., S
n
*
(
a
,
b
) or S
n
+
(
a
,
b
)), we have
(21)
∥
u
-
u
n
∥
m
⩽
C
m
h
t
-
m
+
1
,
where h
=
2
-
n, 0
⩽
t
⩽
N, and C
m is a positive constant depending on m and N.Now, by Proposition  4.3 in [12], we can find a finite dimensional subspace of H
m
+
1
(
a
,
b
) given by
(22)
S
n
(
a
,
b
)
=
span
⁡
{
1
,
Φ
n
,
k
:
Φ
n
,
k
∈
S
n
*
(
a
,
b
)
,
k
∈
I
n
}
,
a basis for which is given by
(23)
{
1
,
Φ
n
,
k
:
Φ
n
,
k
∈
S
n
*
(
a
,
b
)
,
1
-
S
⩽
k
⩽
2
n
S
-
1
}
.Remark 1.
The dimensions of the spaceS
n
0
(
a
,
b
) are 2
n
S
+
S
-
2, that of S
n
*
(
a
,
b
) and S
n
+
(
a
,
b
) are 2
n
S
+
S
-
1 each and of S
n
(
a
,
b
) is 2
n
S
+
S.


## 2.1. Basic Properties of Daubechies Wavelets

For a positive integerN, consider two functions ϕ
,
ψ
∈
L
2
(
ℝ
) defined by
(3)
ϕ
(
x
)
=
∑
k
‍
a
k
ϕ
(
2
x
-
k
)
,
ψ
(
x
)
=
∑
k
‍
b
k
ϕ
(
2
x
-
k
)
,
where {
a
k
}
k
∈
ℤ and {
b
k
}
k
∈
ℤ are two specific sequences [8] such that a
k
=
b
k
=
0 for k
∉
{
0,1
,
…
,
S
}
,
S
=
2
N
-
1. The functions ϕ and ψ are called d
b
N scaling function and d
b
N wavelet function, respectively, and N is called their order. These functions are compactly supported with supp
⁡
(
ϕ
)
=
supp
⁡
(
ψ
)
=
[
0
,
S
] and they are available in wavelet toolbox of MATLAB for 1
⩽
N
⩽
45. They satisfy the following properties:
(4)
∫
ℝ
‍
ϕ
(
x
)
d
x
=
1
,
∑
k
‍
ϕ
(
x
-
k
)
=
1
,
∫
ℝ
‍
x
m
ψ
(
x
)
d
x
=
0
,
0
⩽
m
<
N
.
The integer dilates and translates of ϕ and ψ are defined as
(5)
ϕ
n
,
k
(
x
)
=
2
n
/
2
ϕ
(
2
n
x
-
k
)
,
ψ
n
,
k
(
x
)
=
2
n
/
2
ψ
(
2
n
x
-
k
)
,
for n
,
k
∈
ℤ.Now, for alln
∈
ℤ, define
(6)
V
n
=
L
2
-closure
(
span
⁡
{
ϕ
n
,
k
:
k
∈
ℤ
}
)
.
Then we have
(7)
⋃
n
∈
ℤ
‍
V
n
=
L
2
(
ℝ
)
,
in the sense that
(8)
lim
⁡
n
→
∞
P
n
f
=
f
,
where P
n is the L
2-orthogonal projection of L
2
(
ℝ
) onto V
n.


## 2.2. Approximation of Function Spaces Using Wavelet Integrals

We mention here that for an open interval(
a
,
b
) and for an integer m
⩾
1, the space
(9)
H
m
(
a
,
b
)
=
{
u
∈
H
m
-
1
(
a
,
b
)
:
u
′
∈
H
m
-
1
(
a
,
b
)
}
is called the Sobolev space of order m, which is a Hilbert space with inner product
(10)
〈
u
,
v
〉
m
=
∑
i
=
0
m
‍
∫
a
b
‍
u
(
i
)
v
(
i
)
d
x
and associated norm ∥
·
∥
m, where the derivatives are in weak (or distributional) sense. It is to be noted here that H
0
(
a
,
b
)
=
L
2
(
a
,
b
) is the space of all square (Lebesegue) integrable functions which is also a Hilbert space with inner product 〈
u
,
v
〉
=
∫
a
b
‍
u
v
d
x and associated norm ∥
·
∥.Now, we define the following subspaces ofH
m
(
a
,
b
) which will be used for the solution of PDE problems in this paper:
(11)
H
0
m
(
a
,
b
)
=
{
u
∈
H
m
(
a
,
b
)
:
u
(
a
)
=
u
(
b
)
=
0
}
,
H
*
m
(
a
,
b
)
=
{
u
∈
H
m
(
a
,
b
)
:
u
(
a
)
=
0
}
,
H
+
m
(
a
,
b
)
=
{
u
∈
H
m
(
a
,
b
)
:
u
(
b
)
=
0
}
.LetN be any positive integer and let ϕ and ψ be the d
b
N scaling function and wavelet function, respectively. Then there exists an integer m
,
0
⩽
m
<
N such that ϕ
,
ψ
∈
H
m
(
ℝ
) and the Sobolev space H
m
(
a
,
b
) can be approximated by the restrictions of translates and dilates of ϕ to (
a
,
b
) [4]. Hence H
m
+
1
(
a
,
b
) can be approximated by their integrals.Now we have, supp(
ϕ
)
=
supp
⁡
(
ψ
)
=
[
0
,
S
]. Let n
⩾
0 be any integer and let V
n be defined by ϕ as in (6). Now we define V
n
(
0
,
S
) to be the set of restrictions of all functions in V
n to (
0
,
S
). In fact, we take
(12)
V
n
(
0
,
S
)
=
span
⁡
{
ϕ
n
,
k
|
(
0
,
S
)
:
supp
⁡
(
ϕ
n
,
k
)
∩
(
0
,
S
)
≠
∅
}
.
We shift the support of ϕ from [
0
,
S
] to [
a
,
b
] by using the transformation
(13)
y
=
b
-
a
S
x
+
a
,
where 0
⩽
x
⩽
S and a
⩽
y
⩽
b.Now, consideringx
∈
[
a
,
b
] instead of x
∈
[
0
,
S
], we define the space V
n
(
a
,
b
) as
(14)
V
n
(
a
,
b
)
=
span
⁡
{
ϕ
n
,
k
|
(
a
,
b
)
:
k
∈
I
n
}
,
where
(15)
I
n
=
{
k
∈
ℤ
:
supp
⁡
(
ϕ
n
,
k
)
∩
(
a
,
b
)
≠
∅
}
.
Since (
a
,
b
) is bounded, the space V
n
(
a
,
b
) is finite dimensional and is a closed subspace of H
m
(
a
,
b
). By Proposition  4.1 in [12], dim(V
n
(
a
,
b
)
)
=
2
n
S
+
S
-
1 and a basis for V
n
(
a
,
b
) can be constructed as
(16)
{
ϕ
n
,
k
∈
V
n
(
a
,
b
)
:
1
-
S
⩽
k
⩽
2
n
S
-
1
}
.
Now, for ϕ
n
,
k
∈
V
n
(
a
,
b
), we define the spaces(i)
S
n
0
(
a
,
b
)
=
span
⁡
{
Φ
n
,
k
:
k
∈
I
n
}, where
(17)
Φ
n
,
k
(
x
)
=
∫
a
x
‍
ϕ
n
,
k
(
s
)
d
s
-
x
-
a
b
-
a
∫
a
b
‍
ϕ
n
,
k
(
x
)
d
x
,
(ii)
S
n
*
(
a
,
b
)
=
span
⁡
{
Φ
n
,
k
:
k
∈
I
n
}, where
(18)
Φ
n
,
k
(
x
)
=
∫
a
x
‍
ϕ
n
,
k
(
s
)
d
s
,and(iii)
S
n
+
(
a
,
b
)
=
span
⁡
{
Φ
n
,
k
:
k
∈
I
n
}, where
(19)
Φ
n
,
k
(
x
)
=
∫
x
b
‍
ϕ
n
,
k
(
s
)
d
s
.Now by Proposition  4.2 in [12], the spaces S
n
0
(
a
,
b
), S
n
*
(
a
,
b
), and S
n
+
(
a
,
b
) are finite dimensional subspaces of H
0
m
+
1
(
a
,
b
), H
*
m
+
1
(
a
,
b
), and H
+
m
+
1
(
a
,
b
), respectively, and the sets
(20)
{
Φ
n
,
k
∈
S
n
0
(
a
,
b
)
:
2
-
S
⩽
k
⩽
2
n
S
-
1
}
,
{
Φ
n
,
k
∈
S
n
*
(
a
,
b
)
:
1
-
S
⩽
k
⩽
2
n
S
-
1
}
,
{
Φ
n
,
k
∈
S
n
+
(
a
,
b
)
:
1
-
S
⩽
k
⩽
2
n
S
-
1
}
are bases for S
n
0
(
a
,
b
), S
n
*
(
a
,
b
), and S
n
+
(
a
,
b
), respectively.Figure1 shows the basis functions for S
n
0
(
0,1
), S
n
*
(
0,1
), and S
n
+
(
0,1
) for N
=
2 and n
=
2.Basis functions forS
2
0
(
0,1
), S
2
*
(
0,1
), and S
2
+
(
0,1
) for N
=
2.
(a)
(b)
(c)The derivatives of the above basis functions are the interpolated scaling functions themselves, with a little difference, as follows:(i)
Φ
n
,
k
′
(
x
)
=
ϕ
n
,
k
(
x
), for Φ
n
,
k
∈
S
n
*
(
a
,
b
),
(ii)
Φ
n
,
k
′
(
x
)
=
-
ϕ
n
,
k
(
x
), for Φ
n
,
k
∈
S
n
+
(
a
,
b
),
(iii)
Φ
n
,
k
′
(
x
)
=
ϕ
n
,
k
(
x
)
-
(
1
/
(
b
-
a
)
)
∫
a
b
‍
ϕ
n
,
k
(
x
)
d
x, for Φ
n
,
k
∈
S
n
0
(
a
,
b
).By Theorem  4.2 in [12] and a standard duality argument [13], the error estimates of the above approximations in H
m norm can be obtained as follows.For anyu
∈
H
0
t
+
1
(
a
,
b
) (resp., H
*
t
+
1
(
a
,
b
) or H
+
t
+
1
(
a
,
b
)) and u
n
∈
S
n
0
(
a
,
b
) (resp., S
n
*
(
a
,
b
) or S
n
+
(
a
,
b
)), we have
(21)
∥
u
-
u
n
∥
m
⩽
C
m
h
t
-
m
+
1
,
where h
=
2
-
n, 0
⩽
t
⩽
N, and C
m is a positive constant depending on m and N.Now, by Proposition  4.3 in [12], we can find a finite dimensional subspace of H
m
+
1
(
a
,
b
) given by
(22)
S
n
(
a
,
b
)
=
span
⁡
{
1
,
Φ
n
,
k
:
Φ
n
,
k
∈
S
n
*
(
a
,
b
)
,
k
∈
I
n
}
,
a basis for which is given by
(23)
{
1
,
Φ
n
,
k
:
Φ
n
,
k
∈
S
n
*
(
a
,
b
)
,
1
-
S
⩽
k
⩽
2
n
S
-
1
}
.Remark 1.
The dimensions of the spaceS
n
0
(
a
,
b
) are 2
n
S
+
S
-
2, that of S
n
*
(
a
,
b
) and S
n
+
(
a
,
b
) are 2
n
S
+
S
-
1 each and of S
n
(
a
,
b
) is 2
n
S
+
S.


## 3. Problem Formulation

Here, we consider the variational formulation of the following parabolic IBVP:(24)
∂
u
∂
t
-
∂
∂
x
(
α
∂
u
∂
x
)
+
β
u
+
γ
∂
u
∂
x
=
f
(
x
,
t
)
,
a
<
x
<
b
,
t
>
0
,
with the initial condition:
(25)
u
(
x
,
0
)
=
g
(
x
)
,
a
<
x
<
b
,
and with the boundary conditions:
(26)
(
i
)
u
(
a
,
t
)
=
c
(
t
)
,
u
(
b
,
t
)
=
d
(
t
)
(
DBC
)
,
(
ii
)
u
(
a
,
t
)
=
c
(
t
)
,
∂
u
∂
x
(
b
,
t
)
=
d
(
t
)
(
MBC
1
)
,
(
iii
)
∂
u
∂
x
(
a
,
t
)
=
c
(
t
)
,
u
(
b
,
t
)
=
d
(
t
)
(
MBC
2
)
,
(
iv
)
∂
u
∂
x
(
a
,
t
)
=
c
(
t
)
,
∂
u
∂
x
(
b
,
t
)
=
d
(
t
)
(
NBC
)
,
t
>
0
.
Here, the boundary conditions where the (space) derivative of u is present are natural boundary conditions and the others are essential boundary conditions.Multiplying (24) by a function v
∈
H
0
1
(
a
,
b
), H
*
1
(
a
,
b
), H
+
1
(
a
,
b
), or H
1
(
a
,
b
) (according to boundary conditions (26)(i), (26)(ii), (26)(iii), and (26)(iv), resp.), integrating by parts with respect to x in (
a
,
b
), and using the natural boundary conditions, we get(i)
the weak form of problem (24)–(26)(i) as
(27)
∫
a
b
‍
(
∂
u
∂
t
v
+
α
∂
u
∂
x
d
v
d
x
+
β
u
v
+
γ
∂
u
∂
x
v
)
d
x
=
∫
a
b
‍
f
v
d
x
,
(ii)
the weak form of problem (24)–(26)(ii) as
(28)
∫
a
b
‍
(
∂
u
∂
t
v
+
α
∂
u
∂
x
d
v
d
x
+
β
u
v
+
γ
∂
u
∂
x
v
)
d
x
=
∫
a
b
‍
f
v
d
x
+
d
(
t
)
α
(
b
,
t
)
v
(
b
)
,
(iii)
the weak form of problem (24)–(26)(iii) as
(29)
∫
a
b
‍
(
∂
u
∂
t
v
+
α
∂
u
∂
x
d
v
d
x
+
β
u
v
+
γ
∂
u
∂
x
v
)
d
x
=
∫
a
b
‍
f
v
d
x
-
c
(
t
)
α
(
a
,
t
)
v
(
a
)
,and(iv)
the weak form of problem (24)–(26)(iv) as
(30)
∫
a
b
‍
(
∂
u
∂
t
v
+
α
∂
u
∂
x
d
v
d
x
+
β
u
v
+
γ
∂
u
∂
x
v
)
d
x
=
∫
a
b
‍
f
v
d
x
+
d
(
t
)
α
(
b
,
t
)
v
(
b
)
-
c
(
t
)
α
(
a
,
t
)
v
(
a
)
.


## 4. Spatial Approximation

In [12], it is established that N
⩾
2 is sufficient for the solution of second-order (in space) problems. Let N
⩾
2 be any integer and let ϕ be the d
b
N scaling function. Let H
n, n
⩾
0, be the approximating subspaces of the semidiscrete space approximations of the variational problems formulated above. Considering the bases {
Φ
n
,
j
} of H
n =S
n
0
(
a
,
b
), S
n
*
(
a
,
b
), S
n
+
(
a
,
b
), and S
n
(
a
,
b
) constructed in Section 2, the approximate solution of variational problems (27), (28), (29), and (30) can be taken as
(31)
u
n
(
x
,
t
)
=
u
0
(
x
,
t
)
+
∑
j
‍
z
n
,
j
(
t
)
Φ
n
,
j
(
x
)
,
whereu
0
(
x
,
t
)
=
(
(
b
c
(
t
)
-
a
d
(
t
)
)
/
(
b
-
a
)
) + (
(
d
(
t
)
-
c
(
t
)
)
/
(
b
-
a
)
)
x and H
n
=
S
n
0
(
a
,
b
), for problem (27),
u
0
(
x
,
t
)
=
c
(
t
)
-
a
+
x and H
n
=
S
n
*
(
a
,
b
), for problem (28),
u
0
(
x
,
t
)
=
d
(
t
)
-
b
+
x and H
n
=
S
n
+
(
a
,
b
), for problem (29),
u
0
(
x
,
t
)
=
0 and H
n
=
S
n
(
a
,
b
), for problem (30).Applying Galerkin method to each of the variational problems (27), (28), (29), and (30) with the approximate solution (31), we get a system of first-order ordinary differential equations in z
=
[
z
n
,
j
]:
(32)
M
z
˙
+
A
z
=
F
,
where A
=
[
A
i
j
] is the stiffness matrix given by
(33)
A
i
j
=
∫
a
b
‍
(
α
d
Φ
n
,
j
d
x
d
Φ
n
,
i
d
x
+
β
Φ
n
,
j
Φ
n
,
i
+
γ
d
Φ
n
,
j
d
x
Φ
n
,
i
)
d
x
,
M
=
[
M
i
j
] is the mass matrix given by
(34)
M
i
j
=
∫
a
b
‍
Φ
n
,
j
Φ
n
,
i
d
x
,
and F
=
[
F
i
] is the force vector given by(i)
F
i
=
∫
a
b
‍
f
Φ
n
,
i
d
x
-
K
i, for problem (27),
(ii)
F
i
=
∫
a
b
‍
f
Φ
n
,
i
d
x
+
d
(
t
)
α
(
b
,
t
)
Φ
n
,
i
(
b
)
-
K
i, for problem (28),
(iii)
F
i
=
∫
a
b
‍
f
Φ
n
,
i
d
x
-
c
(
t
)
α
(
a
,
t
)
Φ
n
,
i
(
a
)
-
K
i, for problem (29),
(iv)
F
i
=
∫
a
b
‍
f
Φ
n
,
i
d
x
+
d
(
t
)
α
(
b
,
t
)
Φ
n
,
i
(
b
)
-
c
(
t
)
α
(
a
,
t
)
Φ
n
,
i
(
a
)
-
K
i, for problem (30), where
(35)
K
i
=
∫
a
b
‍
(
∂
u
0
∂
t
Φ
n
,
i
+
α
∂
u
0
∂
x
d
Φ
n
,
i
d
x
+
β
u
0
Φ
n
,
i
+
γ
∂
u
0
∂
x
Φ
n
,
i
)
d
x
.
The solution of the system of ordinary differential equations (32) gives the complete solution of the original problem. The value of the vector z at t
=
0 is determined from the system of linear equations
(36)
M
z
(
0
)
=
G
,
obtained by multiplying (25) by v and integrating in (
a
,
b
) with respect to x, where the vector G
=
[
G
i
] is given by G
i
=
∫
a
b
‍
g
Φ
n
,
i
d
x.Thus, the IBVP (24)–(26) becomes a single initial value problem (IVP) given by (32) and (36).


## 5. Time Discretization

### 5.1. Theθ Scheme

The most commonly used method for the solution of IVP given by (32) and (36) is a θ
-family of approximation [16]:
(37)
θ
z
˙
s
+
1
+
(
1
-
θ
)
z
˙
s
=
z
s
+
1
-
z
s
Δ
t
,
0
⩽
θ
⩽
1
,
where z
s refers to the value of z at time t
=
t
s
=
s
Δ
t, and there is a similar meaning of z
˙
s. (In [16], θ is α.)Using approximation (37) for time t
s and t
s
+
1 in (32), we obtain
(38)
M
z
s
+
1
=
M
z
s
+
θ
Δ
t
(
F
s
+
1
-
A
s
+
1
z
s
+
1
)
+
(
1
-
θ
)
Δ
t
(
F
s
-
A
s
z
s
)
.
Rearranging the terms, we get
(39)
M
^
s
+
1
z
s
+
1
=
A
^
s
z
s
+
F
^
s
,
s
+
1
,
where
(40)
M
^
s
+
1
=
M
+
θ
Δ
t
A
s
+
1
,
A
^
s
=
M
-
(
1
-
θ
)
Δ
t
A
s
,
F
^
s
,
s
+
1
=
Δ
t
[
θ
F
s
+
1
+
(
1
-
θ
)
F
s
]
.
The solution z
s
+
1 at time t
s
+
1 is obtained in terms of the solution z
s at time t
s by inverting the matrix M
^
s
+
1. The solution z
0 at time t
=
t
0
=
0 is z
(
0
), obtained from (36).


### 5.2. Stability and Accuracy

Error is introduced into the solutionz
s
+
1 at each time step when approximation (37) is used. If the error grows unboundedly, the solution scheme is said to be unstable. Otherwise it is said to be stable. Accuracy of a numerical scheme is a measure of the closeness between the approximate and the exact solutions, whereas stability is a measure of the boundedness of the approximate solution with time. The size of the time step Δ
t can influence both the accuracy and the stability. A smaller time step results in an improved accuracy. A numerical scheme is said to be conditionally stable, if it is stable only when certain restrictions on the time step are satisfied. The above θ-family of time approximation schemes, for which θ
<
1
/
2, is stable only if
(41)
Δ
t
<
2
(
1
-
2
θ
)
λ
,
where λ is the largest eigenvalue of the problem associated with the original PDE problem.A numerical scheme is called anexplicit scheme, if the matrix M
^
s
+
1 in (39) is diagonal. Otherwise, it is called an implicit scheme.For different values ofθ, the time discretization scheme (37) is classified as follows.(1)
The forward difference(or Euler)scheme: θ
=
0; explicit and conditionally stable; order of accuracy = O
(
Δ
t
).
(2)
The backward difference scheme: θ
=
1; implicit and unconditionally stable; order of accuracy = O
(
Δ
t
).
(3)
The Galerkin scheme: θ
=
2
/
3; implicit and unconditionally stable; order of accuracy = O
(
(
Δ
t
)
2
).
(4)
The Crank-Nicholson scheme: θ
=
1
/
2; implicit and unconditionally stable; order of accuracy = O
(
(
Δ
t
)
2
).


## 5.1. Theθ Scheme

The most commonly used method for the solution of IVP given by (32) and (36) is a θ
-family of approximation [16]:
(37)
θ
z
˙
s
+
1
+
(
1
-
θ
)
z
˙
s
=
z
s
+
1
-
z
s
Δ
t
,
0
⩽
θ
⩽
1
,
where z
s refers to the value of z at time t
=
t
s
=
s
Δ
t, and there is a similar meaning of z
˙
s. (In [16], θ is α.)Using approximation (37) for time t
s and t
s
+
1 in (32), we obtain
(38)
M
z
s
+
1
=
M
z
s
+
θ
Δ
t
(
F
s
+
1
-
A
s
+
1
z
s
+
1
)
+
(
1
-
θ
)
Δ
t
(
F
s
-
A
s
z
s
)
.
Rearranging the terms, we get
(39)
M
^
s
+
1
z
s
+
1
=
A
^
s
z
s
+
F
^
s
,
s
+
1
,
where
(40)
M
^
s
+
1
=
M
+
θ
Δ
t
A
s
+
1
,
A
^
s
=
M
-
(
1
-
θ
)
Δ
t
A
s
,
F
^
s
,
s
+
1
=
Δ
t
[
θ
F
s
+
1
+
(
1
-
θ
)
F
s
]
.
The solution z
s
+
1 at time t
s
+
1 is obtained in terms of the solution z
s at time t
s by inverting the matrix M
^
s
+
1. The solution z
0 at time t
=
t
0
=
0 is z
(
0
), obtained from (36).


## 5.2. Stability and Accuracy

Error is introduced into the solutionz
s
+
1 at each time step when approximation (37) is used. If the error grows unboundedly, the solution scheme is said to be unstable. Otherwise it is said to be stable. Accuracy of a numerical scheme is a measure of the closeness between the approximate and the exact solutions, whereas stability is a measure of the boundedness of the approximate solution with time. The size of the time step Δ
t can influence both the accuracy and the stability. A smaller time step results in an improved accuracy. A numerical scheme is said to be conditionally stable, if it is stable only when certain restrictions on the time step are satisfied. The above θ-family of time approximation schemes, for which θ
<
1
/
2, is stable only if
(41)
Δ
t
<
2
(
1
-
2
θ
)
λ
,
where λ is the largest eigenvalue of the problem associated with the original PDE problem.A numerical scheme is called anexplicit scheme, if the matrix M
^
s
+
1 in (39) is diagonal. Otherwise, it is called an implicit scheme.For different values ofθ, the time discretization scheme (37) is classified as follows.(1)
The forward difference(or Euler)scheme: θ
=
0; explicit and conditionally stable; order of accuracy = O
(
Δ
t
).
(2)
The backward difference scheme: θ
=
1; implicit and unconditionally stable; order of accuracy = O
(
Δ
t
).
(3)
The Galerkin scheme: θ
=
2
/
3; implicit and unconditionally stable; order of accuracy = O
(
(
Δ
t
)
2
).
(4)
The Crank-Nicholson scheme: θ
=
1
/
2; implicit and unconditionally stable; order of accuracy = O
(
(
Δ
t
)
2
).


## 6. Numerical Experiments

Here the methodology for the solution of IBVP (24)–(26) described above is demonstrated with examples of diffusion, diffusion-reaction, convection-diffusion, and convection-diffusion-reaction equations with all the four types of boundary conditions, using all the four time discretization schemes. We compare the computed solutions with the exact solutions and obtain L
∞
,
L
2 and H
1 norm errors. The stiffness matrices for diffusion and diffusion-reaction equations are symmetric and those for convection-diffusion and convection-diffusion-reaction equations are not symmetric. All the systems of linear equations are solved by using Gaussian elimination method and the integrations are carried out by Simpson’s quadratures. All the computations are performed in MATLAB.


### 6.1. Test Problem (42)

Here we solve a diffusion equation (β
=
γ
=
0) with the first mixed boundary conditions (MBC1). Here we take α
=
1, the space domain [
0,1
], the initial value function g
(
x
)
=
1, and the boundary value functions c
(
t
)
=
0
,
d
(
t
)
=
0. If we solve this problem by the method of separation of variables, we get the exact solution as
(42)
u
(
x
,
t
)
=
2
∑
n
=
0
∞
‍
e
-
μ
n
2
t
sin
μ
n
x
μ
n
,
μ
n
=
(
2
n
+
1
)
π
2
.Forθ
=
0, d
b
2
(
n
=
0
) wavelet is used for spatial discretization. As this scheme is conditionally stable, we have to find an upper limit of the time step Δ
t using the stability condition (41). The largest eigenvalue of the associated problem is 526.47. Therefore, the maximum time step is 2
/
526.47
≈
0.0038. Figure 2 shows the exact solution and d
b
2
(
n
=
0
) wavelet solution for Δ
t
=
0.004 (unstable) and Δ
t
=
0.003 (stable) at x
=
1.Test problem (42) (
θ
=
0
). d
b
2
(
n
=
0
) wavelet solution (
-
) and exact solution (
-
-
) at x
=
1.
(a)
(b)Forθ
=
1, d
b
2
(
n
=
1
) wavelet is used for spatial discretization. As this scheme is unconditionally stable, there is no restriction on Δ
t. However, to obtain a sufficiently accurate solution, Δ
t must be taken small enough. Figure 3 shows the exact solution and the computed solutions at time t
=
1 with Δ
t
= 0.5, 0.25, 0.125, and 0.0625. Here it can be seen that the computed solutions converge to the exact solution as Δ
t
→
0.Figure 3
Test problem (42) (θ
=
1). Convergence of d
b2(n
=
1
) wavelet solutions with decreasing time steps at t
=
1.Forθ
=
2
/
3, d
b
2
(
n
=
2
) wavelet is used for spatial discretization. This scheme is also unconditionally stable. Figure 4 shows the computed solutions at times 5
Δ
t, 10
Δ
t, 15
Δ
t, 20
Δ
t, and 25
Δ
t with Δ
t
=0.05.Figure 4
Test problem (42) (θ
=
2
/
3, Δ
t
=
0.05). d
b
2
(
n
=
2
) wavelet solutions at 5, 10, 15, 20, and 25 time steps.Forθ
=
1
/
2, we use d
b
2 and d
b
3 wavelets for spatial discretization. Table 1 shows L
∞, L
2, and H
1 norm errors at t
=
1 for n
=
0,1
,
2,3 with various values of Δ
t.Table 1
Test problem (42) (θ
=
1
/
2), L
∞
,
L
2, and H
1 errors with d
b2 and d
b3 wavelets at t
=
1.
Wavelets
n
Δ
t
L
∞ error
L
2 error
H
1 error
d
b2
0
0.01
3.7184
×
1
0
-
4
1.9334
×
1
0
-
4
2.6769
×
1
0
-
3
1
0.01
4.9078
×
1
0
-
5
2.7388
×
1
0
-
5
7.9817
×
1
0
-
4
2
0.005
6.7404
×
1
0
-
6
3.9984
×
1
0
-
6
2.1403
×
1
0
-
4
3
0.001
8.4199
×
1
0
-
7
3.8148
×
1
0
-
7
5.5207
×
1
0
-
5
d
b3
0
0.01
2.4536
×
1
0
-
5
1.0938
×
1
0
-
5
1.4991
×
1
0
-
4
1
0.005
4.1251
×
1
0
-
6
2.4336
×
1
0
-
6
2.5853
×
1
0
-
5
2
0.001
1.8832
×
1
0
-
7
1.0155
×
1
0
-
7
2.7747
×
1
0
-
6
3
0.0005
3.8261
×
1
0
-
8
2.4823
×
1
0
-
8
3.5706
×
1
0
-
7


### 6.2. Test Problem (43)

Here we solve a diffusion-reaction equation (γ
=
0) with the second mixed boundary conditions (MBC2). Here we take α
=
1
,
β
=
2
cos
⁡
x, and the space domain [
0
,
π
]. The right-hand function f
(
x
,
t
), the initial value function g
(
x
), and the boundary value functions c
(
t
) and d
(
t
) are obtained as per the exact solution:
(43)
u
(
x
,
t
)
=
e
-
t
sin
x
.Forθ
=
0, d
b
2
(
n
=
0
) wavelet is used for spatial discretization. This scheme is conditionally stable. The largest eigenvalue of the associated problem is 235.57. Therefore, the maximum time step is 2
/
235.57
≈
0.0085. Figure 5 shows the exact solution and d
b
2
(
n
=
0
) wavelet solution for Δ
t
=
0.009 (unstable) and Δ
t
=
0.008 (stable) at x
=
π
/
2.Test problem (43) (θ
=
0). d
b
2
(
n
=
0
) wavelet solution (
-
) and exact solution (
-
-
) at x
=
π
/
2.
(a)
(b)Forθ
=
1, d
b
2
(
n
=
1
) wavelet is used for spatial discretization. Figure 6 shows the exact and the computed solutions at time t
=
1 with Δ
t
= 0.5, 0.25, 0.125, and 0.0625. Here also, it can be seen that the computed solutions converge to the exact solution as Δ
t
→
0.Figure 6
Test problem (43) (θ
=
1). Convergence of d
b
2
(
n
=
1
) wavelet solutions with decreasing time steps at t
=
1.Forθ
=
2
/
3, d
b
2
(
n
=
2
) wavelet is used for spatial discretization. Figure 7 shows the computed solutions at times 5
Δ
t, 10
Δ
t, 15
Δ
t, 20
Δ
t, and 25
Δ
t with Δ
t
= 0.05.Figure 7
Test problem (43) (θ
=
2
/
3, Δ
t
=
0.05). d
b
2
(
n
=
2
) wavelet solutions at 5, 10, 15, 20, and 25 time steps.Forθ
=
1
/
2, we use d
b
2 and d
b
3 wavelets for spatial discretization. Table 2 shows L
∞
,
L
2, and H
1 norm errors at t
=
1 for n
=
0,1
,
2,3 with various values of Δ
t.Table 2
Test problem (43) (θ
=
1
/
2), L
∞
,
L
2, and H
1 errors with d
b2 and d
b3 wavelets at t
=
1.
Wavelets
n
Δ
t
L
∞ error
L
2 error
H
1 error
d
b2
0
0.01
9.6415
×
1
0
-
3
9.0609
×
1
0
-
3
3.4684
×
1
0
-
2
1
0.01
1.6168
×
1
0
-
3
1.2191
×
1
0
-
3
1.1454
×
1
0
-
2
2
0.005
2.1436
×
1
0
-
4
1.4578
×
1
0
-
4
3.1952
×
1
0
-
3
3
0.001
2.7194
×
1
0
-
5
1.7754
×
1
0
-
5
8.3777
×
1
0
-
4
d
b3
0
0.01
6.9866
×
1
0
-
4
5.4571
×
1
0
-
4
5.0691
×
1
0
-
3
1
0.005
4.0701
×
1
0
-
5
3.5001
×
1
0
-
5
6.8968
×
1
0
-
4
2
0.001
2.5157
×
1
0
-
6
2.1889
×
1
0
-
6
8.7327
×
1
0
-
5
3
0.0005
1.6789
×
1
0
-
7
1.3807
×
1
0
-
7
1.0964
×
1
0
-
5


### 6.3. Test Problem (44)

Here we solve a convection-diffusion equation (β
=
0) with Dirichlet boundary conditions (DBC). Here we take α
=
1
/
16
,
β
=
1
/
4, and the space domain [
0,1
]. The functions f
(
x
,
t
)
,
g
(
x
)
,
c
(
t
), and d
(
t
) are obtained as per the exact solution:
(44)
u
(
x
,
t
)
=
2
e
-
(
π
2
/
4
)
t
sin
(
2
π
x
)
+
8
x
2
t
.Forθ
=
0, d
b
2
(
n
=
0
) and d
b
3
(
n
=
0
) wavelets are used for spatial discretization. This scheme is conditionally stable. The largest eigenvalues of the associated problem for d
b
2
(
n
=
0
) and d
b
3
(
n
=
0
) wavelets are 32.64 and 190.68, respectively. Therefore, the maximum time steps for these wavelets are 2
/
32.64
≈
0.0613, and 2
/
190.68
≈
0.0105 respectively. Figure 8 shows the exact, unstable, and stable solutions due to d
b
2
(
n
=
0
) and d
b
3
(
n
=
0
) wavelets at x
=
1
/
2.Test problem (44) (θ
=
0
).  d
b2(n
=
0), d
b
3
(
n
=
0
) wavelet solutions (
-
) and exact solutions (
-
-
) at x
=
1
/
2.
(a)
(b)
(c)
(d)Forθ
=
1, d
b
2
(
n
=
1
) wavelet is used for spatial discretization. Figure 9 shows the exact and the computed solutions at time t
=
0.5 with Δ
t
= 0.5, 0.25, 0.125, and 0.0625. Here also, it can be seen that the computed solutions converge to the exact solution as Δ
t
→
0.Figure 9
Test problem (44) (θ
=
1). Convergence of d
b
2
(
n
=
1
) wavelet solutions with decreasing time steps at t
=
0.5.Forθ
=
2
/
3, d
b
2
(
n
=
2
) wavelet is used for spatial discretization. Figure 10 shows the computed solutions at times 5
Δ
t, 10
Δ
t, 15
Δ
t, 20
Δ
t, and 25
Δ
t with Δ
t
= 0.05.Figure 10
Test problem (44) (θ
=
2
/
3, Δ
t
=
0.05). d
b
2
(
n
=
2
) wavelet solutions at 5, 10, 15, 20, and 25 time steps.Forθ
=
1
/
2, we use d
b
2 and d
b
3 wavelets for spatial discretization. Table 3 shows L
∞
,
L
2, and H
1 norm errors at t
=
1 for n
=
0,1
,
2,3 with various values of Δ
t.Table 3
Test problem (44) (θ
=
1
/
2), L
∞
,
L
2, and H
1 errors with d
b2 and d
b3 wavelets at t
=
1.
Wavelets
n
Δ
t
L
∞ error
L
2 error
H
1 error
d
b2
0
0.01
6.7115
×
1
0
-
2
2.8666
×
1
0
-
2
2.5086
×
1
0
-
1
1
0.01
7.2533
×
1
0
-
3
3.3632
×
1
0
-
3
7.4574
×
1
0
-
2
2
0.005
7.7255
×
1
0
-
4
3.5252
×
1
0
-
4
2.0785
×
1
0
-
2
3
0.001
9.3833
×
1
0
-
5
3.8914
×
1
0
-
5
5.4701
×
1
0
-
3
d
b3
0
0.01
5.5939
×
1
0
-
3
2.3186
×
1
0
-
3
5.1679
×
1
0
-
2
1
0.005
3.6101
×
1
0
-
4
1.4803
×
1
0
-
4
8.5848
×
1
0
-
3
2
0.001
1.8615
×
1
0
-
5
9.1200
×
1
0
-
6
1.1315
×
1
0
-
3
3
0.0005
1.1507
×
1
0
-
6
5.6949
×
1
0
-
7
1.4275
×
1
0
-
4


### 6.4. Test Problem (45)

Here we solve a convection-diffusion-reaction equation with Neumann boundary conditions (NBC). Here we takeα
=
1
/
2
,
β
=
x, γ
=
2
-
x
2, and the space domain [
1,2
]. The functions f
(
x
,
t
)
,
g
(
x
)
,
c
(
t
), and d
(
t
) are obtained as per the exact solution
(45)
u
(
x
,
t
)
=
e
2
x
-
t
.Forθ
=
0, d
b
2
(
n
=
0
) and d
b
3
(
n
=
0
) wavelets are used for spatial discretization. This scheme is conditionally stable. The largest eigenvalues of the associated problem for d
b
2
(
n
=
0
) and d
b
3
(
n
=
0
) wavelets are 1124.50 and 4092.00, respectively. Therefore, the maximum time steps for these wavelets are 2
/
1124.50
≈
0.0018 and 2
/
4092.00
≈
0.00049, respectively. Figure 11 shows the exact, unstable, and stable solutions due to d
b
2
(
n
=
0
) wavelets at x
=
2 and due to d
b
3
(
n
=
0
) wavelets at x
=
1.Test problem (45) (
θ
=
0
).  d
b
2
(
n
=
0
), d
b
3
(
n
=
0
) wavelet solutions (
-
) and exact solutions (
-
-
).
(a)
(b)
(c)
(d)Forθ
=
1, d
b
2
(
n
=
1
) wavelet is used for spatial discretization. Figure 12 shows the exact and the computed solutions at time t
=
2 with Δ
t
= 0.5, 0.25, 0.125, and 0.0625. Here also, it can be seen that the computed solutions converge to the exact solution as Δ
t
→
0.Figure 12
Test problem (45) (
θ
=
1
). Convergence of d
b
2
(
n
=
1
) wavelet solutions with decreasing time steps at t
=
2.Forθ
=
2
/
3, d
b
2
(
n
=
2
) wavelet is used for spatial discretization. Figure 13 shows the computed solutions at times 5
Δ
t, 10
Δ
t, 15
Δ
t, 20
Δ
t, and 25
Δ
t, with Δ
t
= 0.05.Figure 13
Test Problem (45) (θ
=
2
/
3, Δ
t
=
0.05). d
b
2
(
n
=
2
) wavelet solutions at 5, 10, 15, 20, and 25 time steps.Forθ
=
1
/
2, we use d
b
2 and d
b
3 wavelets for spatial discretization. Table 4 shows L
∞
,
L
2, and H
1 norm errors at t
=
1 for n
=
0,1
,
2,3 with various values of Δ
t.Table 4
Test problem (45) (θ
=
1
/
2), L
∞
,
L
2, and H
1 errors with d
b2 and d
b3 wavelets at t
=
1.
Wavelets
n
Δ
t
L
∞ error
L
2 error
H
1 error
d
b2
0
0.01
2.7940
×
1
0
-
1
1.4914
×
1
0
-
1
8.3377
×
1
0
-
1
1
0.01
3.3430
×
1
0
-
2
1.4908
×
1
0
-
2
2.2252
×
1
0
-
1
2
0.005
3.6584
×
1
0
-
3
1.3363
×
1
0
-
3
5.8339
×
1
0
-
2
3
0.001
3.8989
×
1
0
-
4
1.2486
×
1
0
-
4
1.4941
×
1
0
-
2
d
b3
0
0.01
5.7558
×
1
0
-
3
2.2151
×
1
0
-
3
4.4906
×
1
0
-
2
1
0.005
3.7623
×
1
0
-
4
1.2626
×
1
0
-
4
6.6206
×
1
0
-
3
2
0.001
2.5667
×
1
0
-
5
8.2674
×
1
0
-
6
9.0062
×
1
0
-
4
3
0.0005
2.4634
×
1
0
-
6
9.4945
×
1
0
-
7
1.1864
×
1
0
-
4


## 6.1. Test Problem (42)

Here we solve a diffusion equation (β
=
γ
=
0) with the first mixed boundary conditions (MBC1). Here we take α
=
1, the space domain [
0,1
], the initial value function g
(
x
)
=
1, and the boundary value functions c
(
t
)
=
0
,
d
(
t
)
=
0. If we solve this problem by the method of separation of variables, we get the exact solution as
(42)
u
(
x
,
t
)
=
2
∑
n
=
0
∞
‍
e
-
μ
n
2
t
sin
μ
n
x
μ
n
,
μ
n
=
(
2
n
+
1
)
π
2
.Forθ
=
0, d
b
2
(
n
=
0
) wavelet is used for spatial discretization. As this scheme is conditionally stable, we have to find an upper limit of the time step Δ
t using the stability condition (41). The largest eigenvalue of the associated problem is 526.47. Therefore, the maximum time step is 2
/
526.47
≈
0.0038. Figure 2 shows the exact solution and d
b
2
(
n
=
0
) wavelet solution for Δ
t
=
0.004 (unstable) and Δ
t
=
0.003 (stable) at x
=
1.Test problem (42) (
θ
=
0
). d
b
2
(
n
=
0
) wavelet solution (
-
) and exact solution (
-
-
) at x
=
1.
(a)
(b)Forθ
=
1, d
b
2
(
n
=
1
) wavelet is used for spatial discretization. As this scheme is unconditionally stable, there is no restriction on Δ
t. However, to obtain a sufficiently accurate solution, Δ
t must be taken small enough. Figure 3 shows the exact solution and the computed solutions at time t
=
1 with Δ
t
= 0.5, 0.25, 0.125, and 0.0625. Here it can be seen that the computed solutions converge to the exact solution as Δ
t
→
0.Figure 3
Test problem (42) (θ
=
1). Convergence of d
b2(n
=
1
) wavelet solutions with decreasing time steps at t
=
1.Forθ
=
2
/
3, d
b
2
(
n
=
2
) wavelet is used for spatial discretization. This scheme is also unconditionally stable. Figure 4 shows the computed solutions at times 5
Δ
t, 10
Δ
t, 15
Δ
t, 20
Δ
t, and 25
Δ
t with Δ
t
=0.05.Figure 4
Test problem (42) (θ
=
2
/
3, Δ
t
=
0.05). d
b
2
(
n
=
2
) wavelet solutions at 5, 10, 15, 20, and 25 time steps.Forθ
=
1
/
2, we use d
b
2 and d
b
3 wavelets for spatial discretization. Table 1 shows L
∞, L
2, and H
1 norm errors at t
=
1 for n
=
0,1
,
2,3 with various values of Δ
t.Table 1
Test problem (42) (θ
=
1
/
2), L
∞
,
L
2, and H
1 errors with d
b2 and d
b3 wavelets at t
=
1.
Wavelets
n
Δ
t
L
∞ error
L
2 error
H
1 error
d
b2
0
0.01
3.7184
×
1
0
-
4
1.9334
×
1
0
-
4
2.6769
×
1
0
-
3
1
0.01
4.9078
×
1
0
-
5
2.7388
×
1
0
-
5
7.9817
×
1
0
-
4
2
0.005
6.7404
×
1
0
-
6
3.9984
×
1
0
-
6
2.1403
×
1
0
-
4
3
0.001
8.4199
×
1
0
-
7
3.8148
×
1
0
-
7
5.5207
×
1
0
-
5
d
b3
0
0.01
2.4536
×
1
0
-
5
1.0938
×
1
0
-
5
1.4991
×
1
0
-
4
1
0.005
4.1251
×
1
0
-
6
2.4336
×
1
0
-
6
2.5853
×
1
0
-
5
2
0.001
1.8832
×
1
0
-
7
1.0155
×
1
0
-
7
2.7747
×
1
0
-
6
3
0.0005
3.8261
×
1
0
-
8
2.4823
×
1
0
-
8
3.5706
×
1
0
-
7


## 6.2. Test Problem (43)

Here we solve a diffusion-reaction equation (γ
=
0) with the second mixed boundary conditions (MBC2). Here we take α
=
1
,
β
=
2
cos
⁡
x, and the space domain [
0
,
π
]. The right-hand function f
(
x
,
t
), the initial value function g
(
x
), and the boundary value functions c
(
t
) and d
(
t
) are obtained as per the exact solution:
(43)
u
(
x
,
t
)
=
e
-
t
sin
x
.Forθ
=
0, d
b
2
(
n
=
0
) wavelet is used for spatial discretization. This scheme is conditionally stable. The largest eigenvalue of the associated problem is 235.57. Therefore, the maximum time step is 2
/
235.57
≈
0.0085. Figure 5 shows the exact solution and d
b
2
(
n
=
0
) wavelet solution for Δ
t
=
0.009 (unstable) and Δ
t
=
0.008 (stable) at x
=
π
/
2.Test problem (43) (θ
=
0). d
b
2
(
n
=
0
) wavelet solution (
-
) and exact solution (
-
-
) at x
=
π
/
2.
(a)
(b)Forθ
=
1, d
b
2
(
n
=
1
) wavelet is used for spatial discretization. Figure 6 shows the exact and the computed solutions at time t
=
1 with Δ
t
= 0.5, 0.25, 0.125, and 0.0625. Here also, it can be seen that the computed solutions converge to the exact solution as Δ
t
→
0.Figure 6
Test problem (43) (θ
=
1). Convergence of d
b
2
(
n
=
1
) wavelet solutions with decreasing time steps at t
=
1.Forθ
=
2
/
3, d
b
2
(
n
=
2
) wavelet is used for spatial discretization. Figure 7 shows the computed solutions at times 5
Δ
t, 10
Δ
t, 15
Δ
t, 20
Δ
t, and 25
Δ
t with Δ
t
= 0.05.Figure 7
Test problem (43) (θ
=
2
/
3, Δ
t
=
0.05). d
b
2
(
n
=
2
) wavelet solutions at 5, 10, 15, 20, and 25 time steps.Forθ
=
1
/
2, we use d
b
2 and d
b
3 wavelets for spatial discretization. Table 2 shows L
∞
,
L
2, and H
1 norm errors at t
=
1 for n
=
0,1
,
2,3 with various values of Δ
t.Table 2
Test problem (43) (θ
=
1
/
2), L
∞
,
L
2, and H
1 errors with d
b2 and d
b3 wavelets at t
=
1.
Wavelets
n
Δ
t
L
∞ error
L
2 error
H
1 error
d
b2
0
0.01
9.6415
×
1
0
-
3
9.0609
×
1
0
-
3
3.4684
×
1
0
-
2
1
0.01
1.6168
×
1
0
-
3
1.2191
×
1
0
-
3
1.1454
×
1
0
-
2
2
0.005
2.1436
×
1
0
-
4
1.4578
×
1
0
-
4
3.1952
×
1
0
-
3
3
0.001
2.7194
×
1
0
-
5
1.7754
×
1
0
-
5
8.3777
×
1
0
-
4
d
b3
0
0.01
6.9866
×
1
0
-
4
5.4571
×
1
0
-
4
5.0691
×
1
0
-
3
1
0.005
4.0701
×
1
0
-
5
3.5001
×
1
0
-
5
6.8968
×
1
0
-
4
2
0.001
2.5157
×
1
0
-
6
2.1889
×
1
0
-
6
8.7327
×
1
0
-
5
3
0.0005
1.6789
×
1
0
-
7
1.3807
×
1
0
-
7
1.0964
×
1
0
-
5


## 6.3. Test Problem (44)

Here we solve a convection-diffusion equation (β
=
0) with Dirichlet boundary conditions (DBC). Here we take α
=
1
/
16
,
β
=
1
/
4, and the space domain [
0,1
]. The functions f
(
x
,
t
)
,
g
(
x
)
,
c
(
t
), and d
(
t
) are obtained as per the exact solution:
(44)
u
(
x
,
t
)
=
2
e
-
(
π
2
/
4
)
t
sin
(
2
π
x
)
+
8
x
2
t
.Forθ
=
0, d
b
2
(
n
=
0
) and d
b
3
(
n
=
0
) wavelets are used for spatial discretization. This scheme is conditionally stable. The largest eigenvalues of the associated problem for d
b
2
(
n
=
0
) and d
b
3
(
n
=
0
) wavelets are 32.64 and 190.68, respectively. Therefore, the maximum time steps for these wavelets are 2
/
32.64
≈
0.0613, and 2
/
190.68
≈
0.0105 respectively. Figure 8 shows the exact, unstable, and stable solutions due to d
b
2
(
n
=
0
) and d
b
3
(
n
=
0
) wavelets at x
=
1
/
2.Test problem (44) (θ
=
0
).  d
b2(n
=
0), d
b
3
(
n
=
0
) wavelet solutions (
-
) and exact solutions (
-
-
) at x
=
1
/
2.
(a)
(b)
(c)
(d)Forθ
=
1, d
b
2
(
n
=
1
) wavelet is used for spatial discretization. Figure 9 shows the exact and the computed solutions at time t
=
0.5 with Δ
t
= 0.5, 0.25, 0.125, and 0.0625. Here also, it can be seen that the computed solutions converge to the exact solution as Δ
t
→
0.Figure 9
Test problem (44) (θ
=
1). Convergence of d
b
2
(
n
=
1
) wavelet solutions with decreasing time steps at t
=
0.5.Forθ
=
2
/
3, d
b
2
(
n
=
2
) wavelet is used for spatial discretization. Figure 10 shows the computed solutions at times 5
Δ
t, 10
Δ
t, 15
Δ
t, 20
Δ
t, and 25
Δ
t with Δ
t
= 0.05.Figure 10
Test problem (44) (θ
=
2
/
3, Δ
t
=
0.05). d
b
2
(
n
=
2
) wavelet solutions at 5, 10, 15, 20, and 25 time steps.Forθ
=
1
/
2, we use d
b
2 and d
b
3 wavelets for spatial discretization. Table 3 shows L
∞
,
L
2, and H
1 norm errors at t
=
1 for n
=
0,1
,
2,3 with various values of Δ
t.Table 3
Test problem (44) (θ
=
1
/
2), L
∞
,
L
2, and H
1 errors with d
b2 and d
b3 wavelets at t
=
1.
Wavelets
n
Δ
t
L
∞ error
L
2 error
H
1 error
d
b2
0
0.01
6.7115
×
1
0
-
2
2.8666
×
1
0
-
2
2.5086
×
1
0
-
1
1
0.01
7.2533
×
1
0
-
3
3.3632
×
1
0
-
3
7.4574
×
1
0
-
2
2
0.005
7.7255
×
1
0
-
4
3.5252
×
1
0
-
4
2.0785
×
1
0
-
2
3
0.001
9.3833
×
1
0
-
5
3.8914
×
1
0
-
5
5.4701
×
1
0
-
3
d
b3
0
0.01
5.5939
×
1
0
-
3
2.3186
×
1
0
-
3
5.1679
×
1
0
-
2
1
0.005
3.6101
×
1
0
-
4
1.4803
×
1
0
-
4
8.5848
×
1
0
-
3
2
0.001
1.8615
×
1
0
-
5
9.1200
×
1
0
-
6
1.1315
×
1
0
-
3
3
0.0005
1.1507
×
1
0
-
6
5.6949
×
1
0
-
7
1.4275
×
1
0
-
4


## 6.4. Test Problem (45)

Here we solve a convection-diffusion-reaction equation with Neumann boundary conditions (NBC). Here we takeα
=
1
/
2
,
β
=
x, γ
=
2
-
x
2, and the space domain [
1,2
]. The functions f
(
x
,
t
)
,
g
(
x
)
,
c
(
t
), and d
(
t
) are obtained as per the exact solution
(45)
u
(
x
,
t
)
=
e
2
x
-
t
.Forθ
=
0, d
b
2
(
n
=
0
) and d
b
3
(
n
=
0
) wavelets are used for spatial discretization. This scheme is conditionally stable. The largest eigenvalues of the associated problem for d
b
2
(
n
=
0
) and d
b
3
(
n
=
0
) wavelets are 1124.50 and 4092.00, respectively. Therefore, the maximum time steps for these wavelets are 2
/
1124.50
≈
0.0018 and 2
/
4092.00
≈
0.00049, respectively. Figure 11 shows the exact, unstable, and stable solutions due to d
b
2
(
n
=
0
) wavelets at x
=
2 and due to d
b
3
(
n
=
0
) wavelets at x
=
1.Test problem (45) (
θ
=
0
).  d
b
2
(
n
=
0
), d
b
3
(
n
=
0
) wavelet solutions (
-
) and exact solutions (
-
-
).
(a)
(b)
(c)
(d)Forθ
=
1, d
b
2
(
n
=
1
) wavelet is used for spatial discretization. Figure 12 shows the exact and the computed solutions at time t
=
2 with Δ
t
= 0.5, 0.25, 0.125, and 0.0625. Here also, it can be seen that the computed solutions converge to the exact solution as Δ
t
→
0.Figure 12
Test problem (45) (
θ
=
1
). Convergence of d
b
2
(
n
=
1
) wavelet solutions with decreasing time steps at t
=
2.Forθ
=
2
/
3, d
b
2
(
n
=
2
) wavelet is used for spatial discretization. Figure 13 shows the computed solutions at times 5
Δ
t, 10
Δ
t, 15
Δ
t, 20
Δ
t, and 25
Δ
t, with Δ
t
= 0.05.Figure 13
Test Problem (45) (θ
=
2
/
3, Δ
t
=
0.05). d
b
2
(
n
=
2
) wavelet solutions at 5, 10, 15, 20, and 25 time steps.Forθ
=
1
/
2, we use d
b
2 and d
b
3 wavelets for spatial discretization. Table 4 shows L
∞
,
L
2, and H
1 norm errors at t
=
1 for n
=
0,1
,
2,3 with various values of Δ
t.Table 4
Test problem (45) (θ
=
1
/
2), L
∞
,
L
2, and H
1 errors with d
b2 and d
b3 wavelets at t
=
1.
Wavelets
n
Δ
t
L
∞ error
L
2 error
H
1 error
d
b2
0
0.01
2.7940
×
1
0
-
1
1.4914
×
1
0
-
1
8.3377
×
1
0
-
1
1
0.01
3.3430
×
1
0
-
2
1.4908
×
1
0
-
2
2.2252
×
1
0
-
1
2
0.005
3.6584
×
1
0
-
3
1.3363
×
1
0
-
3
5.8339
×
1
0
-
2
3
0.001
3.8989
×
1
0
-
4
1.2486
×
1
0
-
4
1.4941
×
1
0
-
2
d
b3
0
0.01
5.7558
×
1
0
-
3
2.2151
×
1
0
-
3
4.4906
×
1
0
-
2
1
0.005
3.7623
×
1
0
-
4
1.2626
×
1
0
-
4
6.6206
×
1
0
-
3
2
0.001
2.5667
×
1
0
-
5
8.2674
×
1
0
-
6
9.0062
×
1
0
-
4
3
0.0005
2.4634
×
1
0
-
6
9.4945
×
1
0
-
7
1.1864
×
1
0
-
4


## 7. Conclusion

In this paper, we have derived a method for numerical solution of linear parabolic partial differential equations in one space dimension. In this algorithm, we apply wavelet-Galerkin method in the space direction and some classical finite difference schemes in the time direction. The basis functions for wavelet-Galerkin method have been constructed by integrating Daubechies scaling functions as in [12]. Indeed, wavelet functions can also be used instead of scaling functions [13]. The most important merits of this algorithm are that it is useful for problems with mixed and Neumann boundary conditions besides Dirichlet one. Also the boundary conditions and the source term can be time-dependent. The method gives high favourable accuracy when compared to the exact solutions for problems of this types. The matter here left to see is what will happen if wavelets are used to discretize both the space and the time directions. The extension of the algorithm in multidimension is the matter of further research.



---

*Source: 346731-2014-02-27.xml*