# Solving Three Dimensional and Time Depending PDEs by Haar Wavelets Method

Abdeljalil Nachaoui, Ekhlass Al-Rawi, Ahmed Qasim

► To cite this version:

Abdeljalil Nachaoui, Ekhlass Al-Rawi, Ahmed Qasim. Solving Three Dimensional and Time Depending PDEs by Haar Wavelets Method. OALib Journal, 2018, 05 (05), pp.1-18. 10.4236/oalib.1104496. hal-02428861

HAL Id: hal-02428861

https://hal.science/hal-02428861

Submitted on 30 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Solving Three Dimensional and Time Depending PDEs by Haar Wavelets Method

Abdeljalil Nachaoui, Ekhlass Al-Rawi, Ahmed Qasim

► To cite this version:

Abdeljalil Nachaoui, Ekhlass Al-Rawi, Ahmed Qasim. Solving Three Dimensional and Time Depending PDEs by Haar Wavelets Method. OALib Journal, 2018, 05 (05), pp.1-18. 10.4236/oalib.1104496. hal-02428861

HAL Id: hal-02428861

https://hal.science/hal-02428861v1

Submitted on 30 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Solving Three Dimensional and Time Depending PDEs by Haar Wavelets Method

Abdeljalil Nachaoui<fnref n="1" />, Ekhlass S. Al-Rawi<fnref n="2" />, Ahmed F. Qasim<fnref n="2" />*

**How to cite this paper:** Nachaoui, A., Al-Rawi, E.S. and Qasim, A.F. (2018) Solving Three Dimensional and Time Depending PDEs by Haar Wavelets Method. *Open Access Library Journal*, 5: e4496.

https://doi.org/10.4236/oalib.1104496

Received: March 12, 2018

Accepted: May 21, 2018

Published: May 24, 2018

Copyright © 2018 by authors and Open Access Library Inc.
This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

<figure><img src="image_2.png" /></figure>

Open Access

## Abstract

Haar wavelets are applied for solution of three dimensional partial differential equations (PDEs) or time depending two dimensional PDEs. The proposed method is mathematically simple and fast. Two techniques are used in numerical solution, the first based on 2D-Haar wavelets and the second based on 3D-Haar wavelets and we compare them. To demonstrate the efficiency of the method, two test problems (solution of the diffusion and Poisson equations) are discussed. Computer simulation showed that 3D-Haar wavelets are better and closer to the exact solution but it is need to more time from 2D-Haar wavelets.

## Subject Areas

Numerical Mathematics, Partial Differential Equation

## Keywords

Haar Wavelets, Partial Differential Equations, Diffusion Equation, Poisson Equation

## 1. Introduction

As a powerful mathematical tool, Wavelet analysis has been widely used in image digital processing, quantum field theory, numerical analysis and many other fields in recent years.

Haar wavelets have been applied extensively for signal processing in communications and physics research, and more mathematically focused on differential equations and even nonlinear problems. After discrediting the differential equation in a conventional way like the finite difference approximation, wavelets can be used for algebraic manipulations in the system of equations obtained which may lead to better condition number of the resulting system [1].

<footnote n="1">¹Jean Leary Laboratory of Mathematics, University of Nantes, Nantes, France</footnote>

<footnote n="2">²College of Computer Sciences and Mathematics, University of Mosul, Mosul, Republic of Iraq
Email: Abdeljalil.Nachaoui@univ-nantes.fr, ekhlass70@yahoo.com, *ahmednumerical@yahoo.com</footnote>
# Solving Three Dimensional and Time Depending PDEs by Haar Wavelets Method

Abdeljalil Nachaoui<fnref n="1" />, Ekhlass S. Al-Rawi<fnref n="2" />, Ahmed F. Qasim<fnref n="2" />*

**How to cite this paper:** Nachaoui, A., Al-Rawi, E.S. and Qasim, A.F. (2018) Solving Three Dimensional and Time Depending PDEs by Haar Wavelets Method. *Open Access Library Journal*, 5: e4496.

https://doi.org/10.4236/oalib.1104496

Received: March 12, 2018

Accepted: May 21, 2018

Published: May 24, 2018

Copyright © 2018 by authors and Open Access Library Inc.
This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

<figure><img src="image_2.png" /></figure>

Open Access

## Abstract

Haar wavelets are applied for solution of three dimensional partial differential equations (PDEs) or time depending two dimensional PDEs. The proposed method is mathematically simple and fast. Two techniques are used in numerical solution, the first based on 2D-Haar wavelets and the second based on 3D-Haar wavelets and we compare them. To demonstrate the efficiency of the method, two test problems (solution of the diffusion and Poisson equations) are discussed. Computer simulation showed that 3D-Haar wavelets are better and closer to the exact solution but it is need to more time from 2D-Haar wavelets.

## Subject Areas

Numerical Mathematics, Partial Differential Equation

## Keywords

Haar Wavelets, Partial Differential Equations, Diffusion Equation, Poisson Equation

## 1. Introduction

As a powerful mathematical tool, Wavelet analysis has been widely used in image digital processing, quantum field theory, numerical analysis and many other fields in recent years.

Haar wavelets have been applied extensively for signal processing in communications and physics research, and more mathematically focused on differential equations and even nonlinear problems. After discrediting the differential equation in a conventional way like the finite difference approximation, wavelets can be used for algebraic manipulations in the system of equations obtained which may lead to better condition number of the resulting system [1].

<footnote n="1">¹Jean Leary Laboratory of Mathematics, University of Nantes, Nantes, France</footnote>

<footnote n="2">²College of Computer Sciences and Mathematics, University of Mosul, Mosul, Republic of Iraq
Email: Abdeljalil.Nachaoui@univ-nantes.fr, ekhlass70@yahoo.com, *ahmednumerical@yahoo.com</footnote>
Wavelet methods have been applied for solving partial differential equations (PDE-s) from beginning of the early 1990s [2] and [3]. In the last two decades this problem has attracted great attention and numerous papers about these topics have been published. Due to this fact, we must confine somewhat our analysis; in the following only PDEs of mathematical physics (elliptic, parabolic and hyperbolic equations) and of elastostatics are considered. From the first field of investigation, the papers [4] [5] [6] and [7] can be cited. As to the elasticity problems, we refer to the papers [8]-[13]. In all these papers, different wavelet families have been applied.

In most cases, the wavelet coefficients were calculated by the Galerkin or collocation method, by it we had to evaluate integrals of some combinations of the wavelet functions (called also connection coefficients).

Among all the wavelet families, the Haar wavelets deserve special attention. They are made up of pairs of piecewise constant functions and are therefore mathematically the simplest of all the wavelet families. A good feature of the Haar wavelets is also the possibility to integrate these wavelets analytically in arbitrary times. A drawback of these wavelets is their discontinuity; since the derivatives do not exist in the breaking points, it is not possible to apply these wavelets directly for solving PDEs. Ü. Lepik has applied this technique for solving different 1-D problems [14] [15] [16] and [17]. Also he used two dimensional Haar wavelets in solving PDFs which contain two variables [18] [19] [20]. The method is fast and with low error.

The aim of the present paper is to develop the Haar wavelet method for solving three dimensional PDEs, which is fast, mathematically simple and guarantees the necessary accuracy for a relatively small number of grid points. The method is an expansion of the 2D Haar wavelets method which discussed in [18]. We developed 3D and 2D Haar wavelets approximations to the solution of the partial differential equation. We obtain the 2D Haar wavelets method as an approximation of the 3D Haar wavelets.

The paper is organized as follows. In Section 2 formulas for calculating the Haar wavelets and their integrals are reported. In Sections 3, the method of solution is described by using 2D and 3D Haar wavelet respectively. In Sections 4 application of Haar wavelets method is presented in solving two problems (integration of the diffusion and Poisson equations). Conclusions and possible further directions of research are offered in Section 5.

## 2. Haar Wavelets and Their Integrals

The Haar functions are an orthogonal family of switched rectangular waveforms where amplitudes can differ from one function to another. They are defined in the interval [A,B] by [18]:

$$h_l(x) = \begin{cases} 1 & \text{for } x \in [\xi_1(i), \xi_2(i)], \\ -1 & \text{for } x \in [\xi_2(i), \xi_3(i)], \\ 0 & \text{elsewhere} \end{cases} \quad (1)$$
Wavelet methods have been applied for solving partial differential equations (PDE-s) from beginning of the early 1990s [2] and [3]. In the last two decades this problem has attracted great attention and numerous papers about these topics have been published. Due to this fact, we must confine somewhat our analysis; in the following only PDEs of mathematical physics (elliptic, parabolic and hyperbolic equations) and of elastostatics are considered. From the first field of investigation, the papers [4] [5] [6] and [7] can be cited. As to the elasticity problems, we refer to the papers [8]-[13]. In all these papers, different wavelet families have been applied.

In most cases, the wavelet coefficients were calculated by the Galerkin or collocation method, by it we had to evaluate integrals of some combinations of the wavelet functions (called also connection coefficients).

Among all the wavelet families, the Haar wavelets deserve special attention. They are made up of pairs of piecewise constant functions and are therefore mathematically the simplest of all the wavelet families. A good feature of the Haar wavelets is also the possibility to integrate these wavelets analytically in arbitrary times. A drawback of these wavelets is their discontinuity; since the derivatives do not exist in the breaking points, it is not possible to apply these wavelets directly for solving PDEs. Ü. Lepik has applied this technique for solving different 1-D problems [14] [15] [16] and [17]. Also he used two dimensional Haar wavelets in solving PDFs which contain two variables [18] [19] [20]. The method is fast and with low error.

The aim of the present paper is to develop the Haar wavelet method for solving three dimensional PDEs, which is fast, mathematically simple and guarantees the necessary accuracy for a relatively small number of grid points. The method is an expansion of the 2D Haar wavelets method which discussed in [18]. We developed 3D and 2D Haar wavelets approximations to the solution of the partial differential equation. We obtain the 2D Haar wavelets method as an approximation of the 3D Haar wavelets.

The paper is organized as follows. In Section 2 formulas for calculating the Haar wavelets and their integrals are reported. In Sections 3, the method of solution is described by using 2D and 3D Haar wavelet respectively. In Sections 4 application of Haar wavelets method is presented in solving two problems (integration of the diffusion and Poisson equations). Conclusions and possible further directions of research are offered in Section 5.

## 2. Haar Wavelets and Their Integrals

The Haar functions are an orthogonal family of switched rectangular waveforms where amplitudes can differ from one function to another. They are defined in the interval $[A, B]$ by [18]:

$$h_l(x) = \begin{cases} 1 & \text{for } x \in [\xi_1(i), \xi_2(i)], \\ -1 & \text{for } x \in [\xi_2(i), \xi_3(i)], \\ 0 & \text{elsewhere} \end{cases} \quad (1)$$
where

$$ 
\begin{aligned}
\xi_1(i) &= A + 2k\mu\Delta x, \quad \xi_2(i) = A + (2k+1)\mu\Delta x, \\
\xi_3(i) &= A + 2(k+1)\mu\Delta x, \quad \mu = M/m,
\end{aligned}
\qquad (2) $$

The interval $[A,B]$ is participated into 2M subintervals of equal length; the length of each subinterval is $\Delta x = (B-A)/(2M)$. Integer $m = 2^j$ ($j = 0,1,2,\dots,J$) indicates the level of the wavelet; $k = 0,1,2,\dots,m-1$ is the translation parameter. Maximal level of resolution is $J$. The index $i$ is calculated according the formula $i = m + k + 1$; in the case of minimal values. $M = 1, k = 0$ we have $I = 2$, the maximal value of $i$ is $i = 2M = 2^{j+1}$. It is assumed that the value $I = 1$ corresponds to the scaling function for which $h_i(x) = 1$ in $[A,B]$.

The operational matrix of integration P, which is a 2M square matrix, is defined by the equation:

$$ P_{1,i}(x) = \int_{0}^{x} h_{i}(x') dx', $$

In general

$$ P_{v+1,i}(x) = \int_{0}^{x} P_{v,j}(x') dx', \quad v=1,2,\dots \qquad (3) $$

The general form of v-times of integrals [18]:

$$ P_{v,j}(x) = \begin{cases} 0 & \text{for } x < \xi_1(i), \\ \frac{1}{v!} [x - \xi_1(i)]^v & \text{for } x \in [\xi_1(i), \xi_2(i)], \\ \frac{1}{v!} \{[x - \xi_1(i)]^v - 2[x - \xi_2(i)]^v\} & \text{for } x \in [\xi_2(i), \xi_3(i)], \\ \frac{1}{v!} \{[x - \xi_1(i)]^v - 2[x - \xi_2(i)]^v + [x - \xi_3(i)]^v\} & \text{for } x > \xi_3(i). \end{cases} \qquad (4) $$

For solving boundary value problems we need the values $P_{v,j}(B)$, which can be calculated from (4). In special cases $v=1$ or $v=2$, we find

$$ q_1(i) = P_{1,j}(B) = \begin{cases} B-A & \text{for } i=1, \\ 0 & \text{for } i \neq 1, \end{cases} \qquad (5) $$

and

$$ q_2(i) = P_{2,j}(B) = \begin{cases} 0.5(B-A)^2 & \text{for } i=1, \\ 0.25 \frac{(B-A)^2}{m^2} & \text{for } i \neq 1, \end{cases} \qquad (6) $$

In the present paper the collocation method for solving the PDEs is applied. Equations. (1) and (4) are discredited by replacing $x \to x_i$ such that:

$$ x_i = A + \left( l - \frac{1}{2} \right) \Delta x, \quad l = 1, 2, \dots, 2M \qquad (7) $$

It is convenient to introduce the Haar matrices
where

$$ 
\begin{aligned}
\xi_1(i) &= A + 2k\mu\Delta x, \quad \xi_2(i) = A + (2k+1)\mu\Delta x, \\
\xi_3(i) &= A + 2(k+1)\mu\Delta x, \quad \mu = M/m,
\end{aligned}
\qquad (2) $$

The interval $[A,B]$ is participated into 2M subintervals of equal length; the length of each subinterval is $\Delta x = (B-A)/(2M)$. Integer $m = 2^j$ ($j = 0,1,2,\dots,J$) indicates the level of the wavelet; $k = 0,1,2,\dots,m-1$ is the translation parameter. Maximal level of resolution is $J$. The index $i$ is calculated according the formula $i = m + k + 1$; in the case of minimal values. $M = 1, k = 0$ we have $I = 2$, the maximal value of $i$ is $i = 2M = 2^{j+1}$. It is assumed that the value $I = 1$ corresponds to the scaling function for which $h_i(x) = 1$ in $[A,B]$.

The operational matrix of integration P, which is a 2M square matrix, is defined by the equation:

$$ P_{1,i}(x) = \int_{0}^{x} h_{i}(x') dx', $$

In general

$$ P_{v+1,i}(x) = \int_{0}^{x} P_{v,j}(x') dx', \quad v=1,2,\dots \qquad (3) $$

The general form of v-times of integrals [18]:

$$ P_{v,j}(x) = \begin{cases} 0 & \text{for } x < \xi_1(i), \\ \frac{1}{v!} [x - \xi_1(i)]^v & \text{for } x \in [\xi_1(i), \xi_2(i)], \\ \frac{1}{v!} \{[x - \xi_1(i)]^v - 2[x - \xi_2(i)]^v\} & \text{for } x \in [\xi_2(i), \xi_3(i)], \\ \frac{1}{v!} \{[x - \xi_1(i)]^v - 2[x - \xi_2(i)]^v + [x - \xi_3(i)]^v\} & \text{for } x > \xi_3(i). \end{cases} \qquad (4) $$

For solving boundary value problems we need the values $P_{v,j}(B)$, which can be calculated from (4). In special cases $v=1$ or $v=2$, we find

$$ q_1(i) = P_{1,j}(B) = \begin{cases} B-A & \text{for } i=1, \\ 0 & \text{for } i \neq 1, \end{cases} \qquad (5) $$

and

$$ q_2(i) = P_{2,j}(B) = \begin{cases} 0.5(B-A)^2 & \text{for } i=1, \\ 0.25 \frac{(B-A)^2}{m^2} & \text{for } i \neq 1, \end{cases} \qquad (6) $$

In the present paper the collocation method for solving the PDEs is applied. Equations. (1) and (4) are discredited by replacing $x \to x_i$ such that:

$$ x_i = A + \left( l - \frac{1}{2} \right) \Delta x, \quad l = 1, 2, \dots, 2M \qquad (7) $$

It is convenient to introduce the Haar matrices
$H(i,l) = h_i(x_i)$, $P_v(i,l) = P_{v,i}(x_i)$. In the following Sections computer simulations were carried out with the aid of the Matlab programs for which the matrix representation is effective.

### 3. Problem Statement and Method of Solution

Consider two-dimensional partial differential equation of higher order:

$$ 
\begin{gathered}
F(t, x, y, u, Du, D^2 u, \dots, D^{\lambda+\beta+\alpha} u) = f(t, x, y), \\
D^{\lambda+\beta+\alpha} u = \frac{\partial^{(\lambda+\beta+\alpha)} u(t, x, y)}{\partial t^{\lambda} \partial x^{\beta} \partial y^{\alpha}},
\end{gathered}
\quad (8) $$

such that $f(t, x, y)$ is known function or constant.

The independent variables $t, x$ and $y$ belong to a domain $\Omega$ $x \in [A_1, B_1]$, $y \in [A_2, B_2]$, $t \in [A_3, B_3]$, which has the boundary $\partial\Omega$. We have to calculate the function $u(t, x, y)$, which satisfies the required initial and boundary conditions.

#### 3.1. The Solution by the 3D Haar Wavelets

The solution by the 3D Haar wavelets method is started by divides Cuboids $x \in [A_1, B_1]$, $y \in [A_2, B_2]$, $t \in [A_3, B_3]$ into $2M_1$, $2M_2$ and $2M_3$ parts of equal length, respectively.

We assume that the solution is sought in the form:

$$ \frac{\partial^{(\lambda+\beta+\alpha)} u(t, x, y)}{\partial t^{\lambda} \partial x^{\beta} \partial y^{\alpha}} = \sum_{j=1}^{2M_1} \sum_{l=1}^{2M_2} \sum_{i=1}^{2M_3} a_{j,l,j} h_i(t) P_{\beta,j}(x) h_l(y) + \sum_{jj=0}^{\beta-1} \frac{(x)^{jj}}{(jj)!} \frac{\partial^{(\lambda+jj+\alpha)} u(t, 0, y)}{\partial t^{\lambda} \partial x^{jj} \partial y^{\alpha}}, \quad (9) $$

where the elements $a_{j,l,j}$ are constants.

We integration (9) $\beta$-times in regard to $(x)$ from 0 to $x$, we obtain

$$ \frac{\partial^{(\lambda+\alpha)} u(t, x, y)}{\partial t^{\lambda} \partial y^{\alpha}} = \sum_{j=1}^{2M_1} \sum_{l=1}^{2M_2} \sum_{i=1}^{2M_3} a_{j,l,j} h_i(t) P_{\beta,j}(x) h_l(y) + \sum_{jj=0}^{\beta-1} \frac{(x)^{jj}}{(jj)!} \frac{\partial^{(\lambda+jj+\alpha)} u(t, 0, y)}{\partial t^{\lambda} \partial x^{jj} \partial y^{\alpha}}, \quad (10) $$

Now, by integrating (10) $\alpha$-times in regard to $(y)$ from 0 to $y$, we obtain

$$ \frac{\partial^{(\lambda)} u(t, x, y)}{\partial t^{\lambda}} = \sum_{j=1}^{2M_1} \sum_{l=1}^{2M_2} \sum_{i=1}^{2M_3} a_{j,l,j} h_i(t) P_{\beta,j}(x) P_{\alpha,l}(y) + \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(\lambda+ll)} u(t, x, 0)}{\partial t^{\lambda} \partial y^{ll}} \\ + \sum_{jj=0}^{\beta-1} \frac{(x)^{jj}}{(jj)!} \frac{\partial^{(\lambda+jj)} u(t, 0, y)}{\partial t^{\lambda} \partial x^{jj}} - \sum_{jj=0}^{\beta-1} \frac{(x)^{jj}}{(jj)!} \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(\lambda+jj+ll)} u(t, 0, 0)}{\partial t^{\lambda} \partial x^{jj} \partial y^{\alpha}}, \quad (11) $$

Now, by integrating (11) $\lambda$-times in regard to $(t)$ from 0 to $t$, we obtain

$$ (t, x, y) = \sum_{j=1}^{2M_1} \sum_{l=1}^{2M_2} \sum_{i=1}^{2M_3} a_{j,l,j} P_{\lambda,j}(t) P_{\beta,j}(x) P_{\alpha,l}(y) + \sum_{ll=0}^{\lambda-1} \frac{(t)^{ll}}{(ll)!} \frac{\partial^{(ll)} u(0, x, y)}{\partial t^{ll}} \\ + \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(ll)} u(t, x, 0)}{\partial y^{ll}} - \sum_{ll=0}^{\lambda-1} \frac{(t)^{ll}}{(ll)!} \sum_{jj=0}^{\alpha-1} \frac{y^{jj}}{(jj)!} \frac{\partial^{(ll+ll)} u(0, x, 0)}{\partial t^{ll} \partial y^{ll}} \\ + \sum_{jj=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \frac{\partial^{(jj)} u(t, 0, y)}{\partial x^{jj}} - \sum_{ll=0}^{\lambda-1} \frac{(t)^{ll}}{(ll)!} \sum_{jj=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \frac{\partial^{(ll+jj)} u(0, 0, y)}{\partial t^{ll} \partial x^{jj}} $$
$H(i,l) = h_i(x_i)$, $P_v(i,l) = P_{v,i}(x_i)$. In the following Sections computer simulations were carried out with the aid of the Matlab programs for which the matrix representation is effective.

### 3. Problem Statement and Method of Solution

Consider two-dimensional partial differential equation of higher order:

$$ 
\begin{gathered}
F(t, x, y, u, Du, D^2 u, \cdots, D^{\lambda+\beta+\alpha} u) = f(t, x, y), \\
D^{\lambda+\beta+\alpha} u = \frac{\partial^{(\lambda+\beta+\alpha)} u(t, x, y)}{\partial t^{\lambda} \partial x^{\beta} \partial y^{\alpha}},
\end{gathered}
\qquad (8) $$

such that $f(t, x, y)$ is known function or constant.

The independent variables $t$, $x$ and $y$ belong to a domain $\Omega$ $x \in [A_1, B_1]$, $y \in [A_2, B_2]$, $t \in [A_3, B_3]$, which has the boundary $\partial\Omega$. We have to calculate the function $u(t, x, y)$, which satisfies the required initial and boundary conditions.

#### 3.1. The Solution by the 3D Haar Wavelets

The solution by the 3D Haar wavelets method is started by divides Cuboids $x \in [A_1, B_1]$, $y \in [A_2, B_2]$, $t \in [A_3, B_3]$ into $2M_1$, $2M_2$ and $2M_3$ parts of equal length, respectively.

We assume that the solution is sought in the form:

$$ \frac{\partial^{(\lambda+\beta+\alpha)} u(t, x, y)}{\partial t^{\lambda} \partial x^{\beta} \partial y^{\alpha}} = \sum_{j=1}^{2M_1} \sum_{l=1}^{2M_2} \sum_{i=1}^{2M_3} a_{j,l,j} h_i(t) P_{\beta,j}(x) h_l(y) + \sum_{jj=0}^{\beta-1} \frac{(x)^{jj}}{(jj)!} \frac{\partial^{(\lambda+jj+\alpha)} u(t, 0, y)}{\partial t^{\lambda} \partial x^{jj} \partial y^{\alpha}}, \qquad (9) $$

where the elements $a_{j,l,j}$ are constants.

We integration (9) $\beta$-times in regard to $(x)$ from 0 to $x$, we obtain

$$ \frac{\partial^{(\lambda+\alpha)} u(t, x, y)}{\partial t^{\lambda} \partial y^{\alpha}} = \sum_{j=1}^{2M_1} \sum_{l=1}^{2M_2} \sum_{i=1}^{2M_3} a_{j,l,j} h_i(t) P_{\beta,j}(x) h_l(y) + \sum_{jj=0}^{\beta-1} \frac{(x)^{jj}}{(jj)!} \frac{\partial^{(\lambda+jj+\alpha)} u(t, 0, y)}{\partial t^{\lambda} \partial x^{jj} \partial y^{\alpha}}, \qquad (10) $$

Now, by integrating (10) $\alpha$-times in regard to $(y)$ from 0 to $y$, we obtain

$$ \frac{\partial^{(\lambda)} u(t, x, y)}{\partial t^{\lambda}} = \sum_{j=1}^{2M_1} \sum_{l=1}^{2M_2} \sum_{i=1}^{2M_3} a_{j,l,j} h_i(t) P_{\beta,j}(x) P_{\alpha,l}(y) + \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(\lambda+ll)} u(t, x, 0)}{\partial t^{\lambda} \partial y^{ll}} \\ + \sum_{jj=0}^{\beta-1} \frac{(x)^{jj}}{(jj)!} \frac{\partial^{(\lambda+jj)} u(t, 0, y)}{\partial t^{\lambda} \partial x^{jj}} - \sum_{jj=0}^{\beta-1} \frac{(x)^{jj}}{(jj)!} \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(\lambda+jj+ll)} u(t, 0, 0)}{\partial t^{\lambda} \partial x^{jj} \partial y^{\alpha}}, \qquad (11) $$

Now, by integrating (11) $\lambda$-times in regard to $(t)$ from 0 to $t$, we obtain

$$ (t, x, y) = \sum_{j=1}^{2M_1} \sum_{l=1}^{2M_2} \sum_{i=1}^{2M_3} a_{j,l,j} P_{\lambda,i}(t) P_{\beta,j}(x) P_{\alpha,l}(y) + \sum_{ll=0}^{\lambda-1} \frac{(t)^{ll}}{(ll)!} \frac{\partial^{(ii)} u(0, x, y)}{\partial t^{ii}} \\ + \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(ll)} u(t, x, 0)}{\partial y^{ll}} - \sum_{ll=0}^{\lambda-1} \frac{(t)^{ll}}{(ll)!} \sum_{jj=0}^{\alpha-1} \frac{y^{jj}}{(jj)!} \frac{\partial^{(ii+ll)} u(0, x, 0)}{\partial t^{ii} \partial y^{jj}} \\ + \sum_{jj=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \frac{\partial^{(jj)} u(t, 0, y)}{\partial x^{jj}} - \sum_{ll=0}^{\lambda-1} \frac{(t)^{ll}}{(ll)!} \sum_{jj=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \frac{\partial^{(ii+jj)} u(0, 0, y)}{\partial t^{ii} \partial x^{jj}} $$
$$ -\sum_{jj=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(jj+ll)} u(t, 0, 0)}{\partial x^{jj} \partial y^{ll}} + \sum_{jj=0}^{\lambda-1} \frac{t^{jj}}{(jj)!} \sum_{ll=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(jj+ll)} u(0, 0, 0)}{\partial t^{jj} \partial x^{jj} \partial y^{ll}} \quad (12) $$

In this formula, the integrals $P_{\lambda,j}(t), P_{\beta,j}(x)$ and $P_{\alpha,j}(y)$ are calculated according to (4) and the other terms in the Equation (12) are calculated according to the type of the initial and boundary conditions (Dirichlet, Neumann, and Mixed boundary conditions). Details of this method are explained by solving two examples.

## 3.2. The Solution by the 2D Haar Wavelets

When we use the 2D Haar wavelets method, we divide the interval $t \in [A_3, B_3]$ into N equal parts of length $\Delta t = (B_3 - A_3)/N$ and let $t_s = (s-1)\Delta t$ $s = 1, 2, \dots, N$ and:

$$ a_{j,j}(t) = \sum_{i=1}^{2M_3} a_{j,j,i} h_i(t), \quad (13) $$

The mean idea of 2D Haar wavelets is to assume that $a_{j,j}(t)$ are constants in each subinterval $t \in (t_s, t_{s+1}]$, then

For all $t \in (t_s, t_{s+1}]$ the Equation (11) becomes

$$ \frac{\partial^{(\lambda)} u(t, x, y)}{\partial t^{\lambda}} = \sum_{j=1}^{2M_1} \sum_{l=1}^{2M_2} a_{j,j} P_{\beta,j}(x) P_{\alpha,j}(y) + \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(\lambda+ll)} u(t, x, 0)}{\partial t^{\lambda} \partial y^{ll}} + \sum_{jj=0}^{\beta-1} \frac{(x)^{jj}}{(jj)!} \frac{\partial^{(\lambda+jj)} u(t, 0, y)}{\partial t^{\lambda} \partial x^{jj}} - \sum_{jj=0}^{\beta-1} \frac{(x)^{jj}}{(jj)!} \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(\lambda+jj+ll)} u(t, 0, 0)}{\partial t^{\lambda} \partial x^{jj} \partial y^{ll}}, \quad (14) $$

We integration (14) $\lambda$-times in regard to $(t)$ from $t_s$ to $t$, we obtain

$$ \begin{aligned} u(t, x, y) = & \frac{(t-t_s)^{\lambda}}{(\lambda)!} \sum_{j=1}^{2M_1} \sum_{l=1}^{2M_2} a_{j,j} P_{\beta,j}(x) P_{\alpha,j}(y) + \sum_{ll=0}^{\lambda-1} \frac{(t-t_s)^{ll}}{(ll)!} \frac{\partial^{(ll)} u(t_s, x, y)}{\partial t^{ll}} \\ & + \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(ll)} u(t, x, 0)}{\partial y^{ll}} - \sum_{ll=0}^{\lambda-1} \frac{(t-t_s)^{ll}}{(ll)!} \frac{\partial^{(ll+ll)} u(t_s, x, 0)}{\partial t^{ll} \partial y^{ll}} \\ & + \sum_{jj=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \frac{\partial^{(jj)} u(t, 0, y)}{\partial x^{jj}} - \sum_{ll=0}^{\lambda-1} \frac{(t-t_s)^{ll}}{(ll)!} \sum_{jj=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \frac{\partial^{(jj+ll)} u(t_s, 0, y)}{\partial t^{ll} \partial x^{jj}} \\ & - \sum_{jj=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(jj+ll)} u(t, 0, 0)}{\partial x^{jj} \partial y^{ll}} \\ & + \sum_{ll=0}^{\lambda-1} \frac{(t-t_s)^{ll}}{(ll)!} \sum_{jj=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(ll+jj+ll)} u(t_s, 0, 0)}{\partial t^{ll} \partial x^{jj} \partial y^{ll}}. \end{aligned} \quad (15) $$

As in 3D Haar wavelets, integrals $P_{\beta,j}(x)$ and $P_{\alpha,j}(y)$ are calculated according to (4) and the other terms in the equation (15) are calculated according to the type of the initial and boundary conditions (Dirichlet, Neumann, and Mixed boundary conditions).
$$ -\sum_{jj=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(jj+ll)} u(t, 0, 0)}{\partial x^{jj} \partial y^{ll}} + \sum_{jj=0}^{\lambda-1} \frac{t^{jj}}{(jj)!} \sum_{ll=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(jj+ll)} u(0, 0, 0)}{\partial t^{jj} \partial x^{jj} \partial y^{ll}}, \quad (12) $$

In this formula, the integrals $P_{\lambda,j}(t)$, $P_{\beta,j}(x)$ and $P_{\alpha,j}(y)$ are calculated according to (4) and the other terms in the Equation (12) are calculated according to the type of the initial and boundary conditions (Dirichlet, Neumann, and Mixed boundary conditions). Details of this method are explained by solving two examples.

## 3.2. The Solution by the 2D Haar Wavelets

When we use the 2D Haar wavelets method, we divide the interval $t \in [A_3, B_3]$ into $N$ equal parts of length $\Delta t = (B_3 - A_3)/N$ and let $t_s = (s-1)\Delta t$ $s = 1, 2, \dots, N$ and:

$$ a_{j,j}(t) = \sum_{i=1}^{2M_3} a_{j,j,i} h_i(t), \quad (13) $$

The mean idea of 2D Haar wavelets is to assume that $a_{j,j}(t)$ are constants in each subinterval $t \in (t_s, t_{s+1}]$, then

For all $t \in (t_s, t_{s+1}]$ the Equation (11) becomes

$$ \frac{\partial^{(\lambda)} u(t, x, y)}{\partial t^{\lambda}} = \sum_{j=1}^{2M_1} \sum_{l=1}^{2M_2} a_{j,j} P_{\beta,j}(x) P_{\alpha,j}(y) + \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(\lambda+ll)} u(t, x, 0)}{\partial t^{\lambda} \partial y^{ll}} + \sum_{jj=0}^{\beta-1} \frac{(x)^{jj}}{(jj)!} \frac{\partial^{(\lambda+jj)} u(t, 0, y)}{\partial t^{\lambda} \partial x^{jj}} - \sum_{jj=0}^{\beta-1} \frac{(x)^{jj}}{(jj)!} \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(\lambda+jj+ll)} u(t, 0, 0)}{\partial t^{\lambda} \partial x^{jj} \partial y^{ll}}, \quad (14) $$

We integration (14) $\lambda$-times in regard to $(t)$ from $t_s$ to $t$, we obtain

$$ u(t, x, y) = \frac{(t - t_s)^{\lambda}}{(\lambda)!} \sum_{j=1}^{2M_1} \sum_{l=1}^{2M_2} a_{j,j} P_{\beta,j}(x) P_{\alpha,j}(y) + \sum_{ll=0}^{\lambda-1} \frac{(t - t_s)^{ll}}{(ll)!} \frac{\partial^{(ll)} u(t_s, x, y)}{\partial t^{ll}} + \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(ll)} u(t, x, 0)}{\partial y^{ll}} - \sum_{ll=0}^{\lambda-1} \frac{(t - t_s)^{ll}}{(ll)!} \frac{\partial^{(ll+ll)} u(t_s, x, 0)}{\partial t^{ll} \partial y^{ll}} + \sum_{jj=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \frac{\partial^{(jj)} u(t, 0, y)}{\partial x^{jj}} - \sum_{jj=0}^{\beta-1} \frac{(t - t_s)^{jj}}{(jj)!} \frac{\partial^{(jj)} u(t_s, 0, y)}{\partial t^{jj} \partial x^{jj}} - \sum_{jj=0}^{\beta-1} \frac{x^{jj}}{(jj)!} \sum_{ll=0}^{\alpha-1} \frac{y^{ll}}{(ll)!} \frac{\partial^{(jj+ll)} u(t, 0, 0)}{\partial t^{ll} \partial x^{jj} \partial y^{ll}} + \sum_{jj=0}^{\beta-1} \frac{(t - t_s)^{jj}}{(jj)!} \sum_{ll=0}^{\alpha-1} \frac{y^{jj}}{(jj)!} \frac{\partial^{(jj+ll)} u(t_s, 0, 0)}{\partial t^{jj} \partial x^{jj} \partial y^{ll}}. \quad (15) $$

As in 3D Haar wavelets, integrals $P_{\beta,j}(x)$ and $P_{\alpha,j}(y)$ are calculated according to (4) and the other terms in the equation (15) are calculated according to the type of the initial and boundary conditions (Dirichlet, Neumann, and Mixed boundary conditions).
## 4. Application and Numerical Results

We apply the 2D Haar wavelet and 3D Haar wavelet methods in solve two problems (diffusion and Poisson equations) and comparison with the exact solution.

### 4.1. Diffusion Equation

Solve the 2D Heat equation in the domain $(0,1) \times \Omega$ where $\Omega = [0,1] \times [0,1]$ and

$$ 
\left. 
\begin{aligned}
\frac{\partial u}{\partial t} - c^2 \Delta u &= f(t, x, y) && \text{on } (0,1) \times \Omega, \\
u(t, 0, y) &= u(t, 1, y) = 0 && t \in (0,1) \text{ and } y \in [0,1], \\
u(t, x, 0) &= u(t, x, 1) = 0 && t \in (0,1) \text{ and } x \in [0,1], \\
u(0, x, y) &= 0 && \text{on } \Omega.
\end{aligned}
\right\} 
\quad (16) 
$$

Here we have $\lambda = 1, \beta = 2$ and $\alpha = 2$ and suppose that $M_1 = M_2 = M_3 = M$.

#### 4.1.1. The Solution by 3D Haar Wavelet Method

The solution by 3D Haar wavelet is begin using the Equation (12) to approximate problem (16) and considering the initial and boundary conditions at $x = 0$ and $y = 0$, we get

$$ 
\begin{aligned}
u(t,x,y) = u(0,x,y) &+ \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} P_{1,i}(t) P_{2,j}(x) P_{2,l}(y) \\
&+ y \left( \left. \frac{\partial u(t,x,y)}{\partial y} \right|_{y=0} \right) + x \left( \left. \frac{\partial u(t,x,y)}{\partial x} \right|_{x=0} \right) - xy \left( \left. \frac{\partial^2 u(t,x,y)}{\partial x \partial y} \right|_{x=y=0} \right),
\end{aligned} 
\quad (17) 
$$

Taking $x = 1$ and using the boundary conditions in the last equation, we obtain

$$ \left( \left. \frac{\partial u(t, x, y)}{\partial x} \right|_{x=0} \right) - y \left( \left. \frac{\partial^2 u(t, x, y)}{\partial x \partial y} \right|_{x=y=0} \right) = - \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} P_{1,i}(t) q_2(j) P_{2,l}(y), $$

Replacing this result back into (17), we obtain

$$ 
\begin{aligned}
u(t,x,y) = u(0,x,y) &+ \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} P_{1,i}(t) P_{2,j}(x) P_{2,l}(y) \\
&+ y \left( \left. \frac{\partial u(t,x,y)}{\partial y} \right|_{y=0} \right) - x \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} P_{1,i}(t) q_2(j) P_{2,l}(y),
\end{aligned} 
\quad (18) 
$$

Similarly and by using the boundary condition at $y = 1$, we obtain

$$ 
\left( \left. \frac{\partial u(t, x, y)}{\partial y} \right|_{y=0} \right) = - \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} P_{1,i}(t) P_{2,j}(x) q_2(l) \\
+ x \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} P_{1,i}(t) q_2(j) q_2(l),
$$
## 4. Application and Numerical Results

We apply the 2D Haar wavelet and 3D Haar wavelet methods in solve two problems (diffusion and Poisson equations) and comparison with the exact solution.

### 4.1. Diffusion Equation

Solve the 2D Heat equation in the domain $(0,1) \times \Omega$ where $\Omega = [0,1] \times [0,1]$ and

$$ 
\left. 
\begin{aligned}
\frac{\partial u}{\partial t} - c^2 \Delta u &= f(t, x, y) && \text{on } (0,1) \times \Omega, \\
u(t, 0, y) &= u(t, 1, y) = 0 && t \in (0,1) \text{ and } y \in [0,1], \\
u(t, x, 0) &= u(t, x, 1) = 0 && t \in (0,1) \text{ and } x \in [0,1], \\
u(0, x, y) &= 0 && \text{on } \Omega.
\end{aligned}
\right\} \quad (16)
$$

Here we have $\lambda = 1, \beta = 2$ and $\alpha = 2$ and suppose that $M_1 = M_2 = M_3 = M$.

#### 4.1.1. The Solution by 3D Haar Wavelet Method

The solution by 3D Haar wavelet is begin using the Equation (12) to approximate problem (16) and considering the initial and boundary conditions at $x = 0$ and $y = 0$, we get

$$ u(t,x,y) = u(0,x,y) + \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} P_{1,i}(t) P_{2,j}(x) P_{2,l}(y) \\ + y \left( \left. \frac{\partial u(t,x,y)}{\partial y} \right|_{y=0} \right) + x \left( \left. \frac{\partial u(t,x,y)}{\partial x} \right|_{x=0} \right) - xy \left( \left. \frac{\partial^2 u(t,x,y)}{\partial x \partial y} \right|_{x=y=0} \right), \quad (17) $$

Taking $x = 1$ and using the boundary conditions in the last equation, we obtain

$$ \left( \left. \frac{\partial u(t,x,y)}{\partial x} \right|_{x=0} \right) - y \left( \left. \frac{\partial^2 u(t,x,y)}{\partial x \partial y} \right|_{y=0} \right) = - \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} P_{1,i}(t) q_2(j) P_{2,l}(y), $$

Replacing this result back into (17), we obtain

$$ u(t,x,y) = u(0,x,y) + \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} P_{1,i}(t) P_{2,j}(x) P_{2,l}(y) \\ + y \left( \left. \frac{\partial u(t,x,y)}{\partial y} \right|_{y=0} \right) - x \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} P_{1,i}(t) q_2(j) P_{2,l}(y), \quad (18) $$

Similarly and by using the boundary condition at $y = 1$, we obtain

$$ \left( \left. \frac{\partial u(t,x,y)}{\partial y} \right|_{y=0} \right) = - \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} P_{1,i}(t) P_{2,j}(x) q_2(l) \\ + x \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} P_{1,i}(t) q_2(j) q_2(l), $$
Replacing this result back into (18), we obtain

$$ 
\begin{aligned}
u(t,x,y) = & u(0,x,y) + \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} P_{1,i}(t) P_{2,j}(x) P_{2,l}(y) \\
& - x \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} P_{1,i}(t) q_2(j) P_{2,l}(y) \\
& - y \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} P_{1,i}(t) P_{2,j}(x) q_2(l) \\
& + xy \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} P_{1,i}(t) q_2(j) q_2(l),
\end{aligned} 
\quad (19) $$

which can be rewritten as:

$$ u(t,x,y) = u(0,x,y) + \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} \{ P_{1,i}(t) [P_{2,j}(x) - xq_2(j)] [P_{2,l}(y) - yq_2(l)] \}, \quad (20) $$

Derivative the Equation (20), we obtain that

$$ \frac{\partial u}{\partial t} = \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} \{ h_i(t) [P_{2,j}(x) - xq_2(j)] [P_{2,l}(y) - yq_2(l)] \}, \quad (21) $$

$$ \frac{\partial^2 u}{\partial x^2} = \left. \frac{\partial^2 u}{\partial x^2} \right|_{t=0} + \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} \{ P_{1,i}(t) h_j(x) [P_{2,l}(y) - yq_2(l)] \}, \quad (22) $$

$$ \frac{\partial^2 u}{\partial y^2} = \left. \frac{\partial^2 u}{\partial y^2} \right|_{t=0} + \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} \{ P_{1,i}(t) [P_{2,j}(x) - xq_2(j)] h_i(y) \}, \quad (23) $$

Substituting Equations (21)-(23) in (16) for any collocation points $t_s, x_r$ and $y_k$ with $s \in \{1, 2, \dots, 2M\}, r \in \{1, 2, \dots, 2M\}, k \in \{1, 2, \dots, 2M\}$, we get

$$ \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} R_{j,j,l,r,k,s} = g(s,r,k) \quad \text{for } 1 \le s \le 2M, 1 \le r \le 2M, 1 \le k \le 2M \quad (24) $$

where

$$ R_{j,j,l,r,k,s} = H(i,s) [P_2(j,r) - x_r q_2(j)] [P_2(l,k) - y_k q_2(l)] \\ - c^2 P_1(i,s) H(j,r) [P_2(l,k) - y_k q_2(l)] \\ - c^2 P_1(i,s) [P_2(j,r) - x_r q_2(j)] H(l,k), \quad (25) $$

and

$$ g(s,r,k) = f(t_s, x_r, y_k) + c^2 \left. \frac{\partial^2 u}{\partial x^2} \right|_{t=0, x=x_r, y=y_k} + c^2 \left. \frac{\partial^2 u}{\partial y^2} \right|_{t=0, x=x_r, y=y_k}, $$

The terms $\left. \frac{\partial^2 u}{\partial x^2} \right|_{t=0, x=x_r, y=y_k}$ and $\left. \frac{\partial^2 u}{\partial y^2} \right|_{t=0, x=x_r, y=y_k}$ are given by initial condition.

With the following notations:

for $\mathfrak{I}_i \in (0,1)$

such that $\Delta\mathfrak{I} = \frac{1}{2M}$, we write:
Replacing this result back into (18), we obtain

$$
\begin{equation}
\begin{aligned}
u(t,x,y) ={}& u(0,x,y) + \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} P_{1,i}(t) P_{2,j}(x) P_{2,l}(y) \\
& - x \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} P_{1,i}(t) q_2(j) P_{2,l}(y) \\
& - y \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} P_{1,i}(t) P_{2,j}(x) q_2(l) \\
& + xy \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} P_{1,i}(t) q_2(j) q_2(l),
\end{aligned}
\tag{19}
\end{equation}
$$

which can be rewritten as:

$$
\begin{equation}
u(t,x,y) = u(0,x,y) + \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} \{ P_{1,i}(t) [P_{2,j}(x) - xq_2(j)] [P_{2,l}(y) - yq_2(l)] \}, \tag{20}
\end{equation}
$$

Derivative the Equation (20), we obtain that

$$
\frac{\partial u}{\partial t} = \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} \{ h_i(t) [P_{2,j}(x) - xq_2(j)] [P_{2,l}(y) - yq_2(l)] \}, \quad (21)
$$

$$
\frac{\partial^2 u}{\partial x^2} = \left. \frac{\partial^2 u}{\partial x^2} \right|_{t=0} + \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} \{ P_{1,i}(t) h_j(x) [P_{2,l}(y) - yq_2(l)] \}, \quad (22)
$$

$$
\frac{\partial^2 u}{\partial y^2} = \left. \frac{\partial^2 u}{\partial y^2} \right|_{t=0} + \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} \{ P_{1,i}(t) [P_{2,j}(x) - xq_2(j)] h_i(y) \}, \quad (23)
$$

Substituting Equations (21)-(23) in (16) for any collocation points $t_s, x_r$ and $y_k$ with $s \in \{1, 2, \dots, 2M\}$, $r \in \{1, 2, \dots, 2M\}$, $k \in \{1, 2, \dots, 2M\}$, we get

$$
\sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,j,l} R_{j,j,l,r,k,s} = g(s,r,k) \quad \text{for } 1 \le s \le 2M, 1 \le r \le 2M, 1 \le k \le 2M \quad (24)
$$

where

$$
\left.
\begin{aligned}
R_{j,j,l,r,k,s} = H(i,s) [P_2(j,r) - x_r q_2(j)] [P_2(l,k) - y_k q_2(l)] \\
& - c^2 P_1(i,s) H(j,r) [P_2(l,k) - y_k q_2(l)] \\
& - c^2 P_1(i,s) [P_2(j,r) - x_r q_2(j)] H(l,k),
\end{aligned}
\right\}
\tag{25}
$$

and

$$
g(s,r,k) = f(t_s, x_r, y_k) + c^2 \left. \frac{\partial^2 u}{\partial x^2} \right|_{t=0, x=x_r, y=y_k} + c^2 \left. \frac{\partial^2 u}{\partial y^2} \right|_{t=0, x=x_r, y=y_k},
$$

The terms $\left. \frac{\partial^2 u}{\partial x^2} \right|_{t=0, x=x_r, y=y_k}$ and $\left. \frac{\partial^2 u}{\partial y^2} \right|_{t=0, x=x_r, y=y_k}$ are given by initial condition.

With the following notations:

for $\Im \iota \in (0,1)$

such that $\Delta \Im = \frac{1}{2M}$, we write:
$$
\begin{align}
H(j,r) &= h_j(\mathfrak{T}_r) & j &= 1, 2, \dots, 2M; r = 1, 2, \dots, 2M, \nonumber \\
P_2(j,r) &= P_{2,j}(\mathfrak{T}_r) & j &= 1, 2, \dots, 2M; r = 1, 2, \dots, 2M, \tag{26} \\
q_2(j) &= P_{2,j}(1) & j &= 1, 2, \dots, 2M, \nonumber
\end{align}
$$

which calculated from Equations (1), (4) and (6) respectively. $\mathfrak{T}$ can be $t, x$ and $y$.

It is clear that, the wavelet coefficients $a_{j,l,i}$ can be obtained by solving the
linear system (24). For simplify, we transform the system into a form with
second-order matrices using the following:

Let

$$
\eta = (2M)^2 (j-1) + 2M(l-1) + i, \quad \mu = (2M)^2 (r-1) + 2M(k-1) + s, \quad \text{Now}
$$

Equation (24) can be rewritten in the following form

$$
\sum_{\eta=1}^{(2M)^3} S(\mu, \eta) B(\eta) = F(\mu) \text{ for } 1 \le \mu \le (2M)^3, \text{ which give the following sys-}
$$

tem of linear equations

$$
S \cdot B = F \tag{27}
$$

where B and F are (2M)³ vectors and S is a (2M)³ × (2M)³ matrix such that

$$
B = \begin{bmatrix} a_{1,1,1} & a_{1,1,2} & \cdots & a_{1,1,2M} & \cdots & a_{2,1,1} & \cdots & a_{2M,2M,1} & \cdots & a_{2M,2M,2M} \end{bmatrix}_{1 \times (2M)^3}
$$

$$
F = \begin{bmatrix} g_{1,1,1} & g_{1,1,2} & \cdots & g_{1,1,2M} & \cdots & g_{2,1,1} & \cdots & g_{2M,2M,1} & \cdots & g_{2M,2M,2M} \end{bmatrix}_{1 \times (2M)^3}
$$

$$
S = \begin{bmatrix}
R_{1,1,1,1,1,1} & R_{1,1,2,1,1,1} & \cdots & R_{1,1,2M,1,1,1} & R_{1,2,1,1,1,1} & \cdots & R_{2M,2M,2M,1,1,1} \\
R_{1,1,1,1,1,2} & R_{1,1,2,1,1,2} & \cdots & R_{1,1,2M,1,1,2} & R_{1,2,1,1,1,2} & \cdots & R_{2M,2M,2M,1,1,2} \\
\vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
R_{1,1,1,2M,2M,2M} & R_{1,1,2,2M,2M,2M} & \cdots & R_{1,1,2M,2M,2M,2M} & R_{1,2,1,2M,2M,2M} & \cdots & R_{2M,2M,2M,2M,2M,2M}
\end{bmatrix}_{(2M)^3 \times (2M)^3}
$$

where $R_{i,l,i,r,k,s}$ are calculated according to Equation (25). After solving system (27) we obtain the wavelet coefficients $a_{j,l,i}$ and thus for any

$(t, x, y) \in (0,1) \times \Omega$ the solution $u(t, x, y)$ is obtained from Equation (20).

4.1.2. The Solution by 2D Haar Wavelet Method

Now we use Equation (15) to approximate problem (16) and considering the ini-
tial and boundary conditions at x = 0 and y = 0, Equation (15) gives for
$(t,x,y) \in (t_s, t_{s+1}] \times \Omega:$

$$
\begin{align}
u(t,x,y) ={}& (t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} P_{2,j}(x) P_{2,l}(y) \nonumber \\
& + y \left( \left. \frac{\partial u(t,x,y)}{\partial y} \right|_{y=0} - \left. \frac{\partial u(t_s,x,y)}{\partial y} \right|_{y=0} \right) \nonumber \\
& + x \left( \left. \frac{\partial u(t,x,y)}{\partial x} \right|_{x=0} - \left. \frac{\partial u(t_s,x,y)}{\partial x} \right|_{x=0} \right) \nonumber \\
& - xy \left( \left. \frac{\partial^2 u(t,x,y)}{\partial x \partial y} \right|_{x=y=0} - \left. \frac{\partial^2 u(t_s,x,y)}{\partial x \partial y} \right|_{x=y=0} \right), \tag{28}
\end{align}
$$
$$
\begin{align}
H(j,r) &= h_j(\mathfrak{T}_r) & j &= 1, 2, \dots, 2M; r = 1, 2, \dots, 2M, \nonumber \\
P_2(j,r) &= P_{2,j}(\mathfrak{T}_r) & j &= 1, 2, \dots, 2M; r = 1, 2, \dots, 2M, \tag{26} \\
q_2(j) &= P_{2,j}(1) & j &= 1, 2, \dots, 2M, \nonumber
\end{align}
$$

which calculated from Equations (1), (4) and (6) respectively. $\mathfrak{T}$ can be $t$, $x$ and $y$.

It is clear that, the wavelet coefficients $a_{j,l,i}$ can be obtained by solving the
linear system (24). For simplify, we transform the system into a form with
second-order matrices using the following:

Let

$$
\eta = (2M)^2 (j-1) + 2M(l-1) + i, \quad \mu = (2M)^2 (r-1) + 2M(k-1) + s, \quad \text{Now}
$$

Equation (24) can be rewritten in the following form

$$
\sum_{\eta=1}^{(2M)^3} S(\mu, \eta) B(\eta) = F(\mu) \text{ for } 1 \le \mu \le (2M)^3, \text{ which give the following sys-}
$$

tem of linear equations

$$
S \cdot B = F \tag{27}
$$

where B and F are (2M)³ vectors and S is a (2M)³ × (2M)³ matrix such that

$$
B = \begin{bmatrix} a_{1,1,1} & a_{1,1,2} & \cdots & a_{1,1,2M} & \cdots & a_{2,1,1} & \cdots & a_{2M,2M,1} & \cdots & a_{2M,2M,2M} \end{bmatrix}_{1 \times (2M)^3}
$$

$$
F = \begin{bmatrix} g_{1,1,1} & g_{1,1,2} & \cdots & g_{1,1,2M} & \cdots & g_{2,1,1} & \cdots & g_{2M,2M,1} & \cdots & g_{2M,2M,2M} \end{bmatrix}_{1 \times (2M)^3}
$$

$$
S = \begin{bmatrix}
R_{1,1,1,1,1,1} & R_{1,1,2,1,1,1} & \cdots & R_{1,1,2M,1,1,1} & R_{1,2,1,1,1,1} & \cdots & R_{2M,2M,2M,1,1,1} \\
R_{1,1,1,1,1,2} & R_{1,1,2,1,1,2} & \cdots & R_{1,1,2M,1,1,2} & R_{1,2,1,1,1,2} & \cdots & R_{2M,2M,2M,1,1,2} \\
\vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
R_{1,1,1,2M,2M,2M} & R_{1,1,2,2M,2M,2M} & \cdots & R_{1,1,2M,2M,2M,2M} & R_{1,2,1,2M,2M,2M} & \cdots & R_{2M,2M,2M,2M,2M,2M,2M}
\end{bmatrix}_{(2M)^3 \times (2M)^3}
$$

where $R_{i,l,i,r,k,s}$ are calculated according to Equation (25). After solving system (27) we obtain the wavelet coefficients $a_{j,l,i}$ and thus for any

$(t, x, y) \in (0,1) \times \Omega$ the solution $u(t, x, y)$ is obtained from Equation (20).

4.1.2. The Solution by 2D Haar Wavelet Method

Now we use Equation (15) to approximate problem (16) and considering the ini-
tial and boundary conditions at $x = 0$ and $y = 0$, Equation (15) gives for
$(t, x, y) \in (t_s, t_{s+1}] \times \Omega:$

$$
\begin{align}
u(t,x,y) ={}& (t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} P_{2,j}(x) P_{2,l}(y) + u(t_s,x,y) \nonumber \\
& + y \left( \left. \frac{\partial u(t,x,y)}{\partial y} \right|_{y=0} - \left. \frac{\partial u(t_s,x,y)}{\partial y} \right|_{y=0} \right) \nonumber \\
& + x \left( \left. \frac{\partial u(t,x,y)}{\partial x} \right|_{x=0} - \left. \frac{\partial u(t_s,x,y)}{\partial x} \right|_{x=0} \right) \nonumber \\
& - xy \left( \left. \frac{\partial^2 u(t,x,y)}{\partial x \partial y} \right|_{x=y=0} - \left. \frac{\partial^2 u(t_s,x,y)}{\partial x \partial y} \right|_{x=y=0} \right), \tag{28}
\end{align}
$$
Taking $x = 1$ and using the boundary conditions in the last equation, we obtain that for $t \in (t_s, t_{s+1}]$ and $y \in [0,1]$

$$ \left( \left. \frac{\partial u(t, x, y)}{\partial x} \right|_{x=0} - \left. \frac{\partial u(t_s, x, y)}{\partial x} \right|_{x=0} \right) - y \left( \left. \frac{\partial^2 u(t, x, y)}{\partial x \partial y} \right|_{x=y=0} - \left. \frac{\partial u^2(t_s, x, y)}{\partial x \partial y} \right|_{x=y=0} \right) = -(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} q_2(j) P_{2,l}(y), $$

Replacing this result back into (28), we obtain

$$ \begin{aligned} u(t,x,y) = & (t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} P_{2,j}(x) P_{2,l}(y) + u(t_s,x,y) \\ & + y \left( \left. \frac{\partial u(t,x,y)}{\partial y} \right|_{y=0} - \left. \frac{\partial u(t_s,x,y)}{\partial y} \right|_{y=0} \right) \\ & - x(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} q_2(j) P_{2,l}(y), \end{aligned} \quad (29) $$

Similarly and by using the boundary condition at $y = 1$, we obtain for $t \in (t_s, t_{s+1}]$ and $x \in [0,1]$

$$ \left( \left. \frac{\partial u(t, x, y)}{\partial y} \right|_{y=0} - \left. \frac{\partial u(t_s, x, y)}{\partial y} \right|_{y=0} \right) = -(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} P_{2,j}(x) q_2(l) + x(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} q_2(j) q_2(l), $$

Replacing

this result back into (29), we obtain

$$ \begin{aligned} u(t,x,y) = & u(t_s,x,y) + (t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} P_{2,j}(x) P_{2,l}(y) \\ & - x(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} q_2(j) P_{2,l}(y) \\ & - y(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} P_{2,j}(x) q_2(l) \\ & + xy(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} q_2(j) q_2(l), \end{aligned} \quad (30) $$

which can be rewritten as

$$ u(t,x,y) = u(t_s,x,y) + (t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \{ [P_{2,j}(x) - xq_2(j)] [P_{2,l}(y) - yq_2(l)] \}, \quad (31) $$

Derivative the Equation (31), we obtain that

$$ \frac{\partial u}{\partial t} = \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \{ [P_{2,j}(x) - xq_2(j)] [P_{2,l}(y) - yq_2(l)] \}, \quad (32) $$

$$ \left. \frac{\partial^2 u}{\partial x^2} = \frac{\partial^2 u}{\partial x^2} \right|_{t=t_s} + (t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \{ h_j(x) [P_{2,l}(y) - yq_2(l)] \}, \quad (33) $$

$$ \left. \frac{\partial^2 u}{\partial y^2} = \frac{\partial^2 u}{\partial y^2} \right|_{t=t_s} + (t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \{ [P_{2,j}(x) - xq_2(j)] h_l(y) \}, \quad (34) $$
Taking $x = 1$ and using the boundary conditions in the last equation, we obtain that for $t \in (t_s, t_{s+1}]$ and $y \in [0,1]$

$$ \left( \left. \frac{\partial u(t, x, y)}{\partial x} \right|_{x=0} - \left. \frac{\partial u(t_s, x, y)}{\partial x} \right|_{x=0} \right) - y \left( \left. \frac{\partial^2 u(t, x, y)}{\partial x \partial y} \right|_{x=y=0} - \left. \frac{\partial u^2(t_s, x, y)}{\partial x \partial y} \right|_{x=y=0} \right) = -(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} q_2(j) P_{2,l}(y), $$

Replacing this result back into (28), we obtain

$$ \begin{aligned} u(t, x, y) = & (t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} P_{2,j}(x) P_{2,l}(y) + u(t_s, x, y) \\ & + y \left( \left. \frac{\partial u(t, x, y)}{\partial y} \right|_{y=0} - \left. \frac{\partial u(t_s, x, y)}{\partial y} \right|_{y=0} \right) \\ & - x(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} q_2(j) P_{2,l}(y), \end{aligned} \quad (29) $$

Similarly and by using the boundary condition at $y = 1$, we obtain for $t \in (t_s, t_{s+1}]$ and $x \in [0,1]$

$$ \left( \left. \frac{\partial u(t, x, y)}{\partial y} \right|_{y=0} - \left. \frac{\partial u(t_s, x, y)}{\partial y} \right|_{y=0} \right) = -(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} P_{2,j}(x) q_2(l) + x(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} q_2(j) q_2(l), $$

Replacing

this result back into (29), we obtain

$$ \begin{aligned} u(t, x, y) = & u(t_s, x, y) + (t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} P_{2,j}(x) P_{2,l}(y) \\ & - x(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} q_2(j) P_{2,l}(y) \\ & - y(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} P_{2,j}(x) q_2(l) \\ & + xy(t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} q_2(j) q_2(l), \end{aligned} \quad (30) $$

which can be rewritten as

$$ u(t, x, y) = u(t_s, x, y) + (t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \{[P_{2,j}(x) - xq_2(j)][P_{2,l}(y) - yq_2(l)]\}, \quad (31) $$

Derivative the Equation (31), we obtain that

$$ \frac{\partial u}{\partial t} = \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \{[P_{2,j}(x) - xq_2(j)][P_{2,l}(y) - yq_2(l)]\}, \quad (32) $$

$$ \frac{\partial^2 u}{\partial x^2} = \left. \frac{\partial^2 u}{\partial x^2} \right|_{t=t_s} + (t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \{h_j(x)[P_{2,l}(y) - yq_2(l)]\}, \quad (33) $$

$$ \frac{\partial^2 u}{\partial y^2} = \left. \frac{\partial^2 u}{\partial y^2} \right|_{t=t_s} + (t-t_s) \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \{[P_{2,j}(x) - xq_2(j)]h_l(y)\}, \quad (34) $$
Substituting Equation (32)-(34) in (16) for any collocation points $x_r, y_k$ with $r \in \{1, 2, \dots, 2M\}$, $k \in \{1, 2, \dots, 2M\}$ and replacing $t$ by $t_{s+1}$ and $\Delta t$ by the value $\Delta t = t_{s+1} - t_s$, we get

$$ \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} R_{j,l,r,k} = g(s,r,k) \quad \text{for } 1 \le r \le 2M, 1 \le k \le 2M, 1 \le s \le N \quad (35) $$

where

$$ R_{j,l,r,k} = \begin{bmatrix} P_2(j,r) - x_r q_2(j) & [P_2(l,k) - y_k q_2(l)] \\ -c^2 \Delta t H(j,r) & [P_2(l,k) - y_k q_2(l)] \\ -c^2 \Delta t [P_2(j,r) - x_r q_2(j)] & H(l,k) \end{bmatrix}, \quad (36) $$

and

$$ g(s,r,k) = f(t_s, x_r, y_k) + c^2 \left. \frac{\partial^2 u}{\partial x^2} \right|_{t=t_s, x=x_r, y=y_k} + c^2 \left. \frac{\partial^2 u}{\partial y^2} \right|_{t=t_s, x=x_r, y=y_k}, $$

The terms $\frac{\partial^2 u}{\partial x^2}|_{t=t_s, x=x_r, y=y_k}$ and $\frac{\partial^2 u}{\partial y^2}|_{t=t_s, x=x_r, y=y_k}$ are given by initial condition

and after this are calculated from the Equations (33) and (34) respectively.

Here we transform the system from the fourth-order matrices into a second-order matrices by the following:

Let

$$ \eta = 2M(j-1) + l, \quad \mu = 2M(r-1) + k, $$

Now Equation (35) obtains the form:

$$ \sum_{\eta=1}^{(2M)^2} S(\mu, \eta) B(\eta) = F(\mu), \quad \text{for } 1 \le \mu \le (2M)^2, \text{ we get the following system of} $$

linear equations

$$ S \cdot B = F \quad (37) $$

Here $B$ and $F$ are $(2M)^2$ vectors and $S$ is a $(2M)^2 \times (2M)^2$ matrix such that:

$$ B = \begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,2M} & a_{2,1} & \cdots & a_{2,2M} & \cdots & a_{2M,1} & a_{2M,2} & \cdots & a_{2M,2M} \end{bmatrix}_{l \times (2M)^2} $$

$$ F = \begin{bmatrix} g_{s,1,1} & g_{s,1,2} & \cdots & g_{s,1,2M} & g_{s,2,1} & \cdots & g_{s,2,2M} & \cdots & g_{s,2M,1} & \cdots & g_{s,2M,2M} \end{bmatrix}_{l \times (2M)^2} $$

$$ S = \begin{bmatrix} R_{1,1,1,1} & R_{1,2,1,1} & \cdots & R_{1,2M,1,1} & R_{2,1,1,1} & \cdots & R_{2,2M,1,1} & \cdots & R_{2M,2M,1,1} \\ R_{1,1,1,2} & R_{1,2,1,2} & \cdots & R_{1,2M,1,2} & R_{2,1,1,2} & \cdots & R_{2,2M,1,2} & \cdots & R_{2M,2M,1,2} \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \vdots \\ R_{1,1,2M,2M} & R_{1,2,2M,2M} & \cdots & R_{1,2M,2M,2M} & R_{2,1,2M,2M} & \cdots & R_{2M,1,2M,2M} & \cdots & R_{2M,2M,2M,2M} \end{bmatrix}_{(2M)^2 \times (2M)^2} $$

where $R_{i,l,i,r,k,s}$ are calculated according to Equation (36). After solving system

(37) we obtain the wavelet coefficients $a_{j,l,i}$ and thus for any

$(t, x, y) \in (0,1) \times \Omega$ the solution $u(t, x, y)$ is obtained from Equation (31).
Substituting Equation (32)-(34) in (16) for any collocation points $x_r, y_k$ with $r \in \{1, 2, \dots, 2M\}$, $k \in \{1, 2, \dots, 2M\}$ and replacing $t$ by $t_{s+1}$ and $\Delta t$ by the value $\Delta t = t_{s+1} - t_s$, we get

$$ \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} R_{j,l,r,k} = g(s,r,k) \quad \text{for } 1 \le r \le 2M, 1 \le k \le 2M, 1 \le s \le N \quad (35) $$

where

$$ R_{j,l,r,k} = \begin{bmatrix} P_2(j,r) - x_r q_2(j) & [P_2(l,k) - y_k q_2(l)] \\ -c^2 \Delta t H(j,r) & [P_2(l,k) - y_k q_2(l)] \\ -c^2 \Delta t [P_2(j,r) - x_r q_2(j)] & H(l,k) \end{bmatrix} \quad (36) $$

and

$$ g(s,r,k) = f(t_s, x_r, y_k) + c^2 \left. \frac{\partial^2 u}{\partial x^2} \right|_{t=t_s, x=x_r, y=y_k} + c^2 \left. \frac{\partial^2 u}{\partial y^2} \right|_{t=t_s, x=x_r, y=y_k}, $$

The terms $\frac{\partial^2 u}{\partial x^2}|_{t=t_s, x=x_r, y=y_k}$ and $\frac{\partial^2 u}{\partial y^2}|_{t=t_s, x=x_r, y=y_k}$ are given by initial condition

and after this are calculated from the Equations (33) and (34) respectively.

Here we transform the system from the fourth-order matrices into a second-order matrices by the following:

Let

$$ \eta = 2M(j-1) + l, \quad \mu = 2M(r-1) + k, $$

Now Equation (35) obtains the form:

$$ \sum_{\eta=1}^{(2M)^2} S(\mu, \eta) B(\eta) = F(\mu), \quad \text{for } 1 \le \mu \le (2M)^2, \text{ we get the following system of} $$

linear equations

$$ S \cdot B = F \quad (37) $$

Here $B$ and $F$ are $(2M)^2$ vectors and $S$ is a $(2M)^2 \times (2M)^2$ matrix such that:

$$ B = \begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,2M} & a_{2,1} & \cdots & a_{2,2M} & \cdots & a_{2M,1} & a_{2M,2} & \cdots & a_{2M,2M} \end{bmatrix}_{l \times (2M)^2} $$

$$ F = \begin{bmatrix} g_{s,1,1} & g_{s,1,2} & \cdots & g_{s,1,2M} & g_{s,2,1} & \cdots & g_{s,2,2M} & \cdots & g_{s,2M,1} & \cdots & g_{s,2M,2M} \end{bmatrix}_{l \times (2M)^2} $$

$$ S = \begin{bmatrix} R_{1,1,1,1} & R_{1,2,1,1} & \cdots & R_{1,2M,1,1} & R_{2,1,1,1} & \cdots & R_{2,2M,1,1} & \cdots & R_{2M,2M,1,1} \\ R_{1,1,1,2} & R_{1,2,1,2} & \cdots & R_{1,2M,1,2} & R_{2,1,1,2} & \cdots & R_{2,2M,1,2} & \cdots & R_{2M,2M,1,2} \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \vdots \\ R_{1,1,2M,2M} & R_{1,2,2M,2M} & \cdots & R_{1,2M,2M,2M} & R_{2,1,2M,2M} & \cdots & R_{2M,1,2M,2M} & \cdots & R_{2M,2M,2M,2M} \end{bmatrix}_{(2M)^2 \times (2M)^2} $$

where $R_{i,l,i,r,k,s}$ are calculated according to Equation (36). After solving system

(37) we obtain the wavelet coefficients $a_{j,l,i}$ and thus for any

$(t, x, y) \in (0,1) \times \Omega$ the solution $u(t, x, y)$ is obtained from Equation (31).
### 4.1.3. Numerical Results

Taking $c = 1$, and

$$f(t, x, y) = x^2(x-1)y^2(y-1)(3t^2 - 2t) - (6x-2)y^2(y-1)t^2(t-1)$$ such that
$$-x^2(x-1)(6y-2)t^2(t-1),$$

the exact solution for (16) is

$$u(t, x, y) = x^2(x-1)y^2(y-1)t^2(t-1), \quad (38)$$

In the following we use the MATLAB norm error

$$\delta_{2M} = \text{norm}(u - u_{\text{ex}}, 2)/2M,$$

The results for 2D and 3D Haar wavelets method are compared in Table 1 for different values 2M.

We observe from Table 1 that when $2M$ is not sufficiently large value, means that $\Delta t$ is not sufficiently small value then the error is big, we can obtain more precision by the 3D Haar wavelets than 2D Haar wavelets and in less time. We observe that the precision obtained by 3D Haar wavelets in the case $2M = 4$ cannot be obtained for 2D Haar wavelets even one takes $2M = 16$ and spends 100 times.

## 4.2. Poisson Equation

Consider 3D Poisson equation [21]

$$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} = f(x, y, z), \quad 0 \le x, y, z \le 1, \quad (39)$$

$u(x, y, z) = 0$ along the boundaries. Here we have $\lambda = 2, \beta = 2$ and $\alpha = 2$ and suppose that $M_1 = M_2 = M_3 = M$.

### 4.2.1. The Solution by 3D Haar Wavelet Method

Now the solution by 3D Haar wavelet is begin using the Equation (12) to approximate problem (16) and considering the initial and boundary conditions at $x = 0, y = 0$ and $z = 0$, we get

$$\begin{aligned} u(x, y, z) = & \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,j} P_{2,j}(z) P_{2,l}(x) P_{2,l}(y) + y \frac{\partial u(x, y, z)}{\partial y} \Big|_{y=0} \\ & + x \frac{\partial u(x, y, z)}{\partial x} \Big|_{x=0} - xy \frac{\partial^2 u(x, y, z)}{\partial x \partial y} \Big|_{x=y=0} + z \frac{\partial u(x, y, z)}{\partial z} \Big|_{z=0} \\ & - xz \frac{\partial^2 u(x, y, z)}{\partial x \partial z} \Big|_{x=z=0} - yz \frac{\partial^2 u(x, y, z)}{\partial y \partial z} \Big|_{y=z=0} + xyz \frac{\partial^3 u(x, y, z)}{\partial x \partial y \partial z} \Big|_{x=y=z=0} \end{aligned} \quad (40)$$

We use the same technique as in (1) to find the unknown terms in Equation (40).

This done in three steps:

1) Substitute the boundary condition when $x = 1$ in Equation (40) and replace this result back into (40).
### 4.1.3. Numerical Results

Taking $c = 1$, and

$$f(t, x, y) = x^2 (x-1) y^2 (y-1) (3t^2 - 2t) - (6x-2) y^2 (y-1) t^2 (t-1)$$ such that
$$-x^2 (x-1) (6y-2) t^2 (t-1),$$

the exact solution for (16) is

$$u(t, x, y) = x^2 (x-1) y^2 (y-1) t^2 (t-1), \quad (38)$$

In the following we use the MATLAB norm error

$$\delta_{2M} = \text{norm}(u - u_{\text{ex}}, 2) / 2M,$$

The results for 2D and 3D Haar wavelets method are compared in Table 1 for different values 2M.

We observe from Table 1 that when $2M$ is not sufficiently large value, means that $\Delta t$ is not sufficiently small value then the error is big, we can obtain more precision by the 3D Haar wavelets than 2D Haar wavelets and in less time. We observe that the precision obtained by 3D Haar wavelets in the case $2M = 4$ cannot be obtained for 2D Haar wavelets even one takes $2M = 16$ and spends 100 times.

## 4.2. Poisson Equation

Consider 3D Poisson equation [21]

$$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} = f(x, y, z), \quad 0 \le x, y, z \le 1, \quad (39)$$

$u(x, y, z) = 0$ along the boundaries. Here we have $\lambda = 2, \beta = 2$ and $\alpha = 2$ and suppose that $M_1 = M_2 = M_3 = M$.

### 4.2.1. The Solution by 3D Haar Wavelet Method

Now the solution by 3D Haar wavelet is begin using the Equation (12) to approximate problem (16) and considering the initial and boundary conditions at $x = 0, y = 0$ and $z = 0$, we get

$$\begin{aligned} u(x, y, z) = & \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,j} P_{2,j}(z) P_{2,l}(x) P_{2,l}(y) + y \frac{\partial u(x, y, z)}{\partial y} \Big|_{y=0} \\ & + x \frac{\partial u(x, y, z)}{\partial x} \Big|_{x=0} - xy \frac{\partial^2 u(x, y, z)}{\partial x \partial y} \Big|_{x=y=0} + z \frac{\partial u(x, y, z)}{\partial z} \Big|_{z=0} \\ & - xz \frac{\partial^2 u(x, y, z)}{\partial x \partial z} \Big|_{x=z=0} - yz \frac{\partial^2 u(x, y, z)}{\partial y \partial z} \Big|_{y=z=0} + xyz \frac{\partial^3 u(x, y, z)}{\partial x \partial y \partial z} \Big|_{x=y=z=0} \end{aligned} \quad (40)$$

We use the same technique as in (1) to find the unknown terms in Equation (40).

This done in three steps:

1) Substitute the boundary condition when $x = 1$ in Equation (40) and replace this result back into (40).
<table><caption>**Table 1.** Compared the solution of the 2D Heat equation when $\Delta t = 1/2M$ $t = (4M - 1)/(4M)$.</caption><thead><tr><th>2M</th><th>The method</th><th>δ<sub>2M</sub></th><th>Time (s)</th></tr></thead><tbody><tr><td>4</td><td>2D Haar</td><td>3.5900e-004</td><td>1.4700</td></tr><tr><td>4</td><td>3D Haar</td><td>8.7770e-006</td><td>4.2900</td></tr><tr><td>8</td><td>2D Haar</td><td>1.4686e-004</td><td>11.2140</td></tr><tr><td>8</td><td>3D Haar</td><td>1.6790e-006</td><td>33.3030</td></tr><tr><td>16</td><td>2D Haar</td><td>4.0944e-005</td><td>87.1730</td></tr><tr><td>16</td><td>3D Haar</td><td>4.2036e-007</td><td>281.6000</td></tr></tbody></table>

2) Substitute the boundary condition when $y = 1$ in equation resulting 1 and replace this result back into 1.

3) Substitute the boundary condition when $z = 1$ in equation resulting 2 and replace this result back into 2.

Next the boundary conditions are satisfied

$$u(x,y,z) = \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,j} \{[P_{2,i}(z) - zq_2(i)][P_{2,j}(x) - xq_2(j)][P_{2,l}(y) - yq_2(l)]\}, \quad (41)$$

Derivative the Equation (41) we obtain that:

$$\frac{\partial^2 u}{\partial x^2} = \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,j} \{[P_{2,i}(z) - zq_2(i)]h_j(x)[P_{2,j}(y) - yq_2(l)]\},$$

$$\frac{\partial^2 u}{\partial y^2} = \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,j} \{[P_{2,i}(z) - zq_2(i)][P_{2,j}(x) - xq_2(j)]h_l(y)\},$$

$$\frac{\partial^2 u}{\partial z^2} = \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,j} \{h_l(z)[P_{2,j}(x) - xq_2(j)][P_{2,l}(y) - yq_2(l)]\},$$

We substituting above equations in (39) for any collocation points $x_r, y_k$ and $z_s$ with $s \in \{1,2,\dots,2M\}, r \in \{1,2,\dots,2M\}, k \in \{1,2,\dots,2M\}$, we get

$$\sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,j} R_{j,j,l,r,k,s} = g(s,r,k) \quad \text{for } 1 \le s \le 2M, 1 \le r \le 2M, 1 \le k \le 2M, \quad (42)$$

where

$$R_{j,j,j,r,k,s} = [P_2(i,s) - z_s q_2(i)] H(j,r) [P_2(l,k) - y_k q_2(l)] \\ + [P_2(i,s) - z_s q_2(i)] [P_2(j,r) - x_r q_2(j)] H(l,k) \\ + H(i,s) [P_2(j,r) - x_r q_2(j)] [P_2(l,k) - y_k q_2(l)], \quad (43)$$

and

$$g(s,r,k) = f(x_r, y_k, z_s),$$

After solving system (42) we obtain the wavelet coefficients $a_{j,l,i}$ and thus for any $(t,x,y) \in (0,1) \times \Omega$ the solution $u(t,x,y)$ is obtained from Equation (41).
<table><caption>**Table 1.** Compared the solution of the 2D Heat equation when $\Delta t = 1/2M$ $t = (4M - 1)/(4M)$.</caption><thead><tr><th>2M</th><th>The method</th><th>δ<sub>2M</sub></th><th>Time (s)</th></tr></thead><tbody><tr><td>4</td><td>2D Haar</td><td>3.5900e-004</td><td>1.4700</td></tr><tr><td>4</td><td>3D Haar</td><td>8.7770e-006</td><td>4.2900</td></tr><tr><td>8</td><td>2D Haar</td><td>1.4686e-004</td><td>11.2140</td></tr><tr><td>8</td><td>3D Haar</td><td>1.6790e-006</td><td>33.3030</td></tr><tr><td>16</td><td>2D Haar</td><td>4.0944e-005</td><td>87.1730</td></tr><tr><td>16</td><td>3D Haar</td><td>4.2036e-007</td><td>281.6000</td></tr></tbody></table>

2) Substitute the boundary condition when $y = 1$ in equation resulting 1 and replace this result back into 1.

3) Substitute the boundary condition when $z = 1$ in equation resulting 2 and replace this result back into 2.

Next the boundary conditions are satisfied:

$$u(x,y,z) = \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} \{ [P_{2,i}(z) - zq_2(i)] [P_{2,j}(x) - xq_2(j)] [P_{2,l}(y) - yq_2(l)] \}, \quad (41)$$

Derivative the Equation (41) we obtain that:

$$\frac{\partial^2 u}{\partial x^2} = \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} \{ [P_{2,i}(z) - zq_2(i)] h_j(x) [P_{2,j}(y) - yq_2(l)] \},$$

$$\frac{\partial^2 u}{\partial y^2} = \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} \{ [P_{2,i}(z) - zq_2(i)] [P_{2,j}(x) - xq_2(j)] h_l(y) \},$$

$$\frac{\partial^2 u}{\partial z^2} = \sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} \{ h_l(z) [P_{2,j}(x) - xq_2(j)] [P_{2,l}(y) - yq_2(l)] \},$$

We substituting above equations in (39) for any collocation points $x_r, y_k$ and $z_s$ with $s \in \{1, 2, \dots, 2M\}, r \in \{1, 2, \dots, 2M\}, k \in \{1, 2, \dots, 2M\}$, we get:

$$\sum_{j=1}^{2M} \sum_{l=1}^{2M} \sum_{i=1}^{2M} a_{j,l,i} R_{j,j,l,r,k,s} = g(s,r,k) \quad \text{for } 1 \le s \le 2M, 1 \le r \le 2M, 1 \le k \le 2M, \quad (42)$$

where:

$$R_{j,j,j,r,k,s} = \left[ \begin{aligned} & P_2(i,s) - z_s q_2(i) H(j,r) [P_2(l,k) - y_k q_2(l)] \\ & + P_2(i,s) - z_s q_2(i) [P_2(j,r) - x_r q_2(j)] H(l,k) \\ & + H(i,s) [P_2(j,r) - x_r q_2(j)] [P_2(l,k) - y_k q_2(l)] \end{aligned} \right], \quad (43)$$

and:

$$g(s,r,k) = f(x_r, y_k, z_s),$$

After solving system (42) we obtain the wavelet coefficients $a_{j,l,i}$ and thus for any $(t, x, y) \in (0,1) \times \Omega$ the solution $u(t,x,y)$ is obtained from Equation (41).
### 4.2.2. The Solution by 2D Haar Wavelet Method

For the solution by 2D Haar wavelet method we need to divide one of the intervals $x \in [A_1, B_1]$, $y \in [A_2, B_2]$, $z \in [A_3, B_3]$ into N equal parts and we will divide the interval $z \in [A_3, B_3]$ into N equal parts of length $\Delta z = (B_3 - A_3)/N$ and denote to $z_s = (s-1)\Delta z$ $s=1,2,\dots,N$.

By using Equation (15) and considering the initial-boundary conditions when $x=0, y=0$ and $z=0$, we get

$$ 
\begin{aligned}
u(x,y,z) = & \frac{(z-z_s)^2}{2} \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} P_{2,j}(x) P_{2,l}(y) + u(x,y,z_s) \\
& + y \left( \left. \frac{\partial u(x,y,z)}{\partial y} \right|_{y=0} - \left. \frac{\partial u(x,y,z_s)}{\partial y} \right|_{y=0} \right) \\
& + x \left( \left. \frac{\partial u(x,y,z)}{\partial x} \right|_{x=0} - \left. \frac{\partial u(x,y,z_s)}{\partial x} \right|_{x=0} \right) \\
& - xy \left( \left. \frac{\partial^2 u(x,y,z)}{\partial x \partial y} \right|_{x=y=0} - \left. \frac{\partial^2 u(x,y,z_s)}{\partial x \partial y} \right|_{x=y=0} \right) + z \frac{\partial u(x,y,z_s)}{\partial z} \\
& - xz \left. \frac{\partial^2 u(x,y,z_s)}{\partial x \partial z} \right|_{x=0} - yz \left. \frac{\partial^2 u(x,y,z_s)}{\partial y \partial z} \right|_{y=0} + xyz \frac{\partial^3 u(x,y,z_s)}{\partial x \partial y \partial z} \Bigg|_{x=y=0} ,
\end{aligned}
\quad (44) $$

where the element $a_{j,l}$ is constant in the subinterval $z \in (z_s, z_{s+1}]$.

Next the boundary conditions when $x=1, y=1$ and $z=1$ are satisfied

$$ 
\begin{aligned}
u(x,y,z) = & \left[ 1 - \frac{(z-z_s)}{(1-z_s)} \right] u(x,y,z_s) + \left[ \frac{(z-z_s)^2}{2} - \frac{(z-z_s)(1-z_s)}{2} \right] \\
& \times \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \{ [P_{2,l}(x) - xq_2(j)] [P_{2,l}(y) - yq_2(l)] \},
\end{aligned}
\quad (45) $$

from Equation (45), we get

$$ 
\begin{aligned}
\frac{\partial^2 u}{\partial x^2} = & \left[ 1 - \frac{(z-z_s)}{(1-z_s)} \right] \left. \frac{\partial^2 u}{\partial x^2} \right|_{z=z_s} + \left[ \frac{(z-z_s)^2}{2} - \frac{(z-z_s)(1-z_s)}{2} \right] \\
& \times \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \{ h_j(x) [P_{2,l}(y) - yq_2(l)] \},
\end{aligned}
$$

$$ 
\begin{aligned}
\frac{\partial^2 u}{\partial y^2} = & \left[ 1 - \frac{(z-z_s)}{(1-z_s)} \right] \left. \frac{\partial^2 u}{\partial y^2} \right|_{z=z_s} + \left[ \frac{(z-z_s)^2}{2} - \frac{(z-z_s)(1-z_s)}{2} \right] \\
& \times \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \{ [P_{2,l}(x) - xq_2(j)] h_l(y) \},
\end{aligned}
$$

$$ 
\frac{\partial^2 u}{\partial z^2} = \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \{ [P_{2,l}(x) - xq_2(j)] [P_{2,l}(y) - yq_2(l)] \}, $$

We substituting above equation in (39) for any collocation points $x_r, y_k$ with $r \in \{1, 2, \dots, 2M\}, k \in \{1, 2, \dots, 2M\}$ and replacing $z$ by $z_{s+1}$ and $\Delta z$ by the value $\Delta z = z_{s+1} - z_s$, we get
### 4.2.2. The Solution by 2D Haar Wavelet Method

For the solution by 2D Haar wavelet method we need to divide one of the intervals $x \in [A_1, B_1]$, $y \in [A_2, B_2]$, $z \in [A_3, B_3]$ into N equal parts and we will divide the interval $z \in [A_3, B_3]$ into N equal parts of length $\Delta z = (B_3 - A_3)/N$ and denote to $z_s = (s-1)\Delta z$ $s = 1, 2, \dots, N$.

By using Equation (15) and considering the initial-boundary conditions when $x=0, y=0$ and $z=0$, we get

$$ 
\begin{aligned}
u(x, y, z) = & \frac{(z - z_s)^2}{2} \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} P_{2,j}(x) P_{2,l}(y) + u(x, y, z_s) \\
& + y \left( \left. \frac{\partial u(x, y, z)}{\partial y} \right|_{y=0} - \left. \frac{\partial u(x, y, z_s)}{\partial y} \right|_{y=0} \right) \\
& + x \left( \left. \frac{\partial u(x, y, z)}{\partial x} \right|_{x=0} - \left. \frac{\partial u(x, y, z_s)}{\partial x} \right|_{x=0} \right) \\
& - xy \left( \left. \frac{\partial^2 u(x, y, z)}{\partial x \partial y} \right|_{x=y=0} - \left. \frac{\partial^2 u(x, y, z_s)}{\partial x \partial y} \right|_{x=y=0} \right) + z \frac{\partial u(x, y, z_s)}{\partial z} \\
& - xz \left. \frac{\partial^2 u(x, y, z_s)}{\partial x \partial z} \right|_{x=0} - yz \left. \frac{\partial^2 u(x, y, z_s)}{\partial y \partial z} \right|_{y=0} + xyz \frac{\partial^3 u(x, y, z_s)}{\partial x \partial y \partial z} \Bigg|_{x=y=0}
\end{aligned}
\quad (44) $$

where the element $a_{j,l}$ is constant in the subinterval $z \in (z_s, z_{s+1}]$.

Next the boundary conditions when $x=1, y=1$ and $z=1$ are satisfied

$$ 
\begin{aligned}
u(x, y, z) = & \left[ 1 - \frac{(z - z_s)}{(1 - z_s)} \right] u(x, y, z_s) + \left[ \frac{(z - z_s)^2}{2} - \frac{(z - z_s)(1 - z_s)}{2} \right] \\
& \times \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \left\{ [P_{2,l}(x) - xq_2(l)] [P_{2,l}(y) - yq_2(l)] \right\},
\end{aligned}
\quad (45) $$

from Equation (45), we get

$$ 
\begin{aligned}
\frac{\partial^2 u}{\partial x^2} = & \left[ 1 - \frac{(z - z_s)}{(1 - z_s)} \right] \left. \frac{\partial^2 u}{\partial x^2} \right|_{z=z_s} + \left[ \frac{(z - z_s)^2}{2} - \frac{(z - z_s)(1 - z_s)}{2} \right] \\
& \times \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \left\{ h_j(x) [P_{2,l}(y) - yq_2(l)] \right\},
\end{aligned}
$$

$$ 
\begin{aligned}
\frac{\partial^2 u}{\partial y^2} = & \left[ 1 - \frac{(z - z_s)}{(1 - z_s)} \right] \left. \frac{\partial^2 u}{\partial y^2} \right|_{z=z_s} + \left[ \frac{(z - z_s)^2}{2} - \frac{(z - z_s)(1 - z_s)}{2} \right] \\
& \times \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \left\{ [P_{2,l}(x) - xq_2(l)] h_l(y) \right\},
\end{aligned}
$$

$$ 
\frac{\partial^2 u}{\partial z^2} = \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} \left\{ [P_{2,l}(x) - xq_2(l)] [P_{2,l}(y) - yq_2(l)] \right\}, $$

We substituting above equation in (39) for any collocation points $x_r, y_k$ with $r \in \{1, 2, \dots, 2M\}, k \in \{1, 2, \dots, 2M\}$ and replacing $z$ by $z_{s+1}$ and $\Delta z$ by the value $\Delta z = z_{s+1} - z_s$, we get
$$ \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} R_{j,l,r,k} = g(r,k,s) \quad \text{for } 1 \le r \le 2M, 1 \le k \le 2M, 1 \le s \le N, \quad (46) $$

where

$$ R_{j,l,r,k} = \left[ \begin{aligned} & P_2(j,r) - x_r q_2(j) [P_2(l,k) - y_k q_2(l)] \\ & + \left[ \frac{(\Delta z)^2}{2} - \frac{\Delta z (1-z_s)}{2} \right] H(j,r) [P_2(l,k) - y_k q_2(l)] \\ & + \left[ \frac{(\Delta z)^2}{2} - \frac{\Delta z (1-z_s)}{2} \right] [P_2(j,r) - x_r q_2(j)] H(l,k), \end{aligned} \right] \quad (47) $$

and

$$ g(r,k,s) = f(x_r, y_k, z_s) - \left[ 1 - \frac{\Delta z}{(1-z_s)} \left. \frac{\partial^2 u}{\partial x^2} \right|_{z=z_s, x=x_r, y=y_k} - \left[ 1 - \frac{\Delta z}{(1-z_s)} \left. \frac{\partial^2 u}{\partial y^2} \right|_{z=z_s, x=x_r, y=y_k} \right] \right], $$

The solution $u(x, y, z)$ get it from the Equation (45).

### 4.2.3. Numerical Results

Solve (39) for $f(x, y, z) = \sin(\pi x) \sin(\pi y) \sin(\pi z)$ is:

$$ u(x, y, z) = \frac{-1}{3\pi^2} \sin(\pi x) \sin(\pi y) \sin(\pi z), $$

Also we use the MATLAB norm error $\delta_{2M} = \text{norm}(u - u_{ex}, 2)/2M$. We plotted in Figure 1 the error $\delta_{2M}$ for $2M = 16$ near from the last time $z = 1$ to illustrate the impact of error accumulation on the solution by 2D Haar wavelets.

Results obtained using 2D and 3D Haar wavelets method are compared in Table 2 for different values 2M.

In 2D Haar wavelets method, We can reduce the error and increase accuracy by increasing the subdivisions for the time (t) or the interval $z \in [A_3, B_3]$ in examples (1) and (2) and minimize the value of $\Delta t$ or $\Delta z$. With the installation's divisions for $x \in [A_1, B_1]$ and $y \in [A_2, B_2]$ according to formula of Haar wavelets, This helps increase the accuracy of the solution in the 2D Haar wavelets method, But it requires more time to get on the solution as shown in Table 3.

All computation was made by using MATLAB Language, Intel®, core™ I3-2330M CPU, 2.00 GB (Memory), 2.20 GHz (Processor).

## 5. Conclusions

In this paper, we develop an accurate and efficient Haar wavelets method for solving three dimensional PDEs and time depending PDEs. The benefits of the Haar wavelets approach are sparse matrices of representation, fast transformation and possibility of implementation of fast algorithms. It's worth mentioning that the Haar wavelet solution provides excellent results even for small values of
$$ \sum_{j=1}^{2M} \sum_{l=1}^{2M} a_{j,l} R_{j,l,r,k} = g(r,k,s) \quad \text{for } 1 \le r \le 2M, 1 \le k \le 2M, 1 \le s \le N, \quad (46) $$

where

$$ R_{j,l,r,k} = \left[ \begin{aligned} & P_2(j,r) - x_r q_2(j) [P_2(l,k) - y_k q_2(l)] \\ & + \left[ \frac{(\Delta z)^2}{2} - \frac{\Delta z (1-z_s)}{2} \right] H(j,r) [P_2(l,k) - y_k q_2(l)] \\ & + \left[ \frac{(\Delta z)^2}{2} - \frac{\Delta z (1-z_s)}{2} \right] [P_2(j,r) - x_r q_2(j)] H(l,k), \end{aligned} \right] \quad (47) $$

and

$$ g(r,k,s) = f(x_r, y_k, z_s) - \left[ 1 - \frac{\Delta z}{(1-z_s)} \left. \frac{\partial^2 u}{\partial x^2} \right|_{z=z_s, x=x_r, y=y_k} - \left[ 1 - \frac{\Delta z}{(1-z_s)} \left. \frac{\partial^2 u}{\partial y^2} \right|_{z=z_s, x=x_r, y=y_k} \right] \right], $$

The solution $u(x, y, z)$ get it from the Equation (45).

### 4.2.3. Numerical Results

Solve (39) for $f(x, y, z) = \sin(\pi x)\sin(\pi y)\sin(\pi z)$ is:

$$ u(x, y, z) = \frac{-1}{3\pi^2} \sin(\pi x) \sin(\pi y) \sin(\pi z), $$

Also we use the MATLAB norm error $\delta_{2M} = \text{norm}(u - u_{ex}, 2)/2M$. We plotted in Figure 1 the error $\delta_{2M}$ for $2M = 16$ near from the last time $z = 1$ to illustrate the impact of error accumulation on the solution by 2D Haar wavelets.

Results obtained using 2D and 3D Haar wavelets method are compared in Table 2 for different values 2M.

In 2D Haar wavelets method, We can reduce the error and increase accuracy by increasing the subdivisions for the time (t) or the interval $z \in [A_3, B_3]$ in examples (1) and (2) and minimize the value of $\Delta t$ or $\Delta z$. With the installation's divisions for $x \in [A_1, B_1]$ and $y \in [A_2, B_2]$ according to formula of Haar wavelets, This helps increase the accuracy of the solution in the 2D Haar wavelets method, But it requires more time to get on the solution as shown in Table 3.

All computation was made by using MATLAB Language, Intel®, core™ I3-2330M CPU, 2.00 GB (Memory), 2.20 GHz (Processor).

## 5. Conclusions

In this paper, we develop an accurate and efficient Haar wavelets method for solving three dimensional PDEs and time depending PDEs. The benefits of the Haar wavelets approach are sparse matrices of representation, fast transformation and possibility of implementation of fast algorithms. It's worth mentioning that the Haar wavelet solution provides excellent results even for small values of
<figure><img src="image_3.png" /><figcaption>Figure. 1. The error $\delta_{16}$ at $z = 31/32$ by using 3D and 2D Haar wavelets method.</figcaption></figure>

<table><caption>Table 2. Compared the solution of the 3D Poisson equation when $\Delta z = 1/2M$ and $z = (4M - 1)/(4M)$.</caption><thead><tr><th>2M</th><th>The method</th><th>δ<sub>2M</sub></th><th>Time (s)</th></tr></thead><tbody><tr><td>4</td><td>2D Haar</td><td>3.1669e-004</td><td>1.4000</td></tr><tr><td>4</td><td>3D Haar</td><td>7.7797e-005</td><td>3.2800</td></tr><tr><td>8</td><td>2D Haar</td><td>1.8481e-004</td><td>10.6220</td></tr><tr><td>8</td><td>3D Haar</td><td>1.0420e-005</td><td>25.5220</td></tr><tr><td>16</td><td>2D Haar</td><td>5.3391e-005</td><td>84.2600</td></tr><tr><td>16</td><td>3D Haar</td><td>1.3243e-006</td><td>251.632</td></tr></tbody></table>

<table><caption>Table 3. Illustrates the convergence of the solution of the 3D Poisson equation by using the 2D Haar wavelets method with $2M = 8$ and for different values $\Delta z$.</caption><thead><tr><th>Δz</th><th>δ<sub>2M</sub></th><th>Time (s)</th></tr></thead><tbody><tr><td>1/(2M)</td><td>1.8481e-004</td><td>10.6220</td></tr><tr><td>1/(4M)</td><td>5.2623e-005</td><td>21.6400</td></tr><tr><td>1/(8M)</td><td>1.2907e-005</td><td>42.7320</td></tr></tbody></table>

$(2M)$ as note in Table 1 and Table 2. Also when $2M = 32$, $2M = 64$ ..., we can obtain the results closer to the exact values.

Also we compare between 2D Haar wavelets method and 3D Haar wavelets method in numerical solution for 3D PDEs, and we have found that 3D Haar wavelets are better and closer to the exact solution from 2D Haar wavelets as shown in Table 1 and Table 2.
<figure><img src="image_3.png" /><figcaption>Figure. 1. The error $\delta_{16}$ at $z = 31/32$ by using 3D and 2D Haar wavelets method.</figcaption></figure>

<table><caption>Table 2. Compared the solution of the 3D Poisson equation when $\Delta z = 1/2M$ and $z = (4M - 1)/(4M)$.</caption><thead><tr><td>2M</td><td>The method</td><td>δ<sub>2M</sub></td><td>Time (s)</td></tr></thead><tbody><tr><td>4</td><td>2D Haar</td><td>3.1669e-004</td><td>1.4000</td></tr><tr><td>4</td><td>3D Haar</td><td>7.7797e-005</td><td>3.2800</td></tr><tr><td>8</td><td>2D Haar</td><td>1.8481e-004</td><td>10.6220</td></tr><tr><td>8</td><td>3D Haar</td><td>1.0420e-005</td><td>25.5220</td></tr><tr><td>16</td><td>2D Haar</td><td>5.3391e-005</td><td>84.2600</td></tr><tr><td>16</td><td>3D Haar</td><td>1.3243e-006</td><td>251.632</td></tr></tbody></table>

<table><caption>Table 3. Illustrates the convergence of the solution of the 3D Poisson equation by using the 2D Haar wavelets method with $2M = 8$ and for different values $\Delta z$.</caption><thead><tr><td>Δz</td><td>δ<sub>2M</sub></td><td>Time (s)</td></tr></thead><tbody><tr><td>1/(2M)</td><td>1.8481e-004</td><td>10.6220</td></tr><tr><td>1/(4M)</td><td>5.2623e-005</td><td>21.6400</td></tr><tr><td>1/(8M)</td><td>1.2907e-005</td><td>42.7320</td></tr></tbody></table>

$(2M)$ as note in Table 1 and Table 2. Also when $2M = 32, 2M = 64, \dots$, we can obtain the results closer to the exact values.

Also we compare between 2D Haar wavelets method and 3D Haar wavelets method in numerical solution for 3D PDEs, and we have found that 3D Haar wavelets are better and closer to the exact solution from 2D Haar wavelets as shown in Table 1 and Table 2.
The main benefits of the proposed method are its simplicity (already a small number of grid points guarantee the necessary accuracy) and universality (the same approach is applicable for a wide class of PDEs). The method is very convenient for solving boundary value problems, since the boundary conditions are taken into account automatically. For numerical calculations useful are the matrix programs of MATLAB. The most time-consuming procedure is to calculate the integrals (4). In this paper only linear problems were considered, but the method is applicable also for nonlinear PDEs.

## Acknowledgements

This work was carried out when third author was visiting the mathematics department at Nantes University, France, with the support of ministry of higher education and scientific research, Iraq. The authors Ekhlass S. and Qasim F. would like to express their gratitude to Professor Abdeljalil Nachaoui for introducing to the subject and for giving all necessary support to complete this work and for his hospitality.

## References

[1] Zhi, S., Deng, L.-Y. and Qing, J.C. (2007) Numerical Solution of Differential Equations by Using Haar Wavelets. *Proceeding of the International Conference on Wavelet Analysis and pattern Recognition*, 2-4 NovEMBER 2007, Beijing, China, 1039-1044. https://doi.org/10.1109/ICWAPR.2007.4421585

[2] Bertoluzza, S. (1977) An Adaptive Collocation Method Based on Interpolating Wavelets. In: Dahmen, W., Kurdila, A.J. and Oswald, P., Eds., *Multi-Scale Wavelet Methods for Partial Differential Equations*, Academic Press, San Diego, 109-135.

[3] Beylkin, G. and Keiser, J.M. (1977) An Adaptive Pseudo-Wavelet Approach for Solving Nonlinear Partial Differential Equations. In: Dahmen, W., Kurdila, A.J. and Oswald, P., Eds., *Multi-Scale Wavelet Methods for Partial Differential Equations*, Academic Press, San Diego, 137-197.

[4] Chen, X., Xiang, J., Li, B. and He, Z. (2010) A Study of Multiscale Wavelet-Based Elements for Adaptive Finite Element Analysis. *Advances in Engineering Software*, **41**, 196-205. https://doi.org/10.1016/j.advengsoft.2009.09.008

[5] Hariharan, G., Kannan, K. and Sharma, K.R. (2009) Haar Wavelet Method for Solving Fisher's Equation. *Applied Mathematics and Computation*, **211**, 284-292.
https://doi.org/10.1016/j.amc.2008.12.089

[6] Arora, S., Singh, I., Brar, Y.S. and Kumar, S. (2015) Comparative Study of Haar Wavelet with Numerical Methods for Partial Differential Equations. *International Journal of Pure and Applied Mathematics*, **101**, 489-503.

[7] Arbabi, S., Nazari, A. and Darvishi, M.T. (2017) A Two-Dimensional Haar Wavelets Method for Solving Systems of PDEs. *Applied Mathematics and Computation*, **292**, 33-46. https://doi.org/10.1016/j.amc.2016.07.032

[8] Wang, X.J., Nan, B., Zhu, J. and Koeppé, R. (2014) Regularized 3d Functional Regression for Brainimage Data via Haar Wavelets. *The Annals of Applied Statistics*, **8**, 1045-1064. https://doi.org/10.1214/14-AOAS736

[9] Mohammadi, F. (2016) Numerical Solution of Stochastic Volterra Fredholm Integral Equations Using Haar Wavelets. *Scientific Bulletin “Politehnica” University*
The main benefits of the proposed method are its simplicity (already a small number of grid points guarantee the necessary accuracy) and universality (the same approach is applicable for a wide class of PDEs). The method is very convenient for solving boundary value problems, since the boundary conditions are taken into account automatically. For numerical calculations useful are the matrix programs of MATLAB. The most time-consuming procedure is to calculate the integrals (4). In this paper only linear problems were considered, but the method is applicable also for nonlinear PDEs.

## Acknowledgements

This work was carried out when third author was visiting the mathematics department at Nantes University, France, with the support of ministry of higher education and scientific research, Iraq. The authors Ekhlass S. and Qasim F. would like to express their gratitude to Professor Abdeljalil Nachaoui for introducing to the subject and for giving all necessary support to complete this work and for his hospitality.

## References

[1] Zhi, S., Deng, L.-Y. and Qing, J.C. (2007) Numerical Solution of Differential Equations by Using Haar Wavelets. *Proceeding of the International Conference on Wavelet Analysis and pattern Recognition*, 2-4 NovEMBER 2007, Beijing, China, 1039-1044. https://doi.org/10.1109/ICWAPR.2007.4421585

[2] Bertoluzza, S. (1977) An Adaptive Collocation Method Based on Interpolating Wavelets. In: Dahmen, W., Kurdila, A.J. and Oswald, P., Eds., *Multi-Scale Wavelet Methods for Partial Differential Equations*, Academic Press, San Diego, 109-135.

[3] Beylkin, G. and Keiser, J.M. (1977) An Adaptive Pseudo-Wavelet Approach for Solving Nonlinear Partial Differential Equations. In: Dahmen, W., Kurdila, A.J. and Oswald, P., Eds., *Multi-Scale Wavelet Methods for Partial Differential Equations*, Academic Press, San Diego, 137-197.

[4] Chen, X., Xiang, J., Li, B. and He, Z. (2010) A Study of Multiscale Wavelet-Based Elements for Adaptive Finite Element Analysis. *Advances in Engineering Software*, **41**, 196-205. https://doi.org/10.1016/j.advengsoft.2009.09.008

[5] Hariharan, G., Kannan, K. and Sharma, K.R. (2009) Haar Wavelet Method for Solving Fisher's Equation. *Applied Mathematics and Computation*, **211**, 284-292.
https://doi.org/10.1016/j.amc.2008.12.089

[6] Arora, S., Singh, I., Brar, Y.S. and Kumar, S. (2015) Comparative Study of Haar Wavelet with Numerical Methods for Partial Differential Equations. *International Journal of Pure and Applied Mathematics*, **101**, 489-503.

[7] Arbabi, S., Nazari, A. and Darvishi, M.T. (2017) A Two-Dimensional Haar Wavelets Method for Solving Systems of PDEs. *Applied Mathematics and Computation*, **292**, 33-46. https://doi.org/10.1016/j.amc.2016.07.032

[8] Wang, X.J., Nan, B., Zhu, J. and Koeppé, R. (2014) Regularized 3d Functional Regression for Brainimage Data via Haar Wavelets. *The Annals of Applied Statistics*, **8**, 1045-1064. https://doi.org/10.1214/14-AOAS736

[9] Mohammadi, F. (2016) Numerical Solution of Stochastic Volterra Fredholm Integral Equations Using Haar Wavelets. *Scientific Bulletin “Politehnica” University*
of Bucharest. Series A Applied Mathematics and Physics, **78**, 111-126.

[10] Chun, Z. and Zheng, Z. (2007) Three-Dimensional Analysis of Functionally Graded Plate Based on the Haar Wavelet Method. *Acta Mechanica Solida Sinica*, **20**, 95-102.
https://doi.org/10.1007/s10338-007-0711-3

[11] Majak, J., Pohlak, M., Eerme, M. and Lepikult, T. (2009) Weak Formulation Based Haar Wavelet Method for Solving Differential Equations. *Applied Mathematics and Computation*, **211**, 488-494. https://doi.org/10.1016/j.amc.2009.01.089

[12] Castro, L.M.S., Ferreira, A.J.M., Bertoluzza, S., Patra, R.C. and Reddy, J.N. (2010) A Wavelet Collocation Method for the Static Analysis of Sandwich Plates Using a Layerwise Theory. *Composite Structures*, **92**, 1786-1792.
https://doi.org/10.1016/j.compstruct.2010.01.021

[13] Chen, C. and Hsiao, C.H. (1997) Haar Wavelet Method for Solving Lumped and Distributed Parameter Systems. *IEE Proceedings-Control Theory and Applications*, **144**, 87-94. https://doi.org/10.1049/ip-cta:19970702

[14] Lepik, Ü. (2008) Solving Integral and Differential Equations by the Aid of Nonuniform Haar Wavelets. *Applied Mathematics and Computation*, **198**, 326-332.
https://doi.org/10.1016/j.amc.2007.08.036

[15] Lepik, Ü. (2008) Haar Wavelet Method for Solving Higher Order Differential Equations. *International Journal of Mathematics and Computation*, **1**, 84-94.

[16] Lepik, Ü. (2005) Numerical Solutions of Differential Equations Using Haar Wavelets. *Mathematics and Computers in Simulation*, **68**, 127-143.
https://doi.org/10.1016/j.matcom.2004.10.005

[17] Lepik, Ü. (2007) Numerical Solution of Evolution Equations by the Haar Wavelet Method. *Applied Mathematics and Computation*, **185**, 695-704.
https://doi.org/10.1016/j.amc.2006.07.077

[18] Lepik, Ü. (2011) Solving PDFs with the Aid of Two-Dimensional Haar Wavelets. *Computers and Mathematics with Applications*, **61**, 1873-1879.
https://doi.org/10.1016/j.camwa.2011.02.016

[19] AL-Rawi, E.S. and Qasim, A.F. (2014) CAS Wavelets for Solving General Two Dimensional Partial Differential Equations of Higher Order with Application. *International Journal of Enhanced Research in Science Technology & Engineering*, **3**, 496-507.

[20] Shiralashetti, S.C., Angadi, L.M., Deshi, A.B. and Kantli, M.H. (2016) Haar Wavelet Method for the Numerical Solution of Benjamin-Bona-Mahony Equations. *Journal of Information and Computing Science*, **11**, 136-145.

[21] Zhi, S., Cao, Y.-Y. and Qing, J.C. (2012) Solving 2D and 3D Poisson Equations and Biharmonic Equations by the Haar Wavelet Method. *Applied Mathematical Modeling*, **36**, 5143-5161. https://doi.org/10.1016/j.apm.2011.11.078
of Bucharest. Series A Applied Mathematics and Physics, **78**, 111-126.

[10] Chun, Z. and Zheng, Z. (2007) Three-Dimensional Analysis of Functionally Graded Plate Based on the Haar Wavelet Method. *Acta Mechanica Solida Sinica*, **20**, 95-102.
https://doi.org/10.1007/s10338-007-0711-3

[11] Majak, J., Pohlak, M., Eerme, M. and Lepikult, T. (2009) Weak Formulation Based Haar Wavelet Method for Solving Differential Equations. *Applied Mathematics and Computation*, **211**, 488-494. https://doi.org/10.1016/j.amc.2009.01.089

[12] Castro, L.M.S., Ferreira, A.J.M., Bertoluzza, S., Patra, R.C. and Reddy, J.N. (2010) A Wavelet Collocation Method for the Static Analysis of Sandwich Plates Using a Layerwise Theory. *Composite Structures*, **92**, 1786-1792.
https://doi.org/10.1016/j.compstruct.2010.01.021

[13] Chen, C. and Hsiao, C.H. (1997) Haar Wavelet Method for Solving Lumped and Distributed Parameter Systems. *IEE Proceedings-Control Theory and Applications*, **144**, 87-94. https://doi.org/10.1049/ip-cta:19970702

[14] Lepik, Ü. (2008) Solving Integral and Differential Equations by the Aid of Nonuniform Haar Wavelets. *Applied Mathematics and Computation*, **198**, 326-332.
https://doi.org/10.1016/j.amc.2007.08.036

[15] Lepik, Ü. (2008) Haar Wavelet Method for Solving Higher Order Differential Equations. *International Journal of Mathematics and Computation*, **1**, 84-94.

[16] Lepik, Ü. (2005) Numerical Solutions of Differential Equations Using Haar Wavelets. *Mathematics and Computers in Simulation*, **68**, 127-143.
https://doi.org/10.1016/j.matcom.2004.10.005

[17] Lepik, Ü. (2007) Numerical Solution of Evolution Equations by the Haar Wavelet Method. *Applied Mathematics and Computation*, **185**, 695-704.
https://doi.org/10.1016/j.amc.2006.07.077

[18] Lepik, Ü. (2011) Solving PDFs with the Aid of Two-Dimensional Haar Wavelets. *Computers and Mathematics with Applications*, **61**, 1873-1879.
https://doi.org/10.1016/j.camwa.2011.02.016

[19] AL-Rawi, E.S. and Qasim, A.F. (2014) CAS Wavelets for Solving General Two Dimensional Partial Differential Equations of Higher Order with Application. *International Journal of Enhanced Research in Science Technology & Engineering*, **3**, 496-507.

[20] Shiralashetti, S.C., Angadi, L.M., Deshi, A.B. and Kantli, M.H. (2016) Haar Wavelet Method for the Numerical Solution of Benjamin-Bona-Mahony Equations. *Journal of Information and Computing Science*, **11**, 136-145.

[21] Zhi, S., Cao, Y.-Y. and Qing, J.C. (2012) Solving 2D and 3D Poisson Equations and Biharmonic Equations by the Haar Wavelet Method. *Applied Mathematical Modeling*, **36**, 5143-5161. https://doi.org/10.1016/j.apm.2011.11.078