Access the full text.
Sign up today, get DeepDyve free for 14 days.
Loading next page...
/lp/ou_press/an-exponential-integrator-based-discontinuous-galerkin-method-for-xKgejBRPw0
- Publisher
- Oxford University Press
- Copyright
- © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
- ISSN
- 0272-4979
- eISSN
- 1464-3642
- DOI
- 10.1093/imanum/drx056
- Publisher site
- See Article on Publisher Site
Abstract
Abstract The linear complementarity system (LCS) is defined by a linear ordinary differential equation coupled with a finite-dimensional linear complementarity problem (LCP), which has many applications in engineering and economics. In this article, we reformulate the LCS with the boundary condition as an LCP in the Hilbert space of square-integrable functions, and propose a new numerical method for the LCS by using exponential Euler integrator and discontinuous Galerkin approximation. The precision of the proposed method is better than that of the existing time-stepping method in different magnitude of scale. Convergence analysis and numerical experiments are performed to support the arguments. 1. Introduction Given matrices $$A\in\mathbb{R}^{n\times n}$$, $$B\in\mathbb{R}^{n\times m}$$, $$Q\in\mathbb{R}^{m\times n}$$, $$M\in\mathbb{R}^{m\times m}$$ and $$E_0,E_T\in\mathbb{R}^{n\times n}$$, a vector $$b\in\mathbb{R}^n$$ and vector-valued functions $$f:[0,T]\to\mathbb{R}^n$$ and $$g:[0,T]\to\mathbb{R}^m$$, the linear complementarity system (LCS) is to find a state-control pair $$(x,u)$$ of functions: $$x:[0,T]\to \mathbb{R}^n$$, $$u:[0,T]\to\mathbb{R}^m$$, such that \begin{equation}\label{lcs}\left\{\begin{array}{@{}r@{\ }c@{\ }ll}\dot{x}(t) & = & Ax(t)+Bu(t)+f(t)\quad & t\in[0,T]\\ u(t) & \in & \mathrm{SOL}(M,Qx(t)+g(t))\qquad & t\in[0,T]\\ b & = & E_0x(0)+E_Tx(T), \end{array}\right. \end{equation} (1.1) where for a matrix $$M\in\mathbb{R}^{m\times m}$$ and a vector $$q\in\mathbb{R}^m$$, SOL$$(M,q)$$ denotes the solution set of the linear complementarity problem $$({\rm{denoted\; by\; LCP}}(M,q))$$: $$\mathrm{SOL}(M,q)\ :=\ \left\{u\in\mathbb{R}^m|\ 0\ \le \ u\ \bot\ Mu+q\ \ge\ 0\right\}\!. $$ The LCS is a powerful mathematical modeling tool and finds various applications in, e.g., electrical networks with switching structure, contact mechanical systems and dynamical transportation assignment (Brogliato, 2003; Zhong et al., 2011). For many other applications in engineering and economics, refer to Heemels et al. (2000) and Pang & Stewart (2008). The existing numerical methods for LCSs normally utilize the time-stepping scheme (Pang & Stewart, 2008; Han et al., 2012; Chen & Wang, 2013). For a given mesh and step size \begin{equation}\label{mesh}0=t_0<t_1<\cdots<t_N=T,\qquad h=T/N, \end{equation} (1.2) the scheme computes the approximate solution $$(x^h,u^h)$$, where $$x^h$$ is piecewise linear continuous and $$u^h$$ piecewise constant in $$[0,T]$$ with $$x^h(t_k)=x^{h,k}$$ and $$u^h(t_k)=u^{h,k}$$, such that $$\left\{\begin{array}{@{}r@{\ }c@{\ }l}x^{h,k} & = & x^{h,k-1}+h\left(Ax^{h,k}+Bu^{h,k}+f(t_k)\right)\\ u^{h,k} & \in & \ \mathrm{SOL}\left(M,Qx^{h,k}+g(t_k)\right)\\ b & = & E_0x^{h,0}+E_Tx^{h,N}. \end{array}\right. $$ It was shown in Chen & Wang (2011) that if $$M$$ is a P-matrix, $$E_0=I$$ and $$E_T=0$$ then the initial value problem (IVP) of (1.1) has a classic solution $$(x,u)$$, where $$x$$ is continuously differentiable and $$u$$ is continuous on $$[0,T]$$, and the time-stepping method has first-order convergence.1 In general, the LCS (1.1) does not have a classic solution, and one has to seek a weak solution $$(x,u)$$, where $$x$$ is absolutely continuous and $$u$$ is integrable. The pair $$(x,u)$$, besides the boundary/initial condition, satisfies $$x(t)-x(s)=\int^t_s \left[Ax(\tau)+B u(\tau)+f(\tau)\right] \, {\rm{d}}\tau $$ and $$u(t)\ \in\ \mathrm{SOL}(M,Qx(t)+g(t)) $$ for almost every $$0\le s<t\le T$$. Note that the solutions of the LCS usually do not have a good smoothness, and therefore applying an integrator of high order, such as the collocation methods (Kunkel & Stöver, 2002), will not yield fast convergence. Notice that for the boundary value problem (BVP) of the LCS, the initial state is not prescribed, and we need to find the one, which leads to the terminal state such that the boundary condition is fulfilled. The dependence of the terminal state on the initial one is hard to be tracked. The time-stepping method was studied for the BVP of the LCS in Han et al. (2012), but the convergence rate was not established therein. In this article, we propose a new numerical method for solving (1.1) by combining the exponential integrator and the discontinuous Galerkin approximation. Below we summarize the idea and our contributions. At first, we reformulate the LCS (1.1) as an LCP$$(\mathfrak{L},\hat g)$$ in the Hilbert space $$\mathscr{L}^2(0,T;\mathbb{R}^m)$$ of the $$m$$-dimensional vector-valued square integrable functions over the interval $$[0,T]$$: finding $$u\in\mathscr{L}^2(0,T;\mathbb{R}^m)$$ such that the complementarity condition \begin{equation}\label{LCP_L2}0\ \le \ u(t)\ \bot\ (\mathfrak{L}u)(t)+\hat g(t)\ \ge 0 \end{equation} (1.3) holds a.e. in $$[0,T]$$. Here $$\mathfrak{L}$$ is a bounded linear operator on $$\mathscr{L}^2(0,T;\mathbb{R}^m)$$ w.r.t. $$\|\cdot\|_{\mathscr{L}^2}$$. Then we apply the Galerkin approximation to the equivalent variational inequality (VI) formulation of the LCP$$(\mathfrak{L},\hat g)$$ posed in the convex closed function family $$\mathscr{L}^2_+(0,T;\mathbb{R}^m)$$: \begin{equation}\label{domain}\mathscr{L}^2_+(0,T;\mathbb{R}^m):=\left\{u\in\mathscr{L}^2(0,T;\mathbb{R}^m) \ | \ u(t) \ge 0\ \ a.e.\ \mathrm{in}\ [0,T]\right\}\!. \end{equation} (1.4) The resulted discretized problem is a finite-dimensional LCP$$(M^h,q^h)$$, and its solution offers an approximate control $$u^h$$. An approximate state $$x^h$$ is computed by applying the exponential Euler integrator to the ordinary differential equation (ODE) with the approximation $$u^h$$ instead of $$u$$. To our knowledge, the idea of Galerkin approximation and exponential integrator has not been considered for the LCS. Refer to Hochbruck & Ostermann (2010) for comprehensive treatment of various exponential integrators. We establish the error estimate $$\left\|u^h-u\right\|_{\mathscr{L}^2}\le C\sqrt{\left\|\mathfrak{P}_+^hu-u\right\|_{\mathscr{L}^2}}, $$ where $$\mathfrak{P}_+^hu$$ denotes the projection of $$u$$ onto the family $$\mathscr{U}^h_+$$ of piecewise linear functions (normally discontinuous) that are componentwise non-negative over the time interval. The choice of the function family actually yields the discontinuous Galerkin approximation, and gives the pair $$\left(x^h,u^h\right)$$ with $$x^h$$ absolutely continuous and $$u^h\in\mathscr{U}^h_+$$. If $$u\in\mathscr{C}^2(0,T;\mathbb{R}^m)$$, the space of functions that have continuous second-order derivatives, then the above error estimate indicates first-order convergence in $$\|\cdot\|_{\mathscr{L}^2}$$ due to $$\left\|\mathfrak{P}_+^hu-u\right\|_{\mathscr{L}^2}=O(h^2)$$. To our surprise, the numerical results may suggest second-order convergence, both in $$\|\cdot\|_{\mathscr{L}^2}$$ and in $$\|\cdot\|_2$$ at the grid points (see Fig. 2). Fig. 2. View largeDownload slide Errors of $$(x^h_{\rm g},u^h_{\rm g})$$ in $$\|\cdot\|_{\mathscr{L}^2}$$ and $$\|\cdot\|_2$$ at $$T$$ (Example 5.1 with $$T=3\pi$$). Fig. 2. View largeDownload slide Errors of $$(x^h_{\rm g},u^h_{\rm g})$$ in $$\|\cdot\|_{\mathscr{L}^2}$$ and $$\|\cdot\|_2$$ at $$T$$ (Example 5.1 with $$T=3\pi$$). Of course the choice of the family of continuous piecewise linear functions gives the continuous Galerkin approximation. Note that the family of piecewise constant functions is involved in the time-stepping method. The approximate solutions given by our method is better in different magnitude of scale than those offered by the continuous Galerkin method and the time-stepping method, even for the meshes not so refined (Fig. 3). Fig. 3. View largeDownload slide Ratios of error for Galerkin and time-stepping method (Example 5.1 with $$T=3\pi$$). Fig. 3. View largeDownload slide Ratios of error for Galerkin and time-stepping method (Example 5.1 with $$T=3\pi$$). Note that the time-stepping method involves the evaluation of $$(I-hA)^{-1}$$, whereas our method requires the evaluation of $$\varphi_k(hA)$$, which is not necessarily time-consuming than that of the former, where $$\varphi_k(hA)$$ is some matrices related to the matrix exponential $$e^{hA}$$. Summarizing, the computational cost for our Galerkin approximation amounts to that for the time-stepping method, whereas the precision of the proposed method is much better, as illustrated by the numerical results (Fig. 2). This article is organized as follows: in Section 2, we study the reformulation of the LCS into an LCP in the space $$\mathscr{L}^2(0,T;\mathbb{R}^m)$$ and also investigate the solvability of the LCP therein. We present the Galerkin approximation in Section 3 and the convergence analysis in Section 4. In Section 5, we present the algorithmic details and report the numerical results for two examples coming from a switched mechanical system and a generalized Nash equilibrium problem. 2. Reformulation of the LCS as an LCP in the Hilbert space For a given LCS (1.1), we suppose throughout this article that the matrix $$E_0+E_T\,e^{TA}$$ is nonsingular, where $$e^{TA}$$ denotes the matrix exponential. Define matrix-valued kernel functions: \begin{align}K_1(t,s) & = e^{(t-s)A}B\qquad (\mathrm{for}\, \ t\ge s)\notag\\ K_2(t,s) & = e^{tA}(E_0+E_T\,e^{TA})^{-1}E_T \,e^{(T-s)A}B.\label{K12}\end{align} (2.1) Setting $$K_1(t,s)=0$$ for $$0\le t\le s\le T$$, we then have for any $$(t,s)\in[0,T]^2$$ \begin{align*}\|K_1(t,s)\|_F & \le e^{T\|A\|_F}\cdot \|B\|_F,\\ \|K_2(t,s)\|_F & \le e^{2T\|A\|_F}\left\|\left(E_0+E_T\,e^{TA}\right)^{-1}E_T\right\|_F\left\|B\right\|_F\!, \end{align*} where $$\|\cdot\|_F$$ denotes the Frobenius norm. Define the Volterra integral operators \begin{equation}\label{Volterra}(\mathfrak{L}_ju)(t) = \int^T_0K_j(t,s)u(s)\, {\rm{d}}s, \qquad (\,j=1,2). \end{equation} (2.2) Obviously, $$\mathfrak{L}_1$$ and $$\mathfrak{L}_2$$ are compact in $$\mathscr{L}^2(0,T;\mathbb{R}^m)$$, and therefore bounded, their norms are still denoted by $$\|\cdot\|_{\mathscr{L}^2}$$. Actually, from the bound $$\|K_j(t,s)\|_F$$ established above together with the inequality (see Conway, 1985) $$\|\mathfrak{L}_j\|_{\mathscr{L}^2}\ \le \ \left(\iint_{[0,T]^2}\|K_j(t,s)\|_F^2\, {\rm{d}}s\, {\rm{d}}t\right)^{\frac{1}{2}}, $$ it follows the bound $$\|\mathfrak{L}_j\|_{\mathscr{L}^2} \le \kappa_j$$, where \begin{equation}\label{L bound}\kappa_1= T \,e^{T\|A\|_F} \|B\|_F/\sqrt{2},\quad\kappa_2 = T\,e^{2T\|A\|_F}\left\|\left(E_0+E_T\,e^{TA}\right)^{-1}E_T\right\|_F\|B\|_F. \end{equation} (2.3) Let $$b\in\mathbb{R}^n$$, $$f:[0,T]\to\mathbb{R}^n$$ and $$g:[0,T]\to\mathbb{R}^m$$ be square integrable on $$[0,T]$$, and let the LCS (1.1) have a weak solution $$(x,u)$$. Then $$x$$, the solution of the ODE in (1.1) with $$u$$, can be represented by the constant variation formula (Hochbruck & Ostermann, 2010): $$x(t)=e^{tA}x(0)+\int^t_0 e^{(t-s)A}\left[Bu(s)+f(s)\right]\, {\rm{d}}s, $$ where the boundary condition $$E_0x(0)+E_Tx(T)=b$$ can be written as $$b \ = \ E_0x(0)+E_Tx(T)\ = \ \left(E_0+E_T\,e^{TA}\right)x(0)+ E_T \int^T_0 e^{(T-s)A}\left[Bu(s)\,{\rm d}s+f(s)\right]\, {\rm{d}}s. $$ It gives the initial state $$x(0)\ = \ (E_0+E_T\,e^{TA})^{-1}\left(b-E_T\int^T_0 e^{(T-s)A}\left[Bu(s)+f(s)\right]\, {\rm{d}}s\right) $$ and the solution of the BVP of the ODE: \begin{equation}\label{form x}x(t )=\hat f(t)+(\mathfrak{L}_1u)(t)-(\mathfrak{L}_2u)(t), \end{equation} (2.4) where \begin{equation}\label{hatf}\hat f(t)= e^{tA}\left(E_0+E_T\,e^{TA}\right)^{-1}\left(b-E_T\int^T_0 e^{(T-s)A}f(s)\, {\rm{d}}s\right)+ \int^t_0e^{(t-s)A}f(s)\, {\rm{d}}s. \end{equation} (2.5) By plugging (2.4) into the LCP in (1.1), we can see that $$u$$ is a solution of the following linear complementarity problem, denoted by $$\mathrm{LCP}(\mathfrak{L},\hat g)$$, which is to find $$u\in\mathscr{L}^2(0,T;\mathbb{R}^m)$$ such that for almost $$t\in[0,T]$$ it holds \begin{equation}0\ \le \ u(t)\ \bot\ (\mathfrak{L}u)(t)+\hat g(t)\ \ge 0, \end{equation} (2.6) where $$\mathfrak{L}$$ is the linear bounded operator \begin{equation}\label{L}(\mathfrak{L}u)(t):=Q(\mathfrak{L}_1u)(t)-Q(\mathfrak{L}_2u)(t)+Mu(t) \end{equation} (2.7) and \begin{equation}\label{hatg}\hat g(t) = g(t)+Q\,\hat f(t). \end{equation} (2.8) Now we have reformulated the LCS (1.1) into the $$\mathrm{LCP}(\mathfrak{L},\hat g)$$, an LCP posed in $$\mathscr{L}^2(0,T;\mathbb{R}^m)$$. Theorem 2.1 If $$(x,u)$$ is a weak solution of the LCS (1.1) then $$u$$ is a solution of the LCP$$(\mathfrak{L},\hat g)$$. Conversely, if $$u$$ is a solution of the LCP$$(\mathfrak{L},\hat g)$$ then $$(x,u)$$ is a weak solution of (1.1), where $$x$$ is given by (2.4). Proof. The first part has been shown above. Conversely, if $$u$$ is a solution of the LCP$$(\mathfrak{L},\hat g)$$ then $$x$$ is well defined by (2.4), which is absolutely continuous since its weak derivative $$Ax+Bu+f$$ is integrable. □ Remark 2.2 A specific case of the LCP$$(\mathfrak{L},\hat g)$$ is the so-called convolution complementarity problem (CCP), which models some contact mechanical problems (Stewart, 2006): $$0\ \le\ u(t) \ \bot\ \ \int^t_0 K(t-\tau)u(\tau)\, {\rm{d}}\tau+g(t)\ \ge 0, $$ where $$K(\cdot)$$ is a given kernel. The operator involved in CCP is compact, and some properties of compact operators can be utilized in the algorithmic construction and convergence analysis, while the operator $$\mathfrak{L}$$ defined in (2.7) is not compact in general. An obvious advantage of reformulating the LCS as the $$\mathrm{LCP}(\mathfrak{L},\hat g)$$ lies in that we can use rich theory and abundant numerical methods for operator equations (Pang & Qi, 1993; Chen et al., 1997) to treat the LCS, since $$\mathrm{LCP}(\mathfrak{L},\hat g)$$ can be equivalently reformulated as $$\min\{u(t),(\mathfrak{L}u)(t)+\hat g(t)\}\ = \ 0, $$ where ‘min’ is taken componentwise. To develop a Galerkin approximation scheme, we need the variational formulation of the LCP. It is easy to show that $$u$$ is a solution of the LCP$$(\mathfrak{L},\hat g)$$ if and only if $$u\in\mathscr{L}_+^2(0,T;\mathbb{R}^m)$$, and \begin{equation}\label{VI}\langle \mathfrak{L}u+\hat g, v-u\rangle_{\mathscr{L}^2} \ge\ 0,\qquad \forall v\in \mathscr{L}_+^2(0,T;\mathbb{R}^m), \end{equation} (2.9) where the set $$\mathscr{L}^2_+(0,T;\mathbb{R}^m)$$, defined as in (1.4), is convex and closed in $$\|\cdot\|_{\mathscr{L}^2}$$. When $$M$$ is positive semidefinite, $$\mathfrak{L}(\cdot):=Q\mathfrak{L}_1(\cdot)-Q\mathfrak{L}_2(\cdot)+M(\cdot)$$ is pseudo-monotone as $$\mathfrak{L}_1$$ and $$\mathfrak{L}_2$$ are linear compact (Fig. 27.1 of Zeidler, 1990, p. 596). As a direct consequence of Theorem 32.C of Zeidler (1990) (II/B, p. 875), we have the following results on the solvability of the LCS (1.1). Theorem 2.3 Let $$f\in\mathscr{L}^2(0,T;\mathbb{R}^n)$$ and $$g\in\mathscr{L}^2(0,T;\mathbb{R}^m)$$. Suppose that there is a $$u_0\in\mathscr{L}^2_+(0,T;\mathbb{R}^m)$$ such that \begin{equation}\label{coersive}\frac{\langle \mathfrak{L}u,u-u_0\rangle_{\mathscr{L}^2}}{\|u\|_{\mathscr{L}^2}}\ \to \ +\infty\quad \mathrm{as}\ \, \|u\|_{\mathscr{L}^2}\to \infty. \end{equation} (2.10) (1) If $$M$$ is positive semidefinite then the $$\mathrm{LCP}(\mathfrak{L},\hat g)$$ has a solution $$u\in\mathscr{L}^2(0,T;\mathbb{R}^m)$$. (2) If $$\mathfrak{L}$$ is monotone then the solution set of the $$\mathrm{LCP}(\mathfrak{L},\hat g)$$ is convex and closed in $$\|\cdot\|_{\mathscr{L}^2}$$. It is well known that if the operator $$\mathfrak{L}$$ is strongly monotone then the $$\mathrm{LCP}(\mathfrak{L},\hat g)$$ has a unique solution in $$\mathscr{L}^2(0,T;\mathbb{R}^m)$$ (II/B, Zeidler, 1990). From (2.3) we know the norms of $$\mathfrak{L}_j$$ can be bounded by a function of $$T$$ that is decreasing to $$0$$ when $$T\downarrow 0$$. This indicates that $$\mathfrak{L}$$ is strongly monotone if $$M$$ is positive definite and $$T$$ is small enough. We refine it in the following theorem. Theorem 2.4 Let $$M$$ be positive definite. If $$T>0$$ is small enough then $$\mathfrak{L}$$ defined in (2.7) is strongly monotone, and then the $$\mathrm{LCP}(\mathfrak{L},\hat g)$$ has a unique solution in $$\mathscr{L}^2(0,T;\mathbb{R}^m)$$. Proof. Denote by $$\sigma$$ the smallest eigenvalue of $$\frac{1}{2}(M+M^T)$$. Obviously, $$\sigma>0$$ since $$M$$ is positive definite. Let $$\kappa_j$$ be the constants defined in (2.3). It is easy to see $$\|Q\mathfrak{L}_j\|_{\mathscr{L}^2}\le\|Q\|_F\cdot \kappa_j$$ for $$j=1,2$$. Then for any $$u\in\mathscr{L}^2(0,T;\mathbb{R}^m)$$ we have \begin{align*}\langle u, \mathfrak{L}u\rangle & = \langle u, Mu\rangle +\langle u, Q\mathfrak{L}_1u\rangle-\langle u, Q\mathfrak{L}_2u\rangle\\ & \ge \displaystyle \left(\sigma-\|Q\|_F\sum_{j=1,2}\|\mathfrak{L}_j\|_{\mathscr{L}^2}\right)\cdot \|u\|^2_{\mathscr{L}^2} = [\sigma-\|Q\|_F(\kappa_1+\kappa_2) ]\cdot \|u\|^2_{\mathscr{L}^2}, \end{align*} which follows that $$\mathfrak{L}$$ is strongly monotone for $$T>0$$ small enough since $$\kappa_j$$ decreases to 0 as $$T\to 0^+$$. □ Remark 2.5 For the case of the IVP, namely, $$E_0=I$$ and $$E_T=0$$, we have $$\kappa_2=0$$. With the uniform mesh such that the step size $$h$$ fulfills $$\sqrt{2}\sigma > h\,e^{h\|A\|_F}\|Q\|_F\|B\|_F. $$ Theorem 2.4 yields the unique existence of the solution of the $$\mathrm{LCP}(\mathfrak{L},\hat g)$$ in $$[0,h]$$. Repeating the application of the theorem, one can establish the solvability of (1.1) for any $$T>0$$. The LCS is a special case of the so-called differential variational inequality Gwinner (2013), for which sufficient conditions for the existence of a solution (Theorem 2, pp. 392) was given in Chen & Wang (2014). In the current setting, this theorem reads as follows: Theorem 2.6 Let $$E_0+E_T$$ be nonsingular and $$\hat{x}^0=(E_0+E_T)^{-1}b$$. Denote the solution set of the parameterized LCP by $$\mathscr{S}(t,x):=\mathrm{SOL}(M,Qx(t)+g(t))$$. Suppose that $$\mathscr{S}(0, \hat{x}^0)$$ is nonempty and bounded, and $$M$$ is positive semidefinite. Denote $$\mathscr{F}(t,x):=\{Ax+Bu+f(t)| u\in \mathscr{S}(t,x)\}$$. If $$\mathscr{S}(t,x)$$ is lower semicontinuous near $$(0,\hat{x}^0)$$ or $$\mathscr{F}(t,x)$$ is singleton then there exist $$T, \delta_0, \zeta>0$$ such that the LCS (1.1) has a solution $$(x,u)$$, where $$x\in\mathscr{N}(\hat{x}^0,\delta_0+\zeta T)$$ is continuously differentiable with $$x(0)\in\mathscr{N}(\hat{x}^0,\delta_0)$$, and $$u$$ is continuous and is the least norm element of $$\mathscr{S}(t,x(t))$$, namely, $$u$$ is the least norm solution of the $$\mathrm{LCP}(M,Qx(t)+g(t))$$. 3. Galerkin approximation 3.1 Approximation of $$\mathscr{L}_+^2(0,T;\mathbb{R}^m)$$ We need an approximation of the function family $$\mathscr{L}_+^2(0,T;\mathbb{R}^m)$$. Let the mesh (1.2) be given, denote by $$\mathscr{U}^h$$ the space of piecewise linear functions. $$\mathscr{U}^h$$ is a $$(2Nm)$$-dimensional closed subspace of $$\mathscr{L}^2(0,T;\mathbb{R}^m)$$. Introduce \begin{equation}\label{basis}\hat \psi_1(t)=1-t,\qquad\quad \hat \psi_2(t)=t, \end{equation} (3.1) where $$\chi_k(t)$$ is the characteristic function in $$I_k=(t_{k-1},t_k]$$ and $$e_i$$ denotes the $$m$$-dimensional $$i$$th unit coordinate vector. For $$i=1,\dots,m$$, $$j=1,2$$ and $$k=1,\dots,N$$, define \begin{equation}\label{basis aufbau}\psi_{i,j,k}(t) := \hat\psi_j\left(\frac{t-t_{k-1}}{h}\right)\cdot\chi_k(t)\cdot e_i. \end{equation} (3.2) Obviously, the function family $$\{\psi_{i,j,k}: 1\le i\le m, 1\le j\le 2, 1\le k\le N\}$$ spans the subspace $$\mathscr{U}^h$$, and \begin{equation}\label{mass xishu}\langle \psi_{i,j,k},\psi_{i',j',k'}\rangle_{\mathscr{L}^2} = \left\{\begin{array}{@{}ll}h/3\qquad & i= i',\ k= k',\ j=j'\\ h/6 & i= i',\ k= k',\ j\not=j'\\ 0 & \mathrm{otherwise}. \\ \end{array}\right. \end{equation} (3.3) Note that $$u^h(t)=\sum_{i,j,k}c_{i,j,k} \psi_{i,j,k}(t)\in \mathscr{U}^h$$ has the following form in the interval $$I_k$$ $$u^h(t)=\psi_1\left(\frac{t-t_{k-1}}{h}\right)\left(\begin{array}{c}c_{1,1,k}\\ \vdots\\ c_{m,1,k}\end{array}\right)+\psi_2\left(\frac{t-t_{k-1}}{h}\right) \left(\begin{array}{c}c_{1,2,k}\\ \vdots\\ c_{m,2,k}\end{array}\right)\!. $$ The set $$\mathscr{U}^h_+$$ is convex and closed. And since $$u^h(t)\ge 0$$ holds almost everywhere in $$[0,T]$$ if and only if $$c_{i,j,k}\ge 0$$ holds for all the coefficients, it has the following representation: $$\mathscr{U}^h_+\ :=\ \mathscr{U}^h\cap \mathscr{L}^2_+(0,T;\mathbb{R}^m)\ =\ \left\{\sum_{i,j,k}c_{i,j,k}\cdot \psi_{i,j,k} :\ c_{i,j,k}\ge 0\right\}\!. $$ Denote by $$\mathscr{C}^2_+$$ the set of functions in $$\mathscr{L}^2(0,T;\mathbb{R}^m)$$ that are componentwise non-negative on $$[0,T]$$ and have continuous second-order derivatives. Clearly, $$\mathscr{C}^2_+$$ is a dense subset of $$\mathscr{U}^h_+$$. Below we show that $$\mathscr{U}^h_+$$ is an approximation to $$\mathscr{L}^2_+(0,T;\mathbb{R}^m)$$ in a certain sense. Theorem 3.1 Let $$u\in\mathscr{L}^2_+(0,T;\mathbb{R}^m)$$, and $$q_{i,j,k}=\langle u,\psi_{i,j,k}\rangle_{\mathscr{L}^2}/h$$, and denote by ‘mid’ the componentwise median operation. Then (1) $$u$$ has the unique projection $$\mathfrak{P}^+_hu=\sum_{i,j,k}c^*_{i,j,k}\psi_{i,j,k}$$ onto the closed and convex set $$\mathscr{U}^h_+$$ with respect to the norm $$\|\cdot\|_{\mathscr{L}^2}$$, where the coefficients $$c^*_{i,j,k}$$ are given by \begin{align}c^*_{i,1,k} & = mid \left\{0,4q_{i,1,k}-2q_{i,2,k},3q_{i,1,k}\right\}\notag\\ c^*_{i,2,k} & = mid \left\{0,4q_{i,2,k}-2q_{i,1,k},3q_{i,2,k}\right\}\!.\label{proj xishu}\end{align} (3.4) (2) Moreover, if $$u\in\mathscr{C}^2_+$$ then $$\lim_{h\to 0^+}\|u-\mathfrak{P}^+_hu\|_{\mathscr{L}^2}=0$$. Proof. (1) We mention that the projection is well defined because the subset $$\mathscr{U}^h_+$$ is convex and closed. Noting $$\langle u,\psi_{i,j,k}\rangle_{\mathscr{L}^2}=hq_{i,j,k}$$, and by using the inner product of the basic functions given in (3.3), for any $$u^h=\sum_{i,j,k}c_{i,j,k}\psi_{i,j,k}$$, we compute \begin{align*}& \|u-u^h\|^2_{\mathscr{L}^2}= \langle u,u\rangle_{\mathscr{L}^2}-2\langle u,u^h\rangle_{\mathscr{L}^2}+\langle u^h,u^h\rangle_{\mathscr{L}^2}\\ &\quad = \displaystyle \langle u,u\rangle_{\mathscr{L}^2}-2\sum_{i,j,k}c_{i,j,k}u_{i,j,k}+\sum_{i,j,k}\sum_{i',j',k'}\langle \psi_{i,j,k},\psi_{i',j',k'}\rangle_{\mathscr{L}^2} c_{i,j,k}c_{i',j',k'}\\ &\quad = \displaystyle \langle u,u\rangle_{\mathscr{L}^2}+\frac{h}{3}\sum_{i,k}\left(c^2_{i,1,k}+c^2_{i,2,k}+c_{i,1,k}c_{i,2,k}-6q_{i,1,k}c_{i,1,k}-6q_{i,2,k}c_{i,2,k}\right)\!. \end{align*} This objective function is separate in $$i$$ and $$k$$, namely, the minimization of $$\|u-u^h\|^2_{\mathscr{L}^2}$$ in $$\mathscr{U}^h_+$$ is equivalent to solve the $$mN$$ quadratic minimization problems of the form $$\min_{c_{i,1,k},c_{i,2,k}\ge 0}\left\{ c^2_{i,1,k}+c^2_{i,2,k}+c_{i,1,k}c_{i,2,k}-6q_{i,1,k}c_{i,1,k}-6q_{i,2,k}c_{i,2,k}\right\}$$ for $$i=1,\dots, m$$ and $$k=1,\dots,N$$, which can be equivalently reformulated as an LCP \begin{equation}\label{LCP erwei}0\ \le \left(\begin{array}{c}c_{i,1,k}\\ c_{i,2,k}\end{array}\right)\ \bot\ \left(\begin{array}{cc}2 & 1 \\ 1 & 2\end{array}\right)\left(\begin{array}{c}c_{i,1,k}\\ c_{i,2,k}\end{array}\right)-6\left(\begin{array}{c}q_{i,1,k}\\ q_{i,2,k}\end{array}\right)\ \ge 0. \end{equation} (3.5) The LCP (3.5) has a unique solution. Noting that $$u,\psi_{i,j,k}\in\mathscr{U}^h_+$$, therefore $$q_{i,j,k}=\langle u,\psi_{i,j,k}\rangle_{\mathscr{L}^2}/h\ge 0$$. It is easy to test for the case $$q_{i,2,k}\ge 2q_{i,1,k}$$, the LCP (3.5) has the solution $$c^*_{i,1,k}=0,\quad c^*_{i,2,k}=3q_{i,2,k}; $$ for the case $$q_{i,2,k}< 2q_{i,1,k}$$ but $$q_{i,2,k}> \frac{1}{2}q_{i,1,k}$$, the LCP (3.5) has the solution $$c^*_{i,1,k}=4q_{i,1,k}-2q_{i,2,k},\quad c^*_{i,2,k}=4q_{i,2,k}-2q_{i,1,k}; $$ for the case $$q_{i,2,k}\le \frac{1}{2}q_{i,1,k}$$, the LCP (3.5) has the solution $$c^*_{i,1,k}=3q_{i,1,k},\quad c^*_{i,2,k}=0. $$ These three cases can be included in the form (3.4). The solutions $$c^*_{i,j,k}$$ give a global minimizer of $$\|u-u^h\|^2_{\mathscr{L}^2}$$ in $$\mathscr{U}^h_+$$. (2) Denote by $$\tilde u^h$$ the piecewise linear interpolant of $$u$$: $$\tilde u^h(t)=\varphi_1\left(\frac{t_k-t}{h}\right)u(t_{k-1})+\varphi_2\left(\frac{t-t_{k-1}}{h}\right)u(t_{k}) $$ with the estimate $$\|u-\tilde u^h\|_{\mathscr{L}^2}\le C_uh^{2}$$, where $$C_u$$ is a constant independent of $$h$$ (dependent of $$u$$). Obviously, $$\tilde u^h$$ is also non-negative on $$[0,T]$$ and $$\left\|u-\mathfrak{P}^+_hu\right\|_{\mathscr{L}^2}\ \le \ \left\|u-\tilde u^h\right\|_{\mathscr{L}^2}\ \le \ C_uh^{2}. $$ This completes the proof. □ Remark 3.2 Let $$\mathfrak{P}_hu$$ be the orthogonal projection of $$u\in\mathfrak{L}^2_+(0,T;\mathbb{R}^m)$$ onto its subspace $$\mathscr{U}^h$$ with respect to $$\langle\cdot,\cdot\rangle_{\mathscr{L}^2}$$. We can see that $$\left\|u-\max\{\mathfrak{P}_hu,0\}\right\|_{\mathscr{L}^2}\ \le \ \left\|u-\mathfrak{P}^+_hu\right\|_{\mathscr{L}^2}\!. $$ Normally, the function $$\max\{\mathfrak{P}_hu,0\}$$ is not an element of $$\mathscr{U}^h_+$$ but is in a more refined approximation function family. For example, let $$u:[0,1]\to\mathbb{R}$$ with $$u(t)=1$$ when $$t\in\left[0,\frac{1}{2}\right]$$ and $$u(t)=0$$ otherwise. Then one can compute $$\mathfrak{P}_hu=\frac{5}{4}\varphi_1(t)-\frac{1}{4}\varphi_2(t)$$, which takes negative values when $$t\in\left(\frac{5}{6},1\right]$$. Hence, $$\max\{\mathfrak{P}_hu,0\}=\frac{5}{4}\varphi_1(t)-\frac{1}{4}\varphi_2(t)\ge 0$$ in $$\left[0,\frac{5}{6}\right]$$ and $$\max\{\mathfrak{P}_hu,0\}=0$$ in $$\left[\frac{5}{6},1\right]$$. It is not linear on $$[0,1]$$ with the mesh $$0=t_0<t_1=1$$ and $$h=1$$, but is piecewise linear on a refined mesh with $$\frac{5}{6}$$ being a mesh point. Remark 3.3 Suppose that $$u$$ has finitely many discontinuous points (related to mode switches Brogliato (2003)) and suppose that the mesh is so refined that we have just two situations: all the components of $$u$$ do not vanish in the interior of any subintervals except for those, where some components vanish constantly. Then $$\mathfrak{P}^+_hu$$ is the orthogonal projection of $$u\in\mathfrak{L}^2_+(0,T;\mathbb{R}^m)$$ onto the closed subspace $$\mathscr{U}^h$$ with respect to the inner product $$\langle\cdot,\cdot\rangle_{\mathscr{L}^2}$$, and $$\lim_{h\to 0^+}\|u-\mathfrak{P}^+_hu\|_{\mathscr{L}^2}=0$$. 3.2 Discretization of the LCP$$(\mathfrak{L},\hat g)$$ We apply Galerkin approximation to the variational formulation (2.9), which is equivalent to the LCP$$(\mathfrak{L},\hat g)$$. Namely, we find $$u^h=\sum_{l=1}^{d}z_l\varphi_l\in\mathscr{U}^h_+$$ fulfilling (2.9) in $$\mathscr{U}^h$$, where $$\{\varphi_1,\cdots,\varphi_d\}$$ is a basis of $$\mathscr{U}^h$$, $$\varphi_l=\psi_{i,j,k}$$ for $$l=i+m(\,j-1)+2m(k-1)$$, $$1\le i\le m$$, $$j=1,2$$, $$1\le k\le N$$, $$d=2Nm$$, and where $$\psi_{i,j,k}$$ is defined by (3.2). That is, for any $$v^h=\sum_{l=1}^{d}z'_l\varphi_l\in\mathscr{U}^h_+$$, it holds true \begin{align*}0 & \le \displaystyle \frac{1}{h}\left\langle \mathfrak{L}u^h+\hat g,v^h-u^h\right\rangle_{\mathscr{L}^2} \ = \ \frac{1}{h}\left\langle \mathfrak{L}\sum_{i=1}^{d}z_i\varphi_i+\hat g,\sum_{j=1}^{d}(z'_j-z_j)\varphi_j\right\rangle_{\mathscr{L}^2}\\ & = \displaystyle \frac{1}{h}\sum_{i,j=1}^{d} z_i\left(z'_j-z_j\right)\left\langle \mathfrak{L}\varphi_i,\varphi_j\right\rangle_{\mathscr{L}^2}+\frac{1}{h}\sum_{j=1}^{d} \left(z'_j-z_j\right)\left\langle \hat g,\varphi_j \right\rangle_{\mathscr{L}^2}\\ & = \displaystyle \left(z'-z\right)^T\left(M^h z+q^h\right)\!, \end{align*} where $$z=\left(z_i\right)$$, $$z'=\left(z'_i\right)^T$$, $$q^h=\left(q^h_{j}\right)\in\mathbb{R}^d$$ and $$M^h=\left(M^h_{ji}\right)\in\mathbb{R}^{d\times d}$$, \begin{equation}\label{Mh qh}M^h_{ji}:=\frac{1}{h}\left\langle \mathfrak{L}\varphi_i,\varphi_j\right\rangle_{\mathscr{L}^2}, \qquad\quadq^h_{j}:=\frac{1}{h}\left\langle \hat g,\varphi_j\right\rangle_{\mathscr{L}^2}\!. \end{equation} (3.6) We mention that $$u^h\in\mathscr{U}^h_+$$ if and only if $$z\in\mathbb{R}^d_+$$. Therefore, the above variational formulation yields the $$d$$-dimensional LCP$$\left(M^h,q^h\right)$$: find $$z^h\in\mathbb{R}^d$$ such that \begin{equation}\label{LCPh}0\ \le \ z^h\ \bot\ M^hz^h+q^h\ \ge \ 0. \end{equation} (3.7) Suppose that the $$\mathrm{LCP}(\mathfrak{L},\hat g)$$ has a solution $$u\in\mathscr{L}^2(0,T;\mathbb{R}^m)$$. Denote by $$\mathfrak{P}^h_+u$$ the projection of $$u$$ on to $$\mathscr{U}^h_+$$, denote $$r(t)=\min\{u(t),(\mathfrak{L}u)(t)+\hat g(t)\},\quad r^h(t)=\min\left\{(\mathfrak{P}^h_+u)(t),(\mathfrak{L}\mathfrak{P}^h_+u)(t)+\hat g(t)\right\}\!. $$ For fixed $$t$$ we have a matrix $$D={\rm{diag}}(d_i)\in\mathbb{R}^{m\times m}$$ with $$0\le d_i\le 1$$ such that \begin{align*}r(t)-r^h(t) & = \min\{u(t),(\mathfrak{L}u)(t)+\hat g(t)\}-\min\left\{(\mathfrak{P}^h_+u)(t),(\mathfrak{L}\mathfrak{P}^h_+u)(t)+\hat g(t)\right\}\\ & = (I-D)\left(u-\mathfrak{P}^h_+u\right)(t)+D\mathfrak{L}\left(u-\mathfrak{P}^h_+u\right)(t). \end{align*} See Alefeld et al. (1999). Noting $$\|I-D\|_2\le 1$$ and $$\|D\|_2\le 1$$, we have \begin{align*}\left\|r(t)-r^h(t)\right\|^2_2 & = \left\|(I-D)(u(t)-(\mathfrak{P}^h_+u)(t))+D\mathfrak{L}(u-\mathfrak{P}^h_+u)(t)\right\|_2^2\\ & \le \left(\left\|I-D\right\|_2\left\|u(t)-\mathfrak{P}^h_+u(t)\right\|_2+\left\|D\right\|_2\left\|\mathfrak{L}(u-\mathfrak{P}^h_+u)(t)\right\|_2\right)^2\\ & \le 2\left\|u(t)-\mathfrak{P}^h_+u(t)\right\|^2_2+2\left\|\mathfrak{L}(u-\mathfrak{P}^h_+u)(t)\right\|^2_2. \end{align*} Since $$u\in\mathscr{L}^2(0,T;\mathbb{R}^m)$$ solves the $$\mathrm{LCP}(\mathfrak{L},\hat g)$$, it holds true $$r(t)=0$$ a.e. over $$[0,T]$$. Therefore, $$\left\|r^h\right\|^2_{\mathscr{L}^2} \ = \ \int^T_0\left\|r(t)-r^h(t)\right\|^2_2\,{\rm{d}}t\ \le \ 2\left\|u-\mathfrak{P}^h_+u\right\|^2_{\mathscr{L}^2}+2\left\|\mathfrak{L}(u-\mathfrak{P}^h_+u)\right\|^2_{\mathscr{L}^2}\!. $$ Notice that $$\mathfrak{L}$$ is bounded and $$\left\|u-\mathfrak{P}^h_+u\right\|^2_{\mathscr{L}^2}\to 0$$ as $$h\downarrow 0$$. Then $$\left\|r^h\right\|_{\mathscr{L}^2}\to 0$$ as $$h\downarrow 0$$. We know $$r^h=0$$ a.e. over $$[0,T]$$ when $$\left\|r^h\right\|_{\mathscr{L}^2}=0$$. For this reason, we can say $$\mathfrak{P}^h_+u$$ approximately solves the Galerkin approximation problem (3.7), and it justifies in a certain sense the Galerkin approximation in the LCS setting. Of course, a very small measure $$\left\|r^h\right\|_{\mathscr{L}^2}$$ does not imply the solvability of (3.7). We study it below. Theorem 3.4 Let $$M$$ be positive semidefinite, $$f\in\mathscr{L}^2(0,T;\mathbb{R}^n)$$ and $$g\in\mathscr{L}^2(0,T;\mathbb{R}^m)$$. Then under condition (2.10), the LCP (3.7) has a solution. Proof. Fix the mesh and let $$\mathfrak{P}^h_+u_0$$ be the projection of $$u_0$$ onto $$\mathscr{U}_+^h$$. Let $$u\in\mathscr{U}^h_+$$ with $$\|u\|_{\mathscr{L}^2}\to \infty$$. Then from (2.10) it follows that \begin{align*}+\infty\ \leftarrow\ \displaystyle \frac{\left\langle \mathfrak{L}u,u-u_0\right\rangle_{\mathscr{L}^2}}{\|u\|_{\mathscr{L}^2}}& = \displaystyle\frac{\left\langle \mathfrak{L}u,u-\mathfrak{P}^h_+u_0\right\rangle_{\mathscr{L}^2}}{\|u\|_{\mathscr{L}^2}}+\frac{\left\langle \mathfrak{L}u,\mathfrak{P}^h_+u_0-u_0\right\rangle_{\mathscr{L}^2}}{\|u\|_{\mathscr{L}^2}}\\ & \le \displaystyle\frac{\left\langle \mathfrak{L}u,u-\mathfrak{P}^h_+u_0\right\rangle_{\mathscr{L}^2}}{\|u\|_{\mathscr{L}^2}}+\left\|\mathfrak{L}\right\|_{\mathscr{L}^2}\left\|\mathfrak{P}^h_+u_0-u_0\right\|_{\mathscr{L}^2}\!. \end{align*} This follows that condition (2.10) holds true in $$\mathscr{U}^2_+$$ since $$\|\mathfrak{L}\|_{\mathscr{L}^2} \left\|\mathfrak{P}^h_+u_0-u_0\right\|_{\mathscr{L}^2}$$ is finite. The conclusion is drawn again by Theorem 32.C of Zeidler (1990) (II/B, p. 875). □ The solution $$z^h=\left(z^h_i\right)$$ of the $$\mathrm{LCP}(M^h,q^h)$$ gives a piecewise (discontinuous) linear function on the mesh, which has the following form on the subinterval $$(t_{k-1},t_k]$$: \begin{equation}\label{uh}u^h(t)\ = \ \left(\frac{t_k-t}{h}\right)u^h_{k,1}+\left(\frac{t-t_{k-1}}{h}\right) u^h_{k,1}, \end{equation} (3.8) where $$u^h_{k,1} = \left(\begin{array}{c}z^h_{m(2k-2)+1}\\ \vdots \\ z^h_{m(2k-2)+m}\end{array}\right)\!,\quad u^h_{k,1} =\left(\begin{array}{c}z^h_{m(2k-1)+1}\\ \vdots \\ z^h_{m(2k-1)+m}\end{array}\right)\!. $$ In subsequence, $$u^h$$ is just called as the solution of the LCP$$(M^h,q^h)$$ if no ambiguity caused. 4. Convergence analysis We study the convergence of the Galerkin approximation (3.7). We know that problem (2.9) and its Galerkin approximation (3.7) have the unique solution if the operator $$\mathfrak{L}$$ is strongly monotone, which is true if, for example, $$M$$ is positive definite and $$T$$ is small enough. For the case of strong monotonicity, we prove the following result of the convergence rate by adapting the technique of Falk (1974). Theorem 4.1 Let $$\mathfrak{L}$$ be strongly monotone: $$\langle \mathfrak{L}v,v\rangle_{\mathscr{L}^2}\ge\alpha \|v\|^2_{\mathscr{L}^2}$$ holding with a constant $$\alpha>0$$ for any $$v\in\mathscr{L}^2(0,T;\mathbb{R}^m)$$. Let $$u$$ and $$u^h$$ be the unique solutions of the $$\mathrm{LCP}(\mathfrak{L},\hat g)$$ and the $$\mathrm{LCP}(M^h,q^h)$$, respectively. Then \begin{equation}\label{estimate}\|u^h-u\|_{\mathscr{L}^2}\le C\sqrt{ \inf_{v^h\in\mathscr{U}^h_+}\|v^h-u\|_{\mathscr{L}^2}}\ = \ C \sqrt{\left\|\mathfrak{P}_h^+u-u\right\|_{\mathscr{L}^2}}. \end{equation} (4.1) If $$u\in\mathscr{C}^2(0,T;\mathbb{R}^m)$$, the family of functions that have continuous second-order derivatives, then $$\left\|u^h-u\right\|_{\mathscr{L}^2}=O(h)$$. Remark 4.2 The estimate (4.1) holds for any approximation of $$\mathscr{L}^2_+(0,T;\mathbb{R}^m)$$ that is convex and closed, for example, the space of piecewise constant functions which are non-negative over $$[0,T]$$. Proof. Since $$u$$ solves (2.9) and $$u^h$$ solves (3.7), we have $$\langle\mathfrak{L}u+\hat g,u\rangle_{\mathscr{L}^2}=0$$, $$\langle\mathfrak{L}u+\hat g,u^h\rangle_{\mathscr{L}^2}\ge 0$$ and $$\left\langle\mathfrak{L}u^h+\hat g,v^h\right\rangle_{\mathscr{L}^2}\ge 0$$ for any $$v^h\in\mathscr{U}^h_+$$. These yield \begin{align*}\left\langle\mathfrak{L}(u-u^h),u-u^h\right\rangle_{\mathscr{L}^2}& = \left\langle(\mathfrak{L}u+\hat g)-\left(\mathfrak{L}u^h+\hat g\right)\!,u-u^h\right\rangle_{\mathscr{L}^2}\\ & = -\left\langle\mathfrak{L}u+\hat g,u^h\right\rangle_{\mathscr{L}^2}-\left\langle\mathfrak{L}u^h+\hat g,u\right\rangle_{\mathscr{L}^2}\\ &\le \left\langle\mathfrak{L}u^h+\hat g,v^h\right\rangle_{\mathscr{L}^2}-\left\langle\mathfrak{L}u^h+\hat g,u\right\rangle_{\mathscr{L}^2}\\ & = \left\langle\mathfrak{L}\left(u^h-u\right)\!,v^h-u\right\rangle_{\mathscr{L}^2}+\left\langle\mathfrak{L}u+\hat g,v^h-u\right\rangle_{\mathscr{L}^2}\\ & \le \left\|\mathfrak{L}\right\|_{\mathscr{L}^2}\left\|u^h-u\right\|_{\mathscr{L}^2}\left\|v^h-u\right\|_{\mathscr{L}^2}+ \left\|\mathfrak{L}u+\hat g\right\|_{\mathscr{L}^2}\left\|v^h-u\right\|_{\mathscr{L}^2}\!. \end{align*} Using the condition $$\alpha\left\|u-u^h\right\|_{\mathscr{L}^2} \le\left\langle\mathfrak{L}\left(u-u^h\right)\!,u-u^h\right\rangle_{\mathscr{L}^2}$$ and the inequality $$\left\|\mathfrak{L}\right\|_{\mathscr{L}^2}\left\|u^h-u\right\|_{\mathscr{L}^2}\left\|v^h-u\right\|_{\mathscr{L}^2}\ \le \ \frac{\alpha}{2}\left\|u^h-u\right\|^2_{\mathscr{L}^2}+\frac{\left\|\mathfrak{L}\right\|^2_{\mathscr{L}^2}}{2\alpha}\left\|v^h-u\right\|^2_{\mathscr{L}^2}\!, $$ we obtain $$\frac{\alpha}{2}\left\|u^h-u\right\|^2_{\mathscr{L}^2} \ \le \ \frac{\left\|\mathfrak{L}\right\|^2_{\mathscr{L}^2}}{2\alpha}\left\|v^h-u\right\|^2_{\mathscr{L}^2}+\left\|\mathfrak{L}u+\hat g\right\|_{\mathscr{L}^2}\left\|v^h-u\right\|_{\mathscr{L}^2}\!. $$ This follows (4.1) since $$v^h\in\mathscr{U}^h_+$$ is arbitrary. We take $$v^h$$ as the piecewise linear interpolant of $$u$$, which, obviously, is an element of $$\mathscr{U}^h_+$$ with $$\|u-\tilde u^h\|_{\mathscr{L}^2}\le C_uh^{2}$$ if $$u\in\mathscr{C}_+$$, where $$C_u$$ is a constant independent of $$h$$ (dependent of $$u$$). This completes the proof. □ Theorem 4.3 Assume that the operator $$\mathfrak{L}$$ is monotone, and let $$u^h$$ be the solution of the $$\mathrm{LCP}(M^h,q^h)$$. If $$\{u^h\}$$ is uniformly bounded for $$h$$ small enough then $$\{u^h\}$$ has a subsequence, which weakly converges to a solution of the $$\mathrm{LCP}(\mathfrak{L},\hat g)$$. Proof. If $$\{u^h\}$$ is uniformly bounded then $$\{u^h\}$$ has a subsequence (still denoted by $$\{u^h\}$$ for avoiding the cumbersome presentation), converging to $$u\in\mathscr{L}^2_+(0,T;\mathbb{R}^m)$$ weakly. Now we prove that $$u$$ is a solution of (2.9). Note that $$\mathscr{C}^2_+$$ is the set of functions that have continuous second-order derivatives and are non-negative in $$[0,T]$$. Since $$\mathscr{C}^2_+$$ is dense in $$\mathscr{L}^2_+(0,T;\mathbb{R}^m)$$, it is enough to show $$\langle\mathfrak{L}u+\hat g,v-u\rangle_{\mathscr{L}^2}\ge 0$$ for any $$v\in\mathscr{C}^2_+$$. Let $$v\in\mathscr{C}^2_+$$ and take $$v^h$$ as its piecewise linear interpolant. Then $$v^h$$ is strongly convergent to $$v$$ in $$\|\cdot\|_{\mathscr{L}^2}$$ as $$h\downarrow 0$$. Since $$u^h$$ is weakly convergent to $$u$$, we have $$\left \langle\mathfrak{L}u^h,v^h\right\rangle_{\mathscr{L}^2}\to \left\langle\mathfrak{L}u,v\right\rangle_{\mathscr{L}^2},\qquad \left \langle \hat g,v^h-u^h\right\rangle_{\mathscr{L}^2}\to \left\langle \hat g,v-u\right\rangle_{\mathscr{L}^2}\!. $$ Therefore, $$\left\langle\mathfrak{L}u^h,v^h\right\rangle_{\mathscr{L}^2}+\left\langle \hat g,v^h-u^h\right\rangle_{\mathscr{L}^2}\to \left\langle\mathfrak{L}u,v\right\rangle_{\mathscr{L}^2}+\left\langle \hat g,v-u\right\rangle_{\mathscr{L}^2}$$. And we obtain $$0\ \le \ \left\langle \mathfrak{L}u^h,u^h\right\rangle_{\mathscr{L}^2}\le \left\langle \mathfrak{L}u^h,v^h\right\rangle_{\mathscr{L}^2}+ \left\langle \hat g,v^h-u^h\right\rangle_{\mathscr{L}^2}$$ since $$\left\langle \mathfrak{L}u^h+\hat g,v^h-u^h\right\rangle_{\mathscr{L}^2}\ge 0$$. This follows $$\left\langle\mathfrak{L}u+\hat g,v-u\right\rangle_{\mathscr{L}^2}\ge 0$$ for any $$v\in\mathscr{C}^2_+$$. □ 5. Numerical experiment 5.1 Implementation details Let $$M^h$$ and $$q^h$$ be defined as in (3.6). Below we give the representation of the matrix $$M^h$$ and the column vector $$q^h$$, where the linear interpolants $$f^h$$ and $$g^h$$ of $$f$$ and $$g$$ are used as their approximations. Below we denote by $$f_k=f(t_k)$$ and $$g_k=g(t_k)$$ for $$k=0,\dots, N$$. The representations need the following matrices $$G_k=G_k(hA)$$, where $$h>0$$ is the stepsize. Denote $$G_0=e^{hA}$$ and $$G_k=\sum_{n=k}^\infty\frac{1}{n!}(hA)^{n-k}=\frac{1}{k!}I+\frac{1}{(k+1)!}(hA)+\frac{1}{(k+2)!}(hA)^2+\cdots. $$ The matrices $$G_k$$ can be computed efficiently by using, e.g., the Pade and Krylov subspace approximations. Here we only need $$G_k$$ for $$0\le k\le 4$$, which can be evaluated in high precision by the existing software, e.g., the matrix function toolbox (Higham, http://www.ma.man.ac.uk/higham/mftoolbox). Arguments on other software for computing the matrix functions can be found in Higham & Al-Mohy (2010). Denote by $$\otimes$$ the Kronecker tensor product of two matrices. Then $$M^h$$ has the forms: \begin{align*}M^h & = \displaystyle \frac{1}{6} I_N\otimes \left(\begin{array}{@{}rr@{}}2M & M\\ M & 2M \end{array}\right)-h\sum_{k',k=1}^Ne_{k'}\,e_k^T \otimes \left(\begin{array}{@{}cc@{}}(J_{k',k})_{11} & (J_{k',k})_{12}\\ (J_{k',k})_{21} & (J_{k',k})_{22}\end{array}\right)\\ & \quad \displaystyle +\frac{h}{2}\sum_{k'=2}^N \sum_{k=1}^{k'-1}e_{k'}\,e_k^T \otimes\left(\begin{array}{@{}cc@{}}Q(G_1-G_2) B &QG_2B\\ Q(G_1-G_2) B & QG_2B\end{array}\right)\\ & \quad \displaystyle +h I_N\otimes \left(\begin{array}{@{}cc@{}}Q(G_3-G_4)B & QG_4B \\ Q(G_2-2G_3+G_4)B & Q(G_3-G_4)B \end{array}\right)\!, \end{align*} where $$e_k$$ is the $$k$$th $$N$$-dimensional unit coordinate vector, and for $$k', k=1,\dots,N$$: \begin{align*}(J_{k',k})_{11} & = QG^{k'-1}_0G_2\left(E_0+E_T\,e^{TA}\right)^{-1}E_TG_0^{N-k} (G_1-G_2)B \\ (J_{k',k})_{12} & = QG^{k'-1}_0G_2\left(E_0+E_T\,e^{TA}\right)^{-1}E_T G_0^{N-k} G_2B\\ (J_{k',k})_{21} & = QG^{k'-1}_0(G_1-G_2)\left(E_0+E_T\,e^{TA}\right)^{-1}E_T G_0^{N-k}(G_1-G_2)B\\ (J_{k',k})_{22} & = QG^{k'-1}_0(G_1-G_2)\left(E_0+E_T\,e^{TA}\right)^{-1}E_T G_0^{N-k}G_2B. \end{align*} And we compute $$q^h$$ in a naive manner. Given the approximate control $$u^h$$, we present the following formula for approximating the state $$x$$ at the mesh points by $$x^h_k=x^h(t_k)$$ for $$k=0,1,\ldots N$$: \begin{align*}x^h_k & = \displaystyle G_0^{k}\hat x^{0,h}+ \int^{t_k}_0e^{(t_k-s)A}f(s)\,{\rm{d}}s+h\left(G_1B\sum_{j=1}^k u^h_{j,1}+G_2B\sum_{j=1}^k \left(u^h_{j,2}-u^h_{j,1}\right)\right)\\ & \quad \displaystyle -h G_0^k\left(E_0+E_T\,e^{TA}\right)^{-1}E_T \sum_{k=1}^NG_0^{N-k}\left[(G_1-G_2)Bu^h_{k,1}+G_2Bu^h_{k,2}\right]\!, \end{align*} where $$\hat x^{0,h} \ = \ \left(E_0+E_T\,e^{TA}\right)^{-1}\left(b-E_T\int^T_0 e^{(T-s)A}f(s)\,{\rm{d}}s\right)\!. $$ Here we omit the proof of the form of $$M^h$$ and the justification of the approximation of the state, which can be found in the supplementary material (Wang & Chen, 2017). 5.2 Numerical results In this subsection, we apply the Galerkin approximation method and the time-stepping method to two LCSs, which generate, respectively, the numerical solutions $$\left(x^h_{\rm g},u^h_{\rm g}\right)$$ and $$\left(x^h_{\rm e},u^h_{\rm e}\right)$$, where the subscript ‘e’ stands for ‘Euler’ since the time-stepping method actually makes the use of implicit Euler method to discretize the ODEs involved in the LCSs. Here $$x^h_{\rm g}$$ is recovered in the manner stated above by using a solution $$u^h_{\rm g}$$ of the $$\mathrm{LCP}\left(M^h,q^h\right)$$. The exact solution of the LCS is always denoted by $$(x,u)$$. Here we are interested in the errors of state and control, in $$\|\cdot\|_{\mathscr{L}^2}$$ and in $$\|\cdot\|_2$$, namely we compute the values of $$\begin{array}{llll}\left\|x^h_{\rm g}-x\right\|_{\mathscr{L}^2},& \left\|u^h_{\rm g}-u\right\|_{\mathscr{L}^2},& \left\|x^h_{\rm g}(T)-x(T)\right\|_2,& \left\|u^h_{\rm g}(T)-u(T)\right\|_2\!,\\ \left\|x^h_{\rm e}-x\right\|_{\mathscr{L}^2},& \left\|u^h_{\rm e}-u\right\|_{\mathscr{L}^2},& \left\|x^h_{\rm e}(T)-x(T)\right\|_2,& \left\|u^h_{\rm e}(T)-u(T)\right\|_2\!. \end{array}$$ If, for example, $$\log\left(\left\|x^h_{\rm g}-x\right\|_{\mathscr{L}^2}\right)$$ is affine w.r.t. $$\log(h)$$ then the slope gives the order of the state convergence in $$\|\cdot\|_{\mathscr{L}^2}$$. Therefore, for the two examples we report, the logarithms of the errors in different $$h$$. Here we take the uniform grid, the $$\mathrm{LCP}\left(M^h,q^h\right)$$ is solved by using the PATH LCP solver Dirkse et al., numerical methods are coded and performed in the setting of Octave 4.0. Below we present the details of the numerical example and thereafter present the comments on the numerical results. Example 5.1 The collapse of the Tacoma Narrows suspension bridge can be modeled by ODEs of second order with nonsmooth data. The nonsmooth data are reformulated by an LCP, which leads to an IVP of the LCS. The data of the problem and the exact solution can be found in Chen & Mahmoud (2008). This problem can be reformulated as an $$\mathrm{LCP}(\mathfrak{L},\hat g)$$ with a strongly monotone $$\mathfrak{L}$$ and therefore has a unique solution $$(x,u)$$. Below we present some remarks on the numerical results. (1) Take $$h\approx 0.1$$. The first components of the exact state $$x$$ and their approximations $$\left(x^h_{\rm g}\right)_1$$ and $$\left(x^h_{\rm e}\right)_1$$ are plotted in the upper part of Fig. 1, the counterpart results for the control are plotted in the lower part. The step size is not restrictive, while the numerical solutions offered by Galerkin approximation are very close to the exact one and much more precise than the output of the time-stepping method. Fig. 1. View largeDownload slide Exact and numerical solutions of the LCS (Example 5.1 with $$T=3\pi$$). Fig. 1. View largeDownload slide Exact and numerical solutions of the LCS (Example 5.1 with $$T=3\pi$$). (2) The logarithms of the errors for the Galerkin approximation are plotted in Fig. 2, w.r.t. the logarithms of different $$h$$. Note that the graphs are approximately four straight lines with the slopes all about of 2. This strongly indicates that the numerical solutions given by the Galerkin approximation method have a convergence of order 2, both in norm $$\|\cdot\|_{\mathscr{L}^2}$$ and in $$\|\cdot\|_2$$ at the terminal time. Notice that Theorem 4.1 just gives the first-order convergence. (3) The logarithms of the errors $$\left\|x^h_{\rm e}(T)-x(T)\right\|_2$$ are plotted in the upper part of Fig. 3. The graphs are approximately 2 straight lines with the slopes about of 9.2334e$$-$$01. This indicates a convergence of the time-stepping method in $$\|\cdot\|_2$$ at $$T$$ of the order close to 1. As illustrated by the numerical results, the Galerkin approximation has a much better precision than the time-stepping method. For instance, for $$h\approx 10^{-3}$$, we have $$\left\|x^h_{\rm e}(T)-x(T)\right\|_2/\left\|x^h_{\rm g}(T)-x(T)\right\|_2\approx 2.3565\times 10^{3}$$. The overperformance of the Galerkin approximation in $$\|\cdot \|_{\mathscr{L}^2}$$ is more obvious. In the same setting as above, our numerical results show $$\frac{\left\|x^h_{\rm e}-x\right\|_{\mathscr{L}^2}}{\left\|x^h_{\rm g}-x\right\|_{\mathscr{L}^2}}\approx 2.2168\times 10^{7},\qquad \frac{\left\|u^h_{\rm e}-u\right\|_{\mathscr{L}^2}}{\left\|u^h_{\rm g}-u\right\|_{\mathscr{L}^2}}\approx 2.1610\times 10^{6}. $$ (4) Theorem 4.1 indicates that the approximation error for the Galerkin approximation is bounded from above by $$\left\|\mathfrak{P}_h^+u-u\right\|_{\mathscr{L}^2}$$, where $$\mathfrak{P}_h^+u$$ is the projection of the true solution $$u$$ onto the function family $$\mathscr{U}^h_+$$. Here we plot in the lower part of Fig. 3, the logarithms of the errors of the projection onto three different function families: the family of piecewise constant, continuous and discontinuous piecewise linear functions that are non-negative in $$[0,T]$$. The numerical results explain in another manner the reason that the Galerkin approximation numerically overperforms the time-stepping method (approximating $$u$$ by a piecewise constant function) and justifies the application of the discontinuous Galerkin approximation instead of the continuous one (approximating $$u$$ by a continuous piecewise linear function). Note that in the current case, the projection error for discontinuous piecewise linear function family is very small and not sensitive in $$h$$, it means that a high accuracy of our discontinuous Galerkin approximation can be achieved for an $$h$$ not restrictive. We mention that the projection error is dependent of the geometry of the true solution $$u$$. If $$u$$ is piecewise linear and its discontinuities are located at the grid points then the error is zero. Example 5.2 Consider a dynamic Nash equilibrium problem with 2 players and zero-sum cost functionals. Denote by $$y_i\in\mathbb{R}$$ and $$u_i$$ the $$i$$th player’s state and control variables, respectively, $$i=1,2$$. An equilibrium solution of this problem is a state-control pair $$(y^*_1(t),y^*_2(t), u_1^*(t),u_2^*(t))$$, satisfying \begin{align*}\left.\begin{array}{r@{\ }c@{\ }rl}u_1^*(\cdot) & \in & \mathrm{argmin} & \theta(y_1,y^*_2,u_1,u^*_2)\\ & & \mathrm{s.t.} & \dot{y}_1(t)=-2+2y_1+B_1u_1\\ & & & y_1(0)=-1 \\ & & & u_1 \ge 0,\, e^T(u_1+u^*_2)\le 1 \end{array}\right.\\ \left.\begin{array}{r@{\ }c@{\ }rl}u_2^*(\cdot) & \in & \mathrm{argmax} & \theta(y^*_1,y_2,u_1^*,u_2)\\ & & \mathrm{s.t.} & \dot{y}_2(t)=-2t-y_2+B_2u_2\\ & & & y_2(0)=2 \\ & & & u_2\ge 0,\, e^T(u_1^*+u_2)\le 1, \end{array}\right. \end{align*} where $$e=(1,1)^T$$. Denote $$y=(y_1,y_2)^{T}$$ and $$u=(u_1^T,u_2^T)^T$$. The cost functional $$\theta(y_1,y_2,u_1,u_2)=\theta(y,u) $$ reads $$y(T)^TLy(T)+l^Ty(T)+\int^T_0 y^T\left[Py+Su+h(t)\right]+u^T\left[Ru+d(t)\right]\,{\rm{d}}t. $$ Here the matrices $$B_1$$, $$B_2$$, $$L$$, $$l$$, $$P$$, $$S$$ and $$R$$ are taken as in Chen & Wang (2014). The Pontryagin’s minimum/maximum principle for each optimal control problem yields two coupled constrained Hamilton equations. By the vector $$w\in\mathbb{R}^5$$ of multipliers, it gives the following coupled system of ODE and a mixed LCP (a little different to that of (1.1)): \begin{align*}\dot{x}(t) & = Ax(t)+Bu(t)+f(t)\\ 0\,& = Qx(t)+Mu(t)-C^Tw(t)+g(t)\\ 0\, & \le w(t)\, \bot\, Cu(t)-c \, \ge 0\\ b\ & = E_0x(0)+E_Tx(T). \end{align*} Here $$x=(y_1,p_1,y_2,p_2)^T$$, $$p_i$$ is the costate of $$y_i$$. The detailed reformulation, and the data of the matrices $$A$$, $$B$$, $$Q$$, $$M$$, $$C$$, $$E_0$$, $$E_T$$ and the vectors $$b$$ and $$c$$ can be found in Chen & Wang (2014). Here we construct an exact solution $$(x,u)$$ of the LCS by adapting the functions $$f$$ and $$g$$. We mention that both the matrix $$M$$ and the one defining the mixed LCP are positive semidefinite. For this problem, our numerical method can be established in a very similar manner by using the variational formulation $$\left\langle \hat g+\mathfrak{L}u, v-u\right\rangle_{\mathscr{L}^2} \ge\ 0\qquad \forall v\in \mathscr{L}^2(0,T;\mathbb{R}^4)\times \mathscr{L}_+^2(0,T;\mathbb{R}^5). $$ The Galerkin approximation yields a mixed LCP of dimension $$18N$$, where $$N$$ is the number of the subintervals. We have the following observation on the numerical results. (1) The logarithms of the errors for the Galerkin approximation and the time-stepping method are plotted in Figs 4 and 5, respectively. Note that the four graphs in Fig. 4 are approximately straight lines with the slopes all about $$1.0076$$, this indicates the first-order convergence of the Galerkin approximation. For this example, the precision of the Galerkin approximation is still much better than the time-stepping method, whose convergence order is close to 0, as illustrated by the figures. Fig. 4. View largeDownload slide Errors for Galerkin approximation (Example 5.2 with $$T=1$$). Fig. 4. View largeDownload slide Errors for Galerkin approximation (Example 5.2 with $$T=1$$). Fig. 5. View largeDownload slide Errors for time-stepping method (Example 5.2 with $$T=1$$). Fig. 5. View largeDownload slide Errors for time-stepping method (Example 5.2 with $$T=1$$). Note that $$\mathfrak{L}$$ is not monotone, Theorem 4.1 cannot be applied to provide error estimate and convergence order. (2) If $$\mathfrak{L}$$ is not monotone then the $$\mathrm{LCP}(M^h,q^h)$$ may have no solution or have multiple solutions, for which the LCP solvers may not perform well, their output could be far away from the true solution. This leads to the low order of convergence or even divergence. Even when $$M$$ is positive definite, $$\mathfrak{L}$$ is not necessarily monotone. The monotonicity is also dependent of $$T$$. For illustrating this, we consider the regularized matrix $$M+\lambda I$$ with $$\lambda=0.1$$, which is positive definite. It defines the operator $$\mathfrak{L}_\lambda$$ and then gives the matrix $$M^h_\lambda$$. If $$\mathfrak{L}_\lambda$$ is (strongly) monotone then $$M^h_\lambda$$ is positive (definite) semidefinite, and so its symmetric part has a non-negative (positive) smallest eigenvalue $$\sigma(T,h)$$. We plot the values of $$\sigma(T,h)$$ in Fig. 6 for different $$T$$ and $$h$$. Note that for $$T=1, 0.5$$, $$M^h_\lambda$$ is indefinite, while positive definite for a smaller $$T$$. Fig. 6. View largeDownload slide Smallest eigenvalue $$\sigma(h)$$ of the symmetric part of $$M^h$$ (Example 5.2). Fig. 6. View largeDownload slide Smallest eigenvalue $$\sigma(h)$$ of the symmetric part of $$M^h$$ (Example 5.2). 6. Concluding remarks This article reformulates the LCS into an LCP in an Hilbert space and proposes a numerical method by using discontinuous Galerkin approximation. The method solves a finite-dimensional $$\mathrm{LCP}(M^h,q^h)$$ and uses the matrix exponential related matrices $$\varphi_k(hA)$$ to recover the state, and these matrices are also used to construct the data of the $$\mathrm{LCP}\left(M^h,q^h\right)$$, which keeps a quite good fidelity to the LCS. Numerical results show that the accuracy of the proposed method is much better than the time-stepping method. Our method and the time-stepping method need to evaluate $$\varphi_k(hA)y$$ and $$(I-hA)^{-1}y,$$ respectively, for some column vector $$y$$. The evaluation of the former is not necessarily time-consuming than that of the latter (see Higham & Al-Mohy, 2010). Funding Fundamental Research Funds for the Central Universities (No. 14380011), by National Natural Science Foundation of China (No. 11571166), and by the Priority Academic Program Development of Jiangsu Higher Education Institutions to Z.W. in part; Hong Kong Research Grants Council (PolyU153001/14P) to X.C. in part. Footnotes 1 We call the convergence is of order $$p>0$$ if $$\left\|u-u^h\right\|= O(h^p)$$. References Alefeld G. E. , Chen X. & Potra F. A. ( 1999 ) Numerical validation of solutions of linear complementarity problems . Numer. Math. , 83 , 1 – 23 . Google Scholar Crossref Search ADS Brogliato B. ( 2003 ) Some perspective on analysis and control of complementarity systems . IEEE Trans. Autom. Control , 48 , 918 – 935 . Google Scholar Crossref Search ADS Chen X. & Mahmoud S. ( 2008 ) Implicit Runge-Kutta methods for Lipschitz continuous ordinary differential equations . SIAM J. Numer. Anal. , 46 , 266 – 280 . Chen X. , Nashed Z. & Qi L. ( 1997 ) Convergence of Newton’s method for singular smooth and nonsmooth equations using adaptive outer inverses . SIAM J. Optim. , 7 , 445 – 462 . Google Scholar Crossref Search ADS Chen X. & Wang Z. ( 2011 ) Error bounds for a differential linear variational inequality . IMA J. Numer. Anal. , 32 , 957 – 982 . Google Scholar Crossref Search ADS Chen X. & Wang Z. ( 2013 ) Convergence of regularized time-stepping methods for differential variational inequalities . SIAM J. Optim. , 23 , 1647 – 1671 . Google Scholar Crossref Search ADS Chen X. & Wang Z. ( 2014 ) Differential variational inequality approach to dynamic games with shared constraints . Math. Program. , 146 , 379 – 408 . Google Scholar Crossref Search ADS Conway J. B. ( 1985 ) A Course in Functional Analysis . Berlin : Springer . Dirkse S. , Ferris M. & Munson T. ( 2012 ) The PATH solver . Available at http://pages.cs.wisc.edu/ferris/path.html. Accessed 1 September 2017 . Falk R. S. ( 1974 ) Error estimates for the approximation of a class of variational inequalities . Math. Comput. , 28 , 963 – 971 . Google Scholar Crossref Search ADS Gwinner J. ( 2013 ) On a new class of differential variational inequalities and a stability result . Math. Program. , 139 , 205 – 221 . Google Scholar Crossref Search ADS Han L. , Camlibel M. K. , Pang J.-S. & Heemels W. P. M. H. ( 2012 ) A unified numerical scheme for linear-quadratic optimal control problems with joint control and state constraints . Optim. Meth. Softw. , 27 , 761 – 799 . Google Scholar Crossref Search ADS Heemels W. P. M. H. , Schumacher J. M. & Weiland S. ( 2000 ) Linear complementarity systems . SIAM J. Appl. Math. , 60 , 1234 – 1269 . Google Scholar Crossref Search ADS Higham N. J. & Al-Mohy A. H. ( 2010 ) Computing matrix functions . Acta Numer. , 19 , 159 – 208 . Google Scholar Crossref Search ADS Hochbruck M. & Ostermann A. ( 2010 ) Exponential integrators . Acta Numer. , 19 , 209 – 286 . Google Scholar Crossref Search ADS Kunkel P. & Stöver R. ( 2002 ) Symmetric collocation methods for linear differential-algebraic boundary value problems . Numer. Math. , 91 , 475 – 501 . Google Scholar Crossref Search ADS Pang J.-S. & Qi L. ( 1993 ) Nonsmooth equations: motivation and algorithms . SIAM J. Optim. , 3 , 443 – 465 . Google Scholar Crossref Search ADS Pang J.-S. & Stewart D. E. ( 2008 ) Differential variational inequalities . Math. Program. , 113 , 345 – 424 . Google Scholar Crossref Search ADS Stewart D. E. ( 2006 ) Convolution complementarity problems with application to impact problems . IMA J. Appl. Math. , 71 , 92 – 119 . Google Scholar Crossref Search ADS Wang Z. & Chen X. ( 2017 ) Discretized form of the operator $$\mathfrak{L}$$ and state recovery , supplement material available at http://www.imammb.oxfordjournals.org/. Zeidler E. ( 1990 ) Nonlinear Functional Analysis and Its Applications . Leipzig : Teubner . Zhong R. X. , Sumalee A. , Friesz T. L. & Lam W. H. K. ( 2011 ) Dynamic user equilibrium with side constraints for a traffic network: theoretical development and numerical solution algorithm . Transp. Res. B , 45 , 1035 – 1061 . Google Scholar Crossref Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
Journal
IMA Journal of Numerical Analysis – Oxford University Press
Published: Oct 16, 2018