A numerical method for parabolic complementarity problem

1.   Introduction

We are interested in the following parabolic complementarity problem

where $x\in\Omega\subset\mathbb{R}$, $t\in(0, T)$ and $a(x) > 0$. The functions $u^{*}$ and $g$ are known and the initial and boundary conditions for the unknown solution $u$ is

This kind of dynamic problem arises in many real-world applications, such as engineering, stochastic control, mechanics, economics, and risk measures (see, e.g., ^{[1,2,3,4,5,6,7,8]}). In general, it is impossible to obtain an analytical solution of (1.1) except for some special cases and hence we have to rely on numerical methods in practice.

There are many numerical methods for the problem (1.1) and were studied by various authors. A popular approach is to reformulate the problem as a parabolic PDE using the linear penalty term of the form $\mu\left[p-u^{*}\right]_{+}$, where $\mu \gg 1$ is a constant, $p$ denotes the solution to the penalized equation and $[a]_{+}: = \max \{a, 0\}$. This method has been discussed in ^{[2,3,9,10,11,12]}. It has been proved in ^[2] that the solution of the penalized equation converges to the original one at a rate of $\mathcal{O}\left(\mu^{-\frac{1}{2}}\right)$. In ^[10,13,14] power penalty methods have been used for solving constrained optimization problems by using an improved penalty term

where $k\geq1$ is a constant. (This includes the linear penalty $k = 1$ mentioned above as the special case.) The authors proved in ^[10] that $p$ converges to the exact solution $u$ in a proper Sobolev norm at a rate of $\mathcal{O}\left(\mu^{-\frac{k}{2}}\right)$. However, the penalty term (1.3) has an unbounded derivative when $p-u^{*} \rightarrow 0^{+}$, and therefore it needs to be smoothed locally ^[10]. In ^[11], the authors generalized the penalty term in (1.3) to the form of

where $1\gg\epsilon > 0$ is a smoothing parameter. The authors proved that the penalized solution converges to that of the original one at a rate of $\mathcal{O}\left(\left[\mu^{-k}+\epsilon\left(1+\mu \epsilon^{\frac{1}{k}}\right)\right]^{1 / 2}\right)$.

The aim of this paper is to study the numerical solution of the following penalized equation

In the following, we let $f = g+\mu\epsilon^{\frac{1}{k}}$. We will first discretize (1.5) in space by using the centered finite difference method ^[15]. This leads to large-scale nonlinear ordinary differential equations, for which we use the trapezoidal rule ^[16,17] for time discretization. Such a time discretization avoids solving nonlinear equations, but at each time step, we have to handle a large-scale tridiagonal linear system, which is the major computation. By looking for the structure of the coefficient matrix we propose a backward substitution algorithm for such a linear system. Numerically, we find that this algorithm is more efficient than the widely used built-in command '$\texttt{backslash}$' in Matlab in terms of CPU time.

The rest of this paper is organized as follows. In Section 2 we recall some results concerning the unique solvability of (1.5) and convergence of $p$ to $u$. In Section 3 we present the space and time discretization for (1.5). We present the backward substitution algorithm in Section 4 for solving the large-scale tridiagonal linear system at each time step. Numerical results are given in Section 5 and we conclude this paper in Section 6.

2.   Some preliminary results

In this section, we recall some results for the penalized parabolic equation (1.5) at the continuous level. To this end, we introduce some notations used in the following. With $1 \leq s \leq \infty$, we let $\mathbb{L}^s(\Omega) = \left\{v:\left(\int_{\Omega}|v(x)|^s \mathrm{\; d} x\right)^{\frac{1}{k}} < \infty\right\}$ denote the space of all $s$-power integrable functions on $\Omega$. The inner product on $\mathbb{L}^{2}(\Omega)$ is denoted by $(\cdot, \cdot)_{\Omega}$. We use $\|\cdot\|_{\mathbb{L}^s(\Omega)}$ to denote the norm on $\mathbb{L}^s(\Omega)$ and $\mathbb{H}^s(\Omega)$ to denote the Sobolev space with norm $\|\cdot\|_{k, \Omega}$. Let $\mathbb{C}^s(\Omega)$ (respectively, $\mathbb{C}^s(\bar{\Omega})$) be the function set of which a function and its derivatives of up to order $s$ are continuous on $\Omega$ (respectively, $\bar{\Omega}$). Let $\mathbb{H}_{0}^s(\Omega) = \left\{v \in \mathbb{H}^s(\Omega): v(x) = 0, x \in \partial \Omega\right\}$. For any Hilbert space $\mathbb{H}(\Omega)$, we use $\mathbb{L}^s(0, T; \mathbb{H}(\Omega))$ to denote the space

where $1 \leq k \leq \infty$ and $\|\cdot\|_{\mathbb{H}}$ denotes the natural norm on $\mathbb{H}(\Omega)$. The norm in this space is denoted by $\|\cdot\|_{\mathbb{L}^s(0, T; \mathbb{H})}$, i.e.,

We use $\mathbb{H}^{-1}(\Omega)$ to denote the dual space of $\mathbb{H}^{1}(\Omega)$ and use $\langle\cdot, \cdot\rangle$ to denote the duality pair between a Hilbert space and its dual space.

By an integration by parts, it is clear that the variational problem corresponding to (1.5) is the following problem: find $p(t) \in \mathbb{H}_{0}^{1}(\Omega)$ such that for all $v \in\mathbb{H}_{0}^{1}(\Omega)$ it holds

with the initial condition $p(x, 0) = u_{0}(x)$ in $\Omega$ and $f = g+\mu\epsilon^{\frac{1}{k}}$.

Theorem 1. For $a(x) > 0$, problem (2.1) has a unique solution for fixed $\epsilon \geq 0$ and $\mu \geq 1$.

Proof. Clearly, problem (2.1) is equivalent to

where $\Phi_{\epsilon}(w) = (\phi(w)+\epsilon)^{\frac{1}{k}}-(\phi(0)+\epsilon)^{\frac{1}{k}}$. Clearly,

and $\Phi_{\epsilon}(w)$ is a monotonically increasing function of $w$ because of the monotonicity of $\left(\left[w-u^{*}\right]_{+}+\epsilon\right)^{\frac{1}{k}}$ in $w$. This implies

Since $a(x) > 0$, it holds

where we have used the Poincare-Friedrich inequality, i.e., $\|\nabla v\|_{0} \geq C\|v\|_{1} (\forall v \in \mathbb{H}_{0}^{1}(\Omega))$. By integration by parts, we therefore have

which holds for any $v, w \in\mathbb{H}_{0}^{1}(\Omega)$. From (2.2a) and (2.2b) we have

Similarly, by using (2.2b) it holds

because $\left(\Phi_{\epsilon}(v), v\right) = \left(\Phi_{\epsilon}(v)-\Phi_{\epsilon}(0), v-0\right) \geq 0$.

From the definitions of $\Phi_{\epsilon}$ we have

where we have used $\|v\|_{0}^{\frac{1}{k}} \leq \max \left\{1, \|v\|_{0}\right\}$ for $k \geq 1$. Here, $C_{1} \in \mathbb{L}^{2}(0, T)$ is a generic positive function of $t$ and $C_{2} > 0$ a generic constant independent of $\epsilon$ and $k$. This gives

Dividing both sides of this inequality by $\|w\|_{1}$ and taking the supremum with respect to $w$ we obtain

Now, by using (2.2d), (2.3) and (2.4) the unique existence of the solution of (2.1) is guaranteed by the theory established in ^[4,18].

We now present the convergence of the approximate solution $p$ of the penalized problem (1.5) to that of the original problem (1.1).

Theorem 2. Let $a(x) > 0$ and $\frac{\partial u}{\partial t} \in \mathbb{L}^{k+1}(\Omega)$. Then there exists a constant $C > 0$, independent of $u, p$ and $\mu$ such that

where $\mu, k$ and $\epsilon$ are the parameters used in (1.5).

Proof. In the asymptotic sense (i.e., $\epsilon\rightarrow0$) it is clear that (1.5) is equivalent to

(Here we have used $f = g+\mu\epsilon^{\frac{1}{k}}$.) In ^[3,11,12] it was proved that the solution of (2.6) and the solution of the original problem (1.1) satisfy (2.5). Since

it is clear that the solution of (1.5) and the solution of the original problem (1.1) satisfy (2.5) as well.

3.   Space and time discretizations

We now present space and time discretizations for the penalized equation (1.5). To this end, we partition the computation domain $\Omega\times[0, T]$ by mesh sizes $h$ and $\tau$ and denote an arbitrary grid on this domain by $(jh, n\tau)$, where $j = 0, 1, \dots, J$ and $n = 0, 1, \dots, N$. We first discretize the spatial derivative by the centered finite difference formula

where $a_{l} = a(x_l)$, $p_l(t)\approx p(x_l, t)$ (with $l = j-1, j, j+1$), $r_j = \mathcal{O}(h^2)$ is the truncation error and $j = 1, 2, \dots, J-1$. For $j = 0$ and $j = J$ it holds $p_j(t) = 0$ according to the boundary condition in (1.5).

Dropping this truncation error in the above formulas gives the approximations as

Substituting this into (1.5) gives the semi-discrete system as

with initial value condition ${{p\mathit{\boldsymbol{}}}}(0) = {{u\mathit{\boldsymbol{}}}}_{0}$, where

The Matrix $A$ is a tridiagonal matrix

We now introduce time discretization for (3.1). To match the second-order accuracy of the space discretization, we need a second-order time discretization. To this end, we rewrite (3.1) as an integral equation in the interval $[t_n, t_{n+1}]$

i.e.,

Let ${{f\mathit{\boldsymbol{}}}}(t)-A{{p\mathit{\boldsymbol{}}}}(t)-\mu\left(\left[{{p\mathit{\boldsymbol{}}}}(t)-{{u\mathit{\boldsymbol{}}}}^{*}(t)\right]_{+}+{\boldsymbol {\epsilon}}\right)^{\frac{1}{k}} = (q_1(t), q_2(t), \dots, q_{J-1}(t))^\top$. Then,

Let $q(t)$ be an arbitrary component of the vector function $(q_1(t), q_2(t), \dots, q_{J-1}(t))^\top$. The integral $\int_{t_n}^{t_{n+1}}q(s)ds$ is the area of the trapezoid with a curved edge as shown in Figure 1 on the left. In general, it is impossible to get the precise value of such an area. To get an approximation of this area we can consider a regular trapezoid with the curved edge being replaced by a straight edge; see Figure 1 on the right. Then, we have

where $\eta_n = \mathcal{O}(\tau^2)$ denotes the approximation error [15,Chapter 2].

Dropping the approximation error and letting ${{p\mathit{\boldsymbol{}}}}_n\approx {{p\mathit{\boldsymbol{}}}}(t_n)$, we therefore obtain an approximation of the right hand-side of (3.3):

where $\tilde{{f\mathit{\boldsymbol{}}}}_n = \frac{{{f\mathit{\boldsymbol{}}}}(t_n)+{{f\mathit{\boldsymbol{}}}}(t_{n+1})}{2}$. Substituting this into (3.3) gives the full discretization of the penalized equation (1.5):

where $n = 0, 1, \dots, N-1$. This is a nonlinear system due to the term

To handle such a nonlinear system we propose a fixed-point iteration as follows. We can use Newton's method to handle such a nonlinear problem ^[19], but in this case, we have to compute the Jacobian matrix of the whole nonlinear term. First, we rewrite this system as

(The vector ${{b\mathit{\boldsymbol{}}}}_n$ is a known term.) Then, we solve ${{p\mathit{\boldsymbol{}}}}_{n+1}$ via the following iterations

where $l$ is the iteration index and $l_{\max}$ is specified by the prescribed tolerance.

The following is the convergence analysis for the above fixed-point iteration.

Theorem 3. Let $a(x) > 0$ for $x\in\Omega$, $\epsilon > 0$ and $k > 1$. Then the fixed-point iteration (3.5) converges to the unique solution with a rate at least

if $\rho < 1$, where $a_{\max} = \max_{x\in\Omega}a(x)$.

Proof. From (3.5) we have

where $\tilde{{b\mathit{\boldsymbol{}}}}_n = (I_x+\frac{\tau A}{2})^{-1}{{b\mathit{\boldsymbol{}}}}_n$ and $I_x\in\mathbb{R}^{(J-1)\times(J-1)}$ is an identity matrix. Let

Then, it is sufficient to study the contraction property of the function $F({{x\mathit{\boldsymbol{}}}})$, i.e.,

where $\rho$ is the quantity that we need to look insight into. The analysis of $\rho$ is given as follows. It holds

which implies

From the structure of the matrix $A$ (cf. (3.2b)) we have

where

Let $\lambda(\cdot)$ denote an arbitrary eigenvalue of the involved matrix. It it clear

Since $a(x) > 0$ we know that $D_a$ is an invertible diagonal matrix with positive diagonal elements. Hence, by a similarity transform it holds

where $D_a^{\pm\frac{1}{2}} = {\rm diag}(a_1^{\pm\frac{1}{2}}, a_2^{\pm\frac{1}{2}}, \dots, a_{J-1}^{\pm\frac{1}{2}})$. The matrix $D_a^{\frac{1}{2}} \Lambda D_a^{\frac{1}{2}}$ is a symmetric positive definite matrix and thus it holds (see, e.g., [20, Chapter 5])

By letting $\tilde{{z\mathit{\boldsymbol{}}}} = D_a^{\frac{1}{2}}{{z\mathit{\boldsymbol{}}}}$ we have

Since $\Lambda$ is a symmetric positive definite matrix, from [20,Chapter 5] we have

From ^[21,22] we know that the eigenvalues of the tridiagonal matrix $\Lambda$ are given by

Since $J = \frac{1}{h}$, it holds

This gives

Let $a_{\max} = \max_{x\in\Omega}a(x)$ and $a_{\min} = \min_{x\in\Omega}a(x)$. Then, it holds

Substituting (3.12a) and (3.12b) into (3.11) gives

This together with (3.10a) and (3.10b) gives

This indicates that $I_x+\frac{\tau A}{2}$ is an invertible matrix, since $\lambda_{\min}\left(I_x+\frac{\tau A}{2}\right) > 0$.

We next estimate of the infinity-norm of the inverse matrix $\left(I_x+\frac{\tau A}{2}\right)^{-1}$. By the definition of the infinity-norm of a square matrix, it holds

Since $\lambda\left(I_x+\frac{\tau A}{2}\right)$ is invertible, the eigenvectors of this matrix consist of a complete basis of the space $\mathbb{R}^{J-1}$. Then, without loss of generality, we assume that ${{z\mathit{\boldsymbol{}}}}$ is an arbitrary normalized eigenvector (with corresponding eigenvalue $\lambda$), i.e., $\left(I_x+\frac{\tau A}{2}\right){{z\mathit{\boldsymbol{}}}} = \lambda{{z\mathit{\boldsymbol{}}}}$. Hence,

This implies ($\rm {cf.}$ (3.13)).

In summary, we have

We next explore the relationship between $\|Y({{x\mathit{\boldsymbol{}}}}_1)-Y({{x\mathit{\boldsymbol{}}}}_2)\|_\infty$ and $\|{{x\mathit{\boldsymbol{}}}}_1-{{x\mathit{\boldsymbol{}}}}_2\|_\infty$ with ${Y}({{x\mathit{\boldsymbol{}}}})$ being the function defined by (3.7). Since

it is sufficient to study the contraction property of $y(x): = ([x]_++\epsilon)^{\frac{1}{k}}$ with $x\in\mathbb{R}$. We claim

We consider three cases as shown in Figure 2. For the case $x_1\leq0$ and $x_2\leq0$, it is clear that (3.15) holds. For $x_1\leq0$ and $x_2 > 0$ (i.e., the middle case in Figure 2), it holds

where $\xi\in(0, x_2)$. It remains to consider $x_2 > x_1 > 0$ (cf. Figure 2 on the right). We have

In summary, for all the three cases shown in Figure 2 the estimate (3.15) holds. Hence,

Now, substituting (3.14) and (3.16) into (3.9) gives

This proves the desired result as stated by Theorem 3.

4.   Backward substitution algorithm

In the implementation of (3.5) we have to solve a linear system

for each iteration index $l\geq0$. This would be the major computation for solving the penalized equation (1.5) and in this section, we propose an algorithm to handle this problem. (There are also many other efficient algorithms can be used to solve such a banded linear system, such as the hybrid algorithm in ^[23].)

To clearly describe our idea, we let

Then, it holds

From the first equation we have

Substituting this into the middle equation in (4.3) with $j = 2$ gives

i.e.,

In general we have

i.e.,

where $j = 2, 3, \dots, J-2$. In particular for $j = J-1$ we have $y_{J-2} = \tilde{a}_{J-2}y_{J-1}+\tilde{b}_{J-2}$. This together with the last equation in (4.3) gives

Substituting the first equation into the second one gives

With $y_{J-1}$, we can get $y_{J-2}, y_{J-3}, \dots, y_1$ via a backward substitution according to (4.4a) and (4.4b).

In summary, we can solve the linear system (4.1) as follows.

5.   Numerical results

In this section, we present numerical results to study the proposed fixed-point iteration (3.5) and the backward substitution algorithm in Section 4. We use the following data

For the penalty parameters, we use $\epsilon = 1e-4$, $k = 4$ and several values of $\mu$. All numerical results are implemented by Matlab R2017a installed in a desk computer with Mac OS and 2.7 GHz Intel Core i5. The tolerance for the fixed-point iteration (3.5) is set to $\texttt{tol} = \frac{\min\{h^2, \tau^2\}}{10}$, which is sufficient to match the discretization errors.

Let $h = \tau = \frac{1}{100}$. Then, in Figure 3 we plot the profile of $p(x, t)$ and $u^*(x, t)$ for two values of $\mu$: $\mu = 3$ (top row) and $\mu = 15$ (bottom row). We see that the parameter $\mu$ has a remarkable influence on the solution $p$. In particular, for small $\mu$ (e.g., $\mu = 3$) it does not hold $p\leq u^*$ for all $(x, t)\in\Omega\times(0, T)$ as we can see in the top-right subfigure. For large $\mu$, as we can see in the bottom-right subfigure, it holds $p\leq u^*$ uniformly.

In Figure 4 we show the local details of the constraint $u^*$ and the computed solution $p$ at two different time points: $t = 2$ (left) and $t = 4$ (right). Here we consider three values of penalty parameter $\mu$: $\mu = 3, 7, 15$. The results in Figure 4 clearly indicate that as $\mu$ grows the condition $p\leq u^*$ holds uniformly on the space and time domains. This observation confirms the conclusion in ^[3,10,11,12].

In Figure 5 we show the measured iteration number of the fixed-point iteration (3.5) for three values of the penalty parameter $\mu$. (Such an iteration number is the quantity $L_{\max}$ such that the infinity norm of the residual arrives at the prescribed tolerance $\texttt{tol}$ = $\frac{\min\{h^2, \tau^2\}}{10}$.) We see that as $\mu$ grows the required iteration number of the fixed-point iteration (3.5) increases as well. This is an undesirable phenomenon and needs further study.

Let

where the subscript 'fp' denotes fixed-point iteration and ${\rm It}_{\rm fp}^n$ is the iteration number for the $n$-th time point. The quantity ${\rm It}_{\rm fp}^{\max}$ is the maximal iteration number over all the $N$ time points. In Figure 6 we show ${\rm It}_{\rm fp}^{\max}$ as we refine the discretization mesh sizes $h$ (and $\tau$) from $2^{-6}$ to $2^{-11}$. Clearly, the results in Figure 6 indicates that the convergence rate of the fixed-point iteration algorithm is robust with respect to the mesh sizes. Such robustness is a very important property when high accuracy (i.e., $h$ and $\tau$ should be small) is pursued.

We next study the backward substitution algorithm proposed in Section 4. We will compare it with the built-in command $\texttt{backslash}$ in Matlab. Let $p$ be the solution obtained by using the backward substitution algorithm proposed in Section 4 and $p_{\rm inv}$ is the solution obtained by using the $\texttt{backslash}$ command in Matlab. These two linear solvers are used to handle the same linear system in (3.5) for each fixed-point iteration. In Figure 7 we show the error between $p$ and $p_{\rm inv}$ on the space and time domains. We see that both linear solvers leads to the same numerical solution if we neglect the roundoff error due to floating point operations. This implies that the backward substitution algorithm indeed produces a reliable numerical solution.

Using the two linear solvers for implementing the fixed-point iteration (3.5), we show in Figure 8 the measured CPU time for solving the penalized equation (1.5). Clearly, the computation using the backward substitution algorithm needs much less CPU time and this indicates a great advantage for using it in practice.

6.   Conclusions

We have made a numerical study of the parabolic complementarity problems. We first represent the problem via a penalty method following the idea in ^{[2,3,9,10,11,12]} and then we discretize the penalized equation via a Crank-Nicolson method consisting of a centered finite difference formula for the space derivative and a trapezoidal rule for time discretization. In each time step, we have to solve a nonlinear system due to the penalty term and we proposed a fixed-point iteration to handle such a nonlinear system. The major computation of the fixed-point iteration is to solve a tridiagonal linear problem, for which we proposed a backward substitution algorithm, which is a direct linear solver (i.e., it is not iterative). Extensive numerical results are given, which indicate how the discretization mesh sizes and the penalty parameters affect the convergence rate of the fixed-point iteration. In terms of CPU time, the numerical results also indicate an obvious advantage of the backward substitution algorithm compared to the built-in command $\texttt{backslash}$ in Matlab. This would be a very important advantage to deal with large-scale nonlinear/linear systems, such as the one arising from discretizing the backward stochastic differential equation (BSDE) ^[6,8]. The research of BSDE in risk measure is a hot topic at present, which provides a new direction for us to further study the relation between numerical method for parabolic complementarity problem and risk measure. In the forthcoming study, we will further discuss how to use numerical method proposed in this paper to handle risk measure.

Acknowledgements

The authors are very grateful to the anonymous referees for their careful reading of a preliminary version of the manuscript and their valuable suggestions, which greatly improved the quality of this paper. This work is supported by the Guangdong Basic and Applied Basic Research Foundation (2020A1515110671), the Jiangmen Basic and Applied Basic Research Foundation (2021030100070004859).

Conflict of interest

The authors declare that there is no conflicts of interest.