A robust framework for solving PDEs: Biorthogonal spline wavelet methods

Mohammad, Mutaz; Trounev, Alexander

doi:10.2478/ijmce-2026-0005

Full Article

1

Introduction

Biorthogonal spline wavelets have emerged as a powerful tool for numerical analysis and signal processing, offering several key advantages over other wavelet families, such as Haar, Chebyshev, and Legendre wavelets. Notably, the spline wavelets exhibit compact support, making them ideal for representing localized phenomena in the solution space, as they minimize the interaction between distant points. Furthermore, their high regularity ensures smooth approximations of functions and their derivatives, which is crucial in applications requiring accuracy in representing smooth solutions such as in fluid dynamics and heat transfer problems [1–4]. In comparison to Haar wavelets, which have discontinuous derivatives and are often less efficient in representing smooth functions, spline wavelets provide a more refined approximation, especially for higher-order derivatives. While Chebyshev and Legendre polynomials also offer smooth representations, they do not necessarily possess the compact support and duality property inherent in biorthogonal spline wavelets. The biorthogonal wavelet framework allows the construction of two families of wavelets that are dual to each other, facilitating the approximation and reconstruction of functions with higher accuracy, especially in scenarios involving boundary conditions or localized features. This duality enables better preservation of essential features such as locality and smoothness, making spline wavelets a versatile choice for numerical methods in solving differential equations and signal reconstruction. Thus, spline wavelets can significantly enhance the accuracy and computational efficiency of numerical methods, particularly in problems where smoothness and localization are critical. The connection between wavelets and the Fourier transform is noteworthy, as both frameworks analyze functions in terms of their frequency components. The Fourier transform decomposes a function into an infinite series of sinusoidal functions, providing global frequency information. However, it lacks the ability to offer localized insights into discontinuities or rapid changes in the function. In contrast, wavelets provide both time (or space) localization and frequency localization, making them particularly advantageous for handling functions with sharp gradients or localized features [2]. The construction of biorthogonal spline wavelets is typically based on B-splines, which serve as a foundational element in various applications due to their excellent approximation properties [5]. The B-spline of degree n is defined by its control points and a corresponding knot vector, which dictates the intervals over which the spline is defined. The versatility of B-splines is utilized to create wavelets that are not only efficient for representing functions but also adaptable to different smoothness requirements [6].

To construct biorthogonal spline wavelets, one starts by defining a scaling function associated with a given B-spline. The scaling function, often denoted as ϕ(x), acts as a low-pass filter that facilitates the decomposition of functions into different frequency components, analogous to the low-pass characteristics of the Fourier transform. The wavelet function, denoted as ψ(x), serves as a high-pass filter, capturing the detail information necessary for accurate representation [7]. The relationship between the scaling and wavelet functions can be formalized as follows 1 $ϕ (x) = \sum_{k} h_{k} ψ (2 x - k),$ \phi (x) = \sum\limits_k {{h_k}} \psi (2x - k), where h_k are the coefficients that define the scaling operation. The corresponding wavelet function defined by the above equation generates an orthonormal basis ψ_jk(t) for L²(ℝ), represented as 2 $ψ_{j k} (t) = 2^{j / 2} ψ (2^{j} t - k) .$ {\psi _{jk}}(t) = {2^{j/2}}\psi ({2^j}t - k).

This orthonormal basis facilitates an efficient representation and reconstruction of functions in the L² space. In particular, the scaling function ϕ of order 4 plays a crucial role in the wavelet framework, offering the degree of smoothness necessary for high-quality approximations. The associated dual scaling function, of order 2, ensures that the combined system of scaling and wavelet functions maintain the desirable biorthogonality properties. The interplay between these scaling functions enhances the capabilities of wavelet transforms, allowing for improved analysis and processing of various functions, particularly in numerical methods for PDEs.

In the context of solving PDEs, biorthogonal spline wavelets offer several advantages. They facilitate adaptive mesh refinement, allowing for localized increases in resolution where needed, which is particularly beneficial for capturing sharp gradients or discontinuities in the solution. Moreover, the multiscale nature of wavelet transforms enables the extraction of features at different scales, thus improving both accuracy and efficiency in the numerical solution of PDEs [8]. The application of biorthogonal spline wavelets in solving PDEs has been demonstrated in various studies. For instance, wavelet-based methods have been employed for solving elliptic equations, with significant improvements in convergence rates compared to traditional methods [3]. Furthermore, their capability to handle irregular geometries and boundary conditions makes them suitable for complex physical problems [4]. These methods have also been effectively applied in solving incompressible fluid dynamics problems, from Navier-Stokes to Poisson equations [9]. The versatility of biorthogonal spline wavelets extends beyond solving PDEs, with applications in image processing, where their compact support and regularity allow for efficient algorithms in compression and denoising [10]. Additionally, their use in statistical analysis has garnered attention, as wavelets provide localized information for non-stationary processes [11]. Recent advancements have expanded the application of wavelet-based methods to time-fractional PDEs, where Euler wavelets and novel numerical techniques have proven effective in improving computational precision [12]. Similarly, wavelet-based approaches have been applied to the fractional Navier-Stokes equations [13] and neutral delay differential equations [14]. In conclusion, biorthogonal spline wavelets represent a versatile and efficient approach to the numerical solution of PDEs, combining the benefits of spline approximation with the powerful tools provided by wavelet analysis. As research in this area continues to evolve, further developments are expected to enhance the capabilities of wavelet methods in tackling increasingly complex differential equations. Figure 1 shows the scaling function and its corresponding wavelet of order 4 and dual order 2.

In this study, we apply the biorthogonal spline wavelet framework to solve the Diffusion equation, numerically. To demonstrate the utility of this approach, we present a detailed example of the numerical solution to the Diffusion equation using the wavelet collocation method. To fully assess the effectiveness of this method, we include the comparisons with other established techniques from the literature, such as finite difference and finite element methods. This comparison highlights the advantages of wavelet-based methods, including higher accuracy and computational efficiency [3, 4].

This paper is organized as follows: Section 2 presents the PDEs under consideration, including the Diffusion equations. Section 3 introduces the simulation algorithm, which utilizes the biorthogonal spline wavelet (BSW) technique for numerical solutions. Section 4 discusses the numerical results and their accuracy, followed by a conclusion in Section 5, which highlights the potential of the proposed wavelet-based method for solving complex PDEs in various applications.

2

The PDEs problem

In this section, we introduce the PDEs that will be solved numerically by using BSW method. In particular, we consider the Diffusion equation in the general form given by 3 $\frac{\partial u (x, t)}{\partial t} = D \nabla^{2} u (x, t),$ {{\partial u(x,t)} \over {\partial t}} = D{\nabla ^2}u(x,t), where u(x,t) is the unknown function representing the quantity being diffused (e.g., temperature, concentration), D is the diffusion coefficient, and ∇² is the Laplace operator (spatial second derivatives). For a one-dimensional case, the diffusion equation simplifies to 4 $\frac{\partial u (x, t)}{\partial t} = D \frac{\partial^{2} u (x, t)}{\partial x^{2}} .$ {{\partial u(x,t)} \over {\partial t}} = D{{{\partial ^2}u(x,t)} \over {\partial {x^2}}}.

3

Simulation algorithm

We apply the numerical algorithm projected to solve the model studied. The method involves discretizing the equations using specific desensitization, and BSW technique through a robust algorithm. Let V₁(x) and V₂(x) represent the biorthogonal spline scaling functions of orders 2 and 4, respectively. These functions is utilized in the numerical scheme defined by the following integral relationships $\begin{matrix} J_{1}^{1} (x) = \int_{0}^{x} V_{1} (t) d t, & J_{2}^{1} (x) = \int_{0}^{x} V_{2} (t) d t, \\ J_{1}^{2} (x) = \int_{0}^{x} J_{1}^{1} (t) d t, & J_{2}^{2} (x) = \int_{0}^{x} J_{2}^{1} (t) d t . \end{matrix}$ \matrix{ {J_1^1(x) = \int_0^x {{V_1}} (t)dt,} & {J_2^1(x) = \int_0^x {{V_2}} (t)dt,} \cr {J_1^2(x) = \int_0^x {J_1^1} (t)dt,} & {J_2^2(x) = \int_0^x {J_2^1} (t)dt.} \cr }

For example, the elements of integral matrices corresponding to these integrals are given as follows $\begin{array}{l} V_{1} (x) & = & \frac{1}{2} {\begin{array}{l} {(- 1 + x)}^{3} & for x < 0 \\ 0 & for x = 0 \\ {(- 1 + x)}^{3} - 4 {(- \frac{1}{2} + x)}^{3} + 6 x^{3} - 4 {(\frac{1}{2} + x)}^{3} + {(1 + x)}^{3} & for x > 0 \end{array} \\ V_{2} (x) & = & \frac{1}{6} {\begin{array}{l} {(- 1 + x)}^{3} & for x < 0 \\ 0 & for x = 0 \\ {(- 1 + x)}^{3} - 9 x^{3} + 2 {(1 + 2 x)}^{3} - 9 {(1 + x)}^{3} + {(2 + x)}^{3} & for x > 0. \end{array} \end{array}$ \matrix{ {{V_1}(x)} \hfill & = \hfill & {{1 \over 2}\left\{ {\matrix{ {{{( - 1 + x)}^3}} \hfill & {{\rm{for }}x < 0} \hfill \cr 0 \hfill & {{\rm{for }}x = 0} \hfill \cr {{{( - 1 + x)}^3} - 4{{\left( { - {1 \over 2} + x} \right)}^3} + 6{x^3} - 4{{\left( {{1 \over 2} + x} \right)}^3} + {{(1 + x)}^3}} \hfill & {{\rm{for }}x > 0} \hfill \cr } } \right.} \hfill \cr {{V_2}(x)} \hfill & = \hfill & {{1 \over 6}\left\{ {\matrix{ {{{( - 1 + x)}^3}} \hfill & {{\rm{for }}x < 0} \hfill \cr 0 \hfill & {{\rm{for }}x = 0} \hfill \cr {{{( - 1 + x)}^3} - 9{x^3} + 2{{(1 + 2x)}^3} - 9{{(1 + x)}^3} + {{(2 + x)}^3}} \hfill & {{\rm{for }}x > 0.} \hfill \cr } } \right.} \hfill \cr }

The first integrals of V₁(x) and V₂(x) are given by $\begin{array}{l} J_{1}^{1} (x) & = & \int_{0}^{x} V_{1} (t) d t, \\ J_{2}^{1} (x) & = & \int_{0}^{x} V_{2} (t) d t . \end{array}$ \matrix{ {J_1^1(x)} \hfill & = \hfill & {\int_0^x {{V_1}} (t)dt,} \hfill \cr {J_2^1(x)} \hfill & = \hfill & {\int_0^x {{V_2}} (t)dt.} \hfill \cr }

The second integrals of $J_{1}^{1} (x)$ (J_1^1(x) and $J_{2}^{1} (x)$ J_2^1(x) are $\begin{array}{l} J_{1}^{2} (x) & = & \int_{0}^{x} J_{1}^{1} (t) d t, \\ J_{2}^{2} (x) & = & \int_{0}^{x} J_{2}^{1} (t) d t . \end{array}$ \matrix{ {J_1^2(x)} \hfill & = \hfill & {\int_0^x {J_1^1} (t)dt} \hfill \cr {J_2^2(x)} \hfill & = \hfill & {\int_0^x {J_2^1} (t)dt} \hfill \cr }

Let Ω be the set of functions described by these integral forms. For any τ ∈ Ω, a new function ψ(x) is defined as $ψ (t) = τ I_{[0, 1]} (t),$ \psi (t) = \tau {_{[0,1]}}(t), where 𝕀_[0,1] is the indicator function on the interval [0, 1]. Given the following mappings $ψ_{1} = V_{1}, ψ_{2} = V_{2}, ψ_{1, 1} = J_{1}^{1}, ψ_{2, 1} = J_{2}^{1}, ψ_{1, 2} = J_{1}^{2}, ψ_{2, 2} = J_{2}^{2},$ {\psi _1} = {V_1},\quad {\psi _2} = {V_2},\quad {\psi _{1,1}} = J_1^1,\quad {\psi _{2,1}} = J_2^1,\quad {\psi _{1,2}} = J_1^2,\quad {\psi _{2,2}} = J_2^2, we define a family of wavelet functions indexed by j and k ∈ ℤ as follows $\begin{matrix} ψ_{1} (j, k, t) = ψ_{1} (2^{j} t - k), ψ_{2} (j, k, t) = ψ_{2} (2^{j} t - k), \\ ψ (j, k, t) = ψ_{1} (j, k, t) + ψ_{2} (j, k, t), \\ ψ^{1, 1} (j, k, t) = ψ_{1, 1} (2^{j} t - k), ψ^{1, 2} (j, k, t) = ψ_{1, 2} (2^{j} t - k), \\ ψ^{2, 1} (j, k, t) = ψ_{2, 1} (2^{j} t - k), ψ^{2, 2} (j, k, t) = ψ_{2, 2} (2^{j} t - k), \\ ψ^{1} (j, k, t) = \frac{ψ^{1, 1} (j, k, t) + ψ^{2, 1} (j, k, t)}{j}, ψ^{2} (j, k, t) = \frac{ψ^{2, 1} (j, k, t) + ψ^{2, 2} (j, k, t)}{j^{2}} . \end{matrix}$ \matrix{ {{\psi _1}(j,k,t) = {\psi _1}\left( {{2^j}t - k} \right),\quad {\psi _2}(j,k,t) = {\psi _2}\left( {{2^j}t - k} \right),} \cr {\psi (j,k,t) = {\psi _1}(j,k,t) + {\psi _2}(j,k,t),} \cr {{\psi ^{1,1}}(j,k,t) = {\psi _{1,1}}\left( {{2^j}t - k} \right),\quad {\psi ^{1,2}}(j,k,t) = {\psi _{1,2}}\left( {{2^j}t - k} \right),} \cr {{\psi ^{2,1}}(j,k,t) = {\psi _{2,1}}\left( {{2^j}t - k} \right),\quad {\psi ^{2,2}}(j,k,t) = {\psi _{2,2}}\left( {{2^j}t - k} \right),} \cr {{\psi ^1}(j,k,t) = {{{\psi ^{1,1}}(j,k,t) + {\psi ^{2,1}}(j,k,t)} \over j},\quad {\psi ^2}(j,k,t) = {{{\psi ^{2,1}}(j,k,t) + {\psi ^{2,2}}(j,k,t)} \over {{j^2}}}.} \cr }

Let us consider a function ψ ∈ L₂(ℝ) that can be expressed by using the following series expansion $ψ (x) = \sum_{ℓ = 1}^{2} \sum_{j, k \in ℤ} d^{ℓ} (j, k) ψ^{ℓ} (j, k, x),$ \psi (x) = \sum\limits_{\ell = 1}^2 {\sum\limits_{j,k \in } {{d^\ell }} } (j,k){\psi ^\ell }(j,k,x), where the coefficients d^ℓ(j, k) are given by the inner product $d^{ℓ} (j, k) = 〈 ψ, ψ^{ℓ} (j, k, x) 〉 = \int_{ℝ} ψ (x) ψ^{ℓ} (j, k, x) w (x) d x,$ {d^\ell }(j,k) = \left\langle {\psi, {\psi ^\ell }(j,k,x)} \right\rangle = \int_ \psi (x){\psi ^\ell }(j,k,x)w(x){\mkern 1mu} dx, with 〈·, ·〉 representing the standard inner product in L₂(ℝ), and w(x) being a suitable weight function.

To approximate ψ, we truncate the series as follows $ψ_{n, M} (t) = \sum_{ℓ = 1}^{2} \sum_{j = 0}^{n} \sum_{k = 0}^{M - 1} d^{ℓ} (j, k) ψ^{ℓ} (j, k, t) .$ {\psi _{n,M}}(t) = \sum\limits_{\ell = 1}^2 {\sum\limits_{j = 0}^n {\sum\limits_{k = 0}^{M - 1} {{d^\ell }} } } (j,k){\psi ^\ell }(j,k,t).

The error between the original function ψ and the truncated approximation ψ_n,M is given by ${‖ ψ - ψ_{n, M} ‖}_{2}^{2} = {‖ \sum_{ℓ = 1}^{2} \sum_{j \geq n + 1} \sum_{k \geq M + 1} d^{ℓ} (j, k) ψ^{ℓ} (j, k, t) ‖}_{2}^{2} .$ \left\| {\psi - {\psi _{n,M}}} \right\|_2^2 = \left\| {\sum\limits_{\ell = 1}^2 {\sum\limits_{j \ge n + 1} {\sum\limits_{k \ge M + 1} {{d^\ell }} } } (j,k){\psi ^\ell }(j,k,t)} \right\|_2^2.

Furthermore, we observe that $| 〈 d^{ℓ} (j, k) 〉 | \leq \max_{_{t}} | ψ | {‖ ψ^{ℓ} (j, k, t) ‖}_{1} .$ \left| {\left\langle {{d^\ell }(j,k)} \right\rangle } \right| \le \mathop {\max }\limits_{_t} |\psi |{\left\| {{\psi ^\ell }(j,k,t)} \right\|_1}.

Using the Bessel inequality, the error bound can be derived as ${‖ ψ - ψ_{n, M} ‖}_{2}^{2} \leq \max_{t} | ψ | {‖ ψ^{ℓ} (j, k, t) ‖}_{1} \sum_{ℓ = 1}^{2} \sum_{j \geq n + 1} \sum_{k = M} | d^{ℓ} (j, k) | .$ \left\| {\psi - {\psi _{n,M}}} \right\|_2^2 \le {\max _t}|\psi |{\left\| {{\psi ^\ell }(j,k,t)} \right\|_1}\sum\limits_{\ell = 1}^2 {\sum\limits_{j \ge n + 1} {\sum\limits_{k = M} {\left| {{d^\ell }(j,k)} \right|} } } .

To solve the problem, we construct a vector Ξ_ψ of length M = 2ⁿ⁺¹, n ∈ ℕ. The vector is defined as $Ξ_{ψ} = [ψ_{τ}, σ^{κ} (1, 0, x), \dots, σ^{κ} (2^{j}, k, x), \dots, σ^{κ} (2^{n}, 2^{n - 1}, x)], j = 0, 1, 2, \dots, n; k = 0, 1, 2, \dots, 2^{j - 1},$ {\Xi _\psi } = [{\psi _\tau },{\sigma ^\kappa }(1,0,x), \cdots, {\sigma ^\kappa }({2^j},k,x), \cdots, {\sigma ^\kappa }({2^n},{2^{n - 1}},x)],\quad j = 0,1,2, \cdots, n;\;k = 0,1,2, \cdots, {2^{j - 1}}, where the components are assigned as follows ${\begin{array}{l} ψ_{τ} = 1, σ^{κ} = ψ & when τ = V_{1}, V_{2} and κ = 1, \\ ψ_{τ} = x, σ^{κ} = ψ^{1} & when τ = J_{1}^{1}, J_{2}^{1} and κ = j, \\ ψ_{τ} = \frac{x^{2}}{2}, σ^{κ} = ψ^{2} & when τ = J_{1}^{2}, J_{2}^{2} and κ = j^{2} . \end{array}$ \left\{ {\matrix{ {{\psi _\tau } = 1,{\sigma ^\kappa } = \psi } \hfill & {{\rm{when }}\tau = {V_1},{V_2}{\rm{ and }}\kappa = 1{\rm{,}}} \hfill \cr {{\psi _\tau } = x,{\sigma ^\kappa } = {\psi ^1}} \hfill & {{\rm{when }}\tau = J_1^1,J_2^1{\rm{ and }}\kappa = j,} \hfill \cr {{\psi _\tau } = {{{x^2}} \over 2},{\sigma ^\kappa } = {\psi ^2}} \hfill & {{\rm{when }}\tau = J_1^2,J_2^2{\rm{ and }}\kappa = {j^2}.} \hfill \cr } } \right.

These integrals are crucial components in our approach to solve the given equation. By integrating equation (3) step by step, we derive the following 5 $u_{x x t} = Ψ_{k M}^{T} (x) \cdot V \cdot Ψ_{k M} (t),$ {u_{xxt}} = \Psi _{kM}^T(x)\cdotV\cdot{\Psi _{kM}}(t), where V is the coefficient matrix. By further integration, we obtain expressions for the derivatives and the solution as follows 6 $\begin{array}{l} u_{x x} (x, t) & = & Ψ_{k M}^{T} (x) V I_{1} (t) + G_{1}^{''} (x), \\ u_{x} (x, t) & = & I_{1}^{T} (x) V I_{1} (t) + G_{1}^{'} (x) - G_{1}^{'} (0) + F_{2} (t), \\ u (x, t) & = & I_{2}^{T} (x) V I_{1} (t) + G_{1} (x) - G_{1} (0) - x G_{1}^{'} (0) + x G_{2} (t) + G_{3} (t) . \end{array}$ \matrix{ {{u_{xx}}(x,t)} \hfill & = \hfill & {\Psi _{kM}^T(x)V{I_1}(t) + G_1^{\prime \prime }(x),} \hfill \cr {{u_x}(x,t)} \hfill & = \hfill & {I_1^T(x)V{I_1}(t) + G_1^\prime (x) - G_1^\prime (0) + {F_2}(t),} \hfill \cr {u(x,t)} \hfill & = \hfill & {I_2^T(x)V{I_1}(t) + {G_1}(x) - {G_1}(0) - xG_1^\prime (0) + x{G_2}(t) + {G_3}(t).} \hfill \cr }

Here, V is an N × N matrix derived from collocation points, and the functions G₁(x), G₂(t), and G₃(t) represent arbitrary functions of x and t. The terms I₁(t) and I₂(x) denote integration vectors, which are defined over the time and spatial domains, respectively. Finally, the discretizing of the equation at N collocation points and substituting the necessary functions yields a system of linear equations that must be solved to obtain the numerical solution.

4

Illustration to the problem

In this section, we present a numerical example that illustrates the application of our proposed method to the Diffusion equation. This example is designed to demonstrate the effectiveness and accuracy of the numerical solutions obtained through the discretization process at collocation points. We aim to provide insights into the behavior of the solutions and the impact of different parameters on the results.

Let’s consider the Diffusion equation defined on the interval [0,L] with the following initial and boundary conditions $\begin{array}{l} u_{t} (x, t) = D u_{x x} (x, t), 0 \leq x \leq L, \\ u (x, 0) = f (x), 0 \leq x \leq L, \\ u (0, t) = g_{1} (t), t > 0. \\ u (L, t) = g_{2} (t), t \geq 0. \end{array}$ \matrix{ {{u_t}(x,t) = D{u_{xx}}(x,t),\quad 0 \le x \le L,} \hfill \cr {u(x,0) = f(x),\quad 0 \le x \le L,} \hfill \cr {u(0,t) = {g_1}(t),\quad t > 0.} \hfill \cr {u(L,t) = {g_2}(t),\quad t \ge 0.} \hfill \cr }

For illustration, let’s consider as L = 1, f(x) = –x² + x, g₁(t) = 0, and g₂(t) = 0. We apply our numerical method using N collocation points within the interval [0, 1]. We start by discretizing the governing equation. Using the BSW approach, we represent the function u(x,t) as a combination of wavelet basis functions. The first step is to express the solution u(x,t) as $u (x, t) = \sum_{ℓ = 1}^{2} \sum_{j, k \in ℤ} d^{ℓ} (j, k) ψ^{ℓ} (j, k, x, t),$ u(x,t) = \sum\limits_{\ell = 1}^2 {\sum\limits_{j,k \in } {{d^\ell }} } (j,k){\psi ^\ell }(j,k,x,t), where ψ^𝓁(j,k,x,t) are the wavelet basis functions. In the numerical method, we employ spline wavelet-based collocation points distributed within the interval [0, 1]. These points are designed to align with the properties of the wavelet basis functions, ensuring an optimal representation of the solution.

Let $Δ x = \frac{1}{N}$ \Delta x = {1 \over N} represent the grid spacing, where N denotes the number of collocation points. The collocation points x_l are then defined as $x_{l} = l \cdot Δ x, l = 0, 1, 2, \dots, N .$ {x_l} = l\cdot\Delta x,\quad l = 0,1,2, \cdots, N.

Here, x_l denotes the location of the l-th collocation point within the interval [0, 1], and N is the total number of collocation points. The set of collocation points is given by $x_{0}, x_{1}, x_{2}, \dots, x_{N} .$ {x_0},{x_1},{x_2}, \cdots, {x_N}.

This grid of collocation points is used for the discretization of the governing equation, facilitating the application of the wavelet-based collocation method.

We now substitute the discretized version of u(x,t) into the Diffusion equation. This results in a system of linear equations, which is derived by applying the wavelet collocation method at N collocation points. $u_{x x t} = Ψ_{k M}^{T} (x) \cdot V \cdot Ψ_{k M} (t),$ {u_{xxt}} = \Psi _{kM}^T(x)\cdotV\cdot{\Psi _{kM}}(t), where V is the coefficient matrix that is formed by the discretization.

At each time step, the system of linear equations is solved to obtain the numerical approximation of u(x,t). For example, at t = 1, we obtain the following expressions for the derivatives $\begin{matrix} u_{x x} (x, t) = Ψ_{k M}^{T} (x) V I_{1} (t) + G_{1}^{''} (x), \\ u_{x} (x, t) = I_{1}^{T} (x) V I_{1} (t) + G_{1}^{'} (x) - G_{1}^{'} (0) + F_{2} (t), \\ u (x, t) = I_{2}^{T} (x) V I_{1} (t) + G_{1} (x) - G_{1} (0) - x G_{1}^{'} (0) + x G_{2} (t) + G_{3} (t) . \end{matrix}$ \matrix{ {{u_{xx}}(x,t) = \Psi _{kM}^T(x)V{I_1}(t) + G_1^{\prime \prime }(x),} \cr {{u_x}(x,t) = I_1^T(x)V{I_1}(t) + G_1^\prime (x) - G_1^\prime (0) + {F_2}(t),} \cr {u(x,t) = I_2^T(x)V{I_1}(t) + {G_1}(x) - {G_1}(0) - xG_1^\prime (0) + x{G_2}(t) + {G_3}(t).} \cr }

The system is solved at each time step to calculate the numerical approximation of u(x,t). The analytical solution to the Diffusion equation, for comparison, is given by $u (x, t) = \sum_{K = 1}^{\infty} - \frac{4 ({(- 1)}^{K} - 1) e^{- \frac{π^{2} t K^{2}}{5}} \sin (π K x)}{π^{3} K^{3}} .$ u(x,t) = \sum\limits_{K = 1}^\infty - {{4\left( {{{( - 1)}^K} - 1} \right){e^{ - {{{\pi ^2}t{K^2}} \over 5}}}\sin (\pi Kx)} \over {{\pi ^3}{K^3}}}.

Figures 2, 3, 4, and 5 illustrate various graphical results from the numerical simulation for different partial sums. Table 1 presents a numerical comparison of the approximate solution obtained using the biorthogonal spline wavelet of order 2 (BSW-2) and Haar wavelets. The error used in this simulation is the L_∞ error.

Table 1

Numerical values for x_l,t and corresponding u[x_l,t] for two methods.

t=0	t=0	t=0	t=0.3	t=0.3	t=0.6	t=0.6	t=0.9	t=0.9
x_l	BSW-2	Haar	BSW-2	Haar	BSW-2	Haar	BSW-2	Haar
$\frac{1}{16}$ {1 \over {16}}	0	0	0.1867561	0.1850020	0.2158150	0.2140760	0.0752277	0.0745150
$\frac{3}{16}$ {3 \over {16}}	0	0	0.1451483	0.1440110	0.1705720	0.1695220	0.0548761	0.0542550
$\frac{5}{16}$ {5 \over {16}}	0	0	0.1130518	0.1120270	0.1330290	0.1320680	0.0434034	0.0430950
$\frac{7}{16}$ {7 \over {16}}	0	0	0.0883459	0.0875280	0.1037956	0.1030230	0.0336308	0.0332450
$\frac{9}{16}$ {9 \over {16}}	0	0	0.0688853	0.0682610	0.0810188	0.0803560	0.0263730	0.0261660
$\frac{11}{16}$ {{11} \over {16}}	0	0	0.0537768	0.0532920	0.0632059	0.0627320	0.0205080	0.0203022
$\frac{13}{16}$ {{13} \over {16}}	0	0	0.0419592	0.0415060	0.0493254	0.0488210	0.0160445	0.0159001
$\frac{15}{16}$ {{15} \over {16}}	0	0	0.0327424	0.0324060	0.0384910	0.0381090	0.0124974	0.0124034

In order to evaluate the accuracy and effectiveness of the proposed BSW method, we compare error norms based on the numerical results presented in Table 1. This table provides the values of u[x_l,t] at various spatial points x_l and time steps t, computed using the proposed method. By analyzing these values, we assess the method’s accuracy in approximating the solution over time.

To provide context, we also compare these results with those obtained from other commonly used wavelet families, such as Haar and Symlet wavelets. While these methods are widely utilized, each exhibits unique characteristics that lead to varying levels of accuracy depending on the problem being solved.

From Table 1, we observe that the proposed method performs well across different time steps, effectively capturing the temporal evolution of the solution. In comparison to the Haar wavelet and other methods, the BSW method offers superior accuracy, particularly at higher-order time steps. While Haar and other methods yield reasonable results, the BSW method demonstrates better error convergence and stability, especially for finer temporal resolutions, making it a more precise choice for similar problems.

To further assess convergence behavior, we present the convergence rate graphs in Figure 6. The left panel illustrates the convergence rate for the Haar wavelet method, while the right panel shows the BSW-2 method’s convergence. These graphs highlight the error norm reduction over multiple time steps and provide a clear comparison of the two methods.

Additionally, the log-log plots in Figure 7 offer further insight into the numerical performance of both methods. The left plot presents the performance of the Haar wavelet, while the right plot shows the results for the BSW-2 method. These plots visually confirm that the BSW-2 method achieves a superior rate of convergence, with error norms decreasing more rapidly compared to Haar, particularly at finer time resolutions.

In conclusion, the numerical results, convergence rate graphs, and log-log plots all emphasize the advantages of the BSW method. It consistently outperforms Haar wavelets in both accuracy and computational efficiency, especially over multiple time intervals. This makes it a more reliable and precise tool for solving time-dependent problems, surpassing the capabilities of other wavelet families.

5

Conclusion

In this study, we have successfully developed and implemented a numerical algorithm to tackle the challenges of solving the Diffusion equation. By incorporating advanced techniques, such as collocation methods and biorthogonal spline wavelets, we have achieved significant improvements in both solution accuracy and computational efficiency. The results of our numerical experiments validate the proposed method, demonstrating its close alignment with known analytical solutions and its ability to effectively handle complex boundary conditions and initial states.

Looking ahead, future work will focus on extending this framework to address a wider range of PDEs, particularly those with nonlinear characteristics. Furthermore, we plan to investigate adaptive methods that dynamically adjust collocation points based on the solution’s behavior, further enhancing precision while optimizing computational resources.

A robust framework for solving PDEs: Biorthogonal spline wavelet methods

Full Article

Paradigm

My account