Unified analysis of discontinuous Galerkin approximations of flows in fractured porous media on polygonal and polyhedral grids

Paola F. Antonietti; Chiara Facciolà; Marco Verani; Paola F. Antonietti; Chiara Facciolà; Marco Verani

doi:10.3934/mine.2020017

Mathematics in Engineering

2020, Volume 2, Issue 2: 340-385. doi: 10.3934/mine.2020017

Previous Article Next Article

Research article

Unified analysis of discontinuous Galerkin approximations of flows in fractured porous media on polygonal and polyhedral grids

MOX- Laboratory for Modeling and Scientific Computing, Department of Mathematics, Politecnico of Milan, Piazza Leonardo da Vinci 32, 20133 Milano, Italy

Received: 21 March 2019 Accepted: 08 January 2020 Published: 19 February 2020

We propose a unified formulation based on discontinuous Galerkin methods on polygonal/polyhedral grids for the simulation of flows in fractured porous media. We adopt a model for single-phase flows where the fracture is modeled as a (d-1)-dimensional interface in a d-dimensional bulk domain, and model the flow in the porous medium and in the fracture by means of the Darcy's law. The two problems are then coupled through physically consistent conditions. We focus on the numerical approximation of the coupled bulk-fracture problem and present and analyze, in an unified setting, all the possible combinations of primal-primal, mixed-primal, primal-mixed and mixed-mixed formulations for the bulk and fracture problems, respectively. For all the possible combinations, we prove their well-posedness and derive a priori hp-version error estimates in a suitable (mesh-dependent) energy norm. Finally, preliminary numerical experiments assess the theoretical error estimates and accuracy of the proposed formulations.

Keywords:

Citation: Paola F. Antonietti, Chiara Facciolà, Marco Verani. Unified analysis of discontinuous Galerkin approximations of flows in fractured porous media on polygonal and polyhedral grids[J]. Mathematics in Engineering, 2020, 2(2): 340-385. doi: 10.3934/mine.2020017

Related Papers:

[1]	Teh Raihana Nazirah Roslan, Sharmila Karim, Siti Zulaiha Ibrahim, Ali Fareed Jameel, Zainor Ridzuan Yahya . Stochastic pricing formulation for hybrid equity warrants. AIMS Mathematics, 2022, 7(1): 398-424. doi: 10.3934/math.2022027
[2]	Guiwen Lv, Ping Xu, Yanxue Zhang . Pricing of vulnerable options based on an uncertain CIR interest rate model. AIMS Mathematics, 2023, 8(5): 11113-11130. doi: 10.3934/math.2023563
[3]	Din Prathumwan, Thipsuda Khonwai, Narisara Phoochalong, Inthira Chaiya, Kamonchat Trachoo . An improved approximate method for solving two-dimensional time-fractional-order Black-Scholes model: a finite difference approach. AIMS Mathematics, 2024, 9(7): 17205-17233. doi: 10.3934/math.2024836
[4]	Min-Ku Lee, Jeong-Hoon Kim . Closed-form approximate solutions for stop-loss and Russian options with multiscale stochastic volatility. AIMS Mathematics, 2023, 8(10): 25164-25194. doi: 10.3934/math.20231284
[5]	Ashish Awasthi, Riyasudheen TK . An accurate solution for the generalized Black-Scholes equations governing option pricing. AIMS Mathematics, 2020, 5(3): 2226-2243. doi: 10.3934/math.2020147
[6]	Xiuxian Chen, Zhongyang Sun, Dan Zhu . Mean-variance investment and risk control strategies for a dynamic contagion process with diffusion. AIMS Mathematics, 2024, 9(11): 33062-33086. doi: 10.3934/math.20241580
[7]	Kazem Nouri, Milad Fahimi, Leila Torkzadeh, Dumitru Baleanu . Numerical method for pricing discretely monitored double barrier option by orthogonal projection method. AIMS Mathematics, 2021, 6(6): 5750-5761. doi: 10.3934/math.2021339
[8]	Chao Yue, Chuanhe Shen . Fractal barrier option pricing under sub-mixed fractional Brownian motion with jump processes. AIMS Mathematics, 2024, 9(11): 31010-31029. doi: 10.3934/math.20241496
[9]	Xun Lu, Wei Shi, Changhao Yang, Fan Yang . Exponential integral method for European option pricing. AIMS Mathematics, 2025, 10(6): 12900-12918. doi: 10.3934/math.2025580
[10]	Kanagaraj Muthuselvan, Baskar Sundaravadivoo, Kottakkaran Sooppy Nisar, Fahad Sameer Alshammari . New technique for solving the numerical computation of neutral fractional functional integro-differential equation based on the Legendre wavelet method. AIMS Mathematics, 2024, 9(6): 14288-14309. doi: 10.3934/math.2024694

Abstract

1. Introduction

Modeling derivative products in mathematical finance usually starts with a system of stochastic differential equation that correspond to state variables like stock, interest rate and volatility. For financial analysts, economic auditors and risk managers in financial organizations, securities analysis is very important, thus in last decades they have tried to design new models that help to change markets react and solve some problems in the market. Modeling challenges has created a common link between applied mathematicians and financial researchers as the production of knowledge in this area ^[4,15].

As an important model in the original Black–Scholes option valuation theory the volatility and the interest rate are both assumed constant. A popular expansion of the Black–Scholes model has been proposed by Heston. This expansion allows also for correlation between the volatility and the asset price, so that one can capture skewness effects ^[16]. In addition to the volatility the values of many practically important securities are sensitive to interest rates, necessitating the construction of interest rate models ^[17].

Consider an agent who at each time has a portfolio valued at $X(t).$ This portfolio invests in a money market account paying a constant rate of interest $r,$ and in a stock modeled by the geometric Brownian motion,

$\begin{equation} dS(t) = \alpha S(t)dt+\delta S(t)dW(t), \end{equation}$

(1.1)

Suppose at each time $t,$ the investor holds $\Delta(t)$ shares of stock. The position $\Delta(t)$ can be random but must be adapted to the filtration associated with the Brownian motion $W(t), t\geq0.$ The remainder of the portfolio value, $X(t)-\Delta(t) S(t),$ is invested in the money market account ^[12].

In addition to, we consider the following model which presents the system of three SDEs which they have composed from different assets called hybrid model. The asset price model in following can be viewed as an extension of the popular Heston stochastic volatility model where the interest rate is not constant but also follows a stochastic process, described here by the Hull–White model.

$\begin{equation} dS_\tau = R_\tau S_\tau d_\tau+S_\tau\sqrt{V_\tau} dW_{\tau}^1, S(0) \gt 0, \end{equation}$

(1.2)

$\\dV_\tau = k(\eta-V_\tau)d\tau+\delta_1\sqrt{V_\tau}dW_{\tau}^2,V(0) \gt 0,$

$\\dR_\tau = a(b(\tau)-R_\tau)d\tau+\delta_2\sqrt{R_\tau}dW_{\tau}^3,R(0) \gt 0,$

In this model for $0 < \tau\leq T$ with $T > 0$ the given maturity time of an option. The stocks process $S_\tau, V_\tau,$ and $R_\tau$ denote random variables that represent the asset price, its variance, and the interest rate, at time $\tau$ respectively. The variance process $V_\tau$ was proposed by Heston and interest rate process $R_\tau,$ by Cox, Ingersoll and Ross in Eq.(1.2). Parameters $k, \eta, \delta_1, \delta_2,$ and $a,$ are given positive real numbers and $dW_{\tau}^1, dW_{\tau}^2, dW_{\tau}^3$ are Brownian motions with given correlation factors $\rho_1, \rho_2, \rho_3\in[-1, 1].$ Also $b,$ is a given deterministic positive function of time which can be chosen such that to match the current term structure of interest rates ^[12,13].

The introduction of Legendre wavelet method for the variational problems by Razzaghi and Yousefi in 2000 and 2001 ^[3,11,14], several works applying this method were born. The Legendre wavelets method (LWM) is transformed a boundary value problem (BVP) into a system of algebraic equations ^[14]. The unknown parameter of this system is the vector that its components are the decomposition coefficients of the BVP solution into Legendre wavelets basis. In this paper we apply the multi- Legendre wavelet method to solve B-S and HCIR partial differential equations.

This paper is organized as follows: In section 2, PDE models of Black-Scholes equation and Heston-Cox-Ingersoll-Ross equation have derived from SDE models, Section 3 has described properties and formulation of Legendre wavelet and Multi-dimensional Legendre Wavelet, Section4, metodology of multi-dimensional LWM is applied on PDEs that use in section 2. In section 5, numerical results and figures are presented and section 6, Conclusion is used to review.

2. Deriving PDE models of B-S and HCIR

In his section we discuss about the B-S PDE model and HCIR PDE model that they have derived from SDE models ^[1,2,19].

Theorem 2.1. (Ito-Doeblin formula for Brownian motion). Let $f (t, x)$ be a function for which the partial derivatives $f_{t}(t, x), f_{x}(t, x)$ and $f_{xx}(t, x)$ are defined and continues and let $W(t)$ be a Brownian motion. Then for every $T\geq0,$

$\begin{equation} f(t,W(t)) = f(0,W(0))+\int_{0}^{T}f_{t}(t,W(t))dt \end{equation}$

(2.1)

$+\int_{0}^{T}f_{x}(t,W(t))dt+\int_{0}^{T}f_{xx}(t,W(t))dt,$

Proof: see ^[6].

Remark 2.1: In the Ito- Doeblin formula we can write in differential form ^[6],

$\begin{equation} df(t,W(t)) = f_{t}(t,W(t))+f_{t}(t,W(t))dW(t)+\frac{1}{2}f_{t}(t,W(t))dW(t)dW(t),\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+f_{tx}(t,W(t))dtdW(t)+\frac{1}{2}f_{xx}(t,W(t))dtdt, \end{equation}$

(2.2)

But

$\begin{equation} dW(t)dW(t) = d(t), d(t)dW(t) = dW(t)d(t) = 0, d(t)d(t) = 0, \end{equation}$

(2.3)

And the Ito- Doeblin formula in differential form simplifies to

$\begin{equation} df(t,W(t)) = f_{t}(t,W(t))dt+f_{x}(t,W(t))dW(t)+\frac{1}{2}f_{xx}dt. \end{equation}$

(2.4)

Definition 2.1: let $W(t), t\geq0$ be a Brownian motion and let $F(t), t\geq0$ be an associated filtration. An Ito process as the form following

$\begin{equation} X(t) = X(0)+\int_{0}^{t}\bigtriangleup(u)dW(u)+\int_{0}^{t}\theta(u)d(u), \end{equation}$

(2.5)

Where $X(0)$ is nonrandom and $\triangle(u)$ and $\theta(u)$ are adapted stochastic process.

Lemma 2.1: the quadratic variation of the Ito process Eq.(2.5) is

$\begin{equation} [X,X](t) = \int_{0}^{t}\triangle^{2}(u)du. \end{equation}$

(2.6)

Proof: see ^[6].

Definition 2.2: let $X(t), t\geq0$ be an Ito process as described in Definition 2.1 and let $\Gamma(t), t\geq0$ be an adapted process. We define the integral with respect to an Ito process,

$\begin{equation} \int_{0}^{t}\Gamma(u) dX(u) = \int_{0}^{t}\Gamma(u) \triangle(u)dW(u)+\int_{0}^{t}\Gamma(u)\theta(u)du. \end{equation}$

(2.7)

Theorem 2.2. (Ito-Doeblin formula for an Ito process). Let $X(t), t\geq0$ be an Ito process as described in Definition 2.1 and let $f(t, x)$ be a function for which the partial derivatives $f_{t}(t, x), f_{x}(t, x)$ and $f_{xx}(t, x)$ are defined and continues. Then for every $T \geq0,$

$\begin{equation} f(T,X(T)) = f(0,X(0))+\int_{0}(T)f_{t}(t,X(t))dt \end{equation}$

(2.8)

$+\int_{0}^{T}f_{x}(t,X(t))dX(t)+\frac{1}{2}\int_{0}^{T}f_{xx}(t,X(t))d[X,X](t) =$

$f(0,X(0))+\int_{0}^{T}f_{t}(t,X(t))dt+\int_{0}^{T}f_{x}(t,X(t)) \triangle(t)dW(t)dt$

$+\int_{0}^{T}f_{x}(t,X(t)) \theta dt+\frac{1}{2}\int_{0}^{T}f_{xx}(t,X(t))\triangle(t)^{2}(t)$

We may rewrite Eq.(2.8),

$\begin{equation} df(t,x(t)) = f_{t}(t,X(t))dt+f_{x}(t,X(t))dX(t)+\frac{1}{2}f_{xx}(t,X(t))dX(t)dX(t). \end{equation}$

(2.9)

proof: See ^[6].

A SDE has defined as te following

$\begin{equation} dX(u) = \beta(u,X(u))du+\gamma(u,X(u))dW(u). \end{equation}$

(2.10)

The functions $\beta(u, X(u))$ and $\gamma(u, X(u))$ called drift and diffusion respectively. Where $t\geq0$ final goal in these problems is to find a stochastic proces $X(T)$ which is governed by

$\begin{equation} X(t) = x, \end{equation}$

(2.11)

$\begin{equation} X(t) = X(t)+\int_{t}^{T}\beta(u,X(u))du+\int_{t}^{T}\gamma(u,X(u))dW(u), \end{equation}$

(2.12)

Under conditions of the functions $\beta(u, x)$ and $\gamma(u, x)$ there exists a unique process $X(T), T\geq t$ Satisfying Eq.(2.11) and Eq.(2.12) and generally this is a SDE form of following

$\begin{equation} dX(u) = (a(u)+b(u)X(u))du+(\gamma(u)+\delta(u)X(u))dW(u), \end{equation}$

(2.13)

Where $a(u), b(u), \delta(u)$ and $\gamma(u)$ are non-random functions of time and we permit on the right- hand side of Eq.(2.10) is the randomness inherent in the solution $X(u)$ and in the driving Brownian motion $W(u).$ Then the stochastic differential equation for geometric Brownian motion is

$\begin{equation} dS(u) = \alpha S(u) du+\gamma S(u) dW(u), \end{equation}$

(2.14)

The rate of change $\frac{dS}{S},$ can be decomposed in to two parts in Eq.(2.14): One is deterministic, the other one is random.

The deterministic part can be modeled by $\frac{dS}{S} = \mu dt,$ where $\mu$ is a measure of the growth rate of the asset and random part modeled by a Brownian motion $\delta dW(t)$ that $\delta$ is the variance of the return and is called the volatility. The overall asset price model is then given by

$\begin{equation} \frac{dS}{S} = \mu dt+\gamma dW(t), \end{equation}$

(2.15)

Eq.(2.15), is a SDE and it can be expressed as follows

$\begin{equation} dS(t) = \mu (S(t))dt+\delta(S(t))dW_{s}(t). \end{equation}$

(2.16)

An important lemma for finding their solution is the following Ito lemma.

Lemma 2.2: Suppose $S(t)$ satisfies the stochastic differential Eq.(2.16), and $u(S, t)$ is a smooth function. Then $u(S(t), t)$ satisfies the following stochastic differential equation ^[8].

$\begin{equation} du = (u_{t}+\mu u_{s}+\frac{1}{2} \delta^{2} u_{ss})dt+\delta u_{s} dW_{s}. \end{equation}$

(2.17)

By expanding Eq.(2.17) with applying different variables in Ito lemma to earn B-S PDE and HCIR PDE equations, as following Black-Scholes PDE equation ^[1].

Option transaction contract is divided into different types which the most important types are call option and put option. call option contract means that buyer have to buy the certain asset with given price in certain time, as the same way the put option contract leads to sale of asset with given price and certain time by seller. Maturity validity includes American option and European option ^[5]. Where $W(t),$ is a standard Brownian motion. We also assume that interest rate are constant so that one unit of currency invested in the cash account at time $0$ , will be worth

$B(t) = exp(rt),$

at time $t.$ We will value of call option or put option at time $t.$ By Ito Lemma we know that,

$\begin{equation} du(S,t) = (\mu S(t) \frac{\partial u}{\partial S}+ \frac{\partial u}{\partial t} + \frac{1}{2} \delta^{2} S^{2}(t) \frac{\partial^{2}u}{\partial S^2})dt+\delta S(t) \frac{\partial u}{\partial S} dW(t), \end{equation}$

(2.18)

Let us now consider a self-financing trading strategy where at each time $t,$ we hold units $x_{t},$ of the cash account and $y_{t},$ units of the stock. Then $p_{t},$ value of this strategy satisfies,

$\begin{equation} p_{t} = x_{t}B_{t}+y_{t}S_{t}, \end{equation}$

(2.19)

We will choose $x_{t}$ and $y_{t}$ in such a way that the strategy replicates the value of the option the self-financing assumption implies that

$\begin{equation} dp_{t} = x_{t} dB_{t}+y_{t} dS_{t}, \end{equation}$

(2.20)

$\begin{equation} = rx_{t}B_{t} dt+y_{t}(\mu S_{t} dt+\delta S_{t} dWt) \end{equation}$

(2.21)

$= ( rx_{t}B_{t} dt+y_{t}\mu S_{t}) dt+y_{t}\delta S_{t} dWt,$

Note that Eq.(2.20) is consistent with our earlier definition of financing. In particular gains or losses on the portfolio are due entirely to gain or losses in the underlying securities i.e. the cash account and stock and not due to change in holding $x_{t},$ and $y_{t}.$

Returning to our derivation, we can equate times in Eq.(2.18), to obtain

$\begin{equation} y_{t} = \frac{\partial u}{\partial S}, \end{equation}$

(2.22)

$rx_{t}B_{t} = \frac{\partial u}{\partial S}+\frac{1}{2} \delta^{2} S^{2} \frac{\partial^2(u)}{\partial S^2},$

If we set

$C_{0} = P_{0},$

that $C$ and $P$ are call option and put option in finance, the interval value of strategy, then it must be the case that

$C_{t} = P_{t},$

since call option or put option have the same dynamics this is true by construction after we equate terms in Eq.(2.19), with the corresponding terms in Eq.(2.22).

Substituting Eq.(2.22) into Eq.(2.19), we obtain the Black-Scholes PDE

$\begin{equation} rS_{t} \frac{\partial u}{\partial S}+ \frac{\partial u}{\partial t}+\frac{1}{2} \delta^{2} S^{2} \frac{\partial^2(u)}{\partial S^2} -ru = 0. \end{equation}$

(2.23)

In order to solve Eq.(2.23), boundary conditions must be provided in the case of call option or put option.

By applying system of Eq.(1.1) in Eq.(2.17), for construction of PDE model we need part of deterministic section of a stochastic differential equation model as,

$\begin{equation} du(S,v,r,t) = u_{t}(S,v,r,t)dt+u_{s}(S,v,r,t)d(S(t))+u_{r}(S,v,r,t)d(r(t))+u_{v}(S,v,r,t)d(v(t)) \\ \;\; \; +\frac{1}{2} u_{ss} (S,v,r,t)d(S(t))d(S(t))+\frac{1}{2} u_{rr} (S,v,r,t)d(r(t))d(r(t))+\frac{1}{2} u_{vv} (S,v,r,t)d(v(t))d(v(t)) \\\;\; \;\;\;+u_{sv}(S,v,r,t)d(S(t))d(v(t))+u_{sr}(S,v,r,t)d(S(t))d(r(t))+u_{vr}(S,v,r,t)d(v(t))d(r(t)). \end{equation}$

(2.24)

Heston-Cox-Ingersoll-Ross PDE equation ^[6,7,8,9],

$\begin{equation} \frac{du}{dt} = rS\frac{du}{dS}+k(\eta -v)\frac{du}{dv}+a(b(t-T)-r)\frac{du}{dr}+\frac{1}{2}S^{2}v\frac{\partial^{2} u}{\partial S^2}+\frac{1}{2}\delta{1}^{2}v\frac{\partial^{2} u}{\partial v^2} \end{equation}$

(2.25)

$+\frac{1}{2}\delta_{2}^{2}r \frac{\partial^{2} u}{\partial r^2}+\rho_{12}\delta_{1}Sv \frac{\partial^{2} u}{\partial S\partial v}+\rho_{13}S \sqrt{vr}\frac{\partial^{2} u}{\partial S\partial r}+\rho_{23}\delta_{1}\delta_{2}\sqrt{vr} \frac{\partial^{2} u}{\partial v\partial r}-ru.$

By Eq.(2.25) we are obtained PDE model of the HCIR equation.

3. Construction of multi-dimensional Legendre wavelet

In this study Legendre wavelet well-addressed in ^[5,10,18], are employed.The wavelet basis is constructed from a single function, which is called the mother wavelet. These basis functions are called wavelets and they are an orthonormal set. One of the most important wavelets is Legendre wavelets. The Legendre wavelets are obtained from Legendre polynomials. In the past decade special attention has been given to applications of wavelets. The main characteristic of Legendre wavelet is that it reduces to a system of algebraic equation.

The function $\psi (x) \in L^2 (R)$ is a mother wavelet and $\psi (x)_{uv} = \vert u\vert ^\frac{-1}{2}\psi(\frac{x-v}{u})$ which $u, v\in R$ and $u\neq 0$ is a family of continuous wavelets. If we choose the dilation parameter $u = a^{-n}$ and the translation parameter $v = ma^{-n} b$ where $a > 1, b > 0, n$ and $m$ are positive integer, we have the discrete orthogonal wavelets set.

$\lbrace \psi (x)_{n,m} = \vert a\vert ^\frac{n}{2}\psi (a^{n}x-mb): m,n\in Z\rbrace$

The Legendre wavelet is constructed from Legendre function ^[4,15]. One dimension Legendre wavelets over the interval $[0, 1]$ are defined as,

$\begin{equation} \psi_{m,n}(t) = \sqrt{(m+\frac{1}{2})}2^{\frac{k}{2}}P_{m}(2^{k}x-2n+1) , \frac{\widehat{n}-1}{2^{k-1}} 0 , o.w \end{equation}$

(3.1)

With $n = 1, 2..., 2k-1, m = 0, 1, 2, ..., M-1.$ In Eq.(3.1) $P_{m}$ are ordinary Legendre functions of order defined over the interval [-1, 1]. Legendre wavelet is an orthonormal set as,

$\begin{equation} \int_{0}^{1} \psi (x)_{n,m} \psi (x)_{n^\prime, m^\prime}dx = \delta_{n,n^\prime}\delta_{m,m^\prime}, \end{equation}$

(3.2)

The function $f(x) \in L^{2}([0, 1] \times [0, 1]),$ may be expanded as the following,

$\begin{equation} f(x)\cong \sum\limits_{n = 1}^{ 2^{k}-1} \sum\limits_{m = 0}^{ M-1} c_{n,m} \psi (x)_{n,m}, \end{equation}$

(3.3)

Now we defined four dimensions Legendre wavelets in $L^{4} ([0, 1]\times[0, 1]\times[0, 1]\times[0, 1]),$ as the form

$\begin{equation} \psi_{n_{1,2,3,4},m_{1,2,3,4}}(x) = \begin{cases} A_{m_{1,2,3,4}}(\prod\limits_{i = 1}^{4}P_{m_{i}}(2^{k_{i}}x_{i}-2n_{i}+1)),&\frac{{n_{i}}-1}{2^{k_{i}-1}}\leq x_{i} \leq \frac{{n_{i}}}{2^{n_{i}-1}}, i = 1,..,4\\ 0 , & o.w \end{cases} \end{equation}$

(3.4)

that

$A_{m_{1,2,3,4}} = \sqrt{\prod\limits_{i = 1}^{4}(m_{i}+\frac{1}{2})}2^{\frac{\sum\limits_{i = 1}^{4}k_{i}}{2}},\\\\\\\\ m_{i} = 0,1,...,M_{i}-1, n_{i} = 1,2,...,2^{k_{1}-1} {,i = 1,..,4}$

Four dimensions Legendre wavelets are orthonormal set over $[0, 1]\times[0, 1]\times[0, 1]\times[0, 1].$

$\begin{equation} \int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\psi_{n_{1,2,3,4},m_{1,2,3,4}}(x,y,z,q)\psi_{n^{\prime}_{1,2,3,4},m^{\prime}_{1,2,3,4}}(x,y,z,q)dxdydzdq = \end{equation}$

(3.5)

$\delta _{n_{1},n^{\prime}_{1}},\delta _{n_{2},n^{\prime}_{2}},\delta _{n_{3},n^{\prime}_{3}},\delta _{n_{4},n^{\prime}_{4}},\delta _{m_{1},m^{\prime}_{1}},\delta _{m_{2},m^{\prime}_{2}},\delta _{m_{3},m^{\prime}_{3}},\delta _{m_{4},m^{\prime}_{4}}.$

The function $u(x, y, z, q)\in L^{4}([0, 1]\times[0, 1]\times[0, 1]\times[0, 1]),$ can be expanded as fallows

$\begin{equation} u(x,y,z,q) = X(x)Y(y)Z(z)Q(q)\cong \prod\limits_{i = 1}^{4}(\sum\limits_{n_{i} = 1}^{ 2^{k_{i}}-1}\sum\limits_{m_{i} = 1}^{ 2^{M_{i}}-1} c_{n_{i},m_{i}}\psi _{n_{i},m_{i}}(x,y,z,q)). \end{equation}$

(3.6)

Where

$c_{n_{1,2,3,4},m_{1,2,3,4}} = \int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1} X(x)Y(y)Z(z)Q(q) \psi _{n_{1,2,3,4},m_{1,2,3,4}}dxdydzdq.$

The truncated version of Eq.(3.6), can be expressed as the following,

$\begin{equation} u(x,y,z,q) = C^{T}\Psi (x,y,z,q), \end{equation}$

(3.7)

Where $C$ and $\Psi (x, y, z, q),$ are coefficients matrix and wavelets vector matrix, respectively. The dimensions of those are $\lbrace\prod _{i = 1}^{4}(2^{k_{i}-1}M_{i})\rbrace\times1,$ and are given as the form

$C = [ c_{10101010},...,c_{1010101M_{1}-1},c_{10101020},...,c_{1010102M_{2}-1},c_{10101030},...,c_{1010103M_{3}-1},$

$c_{1010102^{k_{1}-1}0},...,c_{1010102^{k_{1}-1}M_{1}-1},...,c_{1110101M_{1}-1},...,c_{1120},...,c_{112 M_{1}-1},$

$...,c_{2^{k_{1}-1} M_{{1}-1}2^{k_{2}-1} 0 2^{k_{3}-1} 0 2^{k_{4}-1} 0},...,c_{2^{k_{1}-1} M_{{1}-1} 2^{k_{2}-1} M_{{2}-1}2^{k_{3}-1} M_{{3}-1}2^{k_{4}-1} M_{{4}-1},...} ]^T,$

and also

$\Psi = [ \psi_{10101010},...,\psi_{1010101M_{1}-1},\psi_{10101020},...,\psi_{1010102M_{2}-1},\psi_{10101030},...,\psi_{1010103M_{3}-1},$

$\psi_{1010102^{k_{1}-1}0},...,\psi_{1010102^{k_{1}-1}M_{1}-1},...,\psi_{1110101M_{1}-1},...,\psi_{1120},...,\psi_{112 M_{1}-1},$

$...,\psi_{2^{k_{1}-1} M_{{1}-1}2^{k_{2}-1}0 2^{k_{3}-1}0 2^{k_{4}-1} 0},...,\psi_{2^{k_{1}-1} M_{{1}-1}2^{k_{2}-1} M_{{2}-1}2^{k_{3}-1} M_{{3}-1}2^{k_{4}-1} M_{{4}-1},... ]^T}.$

The integration of the product of four Legendre wavelet function vectors is obtained as

$\begin{equation} \int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\Psi (x,y,z,q)\Psi^{T} (x,y,z,q) = I. \end{equation}$

(3.8)

Where $I,$ is diagonal unit matrix.(All of procedures in section 3, are established)

Lemma 3.1: if $h\in C ([0, T], L^{2}[0, 1])$ then the functions $C_{n, m}(t)$ are continuous on $[0, 1]$ ^[13].

Lemma 3.2: if $h\in C^{1} ((0, T), L^{2}[0, 1])$ then the function $C_{n, m}(t)$ coefficient belong to $C^{1} (0, T),$ furthermore if $\frac{\partial h}{\partial t} \in L^{2}((0, T), L^{2}[0, 1])$ then,

$\begin{equation} \frac{dC_{n,m}(t)}{dt} = \int_{0}^{1}\frac{\partial h(t,x)}{\partial t}\psi_{n,m} (x)dx. \end{equation}$

(3.9)

Proof: see ^[13].

3.1. Operational matrix of integration

The integration matrix for one variable $x,$ define as follows

$\begin{equation} \int_{0}^{x}\Psi ({x^{\prime}},y,z,q) dx^{\prime} = P_{x} \Psi (x,y,z,q). \end{equation}$

(3.10)

In which

$\\P_{x} = \frac{1} {M_{1}2^{k_{1}}}\begin{bmatrix} \ L & F & F& \cdots &F \\0 &L & F &\cdots &F\\ \vdots & \vdots & \ddots & \ddots & \vdots \\ 0 & 0 & 0& \cdots & L \end{bmatrix},\\$

$P_{x},$ is a $2^{k_{1}-1}2^{k_{2}-1}2^{k_{3}-1}2^{k_{4}-1} M_{1}M_{2}M_{3}M_{4}\times 2^{k_{1}-1}2^{k_{2}-1}2^{k_{3}-1}2^{k_{4}-1} M_{1}M_{2}M_{3}M_{4}$ matrix and $L, F$ and $O,$ are $2^{k_{1}-1} \times M_{1} \times2^{k_{2}-1} \times M_{2} \times 2^{k_{3}-1} \times M_{3}\times 2^{k_{4}-1} \times M_{4}$ matrices ^[14]. (In case of two variables Eq.(3.10), is established clearly).

Theorem 3.1. A function $f(x, y, z, l) \in L^{4}([0, 1)\times [0, 1)\times[0, 1)\times[0, 1))$ is with bounded second derivative say $\vert f^{\prime\prime }(x, y, z, l) \leq M,$ can be expanded as an infinite sum of Legendre wavelets and the series converges uniformly to the function $f(x, y, z, l)$ that is ^[15],

$\begin{equation} f(x,y,z,l) = \sum\limits _{n_{1} = 1}^{\infty}\sum\limits _{m_{1} = 0}^{\infty} \sum\limits _{n_{2} = 1}^{\infty}\sum\limits _{m_{2} = 0}^{\infty}\sum\limits _{n_{3} = 1}^{\infty}\sum\limits _{m_{3} = 0}^{\infty}\sum\limits _{n_{4} = 1}^{\infty}\sum\limits _{m_{4} = 0}^{\infty}c_{n_{1,..,4} ,m_{1,..,4}} \psi _{n_{1,..,4} ,m_{1,..,4}} (x,y,z,l) , \end{equation}$

(3.11)

Proof:

we have proved this theorem for multi-dimension states as following,

$c_{n_{1,..,4}, m_{1,..,4}} = \int _{0}^{1}\int _{0}^{1}\int _{0}^{1}\int _{0}^{1}f(x,y,z,l) dxdydzdl =$

${\int _{\frac{\widehat{n_{1}}-1}{2^{k}}}^{\frac{\widehat{n_{1}}+1}{2^{k}}}{\int _{\frac{\widehat{n_{2}}-1}{2^{k}}}^{\frac{\widehat{n_{2}}+1}{2^{k}}}{\int _{\frac{\widehat{n_{3}}-1}{2^{k}}}^{\frac{\widehat{n_{3}}+1}{2^{k}}}{\int _{\frac{\widehat{n_{4}}-1}{2^{k}}}^{\frac{\widehat{n_{4}}+1}{2^{k}}} f(x,y,z,l) (\frac{2m_{1}+1}{2})^\frac{1}{2} (\frac{2m_{2}+1}{2})^\frac{1}{2} (\frac{2m_{3}+1}{2})^\frac{1}{2} (\frac{2m_{4}+1}{2})^\frac{1}{2}}}}}$

$2^\frac{k_{1}+k_{2}+k_{3}+k_{4}}{2} P(2^{k_{1}}x-\widehat{n})P(2^{k_{2}}x-\widehat{n})P(2^{k_{3}}x-\widehat{n})P(2^{k_{4}}x-\widehat{n})dxdydzdl.$

Now let,

$2^{k_{1}}x-\widehat{n} = t_{1},2^{k_{2}}x-\widehat{n} = t_{2},2^{k_{3}}x-\widehat{n} = t_{3},2^{k_{4}}x-\widehat{n} = t_{4},$

then,

$dx = \frac{1}{2^{k_{1}}} dt_{1},dy = \frac{1}{2^{k_{2}}} dt_{2},dz = \frac{1}{2^{k_{3}}} dt_{3},dl = \frac{1}{2^{k_{4}}} dt_{4}.$

$c_{n_{1,..,4}, m_{1,..,4}} = \int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1} f(\frac{\widehat{n_{1}}+t_{1}}{2^{k_{1}}}) f(\frac{\widehat{n_{2}}+t_{2}}{2^{k_{2}}})f(\frac{\widehat{n_{3}}+t_{3}}{2^{k_{3}}})f(\frac{\widehat{n_{4}}+t_{4}}{2^{k_{4}}}) (\frac{2m_{1}+1}{2})^ \frac{1}{2}$

$(\frac{2m_{2}+1}{2})^ \frac{1}{2}(\frac{2m_{3}+1}{2})^ \frac{1}{2}(\frac{2m_{4}+1}{2})^ \frac{1}{2} 2^{\frac{k_{1}+k_{2}+k_{3}+k_{4}}{2}} P(t_{1})P(t_{2})P(t_{3})P(t_{4}) \frac{dxdydzdl}{2^{k_{1}+k_{2}+k_{3}+k_{4}}} =$

with get out the constants as followes

$(\frac{2m_{1}+1}{2^{k_{1}+1}})^ \frac{1}{2} (\frac{2m_{2}+1}{2^{k_{2}+1}})^ \frac{1}{2} (\frac{2m_{3}+1}{2^{k_{3}+1}})^ \frac{1}{2} (\frac{2m_{4}+1}{2^{k_{4}+1}})^ \frac{1}{2} \int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1} f(\frac{\widehat{n_{1}}+t_{1}}{2^{k_{1}}})$

$f(\frac{\widehat{n_{2}}+t_{2}}{2^{k_{2}}})f(\frac{\widehat{n_{3}}+t_{3}}{2^{k_{3}}})f(\frac{\widehat{n_{4}}+t_{4}}{2^{k_{4}}}) P_{m_{1}}(t_{1})P_{m_{2}}(t_{2})P_{m_{3}}(t_{3}) P_{m_{4}}(t_{4})dt_{1}dt_{2}dt_{3}dt_{4} =$

and polynomials are substituted

$( \frac{1}{2^{k_{1}+1}(2m_{1}+1)})^\frac{1}{2}( \frac{1}{2^{k_{2}+1}(2m_{2}+1)})^\frac{1}{2}( \frac{1}{2^{k_{3}+1}(2m_{3}+1)})^\frac{1}{2}( \frac{1}{2^{k_{4}+1}(2m_{4}+1)})^\frac{1}{2}$

$\int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1} f(\frac{\widehat{n_{1}}+t_{1}}{2^{k_{1}}}) f(\frac{\widehat{n_{2}}+t_{2}}{2^{k_{2}}})f(\frac{\widehat{n_{3}}+t_{3}}{2^{k_{3}}})f(\frac{\widehat{n_{4}}+t_{4}}{2^{k_{4}}}) d(P_{m_{1}+1}(t_{1})-P_{m_{1}-1}(t_{1}))$

$d(P_{m_{2}+1}(t_{2})-P_{m_{2}-1}(t_{2}))d(P_{m_{3}+1}(t_{3})-P_{m_{3}-1}(t_{3}))d(P_{m_{4}+1}(t_{4})-P_{m_{4}-1}(t_{4})) =$

$(\frac{2m_{1}+1}{2^{k_{1}+1}})^ \frac{1}{2} (\frac{2m_{2}+1}{2^{k_{2}+1}})^ \frac{1}{2} (\frac{2m_{3}+1}{2^{k_{3}+1}})^ \frac{1}{2} (\frac{2m_{4}+1}{2^{k_{4}+1}})^ \frac{1}{2} f(\frac{\widehat{n_{1}}+t_{1}}{2^{k_{1}}}) f(\frac{\widehat{n_{2}}+t_{2}}{2^{k_{2}}}) f(\frac{\widehat{n_{3}}+t_{3}}{2^{k_{3}}}) f(\frac{\widehat{n_{4}}+t_{4}}{2^{k_{4}}})$

$(P_{m_{1}+1}(t_{1})-P_{m_{1}-1}(t_{1})) (P_{m_{2}+1}(t_{2})-P_{m_{2}-1}(t_{2})) (P_{m_{3}+1}(t_{3})-P_{m_{3}-1}(t_{3}))$

$(P_{m_{4}+1}(t_{4})-P_{m_{4}-1}(t_{4}))\vert_{-1}^{1} \vert_{-1}^{1}\vert_{-1}^{1}\vert_{-1}^{1} = -(\frac{1}{2^{3k_{1}+1}(2m_{1}+1)})^\frac{1}{2}$

Constants in terms of polynomials are calculated with the others,

$(\frac{1}{2^{3k_{2}+1}(2m_{2}+1)})^\frac{1}{2}(\frac{1}{2^{3k_{3}+1}(2m_{3}+1)})^\frac{1}{2} (\frac{1}{2^{3k_{4}+1}(2m_{4}+1)})^\frac{1}{2}\int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1}f(\frac{\widehat{n_{1}}+t_{1}}{2^{k_{1}}})$

$f(\frac{\widehat{n_{2}}+t_{2}}{2^{k_{2}}})f(\frac{\widehat{n_{3}}+t_{1}}{2^{k_{3}}})f(\frac{\widehat{n_{4}}+t_{4}}{2^{k_{4}}}) d( \frac{P_{m_{1}+2}(t_{1})-P_{m_{1}}(t_{1})} {2m_{1}+3} -\frac{P_{m_{1}}(t_{1})-P_{m_{1}-2}(t_{1})}{2m_{1}-1})$

$d( \frac{P_{m_{2}+2}(t_{2})-P_{m_{2}}(t_{2})} {2m_{2}+3} -\frac{P_{m_{2}}(t_{2})-P_{m_{2}-2}(t_{2})}{2m_{2}-1}) d( \frac{P_{m_{3}+2}(t_{3})-P_{m_{3}}(t_{3})} {2m_{3}+3} -\frac{P_{m_{3}}(t_{3})-P_{m_{3}-2}(t_{3})}{2m_{3}-1})$

$d( \frac{P_{m_{4}+2}(t_{4})-P_{m_{4}}(t_{4})} {2m_{4}+3} -\frac{P_{m_{4}}(t_{4})-P_{m_{4}-2}(t_{4})}{2m_{4}-1})$

that,

$= \vert \int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1} f^{\prime\prime}(\frac{\widehat{n_{1}}+t_{1}}{2^{k_{1}}})f^{\prime\prime}(\frac{\widehat{n_{2}}+t_{2}}{2^{k_{2}}})f ^{\prime\prime} (\frac{\widehat{n_{3}}+t_{1}}{2^{k_{3}}}) f ^{\prime\prime}(\frac{\widehat{n_{4}}+t_{4}}{2^{k_{4}}})$

$( \frac{P_{m_{1}+2}(t_{1})- P_{m_{1}}(t_{1})} {2m_{1}+3} -\frac{P_{m_{1}}(t_{1})-P_{m_{1}-2}(t_{1})}{2m_{1}-1}) ( \frac{P_{m_{2}+2}(t_{2})-P_{m_{2}}(t_{2})} {2m_{2}+3} -\frac{P_{m_{2}}(t_{2})-P_{m_{2}-2}(t_{2})}{2m_{2}-1})$

$( \frac{P_{m_{3}+2}(t_{3})-P_{m_{3}}(t_{3})} {2m_{3}+3} -\frac{P_{m_{3}}(t_{3})-P_{m_{3}-2}(t_{3})}{2m_{3}-1}) ( \frac{P_{m_{4}+2}(t_{4})-P_{m_{4}}(t_{4})} {2m_{4}+3} -\frac{P_{m_{4}}(t_{4})-P_{m_{4}-2}(t_{4})}{2m_{4}-1})\vert ^2$

$\vert \int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1} f^{\prime\prime}(\frac{\widehat{n_{1}}+t_{1}}{2^{k_{1}}})f^{\prime\prime}(\frac{\widehat{n_{2}}+t_{2}}{2^{k_{2}}})f ^{\prime\prime} (\frac{\widehat{n_{3}}+t_{1}}{2^{k_{3}}}) f ^{\prime\prime}(\frac{\widehat{n_{4}}+t_{4}}{2^{k_{4}}})$

$( \frac{(2m_{1}-1)P _{m_{1}+2}(t_{1})-(4m_{1}+2)P_{m_{1}}(t_{1})+(2m_{1}+3)}{2m_{1}+3})...$

Using integration properties as following,

$( \frac{(2m_{4}-1)P _{m_{4}+2}(t_{4})-(4m_{4}+2)P_{m_{4}}(t_{4})+(2m_{1}+3)}{2m_{4}+3})dt_{1} dt_{2}dt_{3}dt_{4} \vert ^2 \leq$

$\int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1} \vert ( \frac{(2m_{1}-1)P _{m_{1}+2}(t_{1})-(4m_{1}+2)P_{m_{1}}(t_{1})+(2m_{1}+3)}{2m_{1}+3})... dt_{1} dt_{2}dt_{3}dt_{4} \vert ^2$

$2M^{2}2M^{2}2M^{2}2M^{2} \int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1}\int _{-1}^{1}( \frac{(2m_{1}-1)P _{m_{1}+2}(t_{1})-(4m_{1}+2)P_{m_{1}}(t_{1})+(2m_{1}+3)}{2m_{1}+3})$

$... dt_{1} dt_{2}dt_{3}dt_{4} =$

Upper bounds are noted as,

$\frac{16 M^{16}}{(2m_{1}+3)^{2}(2m_{1}-1)^{2}...(2m_{4}+3)^{2} (2m_{4}-1)^{2} } ( (2m_{1}-1)^{2} \frac{2}{2m_{1}+5}+$

$(4m_{1}+2)^{2}\frac{2}{2m_{1}+1} + (2m_{1}+3)^{2}\frac{2}{2m_{1}-3})...((2m_{4}-1)^{2}\frac{2}{2m_{4}+5}+$

$(4m_{4}+2)^{2}\frac{2}{2m_{4}+1}+(2m_{4}+3)^{2}\frac{2}{2m_{4}-3}) \lt$

Now $(2m_{1}-1)^{2} < (2m_{1}+3)^{2},$ and $(2m_{1}+1)^{2} < (2m_{1}+3)^{2},$ then

$= \frac{331776M^{16}}{(2m_{1}-1)((2m_{1}-3)...(2m_{4}-1)(2m_{4}-3))^{2}}.$

Therefore we have,

$\vert c_{n_{1,..,4} ,m_{1,..,4}}\vert \lt (576) ^\frac{1}{2} \frac{1 }{(2n)^\frac{1}{2}} \frac{1}{(2m_{1}-3)^{2}...(2m_{4}-3)^{2}}.$

Thus the series $\sum _{n_{1} = 1}^{\infty}\sum _{m_{1} = 0}^{\infty} \sum _{n_{2} = 1}^{\infty}\sum _{m_{2} = 0}^{\infty}\sum _{n_{3} = 1}^{\infty}\sum _{m_{3} = 0}^{\infty}\sum _{n_{4} = 1}^{\infty}\sum _{m_{4} = 0}^{\infty}c_{n_{1, .., 4}, m_{1, .., 4}},$ is ab-solute convergent, it follows that

$\sum\limits _{n_{1} = 1}^{\infty}\sum\limits _{m_{1} = 0}^{\infty} \sum\limits _{n_{2} = 1}^{\infty}\sum\limits _{m_{2} = 0}^{\infty}\sum\limits _{n_{3} = 1}^{\infty}\sum\limits _{m_{3} = 0}^{\infty}\sum\limits _{n_{4} = 1}^{\infty}\sum\limits _{m_{4} = 0}^{\infty}c_{n_{1,..,4} ,m_{1,..,4}} \psi _{n_{1,..,4} ,m_{1,..,4}} (x,y,z,l) ,$

convergence to the function $f(x, y, z, l),$ uniformly.

Notice: This result for two dimensions in theorem $(3.1),$ is calculated as upper bound,

$\vert c_{n_{1}n_{2}m_{1}m_{2}}\vert \lt (288)^{\frac{1}{2}}\frac{1}{2n^{\frac{5}{2}}}\frac{1}{(2m_{1}-3)^{2}(2m_{2}-3)^{2}} .$

4. Methodology for B-S and HCIR PDEs

In this methodology section, the Black-Scholes model for American put problem have used in the form of moving-boundary problem.

Lemma 4.1: For American options, we have ^[5,6]

ⅰ) The optimal exercise time for American call option is $T$ and we have

$C = c,$

ⅱ) The optimal exercise time for American put option as earlier as possible, i.e. $t$ and we have

$P \geq p,$

ⅲ) The put-call parity for American option:

$S-K \lt C-P \lt S-ke^{-t(T-t)},$

As a consequence,

$P\leq K.$

Letting $S$ be the price of an asset $A$ at time $t$ the American early exercise constraint leads to the following model for the value $P(S, t)$ of an American put option to sell the asset $A.$

$\begin{equation} \frac{\partial P}{\partial t}+ \frac{1}{2} \delta^{2} S^{2} \frac{\partial^{2} P} {\partial S^{2}}+rS \frac{\partial P}{\partial S}-rP = 0, 0\leq t \lt T \end{equation}$

(4.1)

$P( S,T) = max ( E-S,0), S\geqslant 0,$

$\frac{\partial P} {\partial S} ( \overline{S}, t) = -1,$

$P(S,t) = E-S, 0 \leq S \lt S(\overline{t}).$

Where $S(\overline{t})$ represent the free boundary. The parameters $\delta, r$ and $E$ represent the volatility of the underlying asset, the interest rate, and the exercise price of the option. And the partial differential equation (PDE) associated with Eq.(2.25) for the fair values of European- style options forms a time- dependent connection- diffusion- reaction equation with mixed spatial- derivative terms.

$\begin{equation} \frac{\partial u}{\partial t} = \frac{1}{2} S^{2} v\frac{\partial^{2} u} {\partial S^{2}}+\frac{1}{2}\delta_{1}^{2} v\frac{\partial^{2} u} {\partial v^{2}} +\frac{1}{2}\delta_{2}^{2} r\frac{\partial^{2} u} {\partial r^{2}}+\rho _{12} \delta_{1}Sv \frac{\partial^{2} u} {\partial S \partial v}+\rho _{13} \delta_{2}S\sqrt{vr} \frac{\partial^{2} u} {\partial S \partial r} \end{equation}$

(4.2)

$\rho _{23} \delta_{1}\delta_{2}\sqrt{vr} \frac{\partial^{2} u} {\partial v \partial r}+rS \frac{\partial u}{\partial S}+k(\eta -v) \frac{\partial u}{\partial v} + a(b(T-t)-r)\frac{\partial u}{\partial r} -ru.$

On the unbounded three- dimensional spatial domain $S > 0, v > 0, r > 0$ with $0 < t\leq T.$ $u(S, v, r, t)$ denotes the fair value of a European- style option if at the time $\tau = T-t$ the asset price, its variance, and the interest rate be equal to $S, v$ and $r$ respectively. Eq.(4.2) is known as the Heston- Cox- Ingersoll- Ross (HCIR) PDE. This equation usually appears with initial and boundary conditions which are determined by the specific option under consideration. Initial conditions in European call options equation is usually as follows,

$\begin{equation} u(S,v,r,T) = max (0,S-k). \end{equation}$

(4.3)

With a given strike price $k > 0.$

In this methodology section the Eq.(4.1) and Eq.(4.2), have been solved using two dimensional and four dimensional Legendre wavelets respectively,

First we assume for B-S PDE,

$\begin{equation} \frac{\partial P}{\partial t} = -\alpha \frac{\partial^{2} P}{\partial S^{2}}-\alpha \frac{\partial P}{\partial S}+\alpha P, \end{equation}$

(4.4)

With boundary condition $P(S, 0) = \beta_{1} (S, t),$ and $\frac{\partial P(S, 0)}{\partial t} = \beta_{2}(S, t).$

Those coefficients of $\alpha_{i}$ the same parameters and variables which be in Eq.(4.1), in practice they will be calculate with this technique in follow.

$\begin{equation} \frac{\partial P}{\partial t} = C_{1}^{T}(S,t) \Psi (S,t), \end{equation}$

(4.5)

Integrating (4.4) with respect to second variable over $[0, t],$ we get

$\begin{equation} P(S,t) = C_{1}^{T}(S,t) P \Psi (S,t)+ \beta_{2} (S,t), \end{equation}$

(4.6)

And also

$\begin{equation} \frac{\partial P}{\partial S} = \frac{d C^{T}(S,t)}{dS} P\Psi (S,t)+\frac{d \beta _{2}(S,t)}{dS}, \end{equation}$

(4.7)

$\frac{\partial ^{2} P}{\partial S^{2}} = \frac{d^{2} C^{T}(S,t)}{dS^{2}} P \Psi (S,t)+\frac{d ^{2}\beta _{2}(S,t)}{dS^{2}},$

We have also

$\begin{equation} 1 = d^{T} \Psi (x), x = e^{T} \Psi (x). \end{equation}$

(4.8)

Substituting Eqs.(4.4) to (4.8) in Eq.(4.1), we obtain

$\begin{equation} C_{1}^{T}(S,t) \Psi = -\alpha_{1}(\frac{d ^{2} C^{T}(S,t)}{dS^{2}} P \Psi +\frac{d ^{2}\beta _{2}(S,t)}{dS^{2}}d^{T}\Psi )-\alpha_{2} (\frac{d C^{T}(S,t)}{dS} P \Psi +\frac{d \beta _{2}(S,t)}{dS}d^{T}\Psi ) \end{equation}$

(4.9)

$+\alpha_{3} (C_{1}^{T}(S,t)P \Psi+\beta_{2}d^{T}\Psi).$

The solution of B-S PDE, $P(S, t)$ can be calculate when Eq.(4.9) can be solved For unknown coefficients of the vector $C(S, t).$

$C = [ 0.12177112177-0.2214022140e-1t, -0.1661746162+0.3451318952t, 0.7030464532e-1-\\0.1278266279e-1t, -0.9594095941e-1+0.1992619926t]^T.$

Let for HCIR PDE in Eq.(4.2), we assume as following,

$\begin{equation} \frac{\partial u}{\partial t} = \alpha _{1}\frac{\partial ^{2} u}{\partial S^{2}} +\alpha _{2}\frac{\partial ^{2} u}{\partial v^{2}}+\alpha _{3}\frac{\partial ^{2} u}{\partial r^{2}}+\alpha _{4}\frac{\partial ^{2} u}{\partial S \partial v}+\alpha _{5}\frac{\partial ^{2} u}{\partial S \partial r}+\alpha _{6}\frac{\partial ^{2} u}{\partial v \partial r}+\alpha _{7}\frac{\partial u}{\partial S } \end{equation}$

(4.10)

$+\alpha _{8}\frac{\partial u}{ \partial v}+\alpha _{9}\frac{\partial u}{\partial r}+\alpha _{10} u.$

Also with boundary condition

$u(S,v,r,t) = \beta_{1} (S,v,r,t),$

and

$\frac{\partial u(S,v,r,0)}{\partial t} = \beta _{2}(S,v,r,t).$

Parameters of Eq.(4.2), by using method as following,

$\begin{equation} \frac{\partial u}{\partial t} = C_{1}^{T}(S,v,r,t) \Psi (S,v,r,t), \end{equation}$

(4.11)

Integrating right hand side in Eq.(4.10), with respect to differential equations over [0, t] we get As well as

$\begin{equation} \frac{\partial u}{\partial S} = \frac {d C ^{T}(S,v,r,t)}{dS} P\Psi (S,v,r,t)+\frac {d \beta ^{2}(S,v,r,t)}{dS}, \end{equation}$

(4.12)

$\frac{\partial u}{\partial v} = \frac {d C ^{T}(S,v,r,t)}{dv} P\Psi (S,v,r,t)+\frac {d \beta ^{2}(S,v,r,t)}{dv},$

$\frac{\partial u}{\partial r} = \frac {d C ^{T}(S,v,r,t)}{dr} P\Psi (S,v,r,t)+\frac {d \beta ^{2}(S,v,r,t)}{dr}.$

We have

$\begin{equation} \frac{\partial^{2} u}{\partial S^{2}} = \frac {d^{2} C ^{T}(S,v,r,t)}{dS^{2}} P\Psi (S,v,r,t)+\frac {d ^{2}\beta ^{2}(S,v,r,t)}{dS^{2}}, \end{equation}$

(4.13)

$\frac{\partial^{2} u}{\partial v^{2}} = \frac {d^{2} C ^{T}(S,v,r,t)}{dv^{2}} P\Psi (S,v,r,t)+\frac {d ^{2}\beta ^{2}(S,v,r,t)}{dv^{2}},$

$\frac{\partial^{2} u}{\partial r^{2}} = \frac {d^{2} C ^{T}(S,v,r,t)}{dr^{2}} P\Psi (S,v,r,t)+\frac {d ^{2}\beta ^{2}(S,v,r,t)}{dr^{2}}.$

And also

$\begin{equation} \frac{\partial^{2} u}{\partial S\partial v} = \frac{d}{dS}\frac {d C ^{T}(S,v,r,t)}{dv} P\Psi (S,v,r,t)+\frac{d}{dS}\frac {d \beta ^{2}(S,v,r,t)}{dv}, \end{equation}$

(4.14)

$\frac{\partial^{2} u}{\partial S\partial r} = \frac{d}{dS}\frac {d C ^{T}(S,v,r,t)}{dr} P\Psi (S,v,r,t)+\frac{d}{dS}\frac {d \beta ^{2}(S,v,r,t)}{dr},$

$\frac{\partial^{2} u}{\partial v\partial r} = \frac{d}{dv}\frac {d C ^{T}(S,v,r,t)}{dr} P\Psi (S,v,r,t)+\frac{d}{dv}\frac {d \beta ^{2}(S,v,r,t)}{dr}.$

Substituting (4.11) to (4.14) in (48), we obtain

$\begin{equation} C^{T}(S,v,r,t) \Psi = \alpha_{1}(\frac{\partial^{2} C^{T}(S,v,r,t)}{\partial S^{2}}P \Psi+ \frac{d^{2} \beta_{2}(S,v,r,t)}{dS^{2}}d^{T}\Psi ) \end{equation}$

(4.15)

$+\alpha_{2}(\frac{\partial^{2} C^{T}(S,v,r,t)}{\partial v^{2}}P \Psi+ \frac{d^{2} \beta_{2}(S,v,r,t)}{dv^{2}}d^{T}\Psi)+\alpha_{3}(\frac{\partial^{2} C^{T}(S,v,r,t)}{\partial r^{2}}P \Psi+ \frac{d^{2} \beta_{2}(S,v,r,t)}{dr^{2}}d^{T}\Psi)$

$+\alpha_{4}(\frac {d}{dS} \frac{dC^{T}(S,v,r,t)}{dv} P \Psi+ \frac{d}{dS} \frac{d\beta_{2}(S,v,r,t)}{dv}d^{T}\Psi)+\alpha_{5}(\frac {d}{dS} \frac{dC^{T}(S,v,r,t)}{dr}P \Psi+ \frac{d}{dS} \frac{d\beta_{2}(S,v,r,t)}{dr}d^{T}\Psi)$

$+\alpha_{6}(\frac {d}{dv} \frac{dC^{T}(S,v,r,t)}{dr}P \Psi+ \frac{d}{dv} \frac{d\beta_{2}(S,v,r,t)}{dr}d^{T}\Psi)+\alpha_{7}(\frac { \partial C^{T}(S,v,r,t) }{\partial S} P \Psi+ \frac{d\beta_{2}(S,v,r,t)}{dS}d^{T}\Psi)$

$+\alpha_{8}(\frac { \partial C^{T}(S,v,r,t) }{\partial v} P \Psi+ \frac{d\beta_{2}(S,v,r,t)}{dv}d^{T}\Psi)+\alpha_{9}(\frac { \partial C^{T}(S,v,r,t) }{\partial r} P \Psi+ \frac{d\beta_{2}(S,v,r,t)}{dr}d^{T}\Psi)$

$-\alpha_{10}( C^{T}(S,v,r,t) P \Psi+ \beta_{2} d^{T} \Psi) .$

This system can be solved for unknown coefficients, which are calculated,

$C[1,1] = 0.1816780437e-2\times sqrt(v) \times s \times t+.2310089049-0.3191757882e-2 \times s \times t^2-\\0.2243167683e-2 \times sqrt(v) \times s-0.3663030704 \times s \times t+0.5319596471e-3 \times t^2+0.4571378943 \times \\s+0.8386506682e-2 \times t,$

$C[1,2] = -0.4720134035e-3 \times sqrt(v) \times s \times t+0.8292430227e-3 \times s \times t^2+0.5827920594e-3 \times\\ sqrt(v) \times s+0.9516832935e-1\times s \times t-0.1382071705e-3 \times t^2-0.1187679088 \times s-\\0.2178878351e-2 \times t+0.1997897471,$

$C [1,3] = 0,C [1,4] = 0, \;\;\;C [1,5] = 0.2097293811,\\ C [1,6] = 0.1705107084,\;\;\;C [1,7] = 0,C [1,8] = 0, \\ C[1,9] = -0.04660652914, \;\;\; C[1,10] = -0.3789126853e-1,\\ C[1,11] = 0, C[1,12] = 0, \;\;\;C[1,13] = -0.4036243822e-1,\\ C[1,14] = -0.3281480113e-1 , \;\;\; C[1,15] = 0,C[1,16] = 0.$

5. Numerical results

In this section, values of parameters are used in B-S PDE model (): $k_{1} = k_{2} = 1, r = 0.2, E = 0.2, \delta = 0.001,$ and HCIR PDE model ( and ): $k = 0.1, \eta = 0.12, \rho_{12} = 0.6, \rho_{13} = 0.2, \rho_{23} = 0.4, \delta_{1} = 0.04, \delta_{2} = 0.03,$ respectively. Appling values of parameters in models, numerical results are calculated. Due to the accuracy and convergence of spectral method which has been reviewed in $[0, 1],$ the numerical results obtained from this method can be assured. Therefore it can be used for future scams to compare with the same methods.

Figure 1. Figure of solution B-S PDE model by multi-LWM.

DownLoad: Full-Size Img PowerPoint

Figure 2. Figure of solution HCIR PDE model by multi-LWM.

DownLoad: Full-Size Img PowerPoint

Figure 3. Figure of solution HCIR PDE model by ADM.

DownLoad: Full-Size Img PowerPoint

By choosing $M = 2, 4$ and $k = 1$ in Eq.(2.23) and $M = 1, 2$ in Eq.(2.25) the models are solved from the coefficients vector of C that was found. Then the results were compared by Adomian decomposition method (ADM), as a semi-analytical method ^[20,21]. Tables 1 and 2 give values of pricing options solutions in B-S and HCIR models.

Table 1. Corresponding solution of Black-Scholes PDE by LWM, ADM.

$t$	$P_{(LWM, M=2)}$	$P_{(LWM, M=4)}$	$P_{(ADM)}$	$Error$
$0.1$	$0.3060753725$	$0.3060954392$	$0.3068021827$	$0.000726810$
$0.2$	$0.02615249781$	$0. 02664386902$	$0.02613346987$	$0.000019028$
$0.4$	$0.04431200702$	$0.04942003420$	$0.04430950409$	$0.000005503$
$0.6$	$0.05847852762$	$0.05882487996$	$0.05862488320$	$0.000146355$
$0.8$	$0.08565205963$	$0.08534345864$	$0.08534345864$	$0.000308601$
$1$	$0.05083260299$	$0.0509281709$	$0.05075733333$	$0.000075268$

| Show Table

DownLoad: CSV

Table 2. Corresponding solution of Heston-Cox-Ingersoll-Ross PDE by LWM, ADM.

$(t, s)$	$U_{(LWM, M=1)}$	$U_{(LWM, M=2)}$	$U_{(ADM)}$	$Error$
$(5, 50)$	$0.1950128206e+3$	$0.1955141009e+3$	$0.1951180661e+3$	$0.0003960348$
$(10,100)$	$0.2721204331e+4$	$0.2725170985e+4$	$0.2789681554e+4$	$0.0064510569$
$(12,150)$	$0.6421235221e+4$	$0.6430612189e+4$	$0.6448380051e+4$	$0.0017767862$
$(15,200)$	$0.1461103897e+5$	$0.1470153697e+5$	$0.1486797294e+5$	$0.0016643597$
$(17,250)$	$0.2471265331e+5$	$0.2474067533e+5$	$0.2471107229e+5$	$0.0002960304$
$(19,300)$	$0.3852234601e+5$	$0.3858603408e+5$	$0.3807240307e+5$	$0.0051363101$
$(22,380)$	$0.6881012512e+5$	$0.6890201411e+5$	$0.6860384740e+5$	$0.0029816671$

| Show Table

DownLoad: CSV

6. Conclusion

The multi-dimensional Legendre wavelet method has been applied to solve the Black-Scholes and Heston-Cox-Ingersoll-Ross models which these PDE equations have been derived from stochastic differential equations by using one of the important stochastic calculous Lemma named Ito. Generally the present method (LWM) reduces the partial differential equation into a set of algebraic equation and it can be useful in high degree derivatives problems, therefore considering the properties and using structure of LWM leads to save the time. Multi- dimensional Legendre wavelet is orthonormal set and its integration of the product is a diagonal matrix as in the case of the one dimension Legendre wavelet. In addition the Theorem 3.1 has been proved in paper which it could be determine an upper bound for coefficients of vector $C$ . Also in this work results of LWM have been compared with the Adomian Decomposition Method (ADM) as a semi-analytical method and they have been shown for different parameters in Tables 1–2 and Figures 1–3, finally.

Conflict of interest

The authors declare that there is no conflicts of interest in this paper.

References

[1]	Hackbusch W, Sauter SA (1997) Composite finite elements for the approximation of PDEs on domains with complicated micro-structures. Numer Math 75: 447-472. doi: 10.1007/s002110050248
[2]	Hackbusch W, Sauter SA (1997) Composite finite elements for problems containing small geometric details. Part II: Implementation and numerical results. Comput Visual Sci 1: 15-25.
[3]	Hyman J, Shashkov M, Steinberg S (1997) The numerical solution of diffusion problems in strongly heterogeneous non-isotropic materials. J Comput Phys 132: 130-148. doi: 10.1006/jcph.1996.5633
[4]	Brezzi F, Lipnikov K, Simoncini V (2005) A family of Mimetic finite difference methods on polygonal and polyhedral meshes. Math Mod Meth Appl S 15: 1533-1551. doi: 10.1142/S0218202505000832
[5]	Brezzi F, Lipnikov K, Shashkov M (2005) Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes. SIAM J Numer Anal 43: 1872-1896. doi: 10.1137/040613950
[6]	Brezzi F, Lipnikov K, Shashkov M (2006) Convergence of mimetic finite difference method for diffusion problems on polyhedral meshes with curved faces. Math Mod Meth Appl S 16: 275-297. doi: 10.1142/S0218202506001157
[7]	Beirão da Veiga L, Lipnikov K, Manzini G (2014) The Mimetic Finite Difference Method for Elliptic Problems, Cham: Springer.
[8]	Antonietti PF, Formaggia L, Scotti A, et al. (2016) Mimetic finite difference approximation of flows in fractured porous media. ESAIM Math Model Numer Anal 50: 809-832. doi: 10.1051/m2an/2015087
[9]	Sukumar N, Tabarraei A (2004) Conforming Polygonal Finite Elements. Int J Numer Meth Eng 61: 2045-2066. doi: 10.1002/nme.1141
[10]	Dolbow J, Moës N, Belytschko T (2001) An extended finite element method for modeling crack growth with frictional contact. Comput Method Appl M 190: 6825-6846. doi: 10.1016/S0045-7825(01)00260-2
[11]	Tabarraei A, Sukumar N (2008) Extended finite element method on polygonal and quadtree meshes. Comput Method Appl M 197: 425-438. doi: 10.1016/j.cma.2007.08.013
[12]	Fries TP, Belytschko T (2010) The extended/generalized finite element method: An overview of the method and its applications. Int J Numer Meth Eng 84: 253-304.
[13]	Beirão da Veiga L, Brezzi F, Cangiani A, et al. (2013) Basic principles of virtual element methods. Math Mod Meth Appl S 23: 199-214. doi: 10.1142/S0218202512500492
[14]	Beirão da Veiga L, Brezzi F, Marini L, et al. (2016) Virtual element method for general second-order elliptic problems on polygonal meshes. Math Mod Meth Appl S 26: 729-750. doi: 10.1142/S0218202516500160
[15]	Antonietti PF, da Veiga LB, Mora D, et al. (2014) A stream virtual element formulation of the Stokes problem on polygonal meshes. SIAM J Numer Anal 52: 386-404. doi: 10.1137/13091141X
[16]	Beirão da Veiga L, Brezzi F, Marini LD, et al. (2016) Mixed virtual element methods for general second order elliptic problems on polygonal meshes. ESAIM Math Model Numer Anal 50: 727-747. doi: 10.1051/m2an/2015067
[17]	Antonietti PF, da Veiga LB, Scacchi S, et al. (2016) A C1 virtual element method for the CahnHilliard equation with polygonal meshes. SIAM J Numer Anal 54: 34-56. doi: 10.1137/15M1008117
[18]	Di Pietro DA, Ern A, Lemaire S (2014) An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators. Comput Meth Appl Math 14: 461-472.
[19]	Di Pietro DA, Ern A (2015) A hybrid high-order locking-free method for linear elasticity on general meshes. Comput Meth Appl Mech Eng 283: 1-21. doi: 10.1016/j.cma.2014.09.009
[20]	Di Pietro DA, Ern A (2015) Hybrid high-order methods for variable-diffusion problems on general meshes. CR Math 353: 31-34.
[21]	Di Pietro DA, Ern A, Lemaire S (2016) A review of hybrid high-order methods: Formulations, computational aspects, comparison with other methods. In: Barrenechea GR, Brezzi F, Cangiani A, et al. Editors, Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, Cham: Springer.
[22]	Cockburn B, Dond B, Guzmán J (2008) A superconvergent LDG-hybridizable Galerkin method for second-order elliptic problems. Math Comput 77: 1887-1916. doi: 10.1090/S0025-5718-08-02123-6
[23]	Cockburn B, Gopalakrishnan J, Lazarov R (2009) Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems. SIAM J Numer Anal 47: 1319-1365. doi: 10.1137/070706616
[24]	Cockburn B, Gopalakrishnan J, Sayas FJ (2010) A projection-based error analysis of HDG methods. Math Comput 79: 1351-1367. doi: 10.1090/S0025-5718-10-02334-3
[25]	Cockburn B, Guzmán J, Wang H (2009) Superconvergent discontinuous Galerkin methods for second-order elliptic problems. Math Comput 78: 1-24. doi: 10.1090/S0025-5718-08-02146-7
[26]	Antonietti PF, Giani S, Houston P (2013) hp-version composite discontinuous Galerkin methods for elliptic problems on complicated domains. SIAM J Sci Comput 35: A1417-A1439. doi: 10.1137/120877246
[27]	Antonietti PF, Giani S, Houston P (2014) Domain decomposition preconditioners for discontinuous Galerkin methods for elliptic problems on complicated domains. J Sci Comput 60: 203-227. doi: 10.1007/s10915-013-9792-y
[28]	Antonietti PF, Manzini G, Verani M (2018) The fully nonconforming virtual element method for biharmonic problems. Math Mod Meth Appl S 28: 387-407. doi: 10.1142/S0218202518500100
[29]	Ayuso de Dios B, Lipnikov K, Manzini G (2016) The nonconforming virtual element method. ESAIM Math Model Numer Anal 50: 879-904. doi: 10.1051/m2an/2015090
[30]	Cangiani A, Manzini G, Sutton OJ (2017) Conforming and nonconforming virtual element methods for elliptic problems. IMA J Numer Anal 37: 1317-1354.
[31]	Droniou J, Eymard R, Herbin R (2016) Gradient schemes: Generic tools for the numerical analysis of diffusion equations. ESAIM Math Model Numer Anal 50: 749-781. doi: 10.1051/m2an/2015079
[32]	Antonietti PF, Brezzi F, Marini LD (2009) Bubble stabilization of discontinuous Galerkin methods. Comput Meth Appl Mech Eng 198: 1651-1659. doi: 10.1016/j.cma.2008.12.033
[33]	Bassi F, Botti L, Colombo A, et al. (2012) On the flexibility of agglomeration based physical space discontinuous Galerkin discretizations. J Comput Phys 231: 45-65. doi: 10.1016/j.jcp.2011.08.018
[34]	Bassi F, Botti L, Colombo A, et al. (2012) Agglomeration based discontinuous Galerkin discretization of the Euler and Navier-Stokes equations. Comput Fluids 61: 77-85. doi: 10.1016/j.compfluid.2011.11.002
[35]	Bassi F, Botti L, Colombo A (2014) Agglomeration-based physical frame dG discretizations: An attempt to be mesh free. Math Mod Meth Appl S 24: 1495-1539. doi: 10.1142/S0218202514400028
[36]	Cangiani A, Georgoulis EH, Houston P (2014) hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes. Math Mod Meth Appl S 24: 2009-2041. doi: 10.1142/S0218202514500146
[37]	Cangiani A, Dong Z, Georgoulis EH, et al. (2016) hp-version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes. ESAIM Math Model Numer Anal 50: 699-725. doi: 10.1051/m2an/2015059
[38]	Cangiani A, Dong Z, Georgoulis EH (2017) hp-version space-time discontinuous Galerkin methods for parabolic problems on prismatic meshes. SIAM J Sci Comput 39: A1251-A1279. doi: 10.1137/16M1073285
[39]	Antonietti PF, Houston P, Hu X, et al. (2017) Multigrid algorithms for hp-version interior penalty discontinuous Galerkin methods on polygonal and polyhedral meshes. Calcolo 54: 1169-1198. doi: 10.1007/s10092-017-0223-6
[40]	Antonietti PF, Cangiani A, Collis J, et al. (2016) Review of discontinuous Galerkin finite element methods for partial differential equations on complicated domains, In: Lecture Notes in Computational Science and Engineering, Springer, 281-310.
[41]	Cangiani A, Dong Z, Georgoulis EH, et al. (2017) hp-Version Discontinuous Galerkin Methods on Polytopic Meshes, Springer International Publishing.
[42]	Antonietti PF, Mazzieri I (2018) High-order discontinuous Galerkin methods for the elastodynamics equation on polygonal and polyhedral meshes. Comput Meth Appl Mech Eng 342: 414-437. doi: 10.1016/j.cma.2018.08.012
[43]	Antonietti PF, Bonaldi F, Mazzieri I (2018) A high-order discontinuous galerkin approach to the elasto-acoustic problem. arXiv preprint arXiv:180301351.
[44]	Antonietti PF, Pennesi G (2019) V-cycle multigrid algorithms for discontinuous Galerkin methods on non-nested polytopic meshes. J Sci Comput 78: 625-652. doi: 10.1007/s10915-018-0783-x
[45]	Alboin C, Jaffré J, Roberts JE, et al. (2000) Domain decomposition for some transmission problems in flow in porous media. In: Numerical Treatment of Multiphase Flows in Porous Media, Berlin: Springer, 22-34.
[46]	Alboin C, Jaffré J, Roberts JE, et al. (2002) Modeling fractures as interfaces for flow and transport in porous media. In: Fluid Flow and Transport in Porous Media: Mathematical and Numerical Treatment (South Hadley, MA, 2001), Providence: Amer. Math. Soc., 13-24.
[47]	Martin V, Jaffré J, Roberts JE (2005) Modeling fractures and barriers as interfaces for flow in porous media. SIAM J Sci Comput 26: 1667-1691. doi: 10.1137/S1064827503429363
[48]	Frih N, Roberts JE, Saada A (2008) Modeling fractures as interfaces: A model for Forchheimer fractures. Comput Geosci 12: 91-104. doi: 10.1007/s10596-007-9062-x
[49]	Fumagalli A, Scotti A (2013) A numerical method for two-phase flow in fractured porous media with non-matching grids. Adv Water Resour 62: 454-464. doi: 10.1016/j.advwatres.2013.04.001
[50]	Jaffré J, Mnejja M, Roberts J (2011) A discrete fracture model for two-phase flow with matrix-fracture interaction. Procedia Comput Sci 4: 967-973. doi: 10.1016/j.procs.2011.04.102
[51]	Angot P, Boyer F, Hubert F (2009) Asymptotic and numerical modelling of flows in fractured porous media. M2AN Math Model Numer Anal 43: 239-275. doi: 10.1051/m2an/2008052
[52]	Formaggia L, Scotti A, Sottocasa F (2018) Analysis of a mimetic finite difference approximation of flows in fractured porous media. ESAIM Math Model Numer Anal 52: 595-630. doi: 10.1051/m2an/2017028
[53]	Benedetto MF, Berrone S, Pieraccini S, et al. (2014) The virtual element method for discrete fracture network simulations. Comput Meths Appl Mech Eng 280: 135-156. doi: 10.1016/j.cma.2014.07.016
[54]	Benedetto MF, Berrone S, Scialò S (2016) A globally conforming method for solving flow in discrete fracture networks using the virtual element method. Finite Elem Anal Des 109: 23-36. doi: 10.1016/j.finel.2015.10.003
[55]	Chave FA, Di Pietro D, Formaggia L (2018) A hybrid high-order method for Darcy flows in fractured porous media. SIAM J Sci Comput 40: A1063-A1094. doi: 10.1137/17M1119500
[56]	D'Angelo C, Scotti A (2012) A mixed finite element method for darcy flow in fractured porous media with non-matching grids. ESAIM Math Mod Numer Anal 46: 465-489. doi: 10.1051/m2an/2011148
[57]	Flemisch B, Fumagalli A, Scotti A (2016) A review of the XFEM-based approximation of flow in fractured porous media. In: Advances in Discretization Methods, Springer, 47-76.
[58]	Burman E, Hansbo P, Larson MG, et al. (2019) Cut finite elements for convection in fractured domains. Comput Fluids 179: 726-734. doi: 10.1016/j.compfluid.2018.07.022
[59]	Antonietti PF, Facciolà C, Russo A, et al. (2019) Discontinuous Galerkin approximation of flows in fractured porous media on polytopic grids. SIAM J Sci Comput 41: A109-A138. doi: 10.1137/17M1138194
[60]	Antonietti PF, Sarti M, Verani M (2015) Multigrid algorithms for hp-discontinuous Galerkin discretizations of elliptic problems. SIAM J Numer Anal 53: 598-618. doi: 10.1137/130947015
[61]	Masud A, Hughes TJ (2002) A stabilized mixed finite element method for darcy flow. Comput Meth Appl Mech Eng 191: 4341-4370. doi: 10.1016/S0045-7825(02)00371-7
[62]	Brezzi F, Hughes TJ, Marini LD, et al. (2005) Mixed discontinuous Galerkin methods for darcy flow. J Sci Comput 22: 119-145.
[63]	Arnold DN, Brezzi F, Cockburn B, et al. (2001/02) Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J Numer Anal 39: 1749-1779.
[64]	Wheeler MF (1978) An elliptic collocation-finite element method with interior penalties. SIAM J Numer Anal 15: 152-161. doi: 10.1137/0715010
[65]	Arnold DN (1982) An interior penalty finite element method with discontinuous elements. SIAM J Numer Anal 19: 742-760. doi: 10.1137/0719052
[66]	Cockburn B, Shu CW (1998) The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM J Numer Anal 35: 2440-2463. doi: 10.1137/S0036142997316712
[67]	Castillo P, Cockburn B, Perugia I, et al. (2000) An a priori error analysis of the local discontinuous Galerkin method for elliptic problems. SIAM J Numer Anal 38: 1676-1706. doi: 10.1137/S0036142900371003
[68]	Perugia I, Schötzau D (2003) The hp-local discontinuous galerkin method for low-frequency time-harmonic maxwell equations. Math Comput 72: 1179-1214.
[69]	Perugia I, Schötzau D (2002) An hp-analysis of the local discontinuous galerkin method for diffusion problems. J Sci Comput 17: 561-571. doi: 10.1023/A:1015118613130
[70]	Strang G, Fix GJ (1973) An analysis of the finite element method, Prentice-Hall Englewood Cliffs.
[71]	Antonietti PF, Houston P (2011) A class of domain decomposition preconditioners for hpdiscontinuous Galerkin finite element methods. J Sci Comput 46: 124-149. doi: 10.1007/s10915-010-9390-1
[72]	Stein EM (1970) Singular integrals and differentiability properties of functions, Princeton university press.
[73]	Talischi C, Paulino GH, Pereira A, et al. (2012) Polymesher: A general-purpose mesh generator for polygonal elements written in matlab. Struct Multidiscip O 45: 309-328. doi: 10.1007/s00158-011-0706-z
[74]	Brenner K, Hennicker J, Masson R, et al. (2016) Gradient discretization of hybrid-dimensional Darcy flow in fractured porous media with discontinuous pressures at matrix-fracture interfaces. IMA J Numer Anal 37: 1551-1585.
[75]	Antonietti PF, Facciolà C, Verani M (2020) Polytopic discontinuous Galerkin methods for the numerical modelling of flow in porous media with networks of intersecting fractures. MOX Rep 08, arXiv:2002.06420.
[76]	Antonietti PF, De Ponti J, Formaggia L, et al. (2019) Preconditioning techniques for the numerical solution of flow in fractured porous media. MOX Rep 17.
[77]	Antonietti PF, Houston P (2013) Preconditioning high-order discontinuous galerkin discretizations of elliptic problems. In: Domain Decomposition Methods in Science and Engineering XX, Springer, 231-238.
[78]	Antonietti PF, Houston P, Pennesi G, et al. (2020) An agglomeration-based massively parallel non-overlapping additive Schwarz preconditioner for high-order discontinuous Galerkin methods on polytopic grids. Math Comput DOI: https://doi.org/10.1090/mcom/3510.

This article has been cited by:

1.	Kazem Nouri, Milad Fahimi, Leila Torkzadeh, Dumitru Baleanu, Numerical method for pricing discretely monitored double barrier option by orthogonal projection method, 2021, 6, 2473-6988, 5750, 10.3934/math.2021339
2.	Lei Fu, Huanhe Dong, Chaudry Masood Khalique, Hongwei Yang, CONSERVATION LAWS, EXACT SOLUTIONS OF TIME-SPACE FRACTIONAL GENERALIZED GINZBURG-LANDAU EQUATION FOR SHALLOW WAKE FLOWS, 2021, 11, 2156-907X, 874, 10.11948/20200053
3.	Abdulaziz Alsenafi, Fares Alazemi, Javad Alavi, An efficient algorithm to solve the geometric Asian power option price PDE under the stochastic volatility model, 2024, 1017-1398, 10.1007/s11075-024-01794-z

Reader Comments

Your name:*

Email:*
© 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)