#  Option Pricing under Multifactor Black–Scholes Model Using Orthogonal
Spline Wavelets

###### Abstract

The paper focuses on pricing European-style options on several underlying
assets under the Black–Scholes model represented by a nonstationary partial
differential equation. The proposed method combines the Galerkin method with
$L^{2}$-orthogonal sparse grid spline wavelets and the Crank–Nicolson scheme
with Rannacher time-stepping. To this end, we construct an orthogonal cubic
spline wavelet basis on the interval satisfying homogeneous Dirichlet boundary
conditions and design a wavelet basis on the unit cube using the sparse tensor
product. The method brings the following advantages. First, the number of
basis functions is significantly smaller than for the full grid, which makes
it possible to overcome the so-called “curse of dimensionality.” Second, some
matrices involved in the computation are identity matrices, which
significantly simplifies and streamlines the algorithm, especially in higher
dimensions. Further, we prove that discretization matrices have uniformly
bounded condition numbers, even without preconditioning, and that the
condition numbers do not depend on the dimension of the problem. Due to the
use of cubic spline wavelets, the method is higher-order convergent. Numerical
experiments are presented for options on the geometric average.

Option Pricing under Multifactor Black–Scholes Model Using Orthogonal Spline
Wavelets

Dana Černá111Department of Mathematics and Didactics of Mathematics, Technical
University of Liberec, Studentská 2, Liberec, 461 17, Czech Republic, Kateřina
Fiňková222Institute of Information Technology and Electronics, Technical
University of Liberec, Studentská 2, Liberec, 461 17, Czech Republic

Keywords: Black–Scholes model, European option, orthogonal spline wavelet,
sparse grid, wavelet-Galerkin method, condition number

##  1 Introduction

Methodology for the valuation of options has undergone rapid development in
recent years; nowadays, various option pricing models and numerical methods
are available to determine the fair option price. This paper focuses on the
famous Black–Scholes model [3], which assumes that prices of underlying assets
follow a log-normal distribution with constant drift and volatility. The
Black–Scholes model is used for various types of options such as vanilla
options, basket options, options on maximum or minimum of several assets,
options on the geometric average, and real options.

Because the Black–Scholes model is a basic model for option pricing, the
numerical valuation of options under this model is an active field of
research, and there are a great many articles concerning this topic, such as
the finite difference method studied in [28], higher-order spline methods
designed in [25, 33, 37], and radial basis function approach examined in [32].
However, it is well known that when the Black–Scholes model is generalized to
higher dimensions, e.g., dimension $d\geq 3$, the so-called “curse of
dimensionality” causes the number of degrees of freedom to grow exponentially
with the number of dimensions. Therefore, the numerical solution of the model
is more complicated, and far less literature has been devoted to this problem,
see, e.g., [1, 23, 29, 27, 30]. Wavelets have also been used for option
pricing under the Black–Scholes model with one or more underlying assets, such
as the B-spline wavelet-Galerkin method combined with discontinuous Galerkin
time discretization [22], an adaptive wavelet method [36, 26], the Legendre
wavelet method [16], a method based on the wavelet approximation and the
characteristic function [17], and the Shannon wavelet inverse Fourier
technique [18].

In this paper, we deal with the Black–Scholes model for pricing multi-asset
European options. For $d$ asset option, this model is represented by a
nonstationary differential equation with $d$ state variables. First, we adjust
the equation using the transformation into logarithmic prices, drift removal,
approximation of the unbounded domain, and setting appropriate boundary
conditions.

Because the properties of Galerkin-type methods strongly depend on the wavelet
basis used, we focus on constructing an appropriate one, namely an
$L^{2}$-orthogonal cubic spline wavelet basis. In [14], it was proved that if
$\tilde{V}_{j}$ is a multiresolution analysis for the space
$L^{2}\left(\mathbb{R}\right)$, then there exists an orthogonal
multiresolution analysis $V_{j}$ for the space $L^{2}\left(\mathbb{R}\right)$
such that

$\tilde{V}_{q}\subset V_{0}\subset\tilde{V}_{q+n}$ (1)

for some $q\geq 0$ and $n\geq 1$, and construction principles of such a
multiresolution analysis were introduced. Based on [14], orthogonal wavelets
were constructed in [13, 39]. However, the general method from [14] does not
lead to uniquely determined scaling functions, because it leads to the system
of quadratic equations with more than one solution and because the part of the
method is the orthogonalization of particular functions, which also does not
lead to unique functions. In [15], the construction of $L^{2}$-orthogonal
spline wavelets corresponding to the $L^{2}$-orthogonal spline functions from
[14] was proposed. However, the resulting wavelet bases are not all equally
suitable for use in numerical mathematics. While the $L^{2}$ condition number
is one due to the $L^{2}$ orthogonality, the $H^{1}$ condition numbers differ
for various wavelet bases. However, the $H^{1}$ condition number is also
important because it influences the condition number of discretization
matrices and the number of iterations needed to resolve the differential
equation with the required accuracy. Furthermore, to be able to use the basis
for the solution of the differential equation on a bounded domain, it is
necessary to adapt the basis of $L^{2}\left(\mathbb{R}\right)$ to a bounded
interval.

In this paper, we present a construction of orthogonal cubic scaling functions
and wavelets on the real line. Then, we construct special boundary scaling
functions and wavelets to adapt the basis from the real line to the interval.
Finally, we use the sparse tensor product to create a wavelet basis. We show
that when normalized in the $H^{1}$-seminorm, the resulting basis is well
conditioned with respect to the $H^{1}$-seminorm, and study the properties of
the basis. The coefficients of the constructed scaling functions and wavelets
are presented in Appendix A.

Then, we propose and study a wavelet-based method for the multidimensional
Black–Scholes model. The method combines the Galerkin method using sparse
tensor product $L^{2}$-orthogonal spline wavelet bases and the Crank-Nicolson
with Rannacher time-stepping. Due to the sparse tensor product structure and
the orthogonality of the bases, the resulting scheme is greatly simplified
compared to the standard Galerkin method as well as many other methods; it
even overcomes the curse of dimensionality. The simplification is attributable
to some of the matrices involved in computations being identity matrices. It
is known that wavelet-based methods for certain types of equations lead to
uniformly conditioned matrices when diagonal preconditioning is used [11]. In
this paper, we prove that the matrices are uniformly conditioned, even without
preconditioning, and that the condition numbers can be bounded independently
on the dimension of the problem. Because the matrices are symmetrical and
positive definite, the resulting system can be solved using the conjugate
gradient method. The uniform boundedness of condition numbers implies that the
numbers of iterations are also bounded independently of the size of the matrix
and the dimension of the problem. Another advantage is high-order convergence
of the proposed method with constructed cubic spline-wavelet basis.

The paper is organized as follows. Section 2 introduces the model and
adjustment of the equation. Section 3 presents the assumption on the wavelet
basis, construction of orthogonal cubic spline wavelet basis on the real line,
its adaptation to a bounded interval and boundary conditions, and construction
of wavelet basis on the unit cube using the sparse tensor product. The
wavelet-based method is proposed and analyzed in Section 4. Numerical results
are provided for options on the geometric average in Section 5.

##  2  Multifactor Black–Scholes model

The multi-dimensional Black–Scholes model is a generalization of the original
Black–Scholes model for one-asset options designed in [3]. For an option on
$d$ assets, the Black–Scholes model assumes that the price $S_{i}$ of the
$i$-th underlying asset at time $\tilde{t}$ follows the diffusion process

$\frac{{\rm d}S_{i}}{S_{i}}=\mu_{i}\,{\rm d}\tilde{t}+\sigma_{i}\,{\rm
d}W_{i},\quad i=1,\ldots,d;$ (2)

see also [1, 2, 22, 24, 29]. The parameter $\mu_{i}$ represents the drift rate
of the price $S_{i}$, and $W_{i}$ is the standard Brownian motion with
volatility $\sigma_{i}$. The Wiener processes $W_{i}$ and $W_{j}$ are
correlated with a correlation coefficient $\rho_{ij}$. Therefore, the
corresponding covariance matrix $\mathbf{Q}$ is a $d\times d$ matrix with
entries $\mathbf{Q}_{ij}=\rho_{ij}\sigma_{i}\sigma_{j}$.

Using the stochastic equation (2) and Itô calculus, a deterministic equation
for the expected value $V$ of the option is obtained. Let $r$ be a risk-free
interest rate, $T$ be the maturity date, and $t=T-\tilde{t}$ be time to
maturity. Then, the fair value of the option
$V\left(S_{1},\ldots,S_{d},t\right)$ corresponding to prices of underlying
assets $S_{i}$ at the time to maturity $t$ can be computed as the solution of
the partial differential equation

$\frac{\partial V}{\partial
t}=\frac{1}{2}\sum_{i=1}^{d}\sum_{j=1}^{d}\rho_{i,j}\sigma_{i}\sigma_{j}S_{i}S_{j}\frac{\partial^{2}V}{\partial
S_{i}\partial S_{j}}+r\sum_{i=1}^{d}S_{i}\frac{\partial V}{\partial
S_{i}}-rV,$ (3)

where $S_{i}>0$ and $t\in\left(0,T\right)$. Equation (3) is degenerate in the
sense that some second-order terms vanish when $S_{i}=0$.

The equation has to be equipped with appropriate initial and boundary
conditions, which depend on the option type. The initial condition represents
the value of the option at maturity and has the form

$V\left(S_{1},\ldots,S_{d},0\right)=V_{0}\left(S_{1},\ldots,S_{d}\right),\quad
S_{i}>0,$ (4)

where $V_{0}$ is the payoff function. Examples of payoff functions are given
in Section 5. Boundary conditions are discussed for transformed and localized
equations further below.

###  2.1 Transformation into log prices and drift removal

Due to the degeneracy of equation (3), it is not possible to derive a
variational formulation and perform the analysis in the standard Sobolev
spaces; special weighted Sobolev spaces must be used instead, see [1].
Therefore, it is beneficial to remove this degeneracy first. It is well-known
[1, 22] that this can be done using substitution into logarithmic prices
$y_{i}=\log S_{i}$, $i=1,\ldots,d$. Then, the function

$U\left(y_{1},\ldots,y_{d},t\right)=V\left(e^{y_{1}},\ldots,e^{y_{d}},t\right)=V\left(S_{1},\ldots,S_{d},t\right)$
(5)

representing the option value in logarithmic prices is the solution to the
transformed equation

$\frac{\partial U}{\partial
t}=\frac{1}{2}\sum_{i=1}^{d}\sum_{j=1}^{d}\rho_{i,j}\sigma_{i}\sigma_{j}\frac{\partial^{2}U}{\partial
y_{i}\partial
y_{j}}+\sum_{i=1}^{d}\left(r-\frac{\sigma_{i}^{2}}{2}\right)\frac{\partial
U}{\partial y_{i}}-rU.$ (6)

Now, the differential operator on the right-hand side of equation (6) is an
elliptic operator with constant coefficients.

First-order terms in equation (6) represent drift and can be removed using
another substitution, $x_{i}=y_{i}-b_{i}t$, $b_{i}=\sigma_{i}^{2}/2-r$,
$i=1,\ldots,d.$ Then, a function

$W\left(x_{1},\ldots,x_{d},t\right)=U\left(x_{1}+b_{1}t,\ldots,x_{d}+b_{d}t,t\right)=U\left(y_{1},\ldots,y_{d},t\right)$
(7)

is a solution to a symmetric problem

$\frac{\partial W}{\partial
t}=\frac{1}{2}\sum_{i=1}^{d}\sum_{j=1}^{d}\rho_{i,j}\sigma_{i}\sigma_{j}\frac{\partial^{2}W}{\partial
x_{i}\partial x_{j}}-rW.$ (8)

###  2.2 Localization

In order to solve equation (8) numerically, it is convenient to approximate
the unbounded domain $\mathbb{R}^{d}$ for $\left(x_{1},\ldots,x_{d}\right)$
with a bounded domain. Therefore, we choose minimal and maximal values for
$S_{i}$, denote them $S_{i}^{min}$ and $S_{i}^{max}$, respectively, and define

$x_{i}^{min}=\ln S_{i}^{min},\quad x_{i}^{max}=\ln S_{i}^{max},\quad
I_{i}=\left(x_{i}^{min},x_{i}^{max}\right).$ (9)

Then, the equation to be solved is

$\frac{\partial W}{\partial
t}-\frac{1}{2}\sum_{i=1}^{d}\sum_{j=1}^{d}\rho_{i,j}\sigma_{i}\sigma_{j}\frac{\partial^{2}W}{\partial
x_{i}\partial x_{j}}+rW=0,\quad x_{i}\in I_{i},\quad t\in\left(0,T\right).$
(10)

In Section 3, a wavelet basis $\Psi$ on $\Omega=\left(0,1\right)^{d}$ is
designed. To solve (10), we can either use a linear transformation to convert
$\Psi$ to the product domain $I_{1}\times I_{2}\times\ldots\times I_{d}$ or
transform (10) to $\Omega$. Here, we use the second approach. Let
$d_{i}=x_{i}^{max}-x_{i}^{min}$ and

$u\left(x_{1},\ldots,x_{d},t\right)=W\left(\frac{x_{1}-x_{1}^{min}}{d_{1}},\ldots,\frac{x_{d}-x_{d}^{min}}{d_{d}}\right).$
(11)

The function $u$ is a solution to

$\frac{\partial u}{\partial
t}-\mathcal{D}\left(u\right)=0,\quad\left(x_{1},\ldots,x_{d}\right)\in\Omega,\quad
t\in\left(0,T\right),$ (12)

where

$\mathcal{D}\left(u\right)=\sum_{i=1}^{d}\sum_{j=1}^{d}P_{i,j}\frac{\partial^{2}u}{\partial
x_{i}\partial x_{j}}-ru,\quad
P_{i,j}=\frac{\rho_{i,j}\sigma_{i}\sigma_{j}}{2\,d_{i}d_{j}}.$ (13)

The initial condition is

$u\left(x_{1},\ldots,x_{d},0\right)=u_{0}\left(x_{1},\ldots,x_{d}\right),\quad\left(x_{1},\ldots,x_{d}\right)\in\Omega,$
(14)

where $u_{0}$ is the transformed payoff function,

$u_{0}\left(x_{1},\ldots,x_{d}\right)=V_{0}\left(\exp\left(\frac{x_{1}-x_{1}^{min}}{d_{1}}\right),\ldots,\exp\left(\frac{x_{d}-x_{d}^{min}}{d_{d}}\right)\right).$
(15)

Prescribing appropriate boundary conditions is a crucial task because exact
values of $u$ on the boundary $\partial\Omega$ of a domain $\Omega$ are
generally not known, and the computation of approximate values on the boundary
can be quite complicated. Therefore, it is convenient to approximate the
boundary conditions by homogeneous Dirichlet boundary conditions

$u\left(x_{1},\ldots,x_{d},t\right)=0,\quad\left(x_{1},\ldots,x_{d}\right)\in\partial\Omega,\quad
t\in\left(0,T\right).$ (16)

However, because these boundary conditions are artificial and are not
satisfied by the original function $u$, the error can be substantial in the
vicinity of $\partial\Omega$. Thus, to obtain a sufficiently accurate solution
in a chosen region of interest in $\Omega$, it is essential to choose
sufficiently large $\Omega$, see [22].

##  3  Orthogonal Spline Wavelets

The efficiency of Galerkin-type methods depends crucially on the used basis
functions. As already mentioned, this paper focuses on the wavelet-Galerkin
method, and thus, the appropriate wavelet basis is needed. In this section, we
make assumptions on the wavelet basis on the bounded interval, and we give a
detailed description of a construction of a wavelet basis meeting all the
assumptions, specifically an orthogonal cubic spline wavelet basis. First,
orthogonal cubic scaling functions and wavelets on the real line are designed.
Special boundary scaling functions and wavelets are then proposed to adapt the
basis from the real line to the bounded interval and to the homogeneous
Dirichlet boundary conditions. Finally, the sparse tensor product is used to
construct a wavelet basis on the product domain $\Omega$. We show that the
resulting basis is well conditioned with respect to the $H^{1}$-seminorm and
study its properties.

Let $\left\langle\cdot,\cdot\right\rangle$ and $\left\|\cdot\right\|$ denote
the inner product and the norm of the space $L^{2}\left(\Omega\right)$,
respectively. Let $\left\|\cdot\right\|_{s}$ be the norm of the Sobolev space
$H^{s}\left(\Omega\right)$, and $\left|\cdot\right|_{1}$ be the seminorm in
$H^{1}\left(\Omega\right)$. The space of functions from
$H^{1}\left(\Omega\right)$ with zero trace on $\partial\Omega$ is denoted
$H_{0}^{1}\left(\Omega\right)$.

Furthermore, let $\mathcal{J}$ be an index set such that each index
$\lambda\in\mathcal{J}$ takes the form $\lambda=\left(j,k\right)$, and let
$\left|\lambda\right|=j$ represent a level. For
$\mathbf{v}=\left\\{v_{\lambda}\right\\}_{\lambda\in\mathcal{J}}$,
$v_{\lambda}\in\mathbb{R}$, define

$\left\|\mathbf{v}\right\|_{2}=\left(\sum_{\lambda\in\mathcal{J}}\left|v_{\lambda}\right|^{2}\right)^{1/2},\quad
l^{2}\left(\mathcal{J}\right)=\left\\{\mathbf{v}:\left\|\mathbf{v}\right\|_{2}<\infty\right\\}.$
(17)

We make assumptions on a family of functions
$\Psi=\left\\{\psi_{\lambda},\lambda\in\mathcal{J}\right\\}$ to be a suitable
orthogonal spline-wavelet basis for $L^{2}\left(I\right)$, where $I$ is a
bounded interval or $I=\mathbb{R}$.

  * $A1)$

$L^{2}$ orthogonality. The set $\Psi$ is an orthogonal basis of
$L^{2}\left(I\right)$, and basis functions are normalized, i.e.,
$\left\langle\psi_{\lambda},\psi_{\lambda}\right\rangle=1$ and
$\left\langle\psi_{\lambda},\psi_{\mu}\right\rangle=0,$ if $\lambda\neq\mu.$

  * $A2)$

Locality. For all $\lambda\in\mathcal{J}$, a diameter of the support of
$\psi_{\lambda}$ is bounded by $C2^{-\left|\lambda\right|}$ with $C$
independent of $\left|\lambda\right|$ and at a given level, the supports of
only finitely many wavelets overlap at any point $x$.

  * $A3)$

Spline basis. Basis functions are continuous piecewise-polynomial functions of
order $p\geq 2$.

  * $A4)$

Hierarchical structure. The set $\Psi$ has a hierarchical structure

$\Psi=\Phi_{j_{0}}\cup\bigcup_{j=j_{0}}^{\infty}\Psi_{j},$ (18)

where

$\Phi_{j_{0}}=\left\\{\phi_{j_{0},k},k\in\mathcal{I}_{j_{0}}\right\\},\quad\Psi_{j}=\left\\{\psi_{j,k},k\in\mathcal{J}_{j}\right\\}.$
(19)

Functions $\phi_{j,k}$ are called scaling functions, and functions
$\psi_{j,k}$ are called wavelets. Furthermore, assume that all the wavelets
$\psi_{j,k}$ are translations and dilations of certain wavelet generators
$\psi_{l}$, i.e., they are of the form

$\psi_{j,k}=2^{j/2}\psi_{l}\left(2^{j}\cdot-m\right),$ (20)

where $l$ and $m$ are integers dependent on $k$.

  * $A5)$

Vanishing moments. We assume that wavelets $\psi_{j,l}$ have $L\geq 1$
vanishing moments, i.e., for any polynomial $p$ of degree less than $L$,

$\int_{I}p\left(x\right)\,\psi_{j,l}\left(x\right)dx=0,\quad
l\in\mathcal{J}_{j},\ j\geq j_{0}.$ (21)

In some constructions, condition (21) cannot be satisfied by so-called
boundary wavelets. In these cases, we admit an exception and we require that
(21) be satisfied only by inner wavelets.

Furthermore, we need to work with spaces $H_{0}^{1}\left(I\right)$, where $I$
is a bounded interval. First, recall a definition of a Riesz basis.

######  Definition 1.

Let $H$ be a separable Hilbert space with the norm $\left\|\cdot\right\|_{H}$.
The set $\Gamma=\left\\{\gamma_{\lambda},\lambda\in\mathcal{J}\right\\}$ is
called a Riesz basis of $H$ if the closure of the span of $\Gamma$ is $H$ and
there exist constants $c,C\in\left(0,\infty\right)$ such that for all
$\mathbf{b}=\left\\{b_{\lambda}\right\\}_{\lambda\in\mathcal{J}}\in
l^{2}\left(\mathcal{J}\right)$, we have

$c\left\|\mathbf{b}\right\|_{2}\leq\left\|\sum_{\lambda\in\mathcal{J}}b_{\lambda}\gamma_{\lambda}\right\|_{H}\leq
C\left\|\mathbf{b}\right\|_{2}.$ (22)

Constant $c_{\psi}:=\text{sup}\left\\{c:c\ \text{satisfies}\
(\ref{riesz})\right\\}$ is called a lower Riesz bound, constant
$C_{\psi}:=\text{inf}\left\\{C:C\ \text{satisfies}\ (\ref{riesz})\right\\}$ is
called an upper Riesz bound, and $cond\ \Psi=C_{\psi}/c_{\psi}$ is called a
condition number of $\Psi$.

If $\Gamma$ satisfies (22), but the closure of its span is not necessarily
dense in $H$, then $\Gamma$ is called a Riesz sequence in $H$.

Now, we require that functions from $\Psi$ normalized in the $H^{1}$-norm be a
Riesz basis of the space $H_{0}^{1}\left(I\right)$. Thus, we formulate another
assumption $A6)$.

  * $A6)$

$H^{1}$ Riesz basis. The set
$\left\\{\psi_{\lambda}/\left\|\psi_{\lambda}\right\|_{1},\psi_{\lambda}\in\Psi\right\\}$
is a Riesz basis of the space $H_{0}^{1}\left(I\right)$.

The construction of wavelets typically starts with a multiresolution analysis
of the space $L^{2}\left(\mathbb{R}\right)$.

######  Definition 2.

The sequence $\left\\{V_{j}\right\\}_{j\in\mathbb{Z}}$ of closed subspaces in
$L^{2}\left(\mathbb{R}\right)$ is called a multiresolution analysis (MRA) of
multiplicity $p$, if the following conditions are satisfied:

  * i)

$V_{j}\subset V_{j+1}$ for all $j\in\mathbb{Z}$.

  * ii)

For all $j\in\mathbb{Z}$, $f\in V_{j}$ if and only if
$f\left(2^{-j}\cdot\right)\in V_{0}$.

  * iii)

$\bigcup\limits_{j\in\mathbb{Z}}V_{j}$ is dense in
$L^{2}\left(\mathbb{R}\right)$.

  * iv)

$\bigcap\limits_{j\in\mathbb{Z}}V_{j}=\left\\{0\right\\}$.

  * v)

There exist functions $\phi_{1},\ldots,\phi_{p}$ such that

$\Phi_{0}=\left\\{\phi_{i}\left(\cdot-k\right),i=1,\ldots,p,k\in\mathbb{Z}\right\\}$
(23)

is a Riesz basis of $V_{0}$.

Functions $\phi_{i}$, $i=1,\ldots,p,$ are said to generate the multiresolution
analysis $\left\\{V_{j}\right\\}_{j\in\mathbb{Z}}$ and are called scaling
generators. If $\Phi_{0}$ is an orthonormal basis of $V_{0}$, then
$\left\\{V_{j}\right\\}_{j\in\mathbb{Z}}$ is called an orthogonal
multiresolution analysis.

###  3.1  Construction of orthogonal cubic spline scaling basis

In this section, we use some of the general design principles from [14] to
create a scaling basis. We start with the well-known scaling generators
already used in [5, 9]. Let $\xi_{1}$ and $\xi_{2}$ be Hermite cubic splines
defined as

$\xi_{1}(x)=\left\\{\begin{array}[]{cl}\left(x+1\right)^{2}\left(1-2x\right),&x\in[-1,0],\\\
\left(1-x\right)^{2}\left(1+2x\right),&x\in[0,1],\\\ 0,&{\rm
otherwise},\end{array}\right.$ (24)

and

$\xi_{2}(x)=\left\\{\begin{array}[]{cl}\left(x+1\right)^{2}x,&x\in[-1,0],\\\
\left(1-x\right)^{2}x,&x\in[0,1],\\\ 0,&{\rm otherwise}.\end{array}\right.$
(25)

Figure 1 displays graphs of functions $\xi_{1}$ and $\xi_{2}$.

Figure 1:  Hermite cubic spline generators $\xi_{1}$ (left) and $\xi_{2}$
(right).

For $j\in\mathbb{N}$ and $k\in\mathbb{Z}$, let $\xi_{j,k}$ be translations and
dilations of these two generators,

$\xi_{j,2k+l-1}\left(x\right)=2^{j/2}\xi_{l}\left(2^{j}x-k\right)\ {\rm
for}\quad l=1,2.$ (26)

It is well known that spaces

$\tilde{V}_{j}=\overline{{\rm span}\left\\{\xi_{j,k},k\in\mathbb{Z}\right\\}}$
(27)

form a multiresolution analysis, which is not orthogonal; see [5, 9]. The aim
of this section is to construct an orthogonal MRA formed by spaces $V_{j}$
satisfying (1). This MRA cannot be generated by only two functions $\xi_{1}$
and $\xi_{2}$ as in (27), but several new generators must be constructed.

As in [14], these generators are constructed to have their supports in
$\left[-1,1\right]$. Because it is required that (1) be satisfied, these
generators must be linear combinations of functions $\xi_{j,k}$ supported in
$\left[-1,1\right]$ for a certain level $j$. For $j\leq 0$, there are not
enough functions $\xi_{j,k}$ supported in $\left[-1,1\right]$; therefore we
start with the level $j=1$.

First, we limit ourselves to the interval $I_{0}=\left[0,1\right]$. For $j=1$,
there are two functions $\xi_{1,2}$ and $\xi_{1,3}$ supported in $I_{0}$. Due
to the symmetry of $\xi_{1}$ and the antisymmetry of $\xi_{2}$, it is easy to
verify that functions

$\phi_{1}=\frac{\xi_{1,2}}{\left\|\xi_{1,2}\right\|},\quad\phi_{2}=\frac{\xi_{1,3}}{\left\|\xi_{1,3}\right\|}$
(28)

are $L^{2}$-orthogonal. The scaling generators $\phi_{1}$ and $\phi_{2}$ are
displayed in Figure 2.

To construct another two generators with supports in $I_{0}$, we take
functions from $\tilde{V}_{2}$ with supports in $I_{0}$, which are functions
$\xi_{2,k}$ for $k=2,\ldots,7$. Let define spaces

$\mathcal{A}_{0}=\text{span}\left\\{\xi_{1,2},\xi_{1,3}\right\\},\quad\mathcal{A}_{1}=\text{span}\left\\{\xi_{2,k},k=2,\ldots,7\right\\}.$
(29)

This means that $\mathcal{A}_{0}$ is the space generated by all functions
$\xi_{j,k}$ on the coarsest level $j=1$ supported in $I_{0}$ and
$\mathcal{A}_{1}$ is the space generated by all functions $\xi_{j,k}$ on level
$j=2$ supported in $I_{0}$. Next, we find an orthogonal basis of
$\mathcal{A}_{1}\ominus\mathcal{A}_{0}$, where $\ominus$ denotes the
orthogonal complement.

Let $\mathbf{C}$ be the matrix with entries

$C_{i,j}=\left<\xi_{1,i+1},\xi_{2,j+1}\right>,\quad i=1,2,\quad j=1,\ldots,6,$
(30)

and $\mathbf{b}^{1},\ldots,\mathbf{b}^{4}$ be vectors forming a basis of the
null space of $\mathbf{C}$. Defining functions

$w_{i}=\sum_{j=1}^{6}b^{i}_{j}\,\xi_{2,j+1},\quad i=1,\ldots,4,$ (31)

and orthogonalizing them, we obtain functions

$v_{1}=\frac{w_{1}}{\left\|w_{1}\right\|},\quad
v_{i}=\frac{w^{i}-\sum_{j=1}^{i-1}\left<w_{i},v_{j}\right>v_{j}}{\left\|w_{i}-\sum_{j=1}^{i-1}\left<w_{i},v_{j}\right>v_{j}\right\|},\quad
i=2,3,4,$ (32)

which are normalized and form an orthogonal basis for
$\mathcal{A}_{1}\ominus\mathcal{A}_{0}$.

The next step of the construction is based on the theory developed in [14].
For the convenience of the reader, similar notations for spaces are used as
therein. Let

$\displaystyle\mathcal{B}_{0}$ $\displaystyle=$
$\displaystyle\text{span}\left(\left\\{\xi_{1,0}\chi_{[0,1]},\xi_{1,1}\chi_{[0,1]}\right\\}\cup\mathcal{A}_{0}\right),$
(33) $\displaystyle\mathcal{B}_{1}$ $\displaystyle=$
$\displaystyle\text{span}\left(\left\\{\xi_{1,4}\chi_{[0,1]},\xi_{1,5}\chi_{[0,1]}\right\\}\cup\mathcal{A}_{0}\right),$

where $\chi_{[0,1]}$ is the characteristic function of the interval
$\left[0,1\right]$. This means that $\mathcal{B}_{0}$ is the space generated
by left boundary functions and $\mathcal{B}_{1}$ is the space generated by
right boundary functions. Furthermore, for $h=0,1,$ let

$\mathcal{C}_{h}=\mathcal{B}_{h}\ominus\mathcal{A}_{0}.$ (34)

The following lemma is crucial for the construction.

######  Lemma 3.

With the above settings, let $W$ be a subspace of
$\mathcal{A}_{1}\ominus\mathcal{A}_{0}$ such that

$\left(I-P_{W}\right)\mathcal{C}_{0}\perp\left(I-P_{W}\right)\mathcal{C}_{1},$
(35)

where $P_{W}$ is the orthogonal projection, $I$ is identity, and $\perp$
denotes orthogonality. If $\phi_{3}$ and $\phi_{4}$ form an orthogonal basis
for $W$, then $V_{j}$ generated by $\phi_{1}$, $\phi_{2}$, $\phi_{3}$,
$\phi_{4}$, $\xi_{1,0}$ and $\xi_{1,1}$ form an $L^{2}$-orthogonal
multiresolution analysis satisfying (1).

###### Proof.

Lemma 3 is a consequence of a more general Lemma 3.1 from [14]. ∎

It follows from (33) and (34) that the space $\mathcal{C}_{0}$ is generated by
functions

$c_{i}=\xi_{1,i-1}\chi_{[0,1]}-\left<\xi_{1,i-1},\phi_{1}\right>\phi_{1}-\left<\xi_{1,i-1},\phi_{2}\right>\phi_{2},\quad
i=1,2.$ (36)

Similarly, the space $\mathcal{C}_{1}$ is generated by functions

$c_{i}=\xi_{1,i+1}\chi_{[0,1]}-\left<\xi_{1,i+1},\phi_{1}\right>\phi_{1}-\left<\xi_{1,i+1},\phi_{2}\right>\phi_{2},\quad
i=3,4.$ (37)

Because $W\subset\mathcal{A}_{1}\ominus\mathcal{A}_{0}$ and
$v_{1},v_{2},v_{3},$ and $v_{4}$ form a basis for
$\mathcal{A}_{1}\ominus\mathcal{A}_{0}$, basis functions $\phi_{3}$ and
$\phi_{4}$ for $W$ can be expressed as

$\phi_{3}=\sum_{j=1}^{4}s_{1,j}v_{j},\quad\phi_{4}=\sum_{j=1}^{4}s_{2,j}v_{j},$
(38)

where the coefficients $s_{i,j}$ are such that condition (35) is satisfied.
This condition can be rewritten as

$0=\left<c_{i}-\sum_{m=1}^{2}\sum_{j=1}^{4}s_{m,j}\left<c_{i},v_{j}\right>v_{j},c_{k}-\sum_{n=1}^{2}\sum_{l=1}^{4}s_{n,l}\left<c_{k},v_{l}\right>v_{l}\right>,$
(39)

for $i=1,2$, and $k=3,4.$ Due to the orthogonality of functions $v_{j}$, we
have

$\left<c_{i},c_{k}\right>=\sum_{m=1}^{2}\sum_{j=1}^{4}s_{i,j}s_{k,j}\left<c_{i},v_{j}\right>\left<c_{k},v_{j}\right>.$
(40)

Furthermore, we require $\phi_{3}$ and $\phi_{4}$ to be orthogonal and
normalized, which yields additional conditions

$\sum_{j=1}^{4}s_{i,j}^{2}=1,\quad i=1,2,\quad\sum_{j=1}^{4}s_{1,j}s_{2,j}=0.$
(41)

Hence, (40) and (41) form a system of $7$ quadratic equations with $8$
unknowns, which cannot be solved explicitly. Therefore, we solve it
numerically by the Levenberg–Marquardt algorithm [31]. According to Lemma 3,
the resulting functions $\phi_{3}$ and $\phi_{4}$ are other two scaling
generators orthogonal to $\phi_{1}$ and $\phi_{2}$. All the scaling generators
supported in $I_{0}$ are displayed in Figure 2.

Figure 2: _Orthogonal cubic spline generators $\phi_{1}$ (top left),
$\phi_{2}$ (top right), $\phi_{3}$ (bottom left), and $\phi_{4}$ (bottom
right) supported in $\left[0,1\right]$. _

Lemma 3 states that functions $\phi_{1}$, $\phi_{2}$, $\phi_{3}$, $\phi_{4}$,
$\xi_{1,0}$ and $\xi_{1,1}$ generate spaces $V_{j}$, which form an orthogonal
multiresolution analysis satisfying (1). However, only functions $\phi_{1}$,
$\phi_{2}$, $\phi_{3}$, and $\phi_{4}$ are orthogonal. Therefore, we perform
orthogonalization and normalization to obtain functions generating an
orthogonal basis. Hence, the two remaining generators $\phi_{5}$ and
$\phi_{6}$ are given by

$g_{5}=\xi_{1,0}-\sum_{j=1}^{4}\left<\xi_{1,0},\phi_{j}\right>\phi_{j}-\sum_{j=1}^{4}\left<\xi_{1,0},\phi_{j}\left(\cdot+1\right)\right>\phi_{j}\left(\cdot+1\right),\quad\phi_{5}=\frac{g_{5}}{\left\|g_{5}\right\|},$
(42)

and

$g_{6}=\xi_{1,1}-\sum_{j=1}^{5}\left<\xi_{1,1},\phi_{j}\right>\phi_{j}-\sum_{j=1}^{5}\left<\xi_{1,1},\phi_{j}\left(\cdot+1\right)\right>\phi_{j}\left(\cdot+1\right),\quad\phi_{6}=\frac{g_{6}}{\left\|g_{6}\right\|}.$
(43)

Generators $\phi_{5}$ and $\phi_{6}$ supported in $\left[-1,1\right]$ are
displayed in Figure 3.

Figure 3: _Orthogonal cubic spline generators $\phi_{5}$ (left) and $\phi_{6}$
(right) with support $\left[-1,1\right]$. _

Using translations and dilations of the generators, we obtain scaling
functions

$\displaystyle\phi_{j,6k+l}\left(x\right)$ $\displaystyle=$ $\displaystyle
2^{j/2}\phi_{l}\left(2^{j}x-k\right),\ \quad l=1,2,3,4,\quad k\in\mathbb{Z},$
(44) $\displaystyle\phi_{j,6k+l}\left(x\right)$ $\displaystyle=$
$\displaystyle 2^{j/2}\phi_{l}\left(2^{j}x-k-1\right),\ \quad l=5,6,\quad
k\in\mathbb{Z},$

and multiresolution spaces

$V_{j}=\overline{{\rm span}\left\\{\phi_{j,k},k\in\mathbb{Z}\right\\}}.$ (45)

The properties of scaling functions and spaces $V_{j}$ are summarized in the
following theorem.

######  Theorem 4.

Let $\phi_{i}$, $\phi_{j,k},$ and $V_{j}$ be defined as above.

  * a)

Generators $\phi_{i}$, $i=1,\ldots,6$, and scaling functions $\phi_{j,k}$,
$j,k\in\mathbb{Z}$, are cubic splines belonging to
$C^{1}\left(\mathbb{R}\right)$.

  * b)

The sequence $\left\\{V_{j}\right\\}_{j\in\mathbb{Z}}$ is an orthogonal
multiresolution analysis satisfying (1).

###### Proof.

  * a)

Because all the generators $\phi_{i}$, $i=1,\ldots,6,$ are constructed as
linear combinations of cubic splines from $C^{1}\left(\mathbb{R}\right)$, they
are also cubic splines from $C^{1}\left(\mathbb{R}\right)$. Definition (44)
implies that $\phi_{j,k}$, $j,k\in\mathbb{Z}$ are also cubic splines and
belong to $C^{1}\left(\mathbb{R}\right)$.

  * b)

The fact that the sequence $\left\\{V_{j}\right\\}_{j\in\mathbb{Z}}$ is an
orthogonal multiresolution analysis satisfying (1) follows from Lemma 3 and
the above construction.

∎

As already mentioned, the system of equations (40) and (41) has more than one
solution. Moreover, instead of $\phi_{1}$ and $\phi_{2}$, we can take their
orthogonalized linear combinations and, similarly, instead of $w_{i}$, we can
take their linear combinations and orthogonalize them to obtain different
functions $v_{i}$. Therefore, there are infinitely many choices of generators
$\phi_{i}$, $i=1,\ldots,6.$ We performed extensive numerical experiments to
find various solutions. In Appendix A, we present a variant of generators that
leads to a small condition number of the resulting normalized wavelet basis
with respect to the $H^{1}$-seminorm because it also plays a significant role
in applications; see [4, 5].

Furthermore, the part of the construction is the Gram–Schmidt
orthogonalization, which can lead to inaccurate results due to round-off
errors. Therefore, we computed all the coefficients using quadruple precision
and then rounded the results to 16 digits. Computing the scalar products of
the resulting functions, we verified numerically that the presented
coefficients are sufficiently accurate. For the computation, Matlab and
Multiprecision Computing Toolbox Advanpix were used.

###  3.2  Construction of orthogonal cubic spline wavelets

Wavelet generators corresponding to the multiresolution analysis
$\left\\{V_{j}\right\\}_{j\in\mathbb{Z}}$ are constructed as generators of the
orthogonal basis of the space $V_{1}\ominus V_{0}$.

The first step is the construction of generators with supports in $I_{0}$.
These generators must be linear combinations of the basis functions of $V_{1}$
with supports in $I_{0}$, which are the functions
$p_{k}=\phi_{k}\left(2\cdot\right)$ for $k=1,\ldots,4$, and
$p_{k+4}=\phi_{k}\left(2\cdot-1\right)$ for $k=1,\ldots,6$. The generators
must be orthogonal to all functions from $V_{0}$ with the part or whole
support in $I_{0}$, i.e., to functions $q_{i}=\phi_{i}$, $i=1,\ldots,6,$ and
$q_{i+2}=\phi_{i}\left(\cdot-1\right)$ for $i=5,6$. Therefore, the
coefficients of the linear combinations can be computed as the solution of the
homogeneous system of linear equations with the matrix $\mathbf{S}$ with
entries $S_{i,j}=\left<p_{i},q_{j}\right>.$ Because functions $p_{i}$, as well
as functions $q_{j}$, are linearly independent, the null space of the matrix
$\mathbf{S}$ is two dimensional. Let vectors $\mathbf{b}^{1}$ and
$\mathbf{b}^{2}$ form its basis. Then

$u_{i}=\sum_{j=1}^{8}b^{i}_{j}\,q_{j},\quad i=1,2,$ (46)

are functions from $V_{1}$ with supports in $I_{0}$. Using the Gram–Schmidt
orthogonalization and $L^{2}$-normalization of $u_{1}$ and $u_{2}$, we obtain
generators

$\psi_{1}=\frac{u_{1}}{\left\|u_{1}\right\|},\quad\psi_{2}=\frac{u_{2}-\left<u_{2},\psi_{1}\right>\psi_{1}}{\left\|u_{2}-\left<u_{2},\psi_{1}\right>\psi_{1}\right\|}.$
(47)

Graphs of functions $\psi_{1}$ and $\psi_{2}$ are displayed in Figure 4.

Figure 4: _Orthogonal wavelet generators $\psi_{1}$ and $\psi_{2}$ with
support $\left[0,1\right]$. _

Now, we construct generators with supports $\left[-1,1\right]$. First,
consider functions from $V_{1}$, which have parts of the support both in
$\left(-\infty,0\right)$ and $\left(0,\infty\right)$, i.e., functions
$\phi_{1,-1}=2^{1/2}\phi_{5}\left(2\cdot\right)$ and
$\phi_{1,0}=2^{1/2}\phi_{6}\left(2\cdot\right)$. Let us denote

$\mathcal{D}_{j}=\overline{{\rm
span}\left(\left\\{\phi_{j,k},k\in\mathbb{Z}\right\\}\backslash\left\\{\phi_{5}\left(2^{j}\cdot\right),\phi_{6}\left(2^{j}\cdot\right)\right\\}\right)}.$
(48)

The functions $\phi_{5}\left(2\cdot\right)$ and $\phi_{6}\left(2\cdot\right)$
are already orthogonal to the functions from $\mathcal{D}_{1}$ because
$\phi_{2,k}$, $k\in\mathbb{Z}$ are orthogonal. Because
$\mathcal{D}_{0}\subset\mathcal{D}_{1}$, these functions are also orthogonal
to the functions from $\mathcal{D}_{0}$. That is, they are orthogonal to all
scaling functions $\phi_{0,k}$, except $\phi_{0,-1}=\phi_{5}$ and
$\phi_{0,0}=\phi_{6}$. Therefore, functions

$h_{i}=\phi_{i+4}\left(2\cdot\right)-\left<\phi_{i+4}\left(2\cdot\right),\phi_{5}\right>\phi_{5}-\left<\phi_{i+4}\left(2\cdot\right),\phi_{6}\right>\phi_{6},\quad
i=1,2,$ (49)

are orthogonal to all scaling functions $\phi_{0,k}$, $k\in\mathbb{Z}.$

Now, orthogonalization and normalization are again performed to obtain wavelet
generators

$\psi_{3}=\frac{h_{1}}{\left\|h_{1}\right\|},\quad\psi_{4}=\frac{h_{2}-\left<h_{2},\psi_{3}\right>\psi_{3}}{\left\|h_{2}-\left<h_{2},\psi_{3}\right>\psi_{3}\right\|}.$
(50)

The remaining two wavelet generators with the support in $\left[-1,1\right]$
are defined as linear combinations of functions from $V_{1}$ with the whole
support in $\left[-1,1\right]$,

$\psi_{i+4}=\sum_{k=-11}^{10}c_{i,k}\phi_{1,k},\quad i=1,2.$ (51)

We need to find coefficients $c_{i,k}$ such that $\psi_{5}$ and $\psi_{6}$ are
orthogonal to functions from $V_{0}$, which have part of the support or the
whole support in $I_{1}$, i.e., to functions $\phi_{0,k}$ for $k=-7,\ldots,6$.
Furthermore, these generators must be orthogonal to already constructed
wavelets $\psi_{i}$, $i=1,\ldots,4,$ and their translations
$\psi_{i}\left(\cdot+1\right)$ and $\psi_{i}\left(\cdot-1\right)$ for $i=3,4$.

According to the following lemma, we can omit some of these conditions.

######  Lemma 5.

If functions $\psi_{5}$ and $\psi_{6}$ are given by (51) and
$\left\langle\psi_{i},\phi_{0,k}\right\rangle$ for $i=5,6$ and
$k=-7,\ldots,6$, then for $i=5,6,$ and $l=3,4$, we have

$\left\langle\psi_{i},\psi_{l}\left(\cdot-1\right)\right\rangle=0,\quad\left\langle\psi_{i},\psi_{l}\left(\cdot+1\right)\right\rangle=0.$
(52)

###### Proof.

According to definitions of generators $\psi_{l}$, $l=3,4$, we can write

$\psi_{l}\left(\cdot-1\right)=d_{l,1}\phi_{5}\left(2\cdot-2\right)+d_{l,2}\phi_{6}\left(2\cdot-2\right)+d_{l,3}\phi_{5}\left(\cdot-1\right)+d_{l,4}\phi_{6}\left(\cdot-1\right),$
(53)

where $d_{l,i}$ are real coefficients. Assumptions of the lemma imply that

$\left\langle\psi_{i},\phi_{m}\left(\cdot-1\right)\right\rangle=0,\quad
i,m=5,6.$ (54)

Due to (51) and orthogonality of functions $\phi_{1,k}$, we obtain

$\left\langle\psi_{i},\phi_{m}\left(2\cdot-2\right)\right\rangle=\frac{1}{\sqrt{2}}\sum_{k=-11}^{10}c_{i,k}\left\langle\phi_{1,k},\phi_{1,m+6}\right\rangle=0,\quad
i,m=5,6.$ (55)

Relations (53), (54) and (55) imply
$\left\langle\psi_{i},\psi_{l}\left(\cdot-1\right)\right\rangle=0$ for
$i=5,6,$ and $l=3,4$. The proof is similar for functions
$\psi_{l}\left(\cdot+1\right)$, $l=3,4$. ∎

Hence, let $\mathbf{T}$ be a matrix with entries

$\displaystyle T_{i,j}$ $\displaystyle=$
$\displaystyle\left\langle\phi_{0,i-8},\phi_{1,j-12}\right\rangle,i=1,\ldots,14,\,j=1,\ldots,22,$
(56) $\displaystyle T_{i,j}$ $\displaystyle=$
$\displaystyle\left\langle\psi_{i-14}\left(\cdot+1\right),\phi_{1,j-12}\right\rangle,i=15,16,\,j=1,\ldots,22,$
$\displaystyle T_{i,j}$ $\displaystyle=$
$\displaystyle\left\langle\psi_{i-16},\phi_{1,j-12}\right\rangle,i=17,\ldots,20,\,j=1,\ldots,22.$

Denote $\mathbf{b}^{1}$ and $\mathbf{b}^{2}$ vectors forming a basis of the
null space of $\mathbf{T}$ and define

$z_{i}=\sum_{k=-11}^{10}b^{i}_{k}\,\phi_{1,k},\quad i=1,2.$ (57)

Then, we can set

$\psi_{5}=\frac{z_{1}}{\left\|z_{1}\right\|},\quad\psi_{6}=\frac{z_{2}-\left<z_{2},\psi_{5}\right>\psi_{5}}{\left\|z_{2}-\left<z_{2},\psi_{5}\right>\psi_{5}\right\|}.$
(58)

Graphs of all wavelet generators with supports $\left[-1,1\right]$ are shown
in Figure 5, and polynomial coefficients of all wavelet generators are given
in Appendix A.

Figure 5: _Orthogonal wavelet generators $\psi_{3}$ (top left), $\psi_{4}$
(top right), $\psi_{5}$ (bottom left), and $\psi_{6}$ (bottom right) with
supports $\left[-1,1\right]$. _

We obtain wavelets $\psi_{j,k}$ using translations and dilations of the
constructed generators,

$\displaystyle\psi_{j,6k+l}\left(x\right)$ $\displaystyle=$ $\displaystyle
2^{j/2}\psi_{l}\left(2^{j}x-k\right),\ \quad l=1,2,\quad j,k\in\mathbb{Z},$
(59) $\displaystyle\psi_{j,6k+l}\left(x\right)$ $\displaystyle=$
$\displaystyle 2^{j/2}\psi_{l}\left(2^{j}x-k-1\right),\ \quad l=3,4,5,6,\quad
j,k\in\mathbb{Z}.$

These wavelets form wavelet spaces

$W_{j}=\overline{{\rm span}\left\\{\psi_{j,k},k\in\mathbb{Z}\right\\}},\quad
j\in\mathbb{Z}.$ (60)

The following theorem summarizes the properties of the constructed wavelets
and wavelet spaces and shows that

$\Psi=\left\\{\phi_{0,k},k\in\mathbb{Z}\right\\}\cup\left\\{\psi_{j,k},j\geq
0,k\in\mathbb{Z}\right\\}$ (61)

is a wavelet basis of the space $L^{2}\left(\mathbb{R}\right).$ The symbol
$\delta_{i,j}$ denotes the Kronecker delta, i.e., $\delta_{ii}=1$ and
$\delta_{ij}=0$ for $i\neq j$.

######  Theorem 6.

Under the above setting, the following statements hold:

  * a)

Generators $\psi_{i}$, $i=1,\ldots,6$, and wavelet $\psi_{j,k}$,
$j,k\in\mathbb{Z}$, are cubic splines belonging to
$C^{1}\left(\mathbb{R}\right)$.

  * b)

For $i=1,\ldots,6,$ scaling generators $\phi_{i}$ and wavelet generators
$\psi_{i}$ satisfy

$\left\langle\phi_{i},\phi_{j}\right\rangle=\delta_{i,j},\quad\left\langle\psi_{i},\psi_{j}\right\rangle=\delta_{i,j},\quad\left\langle\phi_{i},\psi_{j}\right\rangle=0.$
(62)

  * c)

For $i,j,k,l\in\mathbb{Z},$ scaling functions $\phi_{i,k}$ and wavelets
$\psi_{j,l}$ satisfy orthogonality conditions

$\left\langle\phi_{i,k},\phi_{i,l}\right\rangle=\delta_{k,l},\quad\left\langle\psi_{i,k},\psi_{j,l}\right\rangle=\delta_{i,j}\delta_{k,l},$
(63)

and if $i\leq j$, then

$\left\langle\phi_{i,k},\psi_{j,l}\right\rangle=0.$ (64)

  * d)

Wavelets $\psi_{j,k}$, $j,k\in\mathbb{Z}$, have four vanishing moments.

  * e)

The set $\Psi$ is an orthogonal basis of the space
$L^{2}\left(\mathbb{R}\right).$

  * f)

Basis functions from $\Psi$ are local in a sense that ${\rm
diam\,supp}\,\phi_{0,k}\leq 2$ and ${\rm diam\,supp}\,\psi_{j,k}\leq 2^{1-j}$.

###### Proof.

  * a)

Because $\psi_{j,k}$, $j,k\in\mathbb{Z}$ are linear combinations of scaling
functions, the assertion a) is valid.

  * b)

This statement follows directly from the construction.

  * c)

The orthogonality of scaling functions $\phi_{i,k}$ follows from b) and (44).
Similarly, orthogonality of wavelets $\psi_{i,k}$ and $\psi_{i,l}$ follows
from b) and (59). Due to $\left\langle\phi_{i},\psi_{j}\right\rangle=0$ in b),
$V_{j}\perp W_{j}$. Since $V_{j}\perp W_{j}$ and $V_{i}\subset V_{j}$ and
$W_{i}\subset V_{j}$ for $i\leq j$, we have $V_{i}\perp W_{j}$ and $W_{i}\perp
W_{j}$ for $i\leq j$, which implies (63) and (64).

  * d)

Due to c), wavelets $\psi_{j,k}$, $j,k\in\mathbb{Z}$ are orthogonal to all
functions from $V_{j}$, which is a cubic spline space containing all cubic
polynomials. The orthogonality of $\psi_{j,k}$ to all cubic polynomials means
that wavelets have four vanishing moments.

  * e)

The orthogonality of $\Psi$ follows from d). By Theorem 4 b), $\Psi$ generates
$L^{2}\left(\mathbb{R}\right)$.

  * f)

The length of the support of generators $\phi_{i}$ and $\psi_{i}$,
$i=1,\ldots,6,$ is at most $2$. Hence, due to (44) and (59), the assertion is
valid.

∎

Due to Theorem 6, all the requirements $A1)-A5)$ are met and thus $\Psi$ is
indeed an orthogonal wavelet basis of $L^{2}\left(\mathbb{R}\right).$

###  3.3  Orthogonal wavelet basis on the interval

In this section, we adapt the constructed wavelet basis from the real line to
a bounded interval and homogeneous Dirichlet boundary conditions, and study
the properties of the resulting basis.

It is sufficient to focus on the interval $\left[0,1\right]$ because a wavelet
basis on another bounded interval can be obtained by a linear transformation.
The adaptation consists of preserving those basis functions whose whole
supports are contained in $\left[0,1\right]$ and constructing new boundary
functions.

First, we investigate the situation at the point $0$. The functions $\phi_{5}$
and $\phi_{6}$ are the only two scaling functions on the level zero with
nonzero value at $0$. Because we want to obtain a new generator that is
orthogonal to other functions on $\left[0,1\right]$ and is zero at $0$, we
define this boundary scaling generator as

$\phi_{L}=\left(\phi_{6}\left(0\right)\phi_{5}-\phi_{5}\left(0\right)\phi_{6}\right)\chi_{\left[0,1\right]}.$
(65)

Using a similar approach at the point $1$ leads to a boundary scaling function

$\phi_{R}=\left(\phi_{6}\left(0\right)\phi_{5}\left(\cdot-1\right)-\phi_{5}\left(0\right)\phi_{6}\left(\cdot-1\right)\right)\chi_{\left[0,1\right]}.$
(66)

The graphs of boundary scaling generators are displayed in Figure 6.

Figure 6: _Boundary orthogonal scaling generators $\psi_{L}$ (left) and
$\psi_{R}$ (right) with the support $\left[0,1\right]$. _

When it comes to boundary wavelet generators at the point $0$, the situation
is similar in the sense that new boundary wavelets are linear combinations of
wavelets that are nonzero at $0$. In addition, the new wavelet generators must
be orthogonal to $\phi_{L}$. Therefore, we define functions $g_{L_{1}}$ and
$g_{L_{2}}$ as

$g_{L_{i}}=\sum_{j=1}^{4}b_{j}^{i}\,\psi_{j+2},\quad i=1,2.$ (67)

Because we require $g_{L_{i}}\left(0\right)=0$ and
$\left\langle\phi_{L},g_{L_{i}}\right\rangle=0$, the coefficients $b_{j}^{i}$
are elements of vectors $\mathbf{b}^{i}$ forming a basis of the null space of
the matrix $\mathbf{G}$ given by

$G_{1,j}=\psi_{j+2}\left(0\right),\quad
G_{2,j}=\left<\phi_{L},\psi_{j+2}\right>,\quad j=1,\ldots,4.$ (68)

Generators $\psi_{L_{1}}$ and $\psi_{L_{2}}$ are obtained by orthogonalization
and normalization of $g_{L_{1}}$ and $g_{L_{2}}$, similarly to the generators
above.

The construction of right boundary wavelet generators is analogous. The
functions $h_{R_{1}}$ and $h_{R_{2}}$ are defined as

$h_{R_{i}}=\sum_{j=1}^{4}c_{j}^{i}\,\psi_{j}\left(\cdot-1\right),\quad i=1,2,$
(69)

where vectors $\mathbf{c}^{i}$ form the null space of matrix $\mathbf{H}$ with
entries

$H_{1,j}=\psi_{j+2}\left(0\right),\quad
H_{2,j}=\left<\phi_{R},\psi_{j+2}\left(\cdot-1\right)\right>,\quad
j=1,\ldots,4.$ (70)

Generators $\psi_{R_{1}}$ and $\psi_{R_{2}}$ are obtained by orthogonalization
and normalization of $h_{R_{1}}$ and $h_{R_{2}}$. The boundary wavelet
generators are displayed in Figure 7.

Figure 7: _Boundary orthogonal wavelet generators $\psi_{L_{1}}$ (top left),
$\psi_{L_{2}}$ (top right), $\psi_{R_{1}}$ (bottom left), and $\psi_{R_{2}}$
(bottom right) with the support $\left[0,1\right]$. _

Let us define scaling functions on the level $0$ as $\phi_{0,1}=\phi_{L}$,
$\phi_{0,6}=\phi_{R}$, and $\phi_{0,k}=\phi_{k-1}$ for $k=2,\ldots,5$. On each
level $j\geq 0$, we define boundary wavelets

$\psi_{j,i}=\psi_{L_{i}}\left(2^{j}\cdot\right),\quad\psi_{j,6\cdot
2^{j}+i-2}=\psi_{R_{i}}\left(2^{j}\left(1-\cdot\right)\right),\quad i=1,2.$
(71)

and inner wavelets

$\psi_{j,6k+l+2}\left(x\right)=2^{j/2}\psi_{l}\left(2^{j}x-k\right),\quad
l=1,\ldots,6,\quad 3\leq 6k+l\leq 6\cdot 2^{j}-2.$ (72)

We define an index set $\mathcal{J}_{j}=\left\\{k\in\mathbb{Z},1\leq k\leq
6\cdot 2^{j}\right\\}$ and assign the constructed functions to sets

$\Phi_{0}=\left\\{\phi_{0,k},k\in\mathcal{J}_{0}\right\\},\quad\Psi_{j}=\left\\{\psi_{j,k},k\in\mathcal{J}_{j}\right\\},$
(73)

and

$\Psi^{I}=\Phi_{0}\cup\bigcup_{j=0}^{\infty}\Psi_{j},\quad\Psi^{k}=\Phi_{0}\cup\bigcup_{j=0}^{k}\Psi_{j}.$
(74)

Our aim is to show that $\Psi^{I}$ is an orthogonal wavelet basis of the space
$L^{2}\left(0,1\right).$ The properties $A1)-A5)$ are studied in the following
theorem.

######  Theorem 7.

Under the above setting, the following statements hold:

  * a)

All basis functions from $\Psi^{I}$ are cubic splines belonging to
$C^{1}\left(0,1\right)$.

  * b)

Basis functions from $\Psi^{I}$ satisfy orthogonality conditions

$\left\langle\phi_{0,k},\phi_{0,l}\right\rangle=\delta_{k,l},\quad\left\langle\psi_{i,k},\psi_{j,l}\right\rangle=\delta_{i,j}\delta_{k,l},\quad\left\langle\phi_{0,k},\psi_{j,l}\right\rangle=0.$
(75)

  * c)

Inner wavelets $\psi_{j,k}$, $k=3,\ldots,6\cdot 2^{j}-2$ have four vanishing
moments.

  * d)

The set $\Psi^{I}$ is an orthogonal basis of $L^{2}\left(0,1\right).$

  * f)

Basis functions from $\Psi^{I}$ are local in the sense that ${\rm
diam\,supp}\,\phi_{0,k}\leq 2$ and ${\rm diam\,supp}\,\psi_{j,k}\leq 2^{1-j}$.

###### Proof.

The proof is analogous to the proof of Theorem 6. ∎

Due to Theorem 7, we can conclude that $\Psi^{I}$ is an orthogonal wavelet
basis of $L^{2}\left(0,1\right).$ It remains to show that the assumption $A6)$
is also satisfied. To this end, we employ the following theorem. For the
proof, see [4, 8, 12].

######  Theorem 8.

Let $j_{0}\in\mathbb{N}$ and for $j\geq j_{0}$ let $V_{j}$ and $\tilde{V}_{j}$
be subspaces of the space $H\subset L^{2}\left(I\right)$ such that
$V_{j}\subset V_{j+1}$, $\tilde{V}_{j}\subset\tilde{V}_{j+1}$, and ${\rm dim}\
V_{j}={\rm dim}\ \tilde{V}_{j}<\infty$. Let $\Phi_{j}$ be bases of $V_{j}$,
$\tilde{\Phi}_{j}$ be bases of $\tilde{V}_{j}$, and $\Psi_{j}$ be bases of
$\tilde{V}_{j}^{\perp}\cap V_{j+1}$, where $\tilde{V}_{j}^{\perp}$ denotes the
$L^{2}$-orthogonal complement of $\tilde{V}_{j}$ in $H$. Moreover, let the
Riesz bounds with respect to the $L^{2}$-norm of $\Phi_{j}$,
$\tilde{\Phi}_{j}$, and $\Psi_{j}$ be uniformly bounded. Let $\Psi$ be
composed of $\Phi_{j_{0}}$ and $\Psi_{j}$, $j\geq j_{0}$, as in (18).
Furthermore, we assume that

$\mathbf{\Gamma}_{j}=\left\langle\Phi_{j},\tilde{\Phi}_{j}\right\rangle$ (76)

is invertible and that the spectral norm of $\mathbf{\Gamma}_{j}^{-1}$ is
bounded independently on $j$. In addition, for some positive constants $C$,
$\gamma$ and $d$, such that $\gamma<d$, let

$\inf_{v_{j}\in V_{j}}\left\|v-v_{j}\right\|\leq
C2^{-jt}\left\|v\right\|_{t},\ v\in H^{t}\left(I\right)\cap H,\ 0\leq t\leq
d,$ (77)

and

$\left\|v_{j}\right\|_{s}\leq C2^{js}\left\|v_{j}\right\|,\ v_{j}\in V_{j},\
0\leq s<\gamma,$ (78)

and, similarly, let (77) and (78) hold for $\tilde{\gamma}$ and $\tilde{d}$ on
the dual side. Then

$\left\\{\psi_{\lambda}/\left\|\psi_{\lambda}\right\|_{s},\psi_{\lambda}\in\Psi\right\\}$
(79)

is a Riesz sequence in $H^{s}\left(I\right)$ for
$s\in\left(-\tilde{\gamma},\gamma\right)$.

Estimate (77) characterizes the approximation properties of spaces $V_{j}$ and
is often referred to as the Jackson estimate. Estimate (78) relates to the
smoothness properties of $V_{j}$, and is called the Bernstein estimate.

######  Theorem 9.

For $s\in\left(-1.5,1.5\right)$, the set
$\left\\{\psi_{\lambda}/\left\|\psi_{\lambda}\right\|_{s},\psi_{\lambda}\in\Psi^{I}\right\\}$
is a Riesz sequence in $H^{s}$, and the set
$\left\\{\psi_{\lambda}/\left\|\psi_{\lambda}\right\|_{1},\psi_{\lambda}\in\Psi^{I}\right\\}$
is a Riesz basis of $H_{0}^{1}\left(0,1\right).$

###### Proof.

Due to the orthogonality, matrix $\mathbf{\Gamma}_{j}$ defined by (76) is an
identity matrix. The estimates (77) and (78) are satisfied for $d=\tilde{d}=4$
and $\gamma=\tilde{\gamma}=1.5$. These parameters depend on the polynomial
exactness and smoothness of the primal and dual spaces; see [4, 8]. Due to
these facts and Theorem 8, the proof is complete. ∎

As already mentioned, we performed extensive numerical experiments to
construct an $L^{2}$-orthogonal cubic spline wavelet basis on the interval
adapted to homogeneous boundary conditions, which has also small condition
number with respect to the $H^{1}$-seminorm. This basis is presented in A. In
Figure 8, the condition numbers of the set
$\left\\{\psi_{\lambda}/\left|\psi_{\lambda}\right|_{1},\psi_{\lambda}\in\Psi^{k}\right\\}$
with respect to $H^{1}$-seminorm are presented. These condition numbers are
square roots of the condition numbers of stiffness matrices for these sets and
are close to the condition numbers with respect to the $H^{1}$-norm.

Figure 8: Condition numbers with respect to the $H^{1}$-seminorm for the sets
$\Psi^{k}$ normalized in the $H^{1}$-seminorm.

###  3.4  Orthogonal wavelets on sparse tensor product spaces

Wavelet bases on the unit cube $\Omega$ can be created using the tensor
product of wavelet bases on the interval $\left(0,1\right)$. This can be
achieved by several approaches, such as an isotropic approach [7], an
anisotropic approach [6, 7], or the sparse tensor product [22]. When solving
higher-dimensional problems, the isotropic and anisotropic approaches lead to
an excessive number of functions $N\sim\mathcal{O}\left(2^{dk}\right)$ at
higher levels $k$, and $L^{2}$ error of approximation is of order
$\mathcal{O}\left(N^{-p/d}\right)$, where $p$ is order of spline wavelet. This
phenomenon, known as the “curse of dimensionality,” can be circumvented by
using sparse grids [19, 20, 22, 34, 38]. In this case, the number of basis
functions is $\mathcal{O}\left(2^{k}k^{d-1}\right)$ and the $L^{2}$–error of
approximation is the same as in one dimension, up to log terms. Therefore, we
use the sparse tensor product approach here.

Let $\Psi^{I}$ be a wavelet basis on the interval $I=\left[0,1\right]$ given
by (74). To simplify notations redefine
$\mathcal{J}_{0}=\left\\{-5,\ldots,6\right\\}$ and denote
$\psi_{0,k-6}=\phi_{0,k}$ for $k=1,\ldots,6.$ Then the set of all possible
indices for basis functions is

$\mathcal{J}=\left\\{\left(j,k\right),\ j\geq 0,\
k\in\mathcal{J}_{j}\right\\},$ (80)

and the wavelet basis $\Psi^{I}$ can also be expressed as

$\Psi^{I}=\left\\{\psi_{j,k},\ j\geq 0,\
k\in\mathcal{J}_{j}\right\\}=\left\\{\psi_{\lambda},\
\lambda\in\mathcal{J}\right\\}.$ (81)

Now, consider dimension $d\geq 1$. The index set is defined as

$\mathcal{J}=\left\\{\lambda=\left(\lambda_{1},\ldots,\lambda_{d}\right):\lambda_{i}=\left(j_{i},k_{i}\right),j_{i}\geq
0,k_{i}\in\mathcal{J}_{j_{i}}\right\\}.$ (82)

On the unit cube $\Omega$, we consider a set $\Psi$ containing functions

$\psi_{\lambda}=\bigotimes\limits_{i=1}^{d}\psi_{\lambda_{i}},\quad\lambda=\left(\lambda_{1},\ldots,\lambda_{d}\right)\in\mathcal{J},$
(83)

where $\otimes$ denotes the tensor product, i.e., $\left(u\otimes
v\right)\left(x_{1},x_{2}\right)=u\left(x_{1}\right)v\left(x_{2}\right)$.

The set of all such functions is
$\Psi=\left\\{\psi_{\lambda},\lambda\in\mathcal{J}\right\\}$. We show below
that generalizations of A1) – A6) are valid and, therefore, this set is also
referred to as a wavelet basis. Recall that the level of one-dimensional
function $\psi_{\lambda}$ is denoted $\left|\lambda\right|$. For $k\geq 0$,
the so-called sparse tensor product approach uses the set of functions

$\Psi^{k}=\left\\{\psi_{\lambda}:\lambda=\left(\lambda_{1},\ldots,\lambda_{d}\right),\
\sum_{i=1}^{d}\left|\lambda_{i}\right|\leq k\right\\},$ (84)

whereas the so-called anisotropic tensor product approach employs the set

$\Psi^{k}_{ani}=\left\\{\psi_{\lambda}:\lambda=\left(\lambda_{1},\ldots,\lambda_{d}\right),\
\left|\lambda_{i}\right|\leq k\right\\}.$ (85)

These two kinds of sets significantly differ in a number of functions, as
shown in the following lemma. Symbol $\\#$ denotes the cardinality.

######  Lemma 10.

For $d\geq 1$ and $k\geq 1$,

$\\#\Psi^{k}_{ani}=6^{d}2^{dk}=\mathcal{O}\left(2^{dk}\right),\quad\\#\Psi^{k}=\mathcal{O}\left(2^{k}k^{d-1}\right).$
(86)

More precisely,

$\displaystyle\\#\Psi^{k}$ $\displaystyle=$ $\displaystyle
6^{2}\left(k+2\right)2^{k+1},\ \text{if}\ d=2,$ (87)
$\displaystyle\\#\Psi^{k}$ $\displaystyle=$ $\displaystyle
6^{3}\left(k^{2}+7k+8\right)2^{k},\ \text{if}\ d=3,$
$\displaystyle\\#\Psi^{k}$ $\displaystyle=$ $\displaystyle
6^{4}\left(\frac{k^{3}}{6}+\frac{5k^{2}}{2}+\frac{28k}{3}+8\right)2^{k+1},\
\text{if}\ d=4,$ $\displaystyle\\#\Psi^{k}$ $\displaystyle=$ $\displaystyle
6^{5}\left(\frac{k^{4}}{24}+\frac{(13k^{3}}{12}+\frac{203k^{2}}{24}+\frac{269k}{12}+16\right)2^{k+1},\
\text{if}\ d=5.$

###### Proof.

Clearly, $\\#\Psi^{k}_{ani}=6\cdot 2^{k}$ for $d=1$ and $k\geq 1$, and thus
$\\#\Psi^{k}_{ani}=\left(6\cdot 2^{k}\right)^{d}=6^{d}2^{dk}$ for $d\geq 1$
and $k\geq 1$.

For $d=2$,

$\\#\Psi^{k}=\sum\limits_{i=0}^{k}\sum\limits_{j=0}^{k-i}\\#\mathcal{J}_{i}\\#\mathcal{J}_{j}=\\#\mathcal{J}_{0}\\#\mathcal{J}_{0}+2\\#\mathcal{J}_{0}\sum_{j=1}^{k}\\#\mathcal{J}_{j}+\sum\limits_{i=1}^{k-1}\sum\limits_{j=1}^{k-i}\\#\mathcal{J}_{i}\\#\mathcal{J}_{j}.$
(88)

Using $\\#\mathcal{J}_{0}=12$, $\\#\mathcal{J}_{j}=6\cdot 2^{j}$ for $j\geq
1$, and the formula for a sum of geometric series gives the relation for
$\\#\Psi^{k}$ for $d=2$.

A similar approach used for $d>2$ gives the remaining assertions of the lemma.
∎

The resulting multidimensional functions preserve many properties of original
univariate functions.

######  Theorem 11.

Under the above setting, the following statements hold:

  * a)

Functions $\psi_{\lambda}$ are bicubic splines belonging to
$C^{1}\left(\Omega\right)$.

  * b)

The set $\Psi$ is an orthogonal basis of $L^{2}\left(\Omega\right)$.

  * c)

Functions $\psi_{\lambda}$ satisfy $\left\langle
p,\psi_{\lambda}\right\rangle=0$, where
$\lambda=\left(\lambda_{1},\ldots,\lambda_{d}\right)$ such that
$\lambda_{i}=\left(j_{i},k_{i}\right)$, $3\leq k_{i}\leq 6\cdot 2^{j_{i}}-2$,
and $p$ is a polynomial of degree less than $4$, which means that inner
wavelets have four vanishing moments.

  * d)

Basis functions are local in the generalized sense that

${\rm diam}\,{\rm supp}\,\psi_{\lambda}\leq
C\sqrt{d}\,2^{-\left[\lambda\right]},\quad\left[\lambda\right]=\min_{i=1,\ldots,d}\left|\lambda_{i}\right|.$
(89)

###### Proof.

The properties a), b), and c) follow directly from a definition of
$\psi_{\lambda}$ and properties of the corresponding univariate functions
presented in Theorem 7.

Property d) follows from

${\rm diam}\,{\rm
supp}\,\psi_{\lambda}\leq\sqrt{\sum_{i=1}^{d}C_{i}^{2}2^{-2\left|\lambda_{i}\right|}}\leq
C\sqrt{d}2^{-\left[\lambda\right]},\quad C=\max\limits_{i=1,\ldots,d}C_{i}.$
(90)

∎

Hence, conditions $A1)-A5)$ are fulfilled in a generalized sense. The next
theorem studies condition $A6)$.

######  Theorem 12.

The set
$\left\\{\psi_{\lambda}/\left\|\psi_{\lambda}\right\|_{1},\psi_{\lambda}\in\Psi\right\\}$
is a Riesz basis of $H_{0}^{1}\left(\Omega\right)$.

###### Proof.

The space $H_{0}^{1}\left(\Omega\right)$ can be written as

$H_{0}^{1}\left(\Omega\right)=\bigcap\limits_{i=1}^{d}S^{i}\left(\Omega\right),$
(91)

where

$S^{i}\left(\Omega\right)=L^{2}\left(I_{1}\right)\otimes\ldots\otimes
L^{2}\left(I_{i-1}\right)\otimes H_{0}^{1}\left(I_{i}\right)\otimes
L^{2}\left(I_{i+1}\right)\otimes\ldots L^{2}\left(I_{d}\right),$ (92)

see [19]. Thus, Theorem 9 and (91) yield the required result. ∎

Because $\Psi$ is an $L^{2}$–orthogonal basis of $L^{2}\left(\Omega\right)$,
any function $u\in L^{2}\left(\Omega\right)$ has a series representation

$u=\sum_{\lambda\in\mathcal{J}}c_{\lambda}\psi_{\lambda}.$ (93)

Truncating this expansion, we obtain a projection operator $P_{k}$,

$P_{k}\,u=\sum_{\psi_{\lambda}\in\Psi^{k}}c_{\lambda}\psi_{\lambda},\quad u\in
L^{2}\left(\Omega\right).$ (94)

The approximation properties of wavelet sparse grids have been studied in [20,
22, 34, 38]. For the proof of the next result, refer therein.

######  Theorem 13.

Let $\Psi^{k}$ given by (84) be set of function constructed by the tensor
product of univariate spline-wavelet basis of order $p$ satisfying $A1)-A5)$.
Then, for $u\in H^{p}\left(\Omega\right)$ there exists a constant $C$
independent of level $k$ such that

$\left\|u-P_{k}u\right\|\leq C2^{-pk}k^{\frac{d-1}{2}}\left\|u\right\|_{p}.$
(95)

Hence, the approximation order is $\mathcal{O}\left(N^{-p}\right)$, if
dimension $d=1$. Using Lemma 10 and (95), we may conclude that, as mentioned
above, the approximation rate for $d$-dimensional functions is indeed, up to
the log term, the same as in one dimension, and the curse of dimensionality is
overcome.

##  4  Wavelet Method

This section examines a wavelet-based method for solving equation (12)
equipped with conditions (14) and (16). The method is based on the Galerkin
method using sparse tensor product spaces of orthogonal spline wavelets, such
as cubic spline wavelets from Section 3.4, and the Crank–Nicolson scheme with
Rannacher time stepping.

For $u,v\in H_{0}^{1}\left(\Omega\right)$, the bilinear form
$a:H_{0}^{1}\left(\Omega\right)\times
H_{0}^{1}\left(\Omega\right)\rightarrow\mathbb{R}$ corresponding to the
differential operator $\mathcal{D}$ given by (13) is defined as

$a\left(u,v\right)=\sum_{i=1}^{d}\sum_{j=1}^{d}P_{i,j}\left\langle\frac{\partial
u}{\partial x_{i}},\frac{\partial v}{\partial
x_{j}}\right\rangle+r\left\langle u,v\right\rangle.$ (96)

######  Lemma 14.

The bilinear form $a$ defined by (96) is continuous and coercive; that is,
there exists a constant $\beta$ such that

$\left|a\left(u,v\right)\right|\leq\beta\left\|u\right\|_{1}\left\|v\right\|_{1}\quad\forall\,u,v\in
H_{0}^{1}\left(\Omega\right),$ (97)

and a constant $\alpha>0$ such that

$a\left(u,u\right)\geq\alpha\left\|u\right\|^{2}_{1}\quad\forall\,u\in
H_{0}^{1}\left(\Omega\right).$ (98)

Specifically, the form $a$ satisfies the Gårding inequality; that is, there
exist constants $C_{1},C_{2}>0$ such that, for any $v\in
H_{0}^{1}\left(\Omega\right)$,

$a\left(v,v\right)\geq
C_{1}\left\|v\right\|_{1}^{2}-C_{2}\left\|v\right\|^{2}.$ (99)

###### Proof.

The continuity of the bilinear form $a$ is a consequence of definition (96)
and the Cauchy–Schwarz inequality. The coefficient matrix for the second order
terms of differential operator $\mathcal{D}$ is
$\mathbf{D}\mathbf{Q}\mathbf{D}/2$, where $\mathbf{D}$ is a diagonal matrix
with entries $D_{i,i}=1/d_{i}$, $d_{i}=x_{i}^{max}-x_{i}^{min}$, and
$\mathbf{Q}$ is a matrix defined in Section 2. Because $\mathbf{Q}$ is a
variational matrix, it is positive definite, and thus for any
$0\neq\mathbf{x}\in\mathbb{R}^{d}$,

$\mathbf{x}^{T}\mathbf{D}\mathbf{Q}\mathbf{D}\mathbf{x}=\left(\mathbf{D}\mathbf{x}\right)^{T}\mathbf{Q}\left(\mathbf{D}\mathbf{x}\right)>0.$
(100)

This implies that $\mathbf{D}\mathbf{Q}\mathbf{D}/2$ is also positive
definite, operator $\mathcal{D}$ is elliptic, and the bilinear form $a$ is
coercive. The coercivity property (98) implies (99) with $C_{2}=0$. ∎

Let $L^{2}\left(0,T;X\right)$ and $H^{1}\left(0,T;X\right)$ denote Bochner
spaces and $H^{-1}\left(\Omega\right)$ denote the dual space of
$H_{0}^{1}\left(\Omega\right)$. Let $u_{0}\in L^{2}\left(\Omega\right)$ be the
transformed payoff function. The variational formulation of the equation (12)
reads as: Find $u\in L^{2}\left(0,T;H_{0}^{1}\left(\Omega\right)\right)\cap
H^{1}\left(0,T;H^{-1}\left(\Omega\right)\right)$ such that

$\left\langle\frac{\partial U}{\partial
t},v\right\rangle+a\left(U,v\right)=0,\quad\forall v\in V,\,\text{a.e.
in}\left(0,T\right);\,U\left(\cdot,0\right)=u_{0}.$ (101)

######  Theorem 15.

There exists a unique solution $u$ to the problem (101).

###### Proof.

It is well-known that the existence and uniqueness of the solution follows
directly from the Gårding inequality and continuity of the bilinear form $a$
shown in Lemma 14, as detailed [1, 22, 26]. ∎

Let $\Psi$ be an orthonormal wavelet basis for the space
$L^{2}\left(\Omega\right)$ constructed using the tensor product from a
univariate wavelet basis satisfying assumptions $A1)-A6)$, and consider its
subset

$\Psi^{k}=\left\\{\psi_{\lambda}\in\Psi:\lambda\in\mathcal{J}^{k}\right\\},\quad\mathcal{J}^{k}=\left\\{\lambda:\left\|\lambda\right\|_{1}\leq
k\right\\},$ (102)

where $\left\|\cdot\right\|_{1}$ is the $l^{1}$ norm. This means that the sets
$\Psi^{k}$ generate sparse tensor product spaces similarly as in Section 3.4.
Let $X_{k}={\rm span}\,\Psi^{k}$ and denote the dual of $X_{k}$ as
$X_{k}^{\prime}$.

The wavelet-Galerkin method consists of finding $U^{k}\in
L^{2}\left(0,T;X_{k}\right)\cap H^{1}\left(0,T;X_{k}^{\prime}\right)$ such
that for all $v_{k}\in X_{k}$ and a.e. in $\left(0,T\right)$,

$\left\langle\frac{\partial U_{k}}{\partial
t},v_{k}\right\rangle+a\left(U_{k},v_{k}\right)=0,\quad
U_{k}\left(\cdot,0\right)=u_{0}.$ (103)

######  Theorem 16.

There exists a unique solution to the problem (103).

###### Proof.

This theorem also follows from the continuity of the bilinear form $a$ and the
Gårding inequality shown in Lemma 14. ∎

Now, the Crank–Nicolson scheme and Rannacher time-stepping are applied. Let
$M\in\mathbb{N}$ be a number of time steps, $\tau=T/M$ be the step size,
$t_{l}=l\tau$, $l=0,\ldots,M$, be a time level and
$U_{k}^{l}\left(x_{1},\ldots,x_{d}\right)=U_{k}\left(x_{1},\ldots,x_{d},t_{l}\right)$
be the solution at the $l$-th time level. The initial function $u_{0}$
representing the transformed payoff function is typically not smooth. It is
known that using the Crank–Nicolson scheme with non-smooth initial function
can lead to spurious oscillations and reduction of the accuracy of the scheme.
This can be avoided using Rannacher time-stepping, which consists of replacing
two Crank–Nicolson steps with four steps of the implicit Euler scheme with
half step size [35].

The resulting scheme has the following form. For $l=0,0.5,1,1.5$, use the
implicit-Euler scheme with step $\tau/2$,

$\frac{\left\langle
U_{k}^{l+1/2},v_{k}\right\rangle}{\tau/2}-\frac{\left\langle
U_{k}^{l},v_{k}\right\rangle}{\tau/2}+a\left(U_{k}^{l+1/2},v_{k}\right)=0,\quad\forall
v_{k}\in X_{k}.$ (104)

Then, for $l=2,\ldots,M-1$, use the Crank–Nicolson scheme

$\frac{\left\langle U_{k}^{l+1},v_{k}\right\rangle}{\tau}-\frac{\left\langle
U_{k}^{l},v_{k}\right\rangle}{\tau}+\frac{a\left(U_{k}^{l+1},v_{k}\right)}{2}+\frac{a\left(U_{k}^{l},v_{k}\right)}{2}=0,\quad\forall
v_{k}\in X_{k}.$ (105)

In the following, we focus only on the Crank–Nicolson scheme, which forms the
main part of the method. The study of the implicit Euler scheme is similar.
Inserting $v_{k}=\psi_{\mu}$ and using a series representation

$U_{k}^{l}=\sum_{\psi_{\lambda}\in\Psi^{k}}\left(\mathbf{c}_{k}^{l}\right)_{\lambda}\psi_{\lambda}$
(106)

into (105) yields for $l=2,\ldots,M-1,$ a linear system

$\mathbf{A}^{k}\mathbf{c}_{k}^{l+1}=\mathbf{f}_{k}^{l},$ (107)

where for $\psi_{\lambda},\psi_{\mu}\in\Psi^{k}$,

$\mathbf{A}^{k}_{\mu,\lambda}=\frac{\left\langle\psi_{\lambda},\psi_{\mu}\right\rangle}{\tau}+\frac{a\left(\psi_{\lambda},\psi_{\mu}\right)}{2},\quad\left(\mathbf{f}_{k}^{l}\right)_{\mu}=\frac{\left\langle
U_{k}^{l},\psi_{\mu}\right\rangle}{\tau}-\frac{a\left(U_{k}^{l},\psi_{\mu}\right)}{2}.$
(108)

The elements of matrix $\mathbf{A}^{k}$ can also be expressed as

$\mathbf{A}^{k}_{\mu,\lambda}=\left(\frac{1}{\tau}+\frac{r}{2}\right)\left\langle\psi_{\lambda},\psi_{\mu}\right\rangle+\sum_{i=1}^{d}\sum_{j=1}^{d}\frac{P_{i,j}}{2}\left\langle\frac{\partial\psi_{\lambda}}{\partial
x_{i}},\frac{\partial\psi_{\mu}}{\partial x_{j}}\right\rangle.$ (109)

Due to the orthogonality of functions $\psi_{\lambda_{i}}$ and
$\psi_{\mu_{i}}$, we have
$\left\langle\psi_{\lambda_{i}},\psi_{\mu_{j}}\right\rangle=\delta_{\lambda_{i},\mu_{j}}$.
Hence, we obtain

$\displaystyle\left\langle\frac{\partial\psi_{\lambda}}{\partial
x_{i}},\frac{\partial\psi_{\mu}}{\partial x_{j}}\right\rangle$
$\displaystyle=$ $\displaystyle\left\langle\frac{\partial}{\partial
x_{i}}\prod_{k=1}^{d}\psi_{\lambda_{k}}\left(x_{k}\right),\frac{\partial}{\partial
x_{j}}\prod_{l=1}^{d}\psi_{\mu_{l}}\left(x_{l}\right)\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle\psi^{\prime}_{\lambda_{i}},\psi_{\mu_{i}}\right\rangle\left\langle\psi_{\lambda_{j}},\psi^{\prime}_{\mu_{j}}\right\rangle\prod_{\begin{subarray}{c}k=1\\\
k\neq i,j\end{subarray}}^{d}\delta_{\lambda_{k},\mu_{k}},$

for $i\neq j$, and similarly

$\displaystyle\left\langle\frac{\partial\psi_{\lambda}}{\partial
x_{i}},\frac{\partial\psi_{\mu}}{\partial x_{i}}\right\rangle$
$\displaystyle=$ $\displaystyle\left\langle\frac{\partial}{\partial
x_{i}}\prod_{k=1}^{d}\psi_{\lambda_{k}}\left(x_{k}\right),\frac{\partial}{\partial
x_{i}}\prod_{l=1}^{d}\psi_{\mu_{l}}\left(x_{l}\right)\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle\psi^{\prime}_{\lambda_{i}},\psi^{\prime}_{\mu_{i}}\right\rangle\prod_{\begin{subarray}{c}k=1\\\
k\neq i\end{subarray}}^{d}\delta_{\lambda_{k},\mu_{k}}.$

Therefore, to compute the entries of the matrix $\mathbf{A}^{k}$, it is
necessary to first compute the products
$\left\langle\psi_{\lambda_{i}},\psi_{\mu_{i}}\right\rangle$,
$\left\langle\psi^{\prime}_{\lambda_{i}},\psi_{\mu_{i}}\right\rangle$, and
$\left\langle\psi^{\prime}_{\lambda_{i}},\psi^{\prime}_{\mu_{i}}\right\rangle$.

We consider matrices of these products. Let $\mathbf{M}^{k}$ and
$\mathbf{B}^{k}$ be matrices with elements

$\mathbf{M}^{k}_{\lambda,\mu}=\left\langle\psi^{\prime}_{\lambda},\psi^{\prime}_{\mu}\right\rangle,\quad\mathbf{B}^{k}_{\lambda,\mu}=\left\langle\psi^{\prime}_{\lambda},\psi_{\mu}\right\rangle,\quad\left|\lambda\right|,\left|\mu\right|\leq
k.$ (112)

Due to the orthogonality property of the basis used, the matrix
$\mathbf{I}^{k}$ with elements

$\mathbf{I}^{k}_{\lambda,\mu}=\left\langle\psi_{\lambda},\psi_{\mu}\right\rangle,\quad\left|\lambda\right|,\left|\mu\right|\leq
k,$ (113)

is an identity matrix.

The next lemma shows that the condition numbers of matrices $\mathbf{M}^{k}$
and $\mathbf{B}^{k}$ increase exponentially with the parameter $k$
characterizing the level of a univariate basis. Let $\left\|\cdot\right\|$
denote the spectral matrix norm and ${\rm cond}$ denote the condition number
of a matrix.

######  Lemma 17.

There exists a constant $C\in\mathbb{R}$ such that

$\left\|\mathbf{M}^{k}\right\|\leq C2^{2k},\quad{\rm cond}\,\mathbf{M}^{k}\leq
C2^{2k},\quad\left\|\mathbf{B}^{k}\right\|\leq C2^{k}.$ (114)

###### Proof.

By assumption $A6)$, the set $\Psi$, when normalized in the $H^{1}$-norm, is a
Riesz basis of $H_{0}^{1}\left(I\right)$. Therefore, there exist constants
$A_{1},A_{2}>0$ such that

$A_{1}^{2}\sum_{\lambda\in\mathcal{J}^{k}}d_{\lambda}^{2}\leq\left\|\sum_{\lambda\in\mathcal{J}^{k}}d_{\lambda}\frac{\psi_{\lambda}}{\left\|\psi_{\lambda}\right\|_{H^{1}}}\right\|^{2}_{1}\leq
A_{2}^{2}\sum_{\lambda\in\mathcal{J}^{k}}d_{\lambda}^{2}.$ (115)

By (20), we have
$d\,2^{\left|\lambda\right|}\leq\left\|\psi_{\lambda}\right\|_{1}\leq
D\,2^{\left|\lambda\right|}$. Hence, for
$d_{\lambda}=2^{\left|\lambda\right|}\,c_{\lambda}$, we obtain

$B_{1}^{2}\sum_{\lambda\in\mathcal{J}^{k}}2^{2\left|\lambda\right|}c_{\lambda}^{2}\leq\left\|\sum_{\lambda\in\mathcal{J}^{k}}c_{\lambda}\psi_{\lambda}\right\|^{2}_{1}\leq
B_{2}^{2}\sum_{\lambda\in\mathcal{J}^{k}}2^{2\left|\lambda\right|}c_{\lambda}^{2}$
(116)

for some constants $B_{1},B_{2}>0$. Because the $H^{1}$-norm and the
$H^{1}$-seminorm are equivalent in $H_{0}^{1}\left(I\right)$, the relation
(115) remains valid only with different constants when replacing the norm
$\left\|\cdot\right\|_{1}$ with the seminorm $\left|\cdot\right|_{1}$,

$C_{1}^{2}\sum_{\lambda\in\mathcal{J}^{k}}2^{2\left|\lambda\right|}c_{\lambda}^{2}\leq\left|\sum_{\lambda\in\mathcal{J}^{k}}c_{\lambda}\psi_{\lambda}\right|^{2}_{1}\leq
C^{2}_{2}\sum_{\lambda\in\mathcal{J}^{k}}2^{2\left|\lambda\right|}c_{\lambda}^{2}.$
(117)

This implies that

$\displaystyle C_{1}^{2}\sum_{\lambda\in\mathcal{J}^{k}}c_{\lambda}^{2}$
$\displaystyle\leq$ $\displaystyle
C_{1}^{2}\sum_{\lambda\in\mathcal{J}^{k}}2^{2\left|\lambda\right|}c_{\lambda}^{2}\leq\left|\sum_{\lambda\in\mathcal{J}^{k}}c_{\lambda}\psi_{j,l}\right|_{1}^{2}$
$\displaystyle\leq$ $\displaystyle
C_{2}^{2}\sum_{\lambda\in\mathcal{J}^{k}}2^{2\left|\lambda\right|}c_{\lambda}^{2}\leq
C_{2}^{2}\,2^{2k}\sum_{\lambda\in\mathcal{J}^{k}}c_{\lambda}^{2}.$

Relation (4) can be rewritten to matrix form

$C_{1}^{2}\left\|\mathbf{c}\right\|^{2}\leq\left(\mathbf{c}\right)^{T}\mathbf{M}^{k}\mathbf{c}\leq
C_{2}^{2}\,2^{2k}\left\|\mathbf{c}\right\|^{2}.$ (119)

Because matrix $\mathbf{M}^{k}$ is symmetric, using the Rayleigh quotient for
$\mathbf{M}^{k}$ and (119) imply that the maximal and minimal eigenvalues of
$\mathbf{M}$ satisfy

$\lambda_{min}\left(\mathbf{M}^{k}\right)\geq
C_{1}^{2},\quad\left\|\mathbf{M}^{k}\right\|=\lambda_{max}\left(\mathbf{M}^{k}\right)\leq
C_{2}^{2}2^{2k},$ (120)

and thus ${\rm cond}\,\mathbf{M}^{k}\leq C_{2}^{2}\,2^{2k}/C_{1}^{2}$.

Similarly, using Cauchy–Schwarz inequality,

$\displaystyle\left|\left\langle\sum\limits_{j\in\mathcal{J}^{k}}c_{\lambda}\psi_{\lambda},\sum\limits_{j\in\mathcal{J}^{k}}c_{\lambda}\psi^{\prime}_{\lambda}\right\rangle\right|$
$\displaystyle\leq$
$\displaystyle\left\|\sum\limits_{j\in\mathcal{J}^{k}}c_{\lambda}\psi_{\lambda}\right\|\left|\sum\limits_{j\in\mathcal{J}^{k}}c_{\lambda}\psi_{\lambda}\right|_{1}$
$\displaystyle\leq$ $\displaystyle
C_{2}\left(\sum\limits_{j\in\mathcal{J}^{k}}c_{\lambda}^{2}\right)^{1/2}\left(\sum\limits_{j\in\mathcal{J}^{k}}2^{2\left|\lambda\right|}c_{\lambda}^{2}\right)^{1/2}$
$\displaystyle\leq$ $\displaystyle
C_{2}\,2^{k}\sum\limits_{j\in\mathcal{J}^{k}}c_{\lambda}^{2}.$

This is equivalent to

$\left(\mathbf{c}\right)^{T}\mathbf{B}^{k}\mathbf{c}\leq
C_{2}\,2^{k}\left\|\mathbf{c}\right\|^{2}.$ (122)

In contrast to the matrix $\mathbf{M}^{k}$, the matrix $\mathbf{B}^{k}$ is not
symmetric. Therefore, the estimate of its spectral norm using eigenvalues
cannot be used here as above. However, it is known that

$\left\|\mathbf{N}\right\|\leq
2\sup\limits_{0\neq\mathbf{d}}\frac{\mathbf{d}^{T}\mathbf{N}\mathbf{d}}{\left\|\mathbf{d}\right\|^{2}},$
(123)

holds for any square matrix $\mathbf{N}$; see, e.g., Theorem 1.3-1 in [21].
This implies that, indeed, $\left\|\mathbf{B}^{k}\right\|\leq C2^{k}$,
$C=2\,C_{2}$. ∎

To analyze matrices $\mathbf{A}^{k}$ corresponding to bases constructed by
sparse tensor product, let us first consider discretization matrices
$\mathbf{A}^{k}_{ani}$ corresponding to bases constructed using anisotropic
approach. By (4) and (4), these matrices can be expressed as the sum of the
terms, which are tensor products of matrices $\mathbf{I}^{k}$,
$\mathbf{M}^{k}$, and $\mathbf{B}^{k}$,

$\displaystyle\mathbf{A}^{k}_{ani}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{\tau}+\frac{r}{2}\right)\mathbf{I}^{k}+\sum_{i=1}^{d}\frac{P_{i,i}}{2}\left(\bigotimes\limits_{j=1}^{i-1}\mathbf{I}^{k}\right)\otimes\mathbf{M}^{k}\otimes\left(\bigotimes\limits_{j=i+1}^{d}\mathbf{I}^{k}\right)$
$\displaystyle-$
$\displaystyle\sum_{i=1}^{d-1}\sum_{j=i+1}^{d}P_{i,j}\left(\bigotimes\limits_{k=1}^{i-1}\mathbf{I}^{k}\right)\otimes\mathbf{B}^{k}\otimes\left(\bigotimes\limits_{k=i+1}^{j-1}\mathbf{I}^{k}\right)\otimes\mathbf{B}^{k}\otimes\left(\bigotimes\limits_{k=j+1}^{d}\mathbf{I}^{k}\right).$

This is not valid for sparse tensor product bases. However, matrices
$\mathbf{A}^{k}$ can be divided into blocks such that each block can be
computed using the tensor product of the matrices corresponding to one-
dimensional bases similarly as in (4).

More precisely, for $m=\left(m_{1},\ldots,m_{d}\right)$ and
$n=\left(n_{1},\ldots,n_{d}\right)$ such that $\left\|m\right\|_{1}\leq k$ and
$\left\|n\right\|_{1}\leq k$, let $\mathbf{A}^{m,n}$ be a block of
$\mathbf{A}^{k}$ with elements

$\mathbf{A}^{m,n}_{\lambda,\mu}=\left(\frac{1}{\tau}+\frac{r}{2}\right)\left\langle\psi_{\lambda},\psi_{\mu}\right\rangle+\sum_{i=1}^{d}\sum_{j=1}^{d}\frac{P_{i,j}}{2}\left\langle\frac{\partial\psi_{\lambda}}{\partial
x_{i}},\frac{\partial\psi_{\mu}}{\partial x_{j}}\right\rangle,$ (125)

where $\lambda=\left(\lambda_{1},\ldots,\lambda_{d}\right)$ and
$\mu=\left(\mu_{1},\ldots,\mu_{d}\right)$ are such that
$\left|\lambda_{i}\right|=m_{i}$ and $\left|\mu_{i}\right|=n_{i}$ for
$i=1,\ldots d$. This means that $m$ and $n$ characterize levels of univariate
functions $\psi_{\lambda_{i}}$ and $\psi_{\mu_{i}}$, which form basis
functions $\psi_{\lambda}$ and $\psi_{\mu}$ via the tensor product. Thus, the
matrix $\mathbf{A}^{k}$ can be written as

$\mathbf{A}^{k}=\left(\mathbf{A}^{m,n}\right)_{0\leq\left\|m\right\|_{1},\left\|n\right\|_{1}\leq
k}.$ (126)

Similarly to matrix $\mathbf{A}^{k}_{ani}$, each block $\mathbf{A}^{m,n}$ can
be computed as

$\displaystyle\mathbf{A}^{m,n}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{\tau}+\frac{r}{2}\right)\mathbf{I}+\sum_{i=1}^{d}\frac{P_{i,i}}{2}\bigotimes\limits_{j=1}^{i-1}\mathbf{I}\otimes\mathbf{M}^{m_{i},n_{i}}\otimes\bigotimes\limits_{j=i+1}^{d}\mathbf{I}$
$\displaystyle-$
$\displaystyle\sum_{i=1}^{d-1}\sum_{j=i+1}^{d}P_{i,j}\bigotimes\limits_{k=1}^{i-1}\mathbf{I}\otimes\mathbf{B}^{m_{i},n_{i}}\otimes\bigotimes\limits_{k=i+1}^{j-1}\mathbf{I}\otimes\mathbf{B}^{m_{j},n_{j}}\otimes\bigotimes\limits_{k=j+1}^{d}\mathbf{I}.$

Here, $\mathbf{I}$ is an identity matrix of appropriate size and
$\mathbf{M}^{m,n}$ and $\mathbf{B}^{m,n}$ are matrices with elements

$\mathbf{M}^{m,n}_{\lambda,\mu}=\left\langle\psi^{\prime}_{\lambda},\psi^{\prime}_{\mu}\right\rangle,\quad\mathbf{B}^{m,n}_{\lambda,\mu}=\left\langle\psi^{\prime}_{\lambda},\psi_{\mu}\right\rangle,\quad\left|\lambda\right|=m,\quad\left|\mu\right|=n.$
(128)

The fact that many matrices in (4) are identity matrices greatly simplifies
and streamlines numerical computations with these matrices in comparison with
other non-orthogonal tensor product bases.

Due to the orthogonality of the basis and the block tensor product structure
(4), the matrices $\mathbf{A}^{k}$ are uniformly conditioned even without
preconditioning, which is studied in the following theorem.

######  Theorem 18.

Let $C_{2}$ be an upper Riesz bound with respect to the $H^{1}$ seminorm given
by (117). If

$\tau<\frac{\gamma}{A\,4^{k}},\quad
A=2\,C_{2}^{2}\sum\limits_{i=1}^{d}\sum\limits_{j=1}^{d}P_{i,j},\quad\gamma\in\left(0,1\right),$
(129)

then there exists a constant $C$ independent of $k$ such that ${\rm
cond}\,\mathbf{A}^{k}\leq C$.

###### Proof.

The matrix $\mathbf{A}^{k}$ is a submatrix of the matrix
$\mathbf{A}^{k}_{ani}$ and both matrices are symmetric positive definite.
Therefore, ${\rm cond}\,\mathbf{A}^{k}\leq{\rm cond}\,\mathbf{A}^{k}_{ani}$
and it is sufficient to show the uniform boundedness of condition numbers for
$\mathbf{A}^{k}_{ani}$.

For the sake of simplicity, let denote

$\displaystyle\mathbf{C}^{k}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{d}\frac{P_{i,i}}{2}\left(\bigotimes\limits_{j=1}^{i-1}\mathbf{I}\right)\otimes\mathbf{M}^{k}\otimes\left(\bigotimes\limits_{j=i+1}^{d}\mathbf{I}\right)$
$\displaystyle-$
$\displaystyle\sum_{i=1}^{d-1}\sum_{j=i+1}^{d}P_{i,j}\left(\bigotimes\limits_{k=1}^{i-1}\mathbf{I}\right)\otimes\mathbf{B}^{k}\otimes\left(\bigotimes\limits_{k=i+1}^{j-1}\mathbf{I}\right)\otimes\mathbf{B}^{k}\otimes\left(\bigotimes\limits_{k=j+1}^{d}\mathbf{I}\right).$

The proof of Theorem 17 implies that

$\displaystyle\left\|\mathbf{C}^{k}\right\|\leq\sum_{i=1}^{d}\frac{P_{i,i}}{2}\left\|\mathbf{M}^{k}\right\|+\sum_{i=1}^{d-1}\sum_{j=i+1}^{d}P_{i,j}\left\|\mathbf{B}^{k}\right\|^{2}\leq
2\,C_{2}^{2}\,2^{2k}\sum_{i=1}^{d}\sum_{j=1}^{d}P_{i,j}.$ (131)

Thus, we have

$\left\|\mathbf{A}^{k}_{ani}\right\|=\left\|\left(\frac{1}{\tau}+\frac{r}{2}\right)\mathbf{I}+\mathbf{C}^{k}\right\|\leq\left(\frac{1}{\tau}+\frac{r}{2}\right)+A\,2^{2k}$
(132)

with $A$ given by (129).

Now, we aim to estimate the norm of the inverse of $\mathbf{A}^{k}_{ani}$.
Using

$\left(\left(\frac{1}{\tau}+\frac{r}{2}\right)\mathbf{I}+\mathbf{C}^{k}\right)\left(\left(\frac{1}{\tau}+\frac{r}{2}\right)\mathbf{I}+\mathbf{C}^{k}\right)^{-1}=\mathbf{I},$
(133)

we obtain

$\left(\frac{1}{\tau}+\frac{r}{2}\right)\left(\left(\frac{1}{\tau}+\frac{r}{2}\right)\mathbf{I}+\mathbf{C}^{k}\right)^{-1}=\mathbf{I}-\mathbf{C}^{k}\left(\left(\frac{1}{\tau}+\frac{r}{2}\right)\mathbf{I}+\mathbf{C}^{k}\right)^{-1}.$
(134)

To simplify notations, denote

$\mathbf{F}=\left(\left(\frac{1}{\tau}+\frac{r}{2}\right)\mathbf{I}+\mathbf{C}^{k}\right)^{-1}.$
(135)

Taking spectral norm of (134), we have

$\displaystyle\left(\frac{1}{\tau}+\frac{r}{2}\right)\left\|\mathbf{F}\right\|$
$\displaystyle\leq$ $\displaystyle
1+\left\|\mathbf{C}^{k}\right\|\cdot\left\|\mathbf{F}\right\|$ (136)
$\displaystyle\left\|\mathbf{F}\right\|$ $\displaystyle\leq$
$\displaystyle\frac{1}{\left(\frac{1}{\tau}+\frac{r}{2}\right)-\left\|\mathbf{C}^{k}\right\|},$
(137)

because the denominator is positive due to assumption (129) and relation
(131). Hence,

$\left\|\left(\mathbf{A}^{k}_{ani}\right)^{-1}\right\|\leq\frac{1}{\left(\frac{1}{\tau}+\frac{r}{2}\right)-A2^{2k}},$
(138)

and using assumption (129), we have

${\rm cond}\,\mathbf{A}^{k}\leq{\rm
cond}\,\mathbf{A}^{k}_{ani}\leq\frac{1+A\,2^{2k}\left(\frac{1}{\tau}+\frac{r}{2}\right)^{-1}}{1-A\,2^{2k}\left(\frac{1}{\tau}+\frac{r}{2}\right)^{-1}}\leq\frac{1+\gamma}{1-\gamma}.$
(139)

∎

Note that the condition number of $\mathbf{A}^{k}$ is not only uniformly
bounded, but also that the upper bound does not depend on dimension $d$, which
is the significant result because, for many methods, the condition number
grows exponentially with $d$. Note that condition (129) is not very
restrictive. Indeed, if spline wavelet basis of order $p$ is used, then by
Theorem 13 the optimal approximation order with respect to the $L^{2}$ norm is
$\mathcal{O}\left(2^{-pk}\right)$, if we neglect the term
$k^{\left(d-1\right)/2}$. The optimal $L^{2}$ error for the Crank–Nicolson
scheme is $\mathcal{O}\left(\tau^{2}\right)$. Therefore, the optimal choice is
$\tau=C2^{-pk/2}\leq C4^{-k}$ for $p\geq 4$.

Let $\mathbf{c}^{*}$ be the exact solution of equation (107) and
$\mathbf{c}^{k}$ be an approximate solution after $k$ conjugate gradient
iterations. It is well known that the speed of convergence of $\mathbf{c}^{k}$
depends on the condition number of the matrix $\mathbf{A}^{k}$. More
precisely,

$\left\|\mathbf{c}^{k}-\mathbf{c}^{*}\right\|_{A}\leq 2\left(\frac{\sqrt{{\rm
cond}\,\mathbf{A}^{k}}-1}{\sqrt{{\rm
cond}\,\mathbf{A}^{k}}+1}\right)^{k}\left\|\mathbf{c}^{0}-\mathbf{c}^{*}\right\|_{A},$

where $\left\|\cdot\right\|_{A}$ is given by
$\left\|\mathbf{x}\right\|_{A}=\sqrt{\mathbf{x}^{T}\mathbf{A}\mathbf{x}}$.

Hence, the uniformly bounded condition number of matrix $\mathbf{A}^{k}$
implies a uniformly bounded number of conjugate gradient iterations, if the
same stopping criterion for relative residuals is used for all $k$. Moreover,
the upper bound for a number of iterations does not increase with the
dimension $d$.

##  5 Numerical examples

In this section, numerical examples are presented to confirm the theoretical
results and illustrate the efficiency and reliability of the method. The
proposed wavelet scheme is used to evaluate prices of European-style options
on the geometric average of the underlying assets. The advantage of these
types of options is that after a transformation, the option contract can be
valuated using the standard one-dimensional Black–Scholes model [2, 29]. Then,
the resulting errors can be computed comparing numerical and analytic
solutions.

The geometric average of asset prices $S_{1},\ldots,S_{d}$ is defined as

$S=\left(\prod\limits_{i=1}^{d}S_{i}\right)^{1/d}.$ (140)

The payoff function representing the value of the option at maturity is for
the strike $K$ given by

$V\left(S_{1},\ldots,S_{d},0\right)=\max\left(K-S,0\right)$ (141)

for a put option and

$V\left(S_{1},\ldots,S_{d},0\right)=\max\left(S-K,0\right)$ (142)

for a call option. The graphs of payoff functions for $K=10$ and $d=2$ are
displayed in Figure 9.

Figure 9:  Payoff functions for a put option (left) and a call option (right)
on the geometric average of two assets with the strike $K=10$.

The price of a put option on the geometric average of $d$ assets can be
computed using the standard Black–Scholes formula

$V_{put}\left(S_{1},\ldots,S_{d},t\right)=Ke^{-rt}N\left(-d_{2}\right)-Se^{-\delta
t}N\left(-d_{1}\right),$

where $S$ is given by (140), $N$ is the standard normal density function

$N\left(x\right)=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{x}e^{-\frac{t^{2}}{2}}dt,$
(143)

and

$d_{1}=\frac{\ln\frac{S}{K}+\left(r-\delta+\frac{\sigma^{2}}{2}\right)t}{\sigma\sqrt{t}},\quad
d_{2}=\frac{\ln\frac{S}{K}+\left(r-\delta-\frac{\sigma^{2}}{2}\right)t}{\sigma\sqrt{t}}.$
(144)

The parameters $\sigma$ and $\delta$ are given by

$\sigma^{2}=\frac{1}{d^{2}}\sum_{i=1}^{d}\sum_{j=1}^{d}\rho_{ij}\sigma_{i}\sigma_{j},\quad\delta=\frac{1}{2d}\sum\limits_{i=1}^{d}\sigma_{i}^{2}-\frac{\sigma^{2}}{2}.$
(145)

The price of a call option can then be computed using the put-call parity

$V_{call}\left(S_{1},\ldots,S_{d},t\right)=V_{put}\left(S_{1},\ldots,S_{d},t\right)+Se^{-\delta
t}-Ke^{-rt}.$ (146)

Example. Options on the geometric average of $d$ assets

In this example, numerical results are presented for options on the geometric
average with parameters taken from [2]. The strike price is $K=10$ and the
option is maturing in $T=1$ year. The correlation coefficients are
$\rho_{i,j}=0.25$ for $i\neq j$, the risk-free interest rate for a given
horizon is $r=0.06$ p.a., and the corresponding volatilities are
$\sigma_{i}=0.2$. The parameters characterizing the domain for the computation
are set to $S^{min}=0.1$ and $S^{max}=50$. Figure 10 shows the resulting
functions representing prices of put and call options for these parameters
computed using the proposed method with constructed orthogonal cubic spline
wavelet basis. Because artificial boundary conditions are used, the plots are
shown only in the region $\left(1,30\right)^{2}$, to avoid the area near
$S_{1}=0$ and $S_{2}=0$.

Figure 10: _Functions representing values of a put option (left) and a call
option (right) on the geometric average of two assets._

Table 1 presents the results. Parameter $d$ denotes the number of underlying
assets, i.e., the spatial dimension of the problem. The parameter $k$ is a
number characterizing basis $\Psi^{k}$ constructed using sparse tensor product
approach, $N$ is the number of basis functions, and $M$ is the number of time
steps. The resulting system is solved by the conjugate gradient method without
preconditioning. The iterations stop when the relative residual is smaller
than $10^{-10}$ and the resulting number of iterations is denoted $it$. The
pointwise errors are computed for prices of underlying assets
$P_{1}=\left(K/2,\ldots,K/2\right)\in\mathbb{R}^{d}$,
$P_{2}=\left(K,\ldots,K\right)\in\mathbb{R}^{d}$, and
$P_{3}=\left(3K/2,\ldots,3K/2\right)\in\mathbb{R}^{d}$ at time to maturity
$t=T$,

$e_{k}\left(P_{i}\right)=\left|V_{k}\left(P_{i},T\right)-V\left(P_{i},T\right)\right|,\quad
i=1,2,3,$ (147)

where $V$ represents exact value computed by analytic formula (5) and $V_{k}$
stands for the solution computed using the proposed method with basis
$\Psi^{k}$.

Table 1:  Pointwise errors $e_{k}\left(P_{i}\right)$, $i=1,2,3$, and numbers of iterations for options on the geometric average for $d$ assets.  |  |  |  | put option | call option
---|---|---|---|---|---
| $k$ | $N$ | $M$ | $it$ | $e_{k}\left(P_{1}\right)$ | $e_{k}\left(P_{2}\right)$ | $it$ | $e_{k}\left(P_{2}\right)$ | $e_{k}\left(P_{3}\right)$
$d=2$ | 0 | 36 | 1 | 9 | 4.09e-1 | 3.89e-1 | 8 | 4.67e-1 | 1.98e-1
| 1 | 144 | 4 | 9 | 2.16e-3 | 5.29e-3 | 9 | 2.14e-2 | 6.74e-2
| 2 | 432 | 16 | 8 | 1.77e-3 | 1.89e-3 | 8 | 6.58e-3 | 6.69e-3
| 3 | 1 152 | 64 | 7 | 6.42e-4 | 9.11e-4 | 8 | 8.85e-4 | 1.95e-3
| 4 | 2 880 | 256 | 6 | 7.60e-5 | 7.31e-5 | 6 | 6.93e-5 | 1.11e-5
| 5 | 6 912 | 1 024 | 5 | 4.51e-6 | 3.52e-7 | 5 | 5.84e-6 | 6.56e-6
| 6 | 16 128 | 4 096 | 6 | 2.28e-7 | 2.81e-7 | 5 | 2.28e-7 | 7.83e-7
$d=3$ | 0 | 216 | 1 | 9 | 3.62e-1 | 1.83e-1 | 9 | 3.26e-1 | 3.04e-1
| 1 | 1 728 | 4 | 9 | 1.07e-3 | 1.12e-2 | 9 | 3.52e-2 | 6.81e-2
| 2 | 6 912 | 16 | 8 | 9.72e-4 | 6.45e-4 | 8 | 6.63e-3 | 6.77e-3
| 3 | 22 464 | 64 | 6 | 5.38e-5 | 1.13e-3 | 6 | 1.37e-3 | 2.36e-3
| 4 | 65 664 | 256 | 5 | 3.37e-6 | 1.98e-4 | 5 | 1.93e-4 | 1.56e-4
| 5 | 179 712 | 1 024 | 4 | 9.97e-6 | 8.25e-6 | 4 | 8.36e-6 | 3.62e-6
$d=4$ | 0 | 1 296 | 1 | 10 | 3.36e-1 | 1.69e-1 | 10 | 2.34e-1 | 4.18e-1
| 1 | 20 736 | 4 | 10 | 7.32e-3 | 8.71e-3 | 10 | 3.96e-2 | 7.63e-2
| 2 | 103 680 | 16 | 8 | 3.69e-4 | 7.60e-4 | 8 | 8.05e-3 | 5.12e-3
| 3 | 393 984 | 64 | 6 | 1.17e-4 | 2.57e-4 | 6 | 3.71e-4 | 2.19e-3
$d=5$ | 0 | 7 776 | 1 | 10 | 3.21e-1 | 1.58e-1 | 10 | 1.63e-1 | 5.05e-1
| 1 | 248 832 | 4 | 10 | 1.66e-3 | 6.16e-3 | 10 | 4.52e-2 | 8.45e-2
| 2 | 1 492 992 | 16 | 8 | 6.46e-4 | 2.65e-4 | 7 | 1.14e-2 | 2.53e-3

##  6 Conclusion

In this paper, we presented a wavelet-based method for option pricing under
the Black–Scholes model. Because the wavelet basis used is a cubic spline
wavelet basis, the method is high-order accurate with respect to the variables
representing prices, and the sparse tensor product structure of the basis
enables us to overcome the so-called “curse of dimensionality.” Due to the
orthogonality of the basis, the resulting stiffness matrices have uniformly
bounded condition numbers even without preconditioning. Moreover, the
condition number does not increase with the dimension. The numerical
experiments for European-style options on the geometric average confirm the
efficiency and applicability of the scheme; high-order convergence is
achieved; and the numbers of iterations are small, uniformly bounded and not
growing with the dimension. The described method can be used for pricing any
option that is represented by the presented Black–Scholes model such as basket
options, options on maximum or minimum of several underlying assets, and real
options.

## Acknowledgments

This work was supported by grant No. GA22-17028S funded by the Czech Science
Foundation and grant No. PURE-2020-4003 funded by the Technical University of
Liberec.

## References

  * [1] Y. Achdou, O. Pironneau, O., Computational Methods for Option Pricing, Society for Industrial and Applied Mathematics, Philadelphia, 2005.
  * [2] S. Berridge, J.M. Schumacher, An irregular grid method for high-dimensional free-boundary problems in finance, Futur. Gener. Comp. Syst. 20 (2004) 353–362.
  * [3] F. Black, M. Scholes, The pricing of options and corporate liabilities, J. Polit. Econ. 81 (1973) 637–659.
  * [4] A. Cohen, Numerical Analysis of Wavelet methods, Studies in Mathematics and its Applications 32, Elsevier, Amsterdam, 2003.
  * [5] D. Černá, V. Finěk, On a sparse representation of an n-dimensional Laplacian in wavelet coordinates, Result. Math. 69 (2016) 225–243.
  * [6] D. Černá, V. Finěk, Galerkin method with new quadratic spline wavelets for integral and integro-differential equations, J. Comput. Appl. Math. 363 (2020) 426–443.
  * [7] D. Černá, V. Finěk, Wavelet-Galerkin method for second-order integro-differential equations on product domains, in: H. Singh, H. Dutta, M.M. Cavalcanti, M.M. (Eds.), Topics in Integral and Integro-Differential Equations, Springer, Switzerland, 2021, pp. 1–40.
  * [8] W. Dahmen, Stability of multiscale transformations, J. Fourier Anal. Appl. 4 (1996) 341–362.
  * [9] W. Dahmen, B. Han, R.Q. Jia, A. Kunoth, Biorthogonal multiwavelets on the interval: cubic Hermite splines, Constr. Approx. 16 (2000) 221–259.
  * [10] W. Dahmen, Multiscale and wavelet methods for operator equations, Lecture Notes in Mathematics 1825, 2003, pp. 31–96.
  * [11] W. Dahmen, A. Kunoth, Multilevel preconditioning, Numer. Math. 63 (1992) 315–344.
  * [12] T.J. Dijkema, R. Stevenson, A sparse Laplacian in tensor product wavelet coordinates, Numer. Math. 115 (2010) 433–449.
  * [13] T.J. Dijkema, C. Schwab, R. Stevenson, An adaptive wavelet method for solving high-dimensional elliptic PDEs, Constr Approx 30 (2009) 423-–455.
  * [14] G.C. Donovan, J.S. Geronimo, D.P. Hardin, Intertwining multiresolution analyses and the construction of piecewise-polynomial wavelets, SIAM J. Math. Anal. 27 (1996) 1791–1815.
  * [15] G.C. Donovan, J.S. Geronimo, D.P. Hardin, Orthogonal polynomials and the construction of piecewise polynomial smooth wavelets, SIAM J. Math. Anal. 30 (1999) 1029–1056.
  * [16] R. Doostaki, M.M. Hosseini, Option pricing by the Legendre wavelets method, Comput Econ 59 (2022) 749–-773.
  * [17] L. Ortiz-Gracia, C.W. Oosterlee, Robust pricing of European options with wavelets and the characteristic function. SIAM J. Sci. Comput., 35 (2013) B1055–-B1084.
  * [18] L. Ortiz-Gracia, C.W. Oosterlee, A highly efficient Shannon wavelet inverse Fourier technique for pricing European options. SIAM J. Sci. Comput., 38 (2016) B18–-B143.
  * [19] M. Griebel, P. Oswald, Tensor product type subspace splittings and multilevel iterative methods for anisotropic problems, Adv. Comput. Math. 4 (1995) 171–206.
  * [20] M. Griebel, P. Oswald and T. Schiekofer, Sparse grids for boundary integral equations, Numer. Math. 83 (1999) 279–312.
  * [21] K.E. Gustafson, D.K.M. Rao, Numerical Range, The Field of Values of Linear Operators and Matrices, Springer-Verlag, Berlin, 1997.
  * [22] N. Hilber, O. Reichmann, C. Schwab, C. Winter, Computational Methods for Quantitative Finance, Springer, Berlin, 2013.
  * [23] J. Hu, S. Gan, High order method for Black–Scholes PDE, Comput. Math. Appl. 75 (2018) 2259–2270.
  * [24] J. Jo, Y. Kim, Comparison of numerical schemes on multi-dimensional Black–Scholes equations, Bull. Korean Math. Soc. 50 (2013) 2035–2051.
  * [25] M.K. Kadalbajoo, L.P. Tripathi, A. Kumar, A cubic B-spline collocation method for a numerical solution of the generalized Black–Scholes equation Math, Comput. Modelling 55 (2012) 1483–1505.
  * [26] S.R. Kestler, On the Adaptive Tensor Product Wavelet Galerkin Method with Applications in Finance, PhD thesis, Ulm, 2013.
  * [27] J. Kim, T. Kim, J. Jo, Y. Choi, S. Lee, H. Hwang, M. Yoo, D. Jeong, A practical finite difference method for the three-dimensional Black–Scholes equation, European J. Oper. Res. 252 (2016) 183–190.
  * [28] C. Lee, S. Kwak, Y. Hwang, J. Kim, Accurate and efficient finite difference method for the Black–Scholes model with no far-field boundary conditions, Computat. Econ. (2022).
  * [29] C.C.W. Leentvaar, Pricing multi-asset options with sparse grids, Ph.D. thesis, Delft University of Technology, Delft, 2008.
  * [30] J. Lyu, E. Park, S. Kim, W. Lee, C. Lee, S. Yoon, J. Park, J.J. Kim, Optimal non-uniform finite difference grids for the Black–Scholes equations, Math. Comput. Simul. 182 (2021) 690-–704.
  * [31] D. Marquardt, An algorithm for least-squares estimation of nonlinear parameters, SIAM J. Appl. Math. 11 (1963) 431-–441.
  * [32] S. Milovanović, L. von Sydow, A high order method for pricing of financial derivatives using radial basis function generated finite differences, Math. Comput. Simul. 174 (2020) 205–217.
  * [33] R. Mohammadi, Quintic B-spline collocation approach for solving generalized Black–Scholes equation governing option pricing, Comput. Math. Appl. 69 (2015) 777–797.
  * [34] T. von Petersdorff, C. Schwab, Numerical solution of parabolic equations in high dimensions, ESAIM-Math. Model. Numer. Anal. 38 (2004) 93-–127.
  * [35] R. Rannacher. Finite Element solution of diffusion problems with irregular data. Numerische Mathematik, 43 (1984) 309–327.
  * [36] M. Rometsch, A Wavelet Tour of Option Pricing, PhD thesis, Ulm University, 2010.
  * [37] P. Roul, V.M.K. Prasad Goura, A sixth order numerical method and its convergence for generalized Black–Scholes PDE, J. Comput. Appl. Math. 377 (2022) 112881.
  * [38] N. Reich, Wavelet Compression of Anisotropic Integrodifferential Operators on Sparse Tensor Product Spaces, PhD thesis, ETH Zürich, 2008.
  * [39] A.J. Rupp, High Dimensional Wavelet Methods for Structured Financial Products, PhD thesis, Ulm University, 2013.

##  Appendix A  Coefficients of constructed generators

Table 2:  Coefficients of the scaling generators corresponding to the terms $x^{3}$, $x^{2}$, $x$, and $1$ on $\left[a,b\right]$.  | a | b |  | coefficient |  | coefficient
---|---|---|---|---|---|---
$\phi_{1}$ | 0 | 0.5 | $x^{3}$ | -26.25320493255574 | $x^{2}$ | 19.68990369941681
|  |  | $x$ | 0 | $1$ | 0
| 0.5 | 1 | $x^{3}$ | 26.25320493255574 | $x^{2}$ | -59.06971109825042
|  |  | $x$ | 39.37980739883362 | $1$ | -6.563301233138936
$\phi_{2}$ | 0 | 0.5 | $x^{3}$ | 81.97560612767679 | $x^{2}$ | -40.98780306383839
|  |  | $x$ | 0 | $1$ | 0
| 0.5 | 1 | $x^{3}$ | 81.97560612767679 | $x^{2}$ | -204.9390153191920
|  |  | $x$ | 163.9512122553536 | $1$ | -40.98780306383839
$\phi_{3}$ | 0 | 0.25 | $x^{3}$ | 53.41388035788556 | $x^{2}$ | -35.17323440989752
|  |  | $x$ | 0 | $1$ | 0
| 0.25 | 0.5 | $x^{3}$ | -527.8691353655910 | $x^{2}$ | 604.3590514426924
|  |  | $x$ | -210.7755774781431 | $1$ | 21.80567362442823
| 0.5 | 0.75 | $x^{3}$ | 167.3619498129961 | $x^{2}$ | -332.4004596541814
|  |  | $x$ | 204.5606197347904 | $1$ | -38.57643285514344
| 0.75 | 1 | $x^{3}$ | -381.1201871105696 | $x^{2}$ | 1023.674050039347
|  |  | $x$ | -903.9875387469848 | $1$ | 261.4336758182076
$\phi_{4}$ | 0 | 0.25 | $x^{3}$ | 26.62362537753055 | $x^{2}$ | 15.35514686530018
|  |  | $x$ | 0 | $1$ | 0
| 0.25 | 0.5 | $x^{3}$ | 1269.300093977873 | $x^{2}$ | -1351.930373703732
|  |  | $x$ | 450.6409224219518 | $1$ | -46.62170539180379
| 0.5 | 0.75 | $x^{3}$ | 567.2864791328406 | $x^{2}$ | -1200.276995000695
|  |  | $x$ | 825.4977548526893 | $1$ | -184.2117644273027
| 0.75 | 1 | $x^{3}$ | -399.7440135001654 | $x^{2}$ | 1084.571643661303
|  |  | $x$ | -969.9112468221101 | $1$ | 285.0836166609723
$\phi_{5}$ | 0 | 0.25 | $x^{3}$ | 10.19187443147913 | $x^{2}$ | 28.81827269560206
|  |  | $x$ | 27.06092209676673 | $1$ | 8.434523832643798
| 0.25 | 0.5 | $x^{3}$ | 122.4194193918863 | $x^{2}$ | 242.0020294473649
|  |  | $x$ | 157.4525751037239 | $1$ | 33.65839594516685
| 0.5 | 0.75 | $x^{3}$ | 179.9690473081142 | $x^{2}$ | 185.5152420484449
|  |  | $x$ | 57.80356676763296 | $1$ | 5.149292116379875
| 0.75 | 1 | $x^{3}$ | -365.9470852800144 | $x^{2}$ | -134.8104412073693
|  |  | $x$ | 0 | $1$ | 2.188816056270510
| 0 | 0.25 | $x^{3}$ | 79.26943625448920 | $x^{2}$ | -58.48104390125948
|  |  | $x$ | 0 | $1$ | 2.188816056270510
| 0.25 | 0.5 | $x^{3}$ | -354.3621713021780 | $x^{2}$ | 418.6216724392985
|  |  | $x$ | -157.2454317534039 | $1$ | 18.45674809140953
| 0.5 | 0.75 | $x^{3}$ | -50.59040574008613 | $x^{2}$ | 103.1911403389635
|  |  | $x$ | -69.64372382463778 | $1$ | 15.54205645684875
| 0.75 | 1 | $x^{3}$ | -2.187234289273272 | $x^{2}$ | 6.198136014518121
|  |  | $x$ | -5.834569161216426 | $1$ | 1.823667435971577
$\phi_{6}$ | 0 | 0.25 | $x^{3}$ | -15.40091533093236 | $x^{2}$ | -43.54717325197819
|  |  | $x$ | -40.89160051115930 | $1$ | -12.74534259011347
| 0.25 | 0.5 | $x^{3}$ | -184.9074696818534 | $x^{2}$ | -365.5196063064310
|  |  | $x$ | -237.8079396256591 | $1$ | -50.83368094965349
| 0.5 | 0.75 | $x^{3}$ | -265.5803178284557 | $x^{2}$ | -271.7997826856661
|  |  | $x$ | -83.58347989494250 | $1$ | -7.235513007811689
| 0.75 | 1 | $x^{3}$ | 989.1663297899967 | $x^{2}$ | 460.3081592797291
|  |  | $x$ | 47.20549465929524 | $1$ | -0.689599373051132
| 0 | 0.25 | $x^{3}$ | 273.6635501100016 | $x^{2}$ | -249.4001766673951
|  |  | $x$ | 47.20549465929524 | $1$ | -0.689599373051132
| 0.25 | 0.5 | $x^{3}$ | -610.5427721702383 | $x^{2}$ | 723.4466794016099
|  |  | $x$ | -273.4292479476622 | $1$ | 32.48188156000418
| 0.5 | 0.75 | $x^{3}$ | -98.42902348889925 | $x^{2}$ | 199.2567948682677
|  |  | $x$ | -133.3246749253244 | $1$ | 29.46284759700341
| 0.75 | 1 | $x^{3}$ | -5.217694467034758 | $x^{2}$ | 14.77336749684687
|  |  | $x$ | -13.89365159258947 | $1$ | 4.337978562777356
$\phi_{L}$ | 0 | 0.25 | $x^{3}$ | -328.4179396059963 | $x^{2}$ | 294.5324087893769
|  |  | $x$ | -51.91281405451471 | $1$ | 0
| 0.25 | 0.5 | $x^{3}$ | 794.2029204485535 | $x^{2}$ | -940.6297684522534
|  |  | $x$ | 355.1768633060724 | $1$ | -42.11573420089721
| 0.5 | 0.75 | $x^{3}$ | 125.7725578721856 | $x^{2}$ | -254.8795911518400
|  |  | $x$ | 170.7494579379349 | $1$ | -37.78578051988584
| 0.75 | 1 | $x^{3}$ | 6.495819573712430 | $x^{2}$ | -18.39405161997127
|  |  | $x$ | 17.30064451880524 | $1$ | -5.402412472546407
$\phi_{R}$ | 0 | 0.25 | $x^{3}$ | 23.35622528154187 | $x^{2}$ | -4.027321052314352
|  |  | $x$ | 0 | $1$ | 0
| 0.25 | 0.5 | $x^{3}$ | 280.3889896276431 | $x^{2}$ | -286.9065184334287
|  |  | $x$ | 93.24595537566320 | $1$ | -9.647675950503985
| 0.5 | 0.75 | $x^{3}$ | 400.2203835185070 | $x^{2}$ | -791.8720820637919
|  |  | $x$ | 508.3379735878785 | $1$ | -105.9312183853788
| 0.75 | 1 | $x^{3}$ | -1674.366086197562 | $x^{2}$ | 4222.513185710944
|  |  | $x$ | -3512.375260428360 | $1$ | 964.2281609149772
Table 3:  Coefficients of the wavelet generators corresponding to the terms $x^{3}$, $x^{2}$, $x$, and $1$ on $\left[a,b\right]$.  | a | b |  | coefficient |  | coefficient
---|---|---|---|---|---|---
$\psi_{1}$ | 0 | 0.125 | $x^{3}$ | 512.7201177918754 | $x^{2}$ | -43.88706601604656
|  |  | $x$ | 0 | $1$ | 0
| 0.125 | 0.25 | $x^{3}$ | 6423.533783310561 | $x^{2}$ | -3295.876989270898
|  |  | $x$ | 535.9280902425246 | $1$ | -27.72322666992469
| 0.25 | 0.375 | $x^{3}$ | 6249.063735441256 | $x^{2}$ | -6401.764292545093
|  |  | $x$ | 2121.584875855116 | $1$ | -227.2933721204776
| 0.375 | 0.5 | $x^{3}$ | -10261.35309765599 | $x^{2}$ | 13461.12463292785
|  |  | $x$ | -5810.249716786693 | $1$ | 824.5923576584304
| 0.5 | 0.625 | $x^{3}$ | -291.9822638070393 | $x^{2}$ | 622.1814317187671
|  |  | $x$ | -448.3346409643185 | $1$ | 107.1992658183961
| 0.625 | 0.75 | $x^{3}$ | -4060.994255437327 | $x^{2}$ | 8349.124694297900
|  |  | $x$ | -5690.202791496491 | $1$ | 1285.198591225143
| 0.75 | 0.875 | $x^{3}$ | -528.1868779043331 | $x^{2}$ | 1361.218410522837
|  |  | $x$ | -1169.955815420823 | $1$ | 335.3075313951340
| 0.875 | 1 | $x^{3}$ | -130.6441212212436 | $x^{2}$ | 371.4479197042459
|  |  | $x$ | -350.9634757447609 | $1$ | 110.1596772617586
$\psi_{2}$ | 0 | 0.125 | $x^{3}$ | 181.2722522251391 | $x^{2}$ | -14.15090760709985
|  |  | $x$ | 0 | $1$ | 0
| 0.125 | 0.25 | $x^{3}$ | 2419.148109670498 | $x^{2}$ | -1245.470104784877
|  |  | $x$ | 202.9293684766930 | $1$ | -10.49765988775683
| 0.25 | 0.375 | $x^{3}$ | 2576.675107217909 | $x^{2}$ | -2596.000954579617
|  |  | $x$ | 848.6584813339236 | $1$ | -89.98311932657152
| 0.375 | 0.5 | $x^{3}$ | -8471.127674575379 | $x^{2}$ | 10735.23758057596
|  |  | $x$ | -4488.978621463715 | $1$ | 619.5243500374206
| 0.5 | 0.625 | $x^{3}$ | -3086.030528892990 | $x^{2}$ | 6003.939215143344
|  |  | $x$ | -3796.503115292892 | $1$ | 782.9740450998639
| 0.625 | 0.75 | $x^{3}$ | 12658.02295546319 | $x^{2}$ | -26273.32787303766
|  |  | $x$ | 18100.01806795347 | $1$ | -4137.807295812551
| 0.75 | 0.875 | $x^{3}$ | 1896.924386457221 | $x^{2}$ | -4813.188115086536
|  |  | $x$ | 4069.162266224350 | $1$ | -1146.155599563828
| 0.875 | 1 | $x^{3}$ | 253.9044233453815 | $x^{2}$ | -726.4312505262515
|  |  | $x$ | 691.1492310163587 | $1$ | -218.6224038354886
$\psi_{3}$ | -1 | -0.75 | $x^{3}$ | -12.77748177259837 | $x^{2}$ | -36.12926000986063
|  |  | $x$ | -33.92607470192613 | $1$ | -10.57429646466388
| -0.75 | -0.5 | $x^{3}$ | -153.4635771854685 | $x^{2}$ | -303.3692884353265
|  |  | $x$ | -197.3783313309066 | $1$ | -42.19291944937926
| -0.5 | -0.375 | $x^{3}$ | -59.47798764121307 | $x^{2}$ | 2.242882836829032
|  |  | $x$ | 37.74464778305737 | $1$ | 10.71372598259577
| -0.375 | -0.25 | $x^{3}$ | 1758.865605322497 | $x^{2}$ | 1729.275753199708
|  |  | $x$ | 565.9055972736512 | $1$ | 61.79929755698429
| -0.25 | -0.125 | $x^{3}$ | 3444.276587955715 | $x^{2}$ | 1712.829251641442
|  |  | $x$ | 241.6677872507898 | $1$ | 8.102298002304606
| -0.125 | 0 | $x^{3}$ | -5400.816676245629 | $x^{2}$ | -882.1760744370164
|  |  | $x$ | 7.530202490613302 | $1$ | 2.106485345615197
| 0 | 0.125 | $x^{3}$ | 1247.729690468656 | $x^{2}$ | -454.3272932130656
|  |  | $x$ | 7.530202490613302 | $1$ | 2.106485345615197
| 0.125 | 0.25 | $x^{3}$ | -5778.096772948966 | $x^{2}$ | 3410.751764470566
|  |  | $x$ | -629.4039464575936 | $1$ | 35.05371099919686
| 0.25 | 0.375 | $x^{3}$ | -558.2316506051182 | $x^{2}$ | 527.4560525803158
|  |  | $x$ | -166.4808009519398 | $1$ | 17.96851407930143
| 0.375 | 0.5 | $x^{3}$ | 226.0107269480922 | $x^{2}$ | -258.2984965273795
|  |  | $x$ | 91.98285784857104 | $1$ | -9.815156131403054
| 0.5 | 0.75 | $x^{3}$ | 35.52867654721048 | $x^{2}$ | -72.71045185132442
|  |  | $x$ | 49.25635097317727 | $1$ | -11.03865756260972
| 0.75 | 1 | $x^{3}$ | 1.382563578715178 | $x^{2}$ | -3.919859133343993
|  |  | $x$ | 3.692027530542450 | $1$ | -1.15473197591364
$\psi_{4}$ | -1 | -0.75 | $x^{3}$ | 6.603401663338006 | $x^{2}$ | 18.67158529263201
|  |  | $x$ | 17.53296559525000 | $1$ | 5.464781965955998
| -0.75 | -0.5 | $x^{3}$ | 79.28543622294657 | $x^{2}$ | 156.7295257413204
|  |  | $x$ | 101.9689429489431 | $1$ | 21.79690680867349
| -0.5 | -0.375 | $x^{3}$ | -82.45595592132698 | $x^{2}$ | -161.0461033269963
|  |  | $x$ | -94.50064201116839 | $1$ | -17.21165242233730
| -0.375 | -0.25 | $x^{3}$ | -2246.183853666682 | $x^{2}$ | -2216.018451884065
|  |  | $x$ | -722.9071965676484 | $1$ | -77.98646222284469
| -0.25 | -0.125 | $x^{3}$ | -3800.606529757320 | $x^{2}$ | -1924.442967268418
|  |  | $x$ | -285.6652024928302 | $1$ | -11.18728580653431
| -0.125 | 0 | $x^{3}$ | 11937.30882839563 | $x^{2}$ | 2667.744109956829
|  |  | $x$ | 124.6667844000623 | $1$ | -0.910594592674766
| 0 | 0.125 | $x^{3}$ | 3136.965997055351 | $x^{2}$ | -1389.100592703720
|  |  | $x$ | 124.6667844000623 | $1$ | -0.910594592674766
| 0.125 | 0.25 | $x^{3}$ | -7152.501310626270 | $x^{2}$ | 4271.385867623253
|  |  | $x$ | -808.1360506341049 | $1$ | 47.34127467930285
| 0.25 | 0.375 | $x^{3}$ | -927.7525957207491 | $x^{2}$ | 901.3927836102403
|  |  | $x$ | -290.2798926723839 | $1$ | 31.24010426928711
| 0.375 | 0.5 | $x^{3}$ | 154.3100993503755 | $x^{2}$ | -168.8298334848166
|  |  | $x$ | 55.89187066577802 | $1$ | -5.136151388922590
| 0.5 | 0.75 | $x^{3}$ | 35.12445303753813 | $x^{2}$ | -71.05734845883623
|  |  | $x$ | 47.50862037442566 | $1$ | -10.48944171063683
| 0.75 | 1 | $x^{3}$ | 1.892186966470384 | $x^{2}$ | -5.357215161018598
|  |  | $x$ | 5.037869422626043 | $1$ | -1.572841228077830
$\psi_{5}$ | -1 | -0.75 | $x^{3}$ | -0.001049843371405 | $x^{2}$ | -0.003037046895837
|  |  | $x$ | -0.002924563677460 | $1$ | -0.000937360153027
| -0.75 | -0.5 | $x^{3}$ | -0.135429794979682 | $x^{2}$ | -0.283795160078696
|  |  | $x$ | -0.197295565112758 | $1$ | -0.045480714648891
| -0.5 | -0.375 | $x^{3}$ | 64.90438914754031 | $x^{2}$ | 92.49504460139021
|  |  | $x$ | 43.80167998946610 | $1$ | 6.889274490088309
| -0.375 | -0.25 | $x^{3}$ | 887.7169530732866 | $x^{2}$ | 874.0398882738802
|  |  | $x$ | 282.8362623376594 | $1$ | 30.01300552998866
| -0.25 | -0.125 | $x^{3}$ | 689.7712273209479 | $x^{2}$ | 312.8547139314197
|  |  | $x$ | 39.35849874499266 | $1$ | 1.124736063345467
| -0.125 | 0 | $x^{3}$ | -677.4429634928239 | $x^{2}$ | -100.9319418261331
|  |  | $x$ | 0 | $1$ | 0
| 0 | 0.125 | $x^{3}$ | 222.6292417447635 | $x^{2}$ | 78.22658535063456
|  |  | $x$ | 0 | $1$ | 0
| 0.125 | 0.25 | $x^{3}$ | 13436.75267042187 | $x^{2}$ | -7191.190145000426
|  |  | $x$ | 1197.942146868526 | $1$ | -61.96696676846536
| 0.25 | 0.375 | $x^{3}$ | 3309.658731654356 | $x^{2}$ | -3681.684248151221
|  |  | $x$ | 1342.019311962833 | $1$ | -159.0945338018750
| 0.375 | 0.5 | $x^{3}$ | 1623.482149916625 | $x^{2}$ | -1948.015296157701
|  |  | $x$ | 753.1233433882986 | $1$ | -93.13627378293758
| 0.5 | 0.75 | $x^{3}$ | 334.8296335785243 | $x^{2}$ | -668.4089674911632
|  |  | $x$ | 440.0064019753428 | $1$ | -95.39782070083224
| 0.75 | 1 | $x^{3}$ | 23.73541224062620 | $x^{2}$ | -67.14137193635135
|  |  | $x$ | 63.07650715082408 | $1$ | -19.67054745509894
$\psi_{6}$ | -1 | -0.75 | $x^{3}$ | 29.39534138606089 | $x^{2}$ | 83.11643728617156
|  |  | $x$ | 78.04685041416045 | $1$ | 24.32575451404978
| -0.75 | -0.5 | $x^{3}$ | 351.2049102801593 | $x^{2}$ | 694.0258200020245
|  |  | $x$ | 451.3572769791491 | $1$ | 96.43545853732169
| -0.5 | -0.375 | $x^{3}$ | 1275.659605723797 | $x^{2}$ | 1615.831201469838
|  |  | $x$ | 679.8216368642327 | $1$ | 95.77313004336493
| -0.375 | -0.25 | $x^{3}$ | 11255.59466789715 | $x^{2}$ | 11095.38855970805
|  |  | $x$ | 3579.204551188512 | $1$ | 376.2646074570184
| -0.25 | -0.125 | $x^{3}$ | 5327.416704802243 | $x^{2}$ | 2249.311342113261
|  |  | $x$ | 267.6993104714105 | $1$ | 8.640342704059652
| -0.125 | 0 | $x^{3}$ | -2957.628236410947 | $x^{2}$ | -374.9318262498545
|  |  | $x$ | 0 | $1$ | 0
| 0 | 0.125 | $x^{3}$ | -46.42927388024598 | $x^{2}$ | 4.175934965984637
|  |  | $x$ | 0 | $1$ | 0
| 0.125 | 0.25 | $x^{3}$ | -788.4860165790724 | $x^{2}$ | 412.3997958564289
|  |  | $x$ | -67.27205540860357 | $1$ | 3.479838675245901
| 0.25 | 0.375 | $x^{3}$ | -231.1352598754021 | $x^{2}$ | 260.5391696188435
|  |  | $x$ | -95.84500917174907 | $1$ | 11.40576068238651
| 0.375 | 0.5 | $x^{3}$ | -133.3790261589963 | $x^{2}$ | 159.4675576509780
|  |  | $x$ | -61.28221129495857 | $1$ | 7.502793024182585
| 0.5 | 0.75 | $x^{3}$ | -26.46366034181517 | $x^{2}$ | 52.95015549352690
|  |  | $x$ | -34.95133350039340 | $1$ | 7.602283939115114
| 0.75 | 1 | $x^{3}$ | -1.798568145773455 | $x^{2}$ | 5.088295556487342
|  |  | $x$ | -4.780886675654469 | $1$ | 1.491159264940578
$\psi_{L_{1}}$ | 0 | 0.125 | $x^{3}$ | -1425.503021368171 | $x^{2}$ | 772.9341296224186
|  |  | $x$ | -56.55471808022630 | $1$ | 0
| 0.125 | 0.25 | $x^{3}$ | 16500.24187884111 | $x^{2}$ | -9088.459925693413
|  |  | $x$ | 1568.524503551422 | $1$ | -84.06184109786738
| 0.25 | 0.375 | $x^{3}$ | 3527.021590160490 | $x^{2}$ | -3847.552266045169
|  |  | $x$ | 1380.549477854916 | $1$ | -161.9182463911214
| 0.375 | 0.5 | $x^{3}$ | 1363.434660038630 | $x^{2}$ | -1645.715891236687
|  |  | $x$ | 641.9354328937131 | $1$ | -80.47581521996845
| 0.5 | 0.75 | $x^{3}$ | 281.9469114612293 | $x^{2}$ | -562.0738058307491
|  |  | $x$ | 369.4091589208262 | $1$ | -79.93723101283425
| 0.75 | 1 | $x^{3}$ | 20.47469809743168 | $x^{2}$ | -57.91380903680889
|  |  | $x$ | 54.40352378132468 | $1$ | -16.96441284194740
$\psi_{L_{2}}$ | 0 | 0.125 | $x^{3}$ | 3578.516685827048 | $x^{2}$ | -1483.144649002284
|  |  | $x$ | 125.1644746095171 | $1$ | 0
| 0.125 | 0.25 | $x^{3}$ | -3221.524685507020 | $x^{2}$ | 2257.721533142182
|  |  | $x$ | -491.3001316453148 | $1$ | 31.88837248923357
| 0.25 | 0.375 | $x^{3}$ | 253.7811404150473 | $x^{2}$ | -437.2143373073622
|  |  | $x$ | 204.5479612190696 | $1$ | -27.94181235379833
| 0.375 | 0.5 | $x^{3}$ | 859.2458117883292 | $x^{2}$ | -1011.673202197156
|  |  | $x$ | 379.9617016508116 | $1$ | -44.86748817002474
| 0.5 | 0.75 | $x^{3}$ | 179.5608184460906 | $x^{2}$ | -359.6822707980349
|  |  | $x$ | 237.7345152583694 | $1$ | -51.79100365580411
| 0.75 | 1 | $x^{3}$ | 11.94553851645503 | $x^{2}$ | -33.79701898388310
|  |  | $x$ | 31.75742241840201 | $1$ | -9.905941950973904
$\psi_{R1}$ | 0 | 0.25 | $x^{3}$ | 7.391479105551183 | $x^{2}$ | -1.275014819251857
|  |  | $x$ | 0 | $1$ | 0
| 0.25 | 0.5 | $x^{3}$ | 87.83823678183982 | $x^{2}$ | -89.99728641434138
|  |  | $x$ | 29.27736873324082 | $1$ | -3.031180797309012
| 0.5 | 0.625 | $x^{3}$ | 780.1233069697320 | $x^{2}$ | -1281.837786515931
|  |  | $x$ | 701.9040661939105 | $1$ | -127.9200382757332
| 0.625 | 0.75 | $x^{3}$ | 8761.559020997937 | $x^{2}$ | -17645.25953421656
|  |  | $x$ | 11802.93627344290 | $1$ | -2622.696251230958
| 0.75 | 0.875 | $x^{3}$ | 7584.571826637600 | $x^{2}$ | -19150.38737314525
|  |  | $x$ | 16046.79392231900 | $1$ | -4462.413605869881
| 0.875 | 1 | $x^{3}$ | -18377.05522246807 | $x^{2}$ | 51055.94173564137
|  |  | $x$ | -47183.66988964301 | $1$ | 14504.78337646970
$\psi_{R2}$ | 0 | 0.25 | $x^{3}$ | 28.49333825597376 | $x^{2}$ | -4.913915083023494
|  |  | $x$ | 0 | $1$ | 0
| 0.25 | 0.5 | $x^{3}$ | 340.6257513882141 | $x^{2}$ | -348.7311999647503
|  |  | $x$ | 113.3838149785686 | $1$ | -11.73444239472556
| 0.5 | 0.625 | $x^{3}$ | 1001.089085127504 | $x^{2}$ | -1771.030876846587
|  |  | $x$ | 1040.335991555937 | $1$ | -202.1935281803620
| 0.625 | 0.75 | $x^{3}$ | 7909.591723311074 | $x^{2}$ | -15934.19401753424
|  |  | $x$ | 10648.38838829413 | $1$ | -2361.386826210905
| 0.75 | 0.875 | $x^{3}$ | 1810.118176484035 | $x^{2}$ | -4899.457987996012
|  |  | $x$ | 4389.145954257418 | $1$ | -1300.778614730966
| 0.875 | 1 | $x^{3}$ | 7611.345386975026 | $x^{2}$ | -20790.59460977267
|  |  | $x$ | 18873.94129327007 | $1$ | -5694.692070472430