Haar wavelet method for the coupled degenerate reaction-diffusion PDEs and the ODEs having a non-linear source

TL;DR

This paper introduces a Haar wavelet numerical method for solving coupled degenerate reaction-diffusion PDEs and ODEs with nonlinear sources, demonstrating convergence and applying it to medically relevant models.

Contribution

The paper develops a Haar wavelet-based numerical scheme for complex coupled PDE-ODE systems with nonlinear sources, including convergence analysis and practical medical applications.

Findings

Successfully solved model problems of medical significance.

Demonstrated convergence of the numerical scheme.

Used GMRES solver for efficient linear system solutions.

Abstract

In this work, we propose the Haar wavelet method for the coupled degenerate reaction-diffusion PDEs and the ODEs having non-linear a source with Neumann boundary, applicable in various fields of the natural sciences, engineering, and economics, for example in gas dynamics, certain biological models, assets pricing in economics, composite media etc. Convergence analysis of the proposed numerical scheme has been carried out. We use the GMRES solver to solve the linear system of equations. Numerical solutions for the model problems of medical significance have been successfully solved.

Tables

Table 1. Table 1 : Absolute Error for v 𝑣 v at different points of the domain when J = 5 𝐽 5 J=5 and T=0.5

x		absolute error
	$dt={10}^{-2}$	$dt={10}^{-3}$	$dt={10}^{-4}$
0.0234	$2.83\times {10}^{-2}$	$2.0\times {10}^{-3}$	$1.34\times {10}^{-4}$
0.1172	$2.83\times {10}^{-2}$	$2.0\times {10}^{-3}$	$1.34\times {10}^{-4}$
0.2266	$2.82\times {10}^{-2}$	$2.0\times {10}^{-3}$	$1.33\times {10}^{-4}$
0.3828	$2.81\times {10}^{-2}$	$2.0\times {10}^{-3}$	$1.32\times {10}^{-4}$
0.5391	$2.81\times {10}^{-2}$	$2.0\times {10}^{-3}$	$1.32\times {10}^{-4}$
0.7734	$2.81\times {10}^{-2}$	$2.0\times {10}^{-3}$	$1.31\times {10}^{-4}$
0.9297	$2.81\times {10}^{-2}$	$2.0\times {10}^{-3}$	$1.3\times {10}^{-4}$

Table 2. Table 2 : Absolute Error for u e subscript 𝑢 𝑒 u_{e} at different points of the domain when J = 5 𝐽 5 J=5 and T=0.5

x		absolute error
	$dt={10}^{-2}$	$dt={10}^{-3}$	$dt={10}^{-4}$
0.0234	$2.2048\times {10}^{-5}$	$2.952\times {10}^{-6}$	$1.952\times {10}^{-6}$
0.1172	$1.8499\times {10}^{-4}$	$1.011\times {10}^{-6}$	$2.011\times {10}^{-6}$
0.2266	$5.1482\times {10}^{-4}$	$3.822\times {10}^{-6}$	$2.1780\times {10}^{-6}$
0.3828	$1.1\times {10}^{-3}$	$1.4401\times {10}^{-5}$	$1.599\times {10}^{-6}$
0.5391	$1.7\times {10}^{-3}$	$2.5823\times {10}^{-5}$	$1.77\times {10}^{-7}$
0.7734	$2.4\times {10}^{-3}$	$4.8658\times {10}^{-3}$	$6.58\times {10}^{-7}$
0.9297	$2.8\times {10}^{-3}$	$6.3827\times {10}^{-5}$	$8.27\times {10}^{-7}$

Table 3. Table 3 : Error for u e subscript 𝑢 𝑒 u_{e} and v 𝑣 v at T=1 with reference solution when J = 5 , d ​ t = 10 − 5 formulae-sequence 𝐽 5 𝑑 𝑡 superscript 10 5 J=5,dt=10^{-5} .

	$dt={10}^{-2}$	$dt={10}^{-3}$	$dt={10}^{-4}$
$L^{\infty }\left(ue\right)$ error	$1.17\times {10}^{-4}$	$1.5\times {10}^{-5}$	$1.0\times {10}^{-6}$
$L^{\infty }\left(v\right)$ error	$1.17\times {10}^{-4}$	$1.5\times {10}^{-5}$	$1.0\times {10}^{-6}$

Equations

$$h_i(x)=\{\begin{array}{ll}1& \text{when }{\beta }_1(i)\le x<{\beta }_2(i)\\ -1& \text{when }{\beta }_2(i)\le x<{\beta }_3(i)\\ 0& \text{otherwise},\end{array}$$

$${\beta }_1(i)=A+2k\zeta \Delta x,{\beta }_2(i)=A+(2k+1)\zeta \Delta x,$$

$${\beta }_3(i)=A+2(k+1)\zeta \Delta x,\zeta =M/m.$$

$$h_i(x)=\{\begin{array}{ll}1& \text{when }A\le x<B\\ 0& \text{otherwise}.\end{array}$$

$${\int }_A^B{h}_i(x)h_j(x)dx=\{\begin{array}{ll}{2}^{-j}(B-A)& \text{when }i=j\\ 0& \text{when }i\ne j.\end{array}$$

$$p_{\alpha ,i}(x)={\text{idotsint}}_A^B{h}_i(t)d{t}^{\alpha }=\frac{1}{(\alpha -1)!}{\int }_A^x(x-t{)}^{(\alpha -1)}{h}_i(t)dt,$$

$$p_{\alpha ,i}(x)=\{\begin{array}{ll}0& \text{when }x<{\beta }_1(i)\\ \frac{1}{\alpha !}(x-{\beta }_1(i){)}^{\alpha }& \text{when }{\beta }_1(i)\le x<{\beta }_2(i)\\ \frac{1}{\alpha !}[(x-{\beta }_1(i){)}^{\alpha }-2(x-{\beta }_2(i){)}^{\alpha }]& \text{when }{\beta }_2(i)\le x<{\beta }_2(i)\\ \frac{1}{\alpha !}[(x-{\beta }_1(i){)}^{\alpha }-2(x-{\beta }_2(i){)}^{\alpha }+(x-{\beta }_3(i){)}^{\alpha })]& \text{when }x>{\beta }_3(i).\end{array}$$

$$p_{1,i}(x)=\{\begin{array}{ll}(B-A)& \text{when }i=1\\ 0& \text{otherwise}.\end{array}$$

$$p_{2,i}(x)=\{\begin{array}{ll}(B-A{)}^2/2& \text{when }i=1\\ (B-A{)}^2/4m^2& \text{otherwise}.\end{array}$$

$$y_k=\frac{x_{k-1}+x_k}{2},k=1,2,...,2M.$$

$$f(x)=\sum _{i=1}^{\infty }{\alpha }_i{h}_i(x),$$

$${\alpha }_i=2^j{\int }_0^{1}f(x)h_i(x)dx.$$

$$f(x)=\sum _{i=1}^n{\alpha }_i{h}_i(x).$$

$$\frac{\partial v}{\partial t}-div(D_i(x)\nabla u_i)+f(v,w)=I_1,x\in \Omega ,t\in (0,T)$$

$$\frac{\partial v}{\partial t}-div(D_e(x)\nabla u_e)+f(v,w)=I_1,x\in \Omega ,t\in (0,T)$$

$$\frac{\partial w}{\partial t}-g(v,w)=0,x\in \Omega ,t\in (0,T)$$

$$v(x,0)=v_0(x,0),w(x,0)=w_0(x,0)x\in \Omega$$

$$n^T{D}_{i,e}(x)\nabla v=0x\in \partial \Omega ,t\in (0,T).$$

$$C_m\frac{\partial v}{\partial t}+div(D_e(x)\nabla u_e)+f(v,w)=I_2,x\in \Omega ,t\in (0,T)$$

$$-div((D_i+)D_e)(x)\nabla u_e)-div(D_i(x)\nabla v)=I_1-I_2,x\in \Omega ,t\in (0,T)$$

$$\frac{\partial w}{\partial t}-g(v,w)=0,x\in \Omega ,t\in (0,T)$$

$$v(x,0)=v_0(x,0),w(x,0)=w_0(x,0)x\in \Omega$$

$$n^T{D}_{i,e}(x)\nabla v=0x\in \partial \Omega ,t\in (0,T).$$

$$\frac{\partial v}{\partial t}-\frac{d}{dx}(D_e(x)\frac{d}{dx}{u}_e)+f(v,w)=I_2,x\in \Omega ,t\in (0,T)$$

$$-\frac{d}{dx}((D_i+)D_e)(x)\frac{d}{dx}{u}_e)-\frac{d}{dx}(D_i(x)\frac{d}{dx}v)=I_1-I_2,x\in \Omega ,t\in (0,T)$$

$$\frac{\partial w}{\partial t}-g(v,w)=0,x\in \Omega ,t\in (0,T)$$

$$v(x,0)=v_0(x,0),w(x,0)=w_0(x,0),x\in \Omega$$

$$D(x)\frac{d{u}_{i,e}}{dx}(0,t)=0,D(x)\frac{d{u}_{i,e}}{dx}(1,t)=0x\in \partial \Omega ,t\in (0,T).$$

$$\frac{{\partial }^{3}v}{\partial t\partial x^2}(x,t)=\sum _{i=1}^{2\ast M}{\alpha }_i{h}_i(x),t\in [t_s,t_{s+1}),$$

$$\frac{{\partial }^2{u}_e}{\partial x^2}(x,t)=\sum _{j=1}^{2\ast M}{\beta }_j{h}_j(x),t\in [t_s,t_{s+1}),$$

$$\frac{\partial w}{\partial t}(x,t)=\sum _{m=1}^{2\ast M}{\gamma }_m{h}_m(x),t\in [t_s,t_{s+1}).$$

$$\frac{{\partial }^{2}v}{\partial x^2}(x,t)=(t-t_s)\sum _{i=1}^{2M}{\alpha }_i{h}_i(x)+\frac{{\partial }^{2}v}{\partial x^2}(x,t_s).$$

$$\frac{\partial v}{\partial x}(x,t)=(t-t_s)\sum _{i=1}^{2M}{\alpha }_i{P}_{1,i}(x)+\frac{\partial v}{\partial x}(x,t_s)$$

$$v(x,t)=(t-t_s)\sum _{i=1}^{2M}{\alpha }_i{P}_{2,i}(x)+v(x,t_s)-v(0,t_s)+v(0,t).$$

$$\frac{\partial u_e}{\partial x}(x,t)=\sum _{j=1}^{2M}{\beta }_j{P}_{1,j}(x)+\frac{\partial u_e}{\partial x}(0,t_s),$$

$$u_e(x,t)=(t-t_s)\sum _{j=1}^{2M}{\beta }_j{P}_{2,j}(x)+u_e(0,t).$$

$$\frac{{\partial }^{2}v}{\partial t\partial x}(x,t)=(t-t_s)\sum _{i=1}^{2M}{\alpha }_i{P}_{1,i}(x)+\frac{{\partial }^{2}v}{\partial t\partial x}(0,t),$$

$$\frac{\partial v}{\partial t}(x,t)=(t-t_s)\sum _{i=1}^{2M}{\alpha }_i{P}_{2,i}(x)+\frac{\partial v}{\partial x}(0,t).$$

$$\frac{\partial v}{\partial t}(x_k,t_{s+1})-D_e(x)\frac{d^{2}v}{d{x}^2})(x_k,t_{s+1})-D_e^{\prime }(x)\frac{dv}{dx})(x_k,t_{s+1})+f(v(x_k,t_s),$$

Full text

1 Introduction

Degenerate reaction diffusion system arises in the mathematical modeling of the various fields of the natural sciences,engineering and economics, for example in gas dynamics, certain biological models, assets pricing in economics, composite media etc. The degeneracy into the model is corresponding to the interface between the two separate medium in the physical problem. This types of problems are not only important from the application point of view but equally interesting from the analysis point, since it asks for the design of techniques for the existence, uniqueness and stability of the solutions. In sedimentation processes and traffic flow problems, the concentration of the local solids is modeled by a strongly degenerate parabolic equation [1].

Mathematically modeled reaction-diffusion equations describes the variation in the concentration of one or more substances in the separated spaces with the influence of the local chemical reactions and the diffusion. This description implies that this types of systems are applied in chemistry, however, this system can also describes the dynamical processes of the biology, geology, physics and finance. Mathematically, reaction-diffusion systems take the form of semi-linear parabolic partial differential equations [2]. The system corresponding to the population dynamics of the spruce band-worm for a non-degenerate case in a biological setting is discussed in [3]. On the other hand, similar governing equations also arises in mathematical biology as a well-known reaction-diffusion system modeling the interaction between two chemical species. Under certain conditions, it produces stationary solutions with Turing-type spatial patterns [3, 5] and a standard proof for the existence and uniqueness can be found in [4]. The difficulty in direct measurement of the cardiac electric activity motivates for the mathematical modeling and numerical simulations of this phenomena. Hodking and Huxely in 1952 modeled the first mathematical model to calculate the action potential in a squid giant axon which later modified to describe the several biological phenomena. Tung [6] introduced the first mathematical model, called as Bidomain model , for the study of the cardiac electric activity. This model consists of the two degenerate parabolic reaction diffusion system corresponding to the two spaces separated by the interface membrane. This degenerate structure of the bidomain model is essentially due to the differences between the intra- and extracellular anisotropy of the cardiac tissue. Colli Franzone and Savar [7] present a weak formulation for the bidomain model and show that it has a structure suitable to apply the theory of evolution variational inequalities in Hilbert spaces. Bendahmane and Karlsen [9] prove existence and uniqueness for the bidomain equations using, for the existence part, the Faedo-Galerkin method and compactness theory, and Bourgault, Coudi‘ere, and Pierre [9] prove existence and uniqueness for the bidomain equations, first reformulating the problem into a single parabolic PDE and then applying a semigroup approach. In [] Galerkin finite element error analysis for the coupled nonlinear degenerate system of advection - diffusion equations modeling a two-phase immiscible flow through porous media is derived.

From a computational point of view, this space - time bidomain model has been numerically solved via finite difference method, finite volume method, finite element method, adaptive finite element methods using a posteriori error techniques, domain decomposition method using an alternating direction implicit method in []. In [], multiresolution technique is used to solve the degenerate system including the monodomain and bidomain models.

Recently wavelets are getting much attention for their effective use to numerically solve the mathematical models from various field of science, engineering and biology. Wavelets are well known for their inherent nature to adopt to the complexities such as discontinuities, sharp variation etc. Its properties such as orthogonality, compact support, arbitrary regularity and high order vanishing moments are very attractive. Because of these properties, solutions with discontinuities or fast oscillations in a localized region, can be approximated well by using very few wavelets. This method has been used to find the solution of the Integral equations, ordinary differential equations, partial differential equations, and fractional partial differential equations [10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. Different kinds of popular wavelets such as Daubechies [20], Chebyshev [23], Haar wavelets [24], Battle-Lemarie [21], B-spline [22], Legendre wavelets [26] are being used by the researchers for different models. Out of these wavelets Haar wavelets are very popular because of its simplicity and easy implementation in finite domains. Haar wavelets are piecewise constant functions which are orthogonal and have compact support. They also have scaling property. Because of these properties of haar wavelets, Haar wavelet method has become very popular. In this paper, we will discuss the Haar wavelets and the collocation based haar wavelet method. Because of the discontinuity of Haar wavelets, the derivatives does not exist. In this situation it is not possible to calculate the solution of differential equations. To overcome this difficulty, Chen and Hsiao [25] have proposed an idea that the highest order derivative of the differential equation is expanded into the Haar series, not the function itself. Then on integration one can obtain lower order derivatives and the functions too.

Haar wavelet method has been used to solve the linear and non-linear ordinary differential equations of all order in [28], elliptic and parabolic partial differential equations with Dirichlet and Neumann boundary conditions both in [24, 27, 29] and also for the eigenvalue problems in [30]. The basic technique of the collocation Haar wavelet is to convert the continuous problem into a discrete form with finite number of collocation points. Haar wavelet has been used vastly in the field of signal processing communication, Image processing. In [12], fluid flow boundary layer problem is solved via Haar wavelet collocation method. Haar wavelet method has also been used to solve the nonlocal problem in two-dimension. I. Singh and S. Kumar [29] proposed the Haar wavelet collocation method for the solution of three dimensional Poisson and Helmholtz equations. Haar wavelet method has been used to solve the wave-like equations by B. Naresh et.al [31]. In this work, we will develop the haar wavelet method for the coupled non-linear degenerate PDEs-ODEs system with Neumann boundary condition. We will show the advantage of haar wavelet method, like easy implementation and easily extendable to the higher dimension. In the next section, we will introduce the haar wavelet function, properties, and their integration functions. One-dimensional, two dimensional and the three dimensional haar wavelet method for the non-linear coupled non-linear degenerate PDE-PDE system with Neumann boundary which is coupled with the system of ODEs will be developed in section 3. Convergence analysis of the proposed method is conducted in section 4. Numerical results and discussion for problems in one, two and three dimensions has been discussed in the next section.

2 Haar wavelets

Let us consider the interval $x\in[A,B]$, A, B are finite real numbers. Define $M=2^{J}$, $J$ is the maximum level of resolution. This interval $[A,B]$ is equally divided into $2M$ subintervals such that the length of each subinterval is $\Delta x=(B-A)/2M$. Now, define the dilation and translation parameter $j=0,1,...,J$ and $k=0,1,...,m-1$ respectively, where, $m=2^{j}$. The wavelet number is given by $i=m+k+1$. Family of haar wavelets is defined as follows:

For $i\geq 2$

[TABLE]

where

[TABLE]

For $i=1$

[TABLE]

Haar wavelets are orthogonal, since

[TABLE]

For the solution of differential equation, we have to compute the integral

[TABLE]

where $\alpha=1,2,...,n$ and $i=1,2,...,2M$.

For the case $\alpha=0,p_{0,i}(x)=h_{i}(x)$

This integral is calculated with the help of equation (1), which is given by

[TABLE]

When $\alpha=1,2$,

[TABLE]

[TABLE]

For the grid points $x_{k}=A+kh,k=0,1,...,2M$, $h=\Delta x$, collocation points are as follows:

[TABLE]

After this discretization, we define Haar matrix $H$, and Haar Integral matrices $P_{1},P_{2}$ of size $2M\times 2M$ as $H(i,k)=h_{i}(y_{k}),P_{1}(i,k)=p_{1,i}(y_{k}),P_{2}(i,k)=p_{2,i}(y_{k})$.

2.1 Function approximation

Any function $u(x)\in L^{2}[0,1)$ can be approximated in terms of the Haar wavelet series as

[TABLE]

where the wavelet coefficients $\alpha_{i}$ are obtained by

[TABLE]

Since only the finite number of terms are taken for the computational purpose therefore, the function approximation $f$ is given by

[TABLE]

3 Mathematical Model

Consider the degenerate parabolic reaction diffusion system coupled with a system of ODEs of the type

[TABLE]

where, $v=u_{i}-u_{e}$. We can rewrite the above set of equations as follows:

[TABLE]

3.1 Haar wavelet method for the coupled degenerate reaction diffusion PDE and the ODEs

Consider the following system of coupled degenerate PDEs and ODEs in one dimension:

[TABLE]

Let us take

[TABLE]

Integrating equation (11) w.r.t. $t$ from $t_{s}$ to $t$, we obtain

[TABLE]

Integrating equation (14) w.r.t. $x$ from [math] to $x$ twice and using the boundary condition on $v$, we obtain

[TABLE]

Again, Integrating equation (12) w.r.t. $x$ from [math] to $x$ twice and using the boundary condition on $u_{e}$, we obtain

[TABLE]

Now, Integrating equation (11) w.r.t. $x$ from [math] to $x$ twice and using the boundary condition on $v$, we obtain

[TABLE]

To calculate the solution of the system (6)-(8) at the grid points we will write it in the discrete form as follows:

[TABLE]

Now, using (14),(15) in (21), we will get,

[TABLE]

Similarly, using (7), (17) in (23), we obtain,

[TABLE]

Now, from (3.1) and (3.1), we will get the following matrix system,

[TABLE]

Now, Integrate (13) w.r.t. $t$ from $t_{s}$ to $t$, we get

[TABLE]

Using (13) in equation (23), and linearize the non-linear term taking values at the previous time step, we get

[TABLE]

Matrix system of the above equation is given by,

[TABLE]

At each time step, firstly we will calculate the $w$ at the desired time by solving equation (28) and then obtain $v$ using obtained $w$, at the desired time.

3.2 Haar wavelet method for the coupled degenerate reaction diffusion PDE and the ODEs in two dimension

Consider the coupled degenerate reaction diffusion PDEs and the ODEs given as follows:

[TABLE]

Let us write $\frac{\partial^{5}v}{\partial t\partial x^{2}\partial y^{2}}(x,y,t)$, $\frac{\partial^{4}ue}{\partial x^{2}\partial y^{2}}(x,y,t)$ and $\frac{\partial w}{\partial t}(x,y,t)$ in terms of the Haar wavelet as follows :

[TABLE]

Integrating equation (37) w.r.t $t$ from $t_{s}$ to $t$, we will get

[TABLE]

Now, Integrate equation (40) twice w.r.t $x$ from [math] to $x$ also using the boundary conditions, we will obtain the following

[TABLE]

Similarly, Integrate equation (40) twice w.r.t $y$ from [math] to $y$ also using the boundary conditions, we get

[TABLE]

Now, Integrate (3.2) w.r.t $x$ from [math] to $x$ and Integrate (3.2) w.r.t $y$ from [math] to $y$ also using the boundary conditions,we get

[TABLE]

Again, Integrating (3.2) w.r.t $x$ from [math] to $x$, we obtain the following

[TABLE]

Now, Integrate (37) twice w.r.t $x$ from [math] to $x$ and apply the boundary conditions, we get

[TABLE]

Now, Integrate above equation twice w.r.t $y$ from [math] to $y$ and apply the boundary conditions, we obtain

[TABLE]

Now, Integrate equation (38) twice w.r.t $x$ from [math] to $x$ also using the boundary conditions, we will obtain the following

[TABLE]

Similarly, Integrate equation (38) twice w.r.t $y$ from [math] to $y$ also using the boundary conditions, we get

[TABLE]

Now, Integrate (50) w.r.t $x$ from [math] to $x$ and Integrate (49) w.r.t $y$ from [math] to $y$ also using the boundary conditions,we get

[TABLE]

Again, Integrating (51) w.r.t $x$ from [math] to $x$, we obtain the following

[TABLE]

Again, Integrate (39) w.r.t $t$ from $t_{s}$ to $t$, we acquire

[TABLE]

Now, to calculate the solution of the system (29)-(31) at the grid points we will write it in the discrete form as follows:

[TABLE]

Using (39) at the grid points in (57) and linearize the non-linear terms by treating it explicitly, we get the following

[TABLE]

Matrix system of the above equation is given by,

[TABLE]

where $H_{r},H_{s}$ are the Haar matrices and $c^{t}=(c_{k,l})$, which is given by,

[TABLE]

Now, at each time step we will calculate the wavelet coefficient $\gamma$ and then from (54) at the collocation points we will calculate the solution $w$. So, now we will use this $w$ to calculate the solution $v$ and $u_{e}$.

Again, Calculate equations (49), (50), (51), (52) and (3.2) at the collocation points and substitute in (3.2) and treat non-linear terms explicitly in $v$, we get the following

[TABLE]

[TABLE]

The above equation in matrix form at time $t_{s+1}$ can be written as follows:

[TABLE]

where $K=(k_{ij})$ ia matrix of size $8M^{2}\times 8M^{2}$ and $b^{t}=(b_{i})$ is a column vector of size $8M^{2}\times 1$.

Now from the above equation we will calculate the wavelet coefficients $\alpha,\beta$ and the obtain the solution $v$ with the use of calculated $w$, at the desired time step.

4 Analysis of Convergence

In this section we discuss the convergence analysis for the proposed numerical scheme.

Lemma 1

If $v(x,y)$ and $w(x,y)$ are Lipschitz continuous on domain $[0,1]^{2}$, then the wavelet coefficients $a_{i_{1},i_{2}}$, $b_{l_{1},l_{2}}$ corresponding to $v$ and $w$ satisfy the inequality

[TABLE]

Introducing the norm

[TABLE]

where $\parallel.\parallel_{2}$ is the standard $L^{2}$-norm.

Let $u=\begin{bmatrix}v&u_{e}&w\end{bmatrix}$ be the exact solution of the problem and $u_{H}=\begin{bmatrix}v_{H}&u_{e,h}&w_{H}\end{bmatrix}$ be the solution approximated by the Haar wavelets. The error $(E=u-u_{H})$ is given as follows:

Theorem 2

Let $u_{H}$ be the solution approximated by the Haar wavelet, then

[TABLE]

Proof.

[TABLE]

where

[TABLE]

Let $K_{i_{1}}=\sup_{r_{1}}K_{i_{1},r_{1}}$ and $K_{i_{2}}=\sup_{r_{2}}K_{i_{2},r_{2}}$, first term $(I_{1})$ becomes as follows

[TABLE]

Now, let $K_{1}=\sup_{i}K_{i_{1}}$ and $K_{2}=\sup_{n}K_{i_{2}}$, we get

[TABLE]

Similarly, the second term becomes as

[TABLE]

Again, we know by the definition of Haar wavelets ,

[TABLE]

[TABLE]

Therefore,

[TABLE]

5 Haar wavelet method for the coupled degenerate reaction diffusion PDEs and the ODEs in three dimension

[TABLE]

Let us write \frac{\partial^{7}v}{\partial t\partial x^{2}\partial y^{2}\partial z^{2}}\big{(}x,y,z,t\big{)}, \frac{\partial^{6}u_{e}}{\partial x^{2}\partial y^{2}\partial z^{2}}\big{(}x,y,z,t\big{)} and $\frac{\partial w}{\partial t}$ in terms of the Haar wavelet as follows:

[TABLE]

Integrating equation (72) w.r.t $t$ from $t_{s}$ to $t$, we will get

[TABLE]

Now, Integrate equation (75) twice w.r.t $x$ from [math] to $x$ also using the boundary conditions, we will obtain the following

[TABLE]

Now, Integrate equation (5) twice w.r.t $y$ from [math] to $y$ also using the boundary conditions, we get

[TABLE]

Similarly, Integrate equation (5) twice w.r.t $z$ from [math] to $z$ also using the boundary conditions, we get

[TABLE]

Again, Integrating (75) twice w.r.t $z$ from [math] to $z$ and then twice w.r.t $y$ from [math] to $y$ also using the boundary conditions, we get

[TABLE]

Now, Integrate equation (72) twice w.r.t. x , y and then z also using boundary conditions, we will obtain

[TABLE]

Now, Integrating the above equation (5) w.r.t. $t$ from $t_{s}$ to $t$, we will get

[TABLE]

Now, Integrate equation (73) twice w.r.t $x$ from [math] to $x$ also using the Neumann boundary condition on $u_{e}$, we will obtain the following

[TABLE]

Integrate the above equation twice w.r.t $x$ from [math] to $x$ also using the Neumann boundary condition on $u_{e}$, we will get

[TABLE]

Similarly, first Integrate the equation (73) twice w.r.t $z$ from [math] to $z$ and then integrate the obtain equation twice w.r.t $x$ also using the boundary conditions, we get

[TABLE]

Similarly, first Integrate the equation (73) twice w.r.t $z$ from [math] to $z$ and then integrate the obtain equation twice w.r.t $y$ from [math] to $y$ also using the boundary conditions, we get

[TABLE]

Now, integrate the above equation twice w.r.t. $x$ from [math] to $x$ also using the boundary conditions, we get

[TABLE]

Again, Integrate (74) w.r.t $t$ from $t_{s}$ to $t$, we acquire

[TABLE]

To find the solution at the collocation points, we have to discretized the equation (69) - (71) when $t\rightarrow t_{s+1}$. The discrete form is as follows:

[TABLE]

[TABLE]

[TABLE]

Using (74) at the grid points in (89) and linearize the non-linear terms by treating it explicitly, we obtain the following

[TABLE]

Matrix system of the above equation is given by,

[TABLE]

where $H_{l},H_{m},H_{n}$ are the Haar matrices and $c^{t}=(c_{k_{1},k_{2},k_{3}})$, which is given by,

[TABLE]

Now, at each time step we will calculate the wavelet coefficient $\beta$ and then from (86) at the collocation points we will calculate the solution $w$. So, now we will use this $w$ to calculate the solutions $u_{e}$ and $v$.

Again, Calculate equations (49), (50), (51), (52) and (3.2) at the collocation points and substitute in (3.2) and treat non-linear terms explicitly in $v$, we get the following

Again, Calculate equations (5), (5), (LABEL:3dvx1), (LABEL:3dvy1), (5), (83), (84), (LABEL:3duex1), (LABEL:3duey1), at the collocation points and substitute in equations (5) and (88) and treat non-linear terms explicitly in $v$, we will get the following matrix system at time $t_{s+1}$ :

[TABLE]

where $K=(k_{ij})$ ia matrix of size $8M^{3}\times 8M^{3}$ and $b^{t}=(b_{i})$ is a column vector of size $8M^{3}\times 1$.

Now from the above equation we will calculate the wavelet coefficients $\alpha,\beta$ and obtain the solutions $u_{e}$ and $v$ with the use of calculated $w$, at the desired time step.

6 Numerical Result and Discussions

We solve all the examples using above developed haar wavelet method and calculate the absolute error also. Grid validation test or resolution level test has been done for all the problems and here we are presenting for some of the problems. From grid validation we observe that resolution level $J=4$ in two dimension is good enough to calculate the solution. We use the GMRES solver to solve the linear system of equations.

Example $1$. We consider the one dimensional degenerate coupled PDEs and the ODE having homogeneous Neumann boundary as follows:

[TABLE]

where, Resolution level test for the proposed Haar wavelet method has been presented in Fig. 1.

Pointwise absolute error for different time steps is shown in Table 1 and 2. Solution at $dt=10^{-5},J=7$ is taken as the reference solution. From the Table 1 and 2. it can be seen clearly that absolute error for $v$ and $u_{e}$ decreases significantly with the smaller time step size. The Haar wavelet solution for $v$ and $u_{e}$ at the grid points are given in Fig. 2 and 3.

Example $2$.

We consider the two dimensional degenerate coupled PDEs and the ODE having homogeneous Neumann boundary as follows:

[TABLE]

First of all the grid validation of the proposed algorithm for this problem is presented in Fig. 4 which clearly shows the accuracy of the solution at the different resolution level. So, resolution level $J=4$ is good enough to calculate the results.

Error for different time steps is shown in Table 3. Solution at $dt=10^{-5},J=5$ is taken as the reference solution. From the Table 3, it can be seen clearly that error decreases seriously with the smaller time step size.

Example $3$

[TABLE]

with Neumann boundary conditions on $v$ and $u_{e}$.

where $d_{11}=1.2\times 10^{-3},d_{22}=2.5562\times 10^{-4},d_{33}=2.5562\times 10^{-4}$.

The solution of this system for $v$ and $u_{e}$ using Haar wavelet method is presented in Fig. 7.

Remark: The examples presented above are vastly applicable to the field of cardiac electrophysiology.

7 Conclusion

A Haar Wavelet Method for a class of for the coupled degenerate reaction diffusion PDEs and the ODEs having non-linear source with Neumann boundary has been proposed. The method is both simple and easy to implement in two and three dimensions. Convergence analysis has also been done to ensure the stability and accuracy. Model problems have been successfully solved. Numerical error reduces with the increase in time step size or resolution level. Problems with clinical relevance have also been successfully dealt with.