Exponential integrators for large-scale stiff matrix Riccati   differential equations

Dongping Li

arXiv:1907.12971·math.NA·August 20, 2019

Exponential integrators for large-scale stiff matrix Riccati differential equations

Dongping Li

PDF

Open Access

TL;DR

This paper introduces exponential Rosenbrock-type integrators for efficiently solving large-scale stiff matrix Riccati differential equations, demonstrating high accuracy and computational efficiency through numerical comparisons.

Contribution

It presents novel application of exponential integrators to large-scale stiff matrix Riccati equations, including implementation strategies and low-rank approximation techniques.

Findings

01

Exponential integrators achieve high accuracy for large-scale problems.

02

The methods are computationally efficient compared to traditional approaches.

03

Numerical results confirm the effectiveness of the proposed schemes.

Abstract

Matrix Riccati differential equations arise in many different areas and are particular important within the field of control theory. In this paper we consider numerical integration for large-scale systems of stiff matrix Riccati differential equations. We show how to apply exponential Rosenbrock-type integrators to get approximate solutions. Two typical exponential integration schemes are considered. The implementation issues are addressed and some low-rank approximations are exploited based on high quality numerical algebra codes. Numerical comparisons demonstrate that the exponential integrators can obtain high accuracy and efficiency for solving large-scale systems of stiff matrix Riccati differential equations.

Tables3

Table 1. Table 1: CPU time in seconds and relative errors for each of the methods over the time interval [0, 1]

$m a t r i x$	size(A)	GExpEuler		LrExpEuler		BrExpEuler		Erow3		Additive4		ode15s
		Error	Time	Error	Time	Error	Time	Error	Time	Error	Time	Error	Time
fdm-sym	64 $\times$ 64	1.22e-14	0.95	1.31e-14	3.36	4.58e-14	1.56	1.30e-14	4.25	2.09e-04	2.67	1.78e-14	334.81
fdm-sym	100 $\times$ 100	1.57e-14	1.81	1.73e-14	6.28	4.46e-13	1.82	1.77e-14	7.51	8.53e-04	3.07	2.29e-14	3708.20
fdm-nonsym	64 $\times$ 64	2.01e-14	0.96	2.16e-14	2.12	8.61e-14	2.12	2.15e-14	3.11	5.17e-04	4.02	2.89e-14	423.85
fdm-nonsym	100 $\times$ 100	2.26e-14	2.16	2.78e-14	12.67	3.21e-14	2.58	2.79e-14	15.37	1.89e-03	5.73	3.19e-14	4773.10

Table 2. Table 2: CPU time in seconds and relative errors for each of the methods over the time interval [0, 45] with h = 4.5 ℎ 4.5 h=4.5

$m a t r i x$	$L_{0}$	GExpEuler		LrExpEuler		BrExpEuler		Erow3		Additive4
		Error	Time	Error	Time	Error	Time	Error	Time	Error	Time
rail371	$0_{371 \times 1}$	8.92e-14	0.75	7.19e-12	0.11	7.07e-11	1.10	7.19e-12	0.16	3.38e-11	0.33
rail371	rand(371,2)	8.61e-14	0.99	4.23e-12	0.13	7.25e-11	1.18	4.22e-12	0.19	3.26e-11	0.42
rail1357	$0_{1357 \times 1}$	2.23e-14	21.29	7.66e-12	0.56	2.48e-10	14.55	7.66e-12	0.80	3.60e-11	0.52
rail1357	rand(1357,2)	1.59e-14	29.54	4.18e-12	0.67	1.78e-10	16.11	4.27e-12	0.97	2.57e-11	0.56

Table 3. Table 3: The CPU and the relative errors with respect to the reference solution for N = 1000 . 𝑁 1000 N=1000.

Matrix	size(A)	LrExpEuler		BrExpEuler		Erow3		Additive4		Additive6
		Error	Time	Error	Time	Error	Time	Error	Time	Error	Time
fdm-sym	400 $\times$ 400	8.21e-07	16.64	1.46e-08	10.79	8.21e-07	41.31	3.68e-05	64.69	8.03e-07	86.99
	900 $\times$ 900	7.67e-05	110.77	3.06e-06	77.50	7.67e-05	191.85	5.39e-04	148.54	4.18e-05	189.66
	1600 $\times$ 1600	7.84e-04	500.91	6.21e-05	470.55	7.84e-04	711.84	2.53e-03	293.53	3.42e-04	348.13
	2500 $\times$ 2500	3.14e-03	1641.29	3.70e-04	1953.69	3.14e-03	2067.08	7.51e-03	541.69	1.15e-03	957.70
fdm-nonsym	400 $\times$ 400	1.18e-06	28.03	1.96e-08	19.37	1.18e-06	87.31	4.40e-05	116.79	1.03e-06	170.70
	900 $\times$ 900	7.80e-05	112.67	3.37e-06	116.65	7.80e-05	263.88	5.64e-04	392.62	4.51e-05	491.76
	1600 $\times$ 1600	8.16e-04	558.12	6.52e-05	702.36	8.16e-04	925.01	2.56e-03	917.33	3.58e-04	990.88
	2500 $\times$ 2500	3.22e-03	2013.00	3.92e-04	2822.60	3.22e-03	2629.20	8.00e-03	1671.76	1.22e-03	2120.35

Equations92

\left\{\begin{array}[]{l}X^{\prime}(t)=AX(t)+X(t)D+Q-X(t)GX(t)=:F(X(t)),\\ X(0)=X_{0},\end{array}\right.

\left\{\begin{array}[]{l}X^{\prime}(t)=AX(t)+X(t)D+Q-X(t)GX(t)=:F(X(t)),\\ X(0)=X_{0},\end{array}\right.

\left\{\begin{array}[]{l}X^{\prime}(t)=AX(t)+X(t)A^{T}+Q-X(t)GX(t),\\ X(0)=X_{0},\end{array}\right.

\left\{\begin{array}[]{l}X^{\prime}(t)=AX(t)+X(t)A^{T}+Q-X(t)GX(t),\\ X(0)=X_{0},\end{array}\right.

\left\{\begin{array}[]{l}y^{\prime}=Ly+f(y),~{}y\in R^{N},~{}L\in R^{N\times N},\\ y(0)=y_{0}.\end{array}\right.

\left\{\begin{array}[]{l}y^{\prime}=Ly+f(y),~{}y\in R^{N},~{}L\in R^{N\times N},\\ y(0)=y_{0}.\end{array}\right.

\displaystyle y(t_{n}+h_{n})=e^{h_{n}L}y(t_{n})+\int_{0}^{h_{n}}e^{(h_{n}-s)L}f\big{(}y(t_{n}+s)\big{)}ds.

\displaystyle y(t_{n}+h_{n})=e^{h_{n}L}y(t_{n})+\int_{0}^{h_{n}}e^{(h_{n}-s)L}f\big{(}y(t_{n}+s)\big{)}ds.

X (t) = e^{t A} X_{0} e^{t D} + \int_{0}^{t} e^{(t - τ) A} Q e^{(t - τ) D} d τ - \int_{0}^{t} e^{(t - τ) A} X (τ) GX (τ) e^{(t - τ) D} d τ .

X (t) = e^{t A} X_{0} e^{t D} + \int_{0}^{t} e^{(t - τ) A} Q e^{(t - τ) D} d τ - \int_{0}^{t} e^{(t - τ) A} X (τ) GX (τ) e^{(t - τ) D} d τ .

X (t_{n} + h_{n}) = e^{h_{n} A} X (t_{n}) e^{h_{n} D} + h_{n} \int_{0}^{1} e^{(1 - s) h_{n} A} Q e^{(1 - s) h_{n} D} d s - h_{n} \int_{0}^{1} e^{(1 - s) h_{n} A} X (t_{n} + s h_{n}) GX (t_{n} + s h_{n}) e^{(1 - s) h_{n} D} d s .

X (t_{n} + h_{n}) = e^{h_{n} A} X (t_{n}) e^{h_{n} D} + h_{n} \int_{0}^{1} e^{(1 - s) h_{n} A} Q e^{(1 - s) h_{n} D} d s - h_{n} \int_{0}^{1} e^{(1 - s) h_{n} A} X (t_{n} + s h_{n}) GX (t_{n} + s h_{n}) e^{(1 - s) h_{n} D} d s .

S (X) := A X + X D

S (X) := A X + X D

e^{S} (X) = e^{H} e^{P} (X) = e^{A} X e^{D},

e^{S} (X) = e^{H} e^{P} (X) = e^{A} X e^{D},

X (t_{n} + h_{n}) = e^{h_{n} S} (X (t_{n})) + h_{n} \int_{0}^{1} e^{(1 - s) h_{n} S} (Q) d s - h_{n} \int_{0}^{1} e^{(1 - s) h_{n} S} (X (t_{n} + s h_{n}) GX (t_{n} + s h_{n})) d s .

X (t_{n} + h_{n}) = e^{h_{n} S} (X (t_{n})) + h_{n} \int_{0}^{1} e^{(1 - s) h_{n} S} (Q) d s - h_{n} \int_{0}^{1} e^{(1 - s) h_{n} S} (X (t_{n} + s h_{n}) GX (t_{n} + s h_{n})) d s .

X^{'} (t) = A X (t) + X (t) D + Q (t), X (t_{n}) = X_{n}

X^{'} (t) = A X (t) + X (t) D + Q (t), X (t_{n}) = X_{n}

X (t_{n} + h_{n}) = e^{h_{n} S} (X (t_{n})) + j = 1 \sum m h_{n}^{j} φ_{j} (h_{n} S) (Q^{(j - 1)} (t_{n})) d s + R_{m} (h_{n}),

X (t_{n} + h_{n}) = e^{h_{n} S} (X (t_{n})) + j = 1 \sum m h_{n}^{j} φ_{j} (h_{n} S) (Q^{(j - 1)} (t_{n})) d s + R_{m} (h_{n}),

R_{m} (h_{n}) = \frac{h _{n}^{m + 1}}{( m - 1 )!} \int_{0}^{1} \int_{0}^{1} s^{m} (1 - θ)^{m - 1} e^{(1 - s) h_{n} S} (Q^{(m)} (t_{n} + θ s h_{n})) d θ d s,

R_{m} (h_{n}) = \frac{h _{n}^{m + 1}}{( m - 1 )!} \int_{0}^{1} \int_{0}^{1} s^{m} (1 - θ)^{m - 1} e^{(1 - s) h_{n} S} (Q^{(m)} (t_{n} + θ s h_{n})) d θ d s,

φ_{j} (z) = \int_{0}^{1} e^{(1 - θ) z} \frac{θ ^{j - 1}}{( j - 1 )!} d θ, j \geq 1.

φ_{j} (z) = \int_{0}^{1} e^{(1 - θ) z} \frac{θ ^{j - 1}}{( j - 1 )!} d θ, j \geq 1.

X (t_{n} + h_{n}) = e^{h_{n} S} X_{n} + h_{n} \int_{0}^{1} e^{(1 - s) h_{n} S} (Q (t_{n} + s h_{n})) d s .

X (t_{n} + h_{n}) = e^{h_{n} S} X_{n} + h_{n} \int_{0}^{1} e^{(1 - s) h_{n} S} (Q (t_{n} + s h_{n})) d s .

Q (t_{n} + s h_{n}) = j = 1 \sum m \frac{( s h _{n} ) ^{j - 1}}{( j - 1 )!} Q^{(j - 1)} (t_{n}) + \int_{0}^{1} \frac{( s h _{n} ) ^{m}}{( m - 1 )!} (1 - θ)^{m - 1} Q^{(m)} (t_{n} + θ s h_{n})) d θ

Q (t_{n} + s h_{n}) = j = 1 \sum m \frac{( s h _{n} ) ^{j - 1}}{( j - 1 )!} Q^{(j - 1)} (t_{n}) + \int_{0}^{1} \frac{( s h _{n} ) ^{m}}{( m - 1 )!} (1 - θ)^{m - 1} Q^{(m)} (t_{n} + θ s h_{n})) d θ

φ_{j + 1} (z) = \frac{φ _{j} ( z ) - \frac{1}{j !}}{z}, φ_{0} (z) = e^{z} .

φ_{j + 1} (z) = \frac{φ _{j} ( z ) - \frac{1}{j !}}{z}, φ_{0} (z) = e^{z} .

X (t) = e^{t S} N_{0} + j = 1 \sum m t^{j} φ_{j} (t S) N_{j} .

X (t) = e^{t S} N_{0} + j = 1 \sum m t^{j} φ_{j} (t S) N_{j} .

X^{'} (t) = S_{n} (X) + G_{n} (X),

X^{'} (t) = S_{n} (X) + G_{n} (X),

S_{n} (X) = A_{n} X + X D_{n}, G_{n} (X) = F (X) - S_{n} (X)

S_{n} (X) = A_{n} X + X D_{n}, G_{n} (X) = F (X) - S_{n} (X)

X (t_{n} + h_{n}) = e^{h_{n} S_{n}} (X (t_{n})) + h_{n} \int_{0}^{1} e^{(1 - s) h_{n} S_{n}} (G_{n} (X (t_{n} + s h_{n}))) d s .

X (t_{n} + h_{n}) = e^{h_{n} S_{n}} (X (t_{n})) + h_{n} \int_{0}^{1} e^{(1 - s) h_{n} S_{n}} (G_{n} (X (t_{n} + s h_{n}))) d s .

\left\{\begin{array}[]{l}X_{ni}=e^{c_{i}h_{n}\mathcal{S}_{n}}(X_{n})+h_{n}\sum\limits^{i-1}_{j=1}a_{ij}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{nj})),~{}~{}1\leq i\leq s,\\ X_{n+1}=e^{h_{n}\mathcal{S}_{n}}(X_{n})+h_{n}\sum\limits^{s}_{i=1}b_{i}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{ni})).\end{array}\right.

\left\{\begin{array}[]{l}X_{ni}=e^{c_{i}h_{n}\mathcal{S}_{n}}(X_{n})+h_{n}\sum\limits^{i-1}_{j=1}a_{ij}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{nj})),~{}~{}1\leq i\leq s,\\ X_{n+1}=e^{h_{n}\mathcal{S}_{n}}(X_{n})+h_{n}\sum\limits^{s}_{i=1}b_{i}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{ni})).\end{array}\right.

\displaystyle\begin{array}[]{ll}X_{n+1}&=e^{h_{n}\mathcal{S}_{n}}(X_{n})+h_{n}\varphi_{1}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{n}))\\ &=X_{n}+h_{n}\varphi_{1}(h_{n}\mathcal{S}_{n})(F(X_{n})).\end{array}

\displaystyle\begin{array}[]{ll}X_{n+1}&=e^{h_{n}\mathcal{S}_{n}}(X_{n})+h_{n}\varphi_{1}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{n}))\\ &=X_{n}+h_{n}\varphi_{1}(h_{n}\mathcal{S}_{n})(F(X_{n})).\end{array}

\begin{array}[]{lll}X_{n,2}&=e^{h_{n}\mathcal{S}_{n}}(X_{n})+h_{n}\varphi_{1}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{n}))\\ &=X_{n}+h_{n}\varphi_{1}(h_{n}\mathcal{S}_{n})(F(X_{n})),\\ X_{n+1}&=e^{h_{n}\mathcal{S}_{n}}(X_{n})+h_{n}(\varphi_{1}(h_{n}\mathcal{S}_{n})-2\varphi_{3}(h_{n}\mathcal{S}_{n}))(\mathcal{G}_{n}(X_{n}))+2h_{n}\varphi_{3}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{n,2}))\\ &=X_{n}+h_{n}\varphi_{1}(h_{n}\mathcal{S}_{n})(F(X_{n}))+2h_{n}\varphi_{3}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{n,2})-\mathcal{G}_{n}(X_{n})).\end{array}

\begin{array}[]{lll}X_{n,2}&=e^{h_{n}\mathcal{S}_{n}}(X_{n})+h_{n}\varphi_{1}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{n}))\\ &=X_{n}+h_{n}\varphi_{1}(h_{n}\mathcal{S}_{n})(F(X_{n})),\\ X_{n+1}&=e^{h_{n}\mathcal{S}_{n}}(X_{n})+h_{n}(\varphi_{1}(h_{n}\mathcal{S}_{n})-2\varphi_{3}(h_{n}\mathcal{S}_{n}))(\mathcal{G}_{n}(X_{n}))+2h_{n}\varphi_{3}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{n,2}))\\ &=X_{n}+h_{n}\varphi_{1}(h_{n}\mathcal{S}_{n})(F(X_{n}))+2h_{n}\varphi_{3}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{n,2})-\mathcal{G}_{n}(X_{n})).\end{array}

\displaystyle\mathbb{A}_{n}=\left(\begin{tabular}[]{cccccc}$A_{n}$&$\mathcal{G}_{n}$\\ $0$&$-D_{n}$\end{tabular}\right)\in\mathbb{C}^{(M+N)\times(M+N)}.

\displaystyle\mathbb{A}_{n}=\left(\begin{tabular}[]{cccccc}$A_{n}$&$\mathcal{G}_{n}$\\ $0$&$-D_{n}$\end{tabular}\right)\in\mathbb{C}^{(M+N)\times(M+N)}.

\displaystyle e^{\mathbb{A}_{n}}=\left(\begin{tabular}[]{cccccc}$e^{A_{n}}$&$\int_{0}^{1}e^{(1-s)A_{n}}\mathcal{G}_{n}e^{-sD_{n}}ds$\\ $0$&$e^{-D_{n}}$\end{tabular}\right)\in\mathbb{C}^{(M+N)\times(M+N)}.

\displaystyle e^{\mathbb{A}_{n}}=\left(\begin{tabular}[]{cccccc}$e^{A_{n}}$&$\int_{0}^{1}e^{(1-s)A_{n}}\mathcal{G}_{n}e^{-sD_{n}}ds$\\ $0$&$e^{-D_{n}}$\end{tabular}\right)\in\mathbb{C}^{(M+N)\times(M+N)}.

\displaystyle X_{n+1}=\left((I_{M},0)e^{h_{n}\mathbb{A}_{n}}\left(\begin{tabular}[]{cccccc}$X_{n}$\\ $I_{M}$\end{tabular}\right)\right)e^{h_{n}D_{n}}.

\displaystyle X_{n+1}=\left((I_{M},0)e^{h_{n}\mathbb{A}_{n}}\left(\begin{tabular}[]{cccccc}$X_{n}$\\ $I_{M}$\end{tabular}\right)\right)e^{h_{n}D_{n}}.

φ_{k} (h_{n} S_{n}) N_{k} \approx \frac{1}{( k - 1 )!} j = 0 \sum p ω_{j} s_{j}^{k - 1} e^{(1 - s_{j}) h_{n} S_{n}} (N_{k}) .

φ_{k} (h_{n} S_{n}) N_{k} \approx \frac{1}{( k - 1 )!} j = 0 \sum p ω_{j} s_{j}^{k - 1} e^{(1 - s_{j}) h_{n} S_{n}} (N_{k}) .

φ_{0} (\hat{S}_{n}) N_{0} + φ_{1} (\hat{S}_{n}) N_{1} + φ_{2} (\hat{S}_{n}) N_{2} + \dots + φ_{k} (\hat{S}_{n}) N_{k},

φ_{0} (\hat{S}_{n}) N_{0} + φ_{1} (\hat{S}_{n}) N_{1} + φ_{2} (\hat{S}_{n}) N_{2} + \dots + φ_{k} (\hat{S}_{n}) N_{k},

\left\{\begin{array}[]{l}W_{0}=N_{0},\\ W_{j}=\hat{\mathcal{S}}_{n}(W_{j-1})+N_{j},~{}~{}j=1,\cdots,k,\end{array}\right.

\left\{\begin{array}[]{l}W_{0}=N_{0},\\ W_{j}=\hat{\mathcal{S}}_{n}(W_{j-1})+N_{j},~{}~{}j=1,\cdots,k,\end{array}\right.

j = 0 \sum k φ_{j} (\hat{S}_{n}) N_{j} = φ_{k} (\hat{S}_{n}) (W_{k}) + j = 0 \sum k - 1 \frac{1}{j !} W_{j} .

j = 0 \sum k φ_{j} (\hat{S}_{n}) N_{j} = φ_{k} (\hat{S}_{n}) (W_{k}) + j = 0 \sum k - 1 \frac{1}{j !} W_{j} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsNumerical methods for differential equations · Model Reduction and Neural Networks · Matrix Theory and Algorithms

Full text

Exponential integrators for large-scale stiff matrix Riccati differential equations111The work was supported by the Natural Science Foundation of Jilin Province of China

(20180101224JC)

Dongping Li

[email protected]

a Department of Mathematics, Jilin University, Changchun 130012, PR China

bDepartment of Mathematics, Changchun Normal University, Changchun 130032, PR China

Abstract

Matrix Riccati differential equations arise in many different areas and are particular important within the field of control theory. In this paper we consider numerical integration for large-scale systems of stiff matrix Riccati differential equations. We show how to apply exponential Rosenbrock-type integrators to get approximate solutions. Two typical exponential integration schemes are considered. The implementation issues are addressed and some low-rank approximations are exploited based on high quality numerical algebra codes. Numerical comparisons demonstrate that the exponential integrators can obtain high accuracy and efficiency for solving large-scale systems of stiff matrix Riccati differential equations.

keywords:

Matrix Riccati differential equations, Exponential integrators, $\varphi$ -functions, Low-rank approximation

MSC:

[2010] 65L05 , 65F10, 65F30

\newcaptionstyle

left \usecaptionmargin

\captionlabelfont\captionlabel

\onelinecaption\captiontext\captiontext

1 Introduction

In this paper we are concerned with numerical methods for large-scale systems of stiff matrix Riccati differential equations (MRDEs) of the following form

[TABLE]

where $A\in\mathbb{R}^{M\times M},~{}D\in\mathbb{R}^{N\times N},~{}Q\in\mathbb{R}^{M\times N},~{}G\in\mathbb{R}^{N\times M}$ are given matrices and $X(t)\in\mathbb{R}^{M\times N}$ is the unknown matrix-valued function. MRDEs of this form occur in many important applications such as optimal control, $H_{\infty}$ -control, filtering, boundary value problems for systems of ODEs and many others (see e.g. [1, 4, 20, 22]). In most control problems, the coefficient matrices $A$ and $D$ are obtained from the discretization of operators defined on infinite dimensional spaces, and the fast and slow modes exist. This means that the associated MRDEs will be fairly large and stiff.

An important special case of (1) is the symmetric MRDEs

[TABLE]

here, $Q=Q^{T},G=G^{T}$ and $X_{0}=X_{0}^{T}.$ It is obvious that the solution of symmetric MRDEs is symmetric as $X(t)^{T}$ is also a solution. The symmetric MRDEs is possibly the most widely studied equations due to its importance in linear-quadratic optimal control problems. Another special mention should be paid to MRDEs (1) with $G=0,$ yieiding the so-called matrix Sylvester differential equations (MSDEs). For a thorough description of these equations and some qualitative issues, we refer the reader to [1, 12, 22, 23, 32] and the references appearing therein.

Many numerical methods have been developed in the past for solving MRDEs. Perhaps the most natural numerical technique is to rewrite (1) as an $MN$ -vector ODEs based on Kronecker product, and then to use a standard numerical integrators such as Runge-Kutta or linear multi-step solvers [9]. However, these approaches are not suitable for solution of large stiff MRDEs. They are generally computationally expensive and hard to exploit the structure inherited in some large practical problems. For stiff MRDEs, some matrix-valued versions of implicit time integration schemes, such as the BDF, Rosenbrock methods have been explored through a direct time discretization of (1), see e.g. [7, 10, 11]. Recently, some other unconventional numerical methods have been also developed for MRDEs and related problems, including splitting methods [29, 35, 36] and projection methods [13, 14, 25], etc.

The aim of this paper is to introduce exponential integrators for large-scale stiff problems of the forms 1 and 2. In the past two decades exponential integrators have become a popular tool for solving large-scale stiff semi-linear systems of ODEs

[TABLE]

A general derivation of exponential integrators is based on the variation-of-constants formula

[TABLE]

By approximating the nonlinear terms $f\big{(}y(t_{n}+s)\big{)}$ by an appropriate algebra polynomial, various type of exponential integrators can be exploited. Different approximations result in different types of exponential integrators of either multi-step type or Runge-Kutta type, see e.g. [16, 17, 19, 24]. A main advantage of exponential integrators with stiff order conditions don’t suffer from an order reduction even if the matrix $L$ is a discretization of a unbounded linear operator. For a full overview of exponential integrators and associated software, we refer the readers to the reviews [18, 30] and references therein. Although the matrix differential equations (1) can be reformulated as the form (3) and solved by an exponential integrator, this approach will generate very large $L$ and not be appropriate.

In the present paper we propose matrix-valued versions exponential integrators for stiff MRDEs (1). The methods provides an efficient alternative to implicit integrators for computing solutions of MRDEs. For large-scale systems, in many applications it is often observed, both practical and theoretical, the solution has low numerical rank and can be approximated by products of low-rank matrices [27, 37]. To utilize such structure, we introduce how the low-rank implementation can be applied to exponential integrators. Thus we are able to save computational and memory storage requirements compared to the simple application of exponential integrators.

The remainder of the paper is organized as follows. In Section 2, we give some basic results and properties of MRDEs. In Section 3, the exponential Rosenbroc-type methods are introduced for the application to the MRDEs. Section 4 we show some issues of implementation and the low-rank approximations for both typical exponential integration schemes are exploited. Section 5 is devoted to some numerical examples and comparisons with splitting methods of similar orders. Finally, we draw some conclusions in Section 6.

2 Preliminaries

We start with recalling a general result on the solution of the MRDEs. The following result shows that the MRDEs (1) can be equivalently written in an integral form (see e.g. [26]).

Theorem 1.

The exact solution of the MRDEs (1) is given by

[TABLE]

Proof.

The proof can be done directly by differentiating both sides. ∎

Specifically, over the time interval $[t_{n},~{}t_{n}+h_{n}],$ by using a change of variables $t=t_{n}+sh_{n}$ in (5) to give

[TABLE]

The formula (6) also holds for time-varying coefficient matrices $Q=Q(t)$ and $G=G(t).$ To make constructing methods for MRDEs easier, we use

[TABLE]

to denote the linear operator from the right-hand side of MRDEs (1), which is called Sylvester operator. The operator exponential satisfy the following relation (see [5]):

[TABLE]

here $\mathcal{H}(X)=AX$ and $\mathcal{P}(X)=XD.$ Then, expression (6) has the simplified form

[TABLE]

The first term from the right-hand side of (9) involves operator exponential and represents the homogenous part of the solution, whereas the other two terms consist of integrals, again involving operator exponential. A natural idea to construct exponential integrators is to approximate the integrals on the right-hand side of (9) by a quadrature formula, in which only the nonlinearity term $XGX$ are approximated but the operator exponential are treated exactly. In particular, for MSDEs, we have the following result.

Lemma 1.

Let $A\in\mathbb{R}^{M\times M},~{}D\in\mathbb{R}^{N\times N},$ and $Q(t)\in\mathbb{R}^{M\times N}$ be a sufficiently differential matrix-value function, then the exact solution of the matrix differential equations

[TABLE]

can be represented by the expansion

[TABLE]

where

[TABLE]

Proof.

By formula (8), the solution of equations (10) can be written

[TABLE]

Inserting the Taylor series expansion of $Q(t_{n}+sh_{n})$

[TABLE]

into the formula (14) and applying the definition (13) we arrive at the required result. ∎

The functions $\varphi_{j}(z)$ defined in (13) satisfy the following recurrence relations

[TABLE]

A special case of the nonhomogeneous term in equations (10) is an matrix polynomials, i.e., $Q(t)=\sum\limits_{j=0}^{m-1}\frac{t^{j-1}}{(j-1)!}N_{j},$ $N_{j}\in\mathbb{R}^{M\times N},~{}~{}j=0,1,\cdots,m-1,$ in this case, the exact solution of equations (10) can be represented by the expansion

[TABLE]

3 Exponential Rosenbrock-type integrators for MRDEs

In this section, we consider the time discretization of MRDEs (1). Rewrite equations (1) as

[TABLE]

where $\mathcal{S}_{n}$ denotes the Fréchet derivative of $F$ and $\mathcal{G}_{n}$ the nonlinear remainder at $X_{n},$ respectively:

[TABLE]

with $A_{n}=(A-X_{n}G)$ and $D_{n}=(D-GX_{n}).$

It is obvious that $\mathcal{S}_{n}$ is a Sylvester operator. Formally, by the variation of constants formula (9), the exact solution of (18) can be written as follows:

[TABLE]

The above expression has a similar structure with (4) but Sylvester operator exponential instead of matrix exponential. Thus we can apply various existing exponential integrators to (18). The application of the general exponential Runge-Kutta type methods [16], to the MRDEs (18) yields

[TABLE]

Here, $c_{i}$ is the nodes, and the coefficients $a_{ij}(z),$ $b_{i}(z)$ are linear combinations of the $\varphi_{j}(c_{i}z),$ $\varphi_{j}(z)$ respectively. These coefficients can be determined by a stiff error analysis which can be adapted from the stiff order theory presented in [18, 28]. The process is highly sophisticated and is omitted here. In our context, we only consider two specific exponential integration schemes. They will be used in our numerical experiments in Section 5. The first and simplest exponential integration scheme is the exponential Rosenbrock-type Euler scheme

[TABLE]

The scheme is computationally attractive since it is second order with only one $\varphi$ -function. The second scheme is order three (denoted Erow3), which can be regarded as a modification of exponential Rosenbrock-type Euler scheme

[TABLE]

The internal stage has the same structure as the exponential Rosenbrock-type Euler scheme (24), and the external stage is a perturbation of the internal stage. The above two schemes are usually embedded to create an adaptive time stepping method.

4 Implementation issues

For exponential integrators, the main computational cost is to approximate the exponential and exponential-type functions at each time-step. To our knowledge, there is no explicit method to evaluate the functions of a Sylvester operator in the literatures. For the computation of the first $\varphi$ -function, the following formula gives an indirect way.

Define the augmented matrix $\mathbb{A}_{n}$ by

[TABLE]

Using the formula (10.40) arising in ([15], we have

[TABLE]

Then the scheme (24) can be rewritten as

[TABLE]

For the computation of a single matrix exponential or its action on a thin matrix, a number of methods have been proposed in the literatures for carrying out this task, see e.g. [2, 3, 34] and the review [21]. This approach has the major advantage of simplicity but is likely to be too expensive for large $M$ and $N.$

A more general strategy for approximating $\varphi$ -functions is to apply a numerical quadrature scheme. For a given function $\varphi_{k}$ , the Sylvester operator $\mathcal{S}_{n}$ and an matrix $N_{k},$ we approximate $\varphi_{k}(h_{n}\mathcal{S}_{n})(N_{k})$ by a quadrature approximation with quadrature nodes $s_{j}$ and weights $\omega_{j}:$

[TABLE]

Thus to evaluate $\varphi_{k}(h_{n}\mathcal{S}_{n})N_{k}$ we need to compute $p+1$ operator exponential acting on the same matrix. In practical application if the matrix $N_{k}$ has a low-rank factorization $N_{k}=L_{k}D_{k}U_{k}^{T}$ where both $L_{k}$ and $U_{k}$ are full column rank and $D_{k}$ is nonsingular, the block Krylov subspace method can be applied to the computation of operator exponential involved.

In fact, as shown in the schemes (24) and (25), every stage in an exponential integrator can be expressed as a linear combination of the form

[TABLE]

here $\hat{\mathcal{S}}_{n}=h_{n}\mathcal{S}_{n},$ $N_{i}\in R^{M\times N},i=0,1,\cdots,k.$ Using the recurrence relation (16) we can calculate (30) recursively. Two alternatives are available. The first one is a forward recursion, i.e.,

[TABLE]

then

[TABLE]

The main computational cost of this process includes the action of $k$ Sylvester operator and a $\varphi_{k}(\hat{\mathcal{S}}_{n})(W_{k})$ . Another approach is a backward recursion

[TABLE]

This process requires the computation of $k$ algebra Sylvester equations and a Sylvester operator exponential acting on matrix.

In many practical applications the MRDEs have an low-rank structure and the solution has the low-rank property. In such cases it is necessary to avoid forming the matrices $X_{n}$ explicitly, because this in general leads to dense computations. In the remainder of the section we briefly introduce how to implement the above mentioned two schemes in a low-rank fashion to the symmetric MRDEs. For simplicity let us consider the symmetric MRDEs (2).

Provided $Q,$ $G$ and $X_{0}$ in (2) are symmetric positive semi-definite, and are given in the low-rank form

[TABLE]

with $C\in\mathbb{R}^{l\times N},$ $B\in\mathbb{R}^{N\times q},$ $L_{0}\in\mathbb{R}^{N\times r},$ and $D_{0}\in\mathbb{R}^{r\times r},$ $l,q,r\ll N.$ This implies that the solution $X(t)$ to the MRDEs (2) is also symmetric positive semi-definite for all $t>0.$ First, we consider the exponential Rosenbrock-type Euler scheme (24). Assume that the previous solution approximations $X_{n}$ admit a decomposition of the form $X_{n}=L_{n}D_{n}L_{n}^{T}$ with $L_{n}\in R^{N\times r_{n}},D_{n}\in R^{r_{n}\times r_{n}}.$ $F_{n}$ in scheme (24) can be written the form of $LDL^{T}:$

[TABLE]

The new matrix $\tilde{L}_{n}$ has more columns than $L_{n},$ and also more than their rank. As the number of columns in the decomposition increases, the computation cost will become prohibitively expensive. To overcome this difficulty we can incorporate the column compression strategy [27] to $\tilde{L}_{n},~{}\tilde{D}_{n},$ and find more suitable low-rank factors. Then, we apply the numerical quadrature formula (29) to approximate $h_{n}\varphi_{1}(h_{n}\mathcal{S}_{n})F_{n},$ and the decomposition $Y_{n}T_{n}Y_{n}^{T}$ is given by the factors

[TABLE]

and

[TABLE]

Note that the evaluation of $Y_{n}$ requires computation of $p+1$ products between matrix exponential and a thin matrix. For large matrix $A_{n},$ the block Krylov projection algorithm is a popular choice [33]. An advantage of this computation is that one can project the $p+1$ operator exponential in the same search subspace $\mathcal{K}_{m}(A_{n},\tilde{L}_{n})$ and evaluate them simultaneously.

Now, using the splitting of $h_{n}\varphi_{1}(h_{n}\mathcal{S}_{n})F_{n}$ and of the solution $X_{n}=L_{n}D_{n}L_{n}^{T},$ the approximation $L_{n+1}D_{n+1}L_{n+1}^{T}$ to $X_{n+1}$ is given by

[TABLE]

Again, we can employ column compression strategy to obtain low-rank splitting factors.

We now describe an alternative low-rank implementation of the exponential Rosenbrock-type Euler scheme. Apply the backward recursion (33)-(34) to the scheme (24), giving

[TABLE]

This results in solving one algebraic Lyapunov equation (ALE) which the right hand side has low-rank form in each time-step. There are many methods for solving Lyapunov equations where the right-hand side is of low rank, for instance by a low-rank ADI iteration [6] or Krylov subspace based methods [38]. Due to the availability of low-rank ADI iteration based codes, here we limit ourselves to this procedure.

For order-third exponential integration scheme (25), again, the previous solution approximation $X_{n}=L_{n}D_{n}L_{n}^{T}$ is assumed to be given in low-rank format. Note that the interval stage $X_{n,2}$ is the same as exponential Rosenbrock-type Euler scheme, thus $X_{n,2}$ can be written in the low-rank form $L_{n,2}D_{n,2}L_{n,2}^{T}.$ The external stage $X_{n+1}$ is a perturbation of the matrix $X_{n,2}$ by $h_{n}\varphi_{3}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{n,2})-\mathcal{G}_{n}(X_{n})).$ In order to find a low-rank factorization of the entire right hand side, we first consider the $LDL^{T}$ -type splitting of $\mathcal{G}_{n}(X_{n,2})-\mathcal{G}_{n}(X_{n}).$ Direct calculation shows that

[TABLE]

Inserting the splitting factors of $X_{n}$ and $X_{n,2}$ into (53) finally gives the splitting $\bar{L}_{n}\bar{D}_{n}\bar{L}_{n}^{T}$ with

[TABLE]

Again, using (29), the splitting factors $\bar{Y}_{n},\bar{T}_{n}$ of $h_{n}\varphi_{3}(h_{n}\mathcal{S}_{n})(\mathcal{G}_{n}(X_{n,2})-\mathcal{G}_{n}(X_{n}))$ can be computed as follows:

[TABLE]

and

[TABLE]

Now, using the $LDL^{T}$ -type splitting with $X_{n,2}=L_{n,2}D_{n,2}L_{n,2}^{T}$ we obtain

[TABLE]

with

[TABLE]

and

[TABLE]

In actual implementation, once the new splitting factors are formed, column compression strategy should be performed to eliminate the redundant information.

5 Numerical experiments

In this section, we present some numerical experiments to illustrate the behaviour of exponential integration methods. We compare the numerical performance of the exponential Rosenbrock-type Euler scheme (24) (denoted ExpEuler) and the third order exponential integration scheme Erow3 (25) with the splitting schemes in [36]. For ExpEuler, we consider all the above mentioned three different implementations. They are marked as follows: the general implementation (28) (denoted by GExpEuler), the low-rank implementation (denoted by LrExpEuler) and the backward recursion implementation (52) (denoted by BrExpEuler). For low-rank implementations, the tolerance for column compression strategies are set to $n\cdot\epsilon,$ where $n$ is the system dimension and $\epsilon$ denotes the machine precision. All experiments are performed under Windows 10 and MATLAB R2018b running on a laptop with an Intel Core i7 processor with 1.8 GHz and RAM 8 GB. Unless otherwise stated, we use the relative errors at the final time, measured in the Frobenius norm.

Experiment 1. As the first test, we consider the matrix $A\in\mathbb{R}^{n\times n}$ obtained from the standard 5-point difference discretization of the two-dimensional PDE

[TABLE]

on the domain $\Omega=[0,1]^{2}$ with homogeneous Dirichlet boundary conditions, and $B,C\in\mathbb{R}^{n\times 2}$ with entries chosen randomly from [0; 1]. The matrix $A$ is a negative stiffness matrix which can be generated by MATLAB function fdm2dmatrix from LYAPACK toolbox [31]. We consider the discretization of (65) for two different values of functions $f_{1},f_{2},f_{3},$ namely $f_{1}=f_{2}=f_{3}=0$ and $f_{1}=10x,$ $f_{2}=100y,$ $f_{3}=0,$ respectively. The former generate a symmetric matrix (denoted fdm-sym), while the latter is unsymmetric (denoted fdm-nonsym). The corresponding initial values $X_{0}$ are choosen as a low-rank product $X_{0}=L_{0}L_{0}^{\prime},$ where $L_{0}\in\mathbb{R}^{n\times 2}$ are randomly generated. To ensure the availability of a reference solution, we perform two sizes on the time interval [0, 1], one with $n=64$ and the other with $n=100.$ The reference solutions are obtained by MATLAB built-in function ode15s with an absolute tolerance of $10^{-20}$ and a relative tolerance of $2.22045\cdot 10^{-14}.$ This has been done by vectorizing the MRDEs into a vector-valued ODEs with $n^{2}$ unknowns. We use the above mentioned methods to integrate the four systems over the time interval [0, 1] with time step size $h=0.01.$

Table 1 lists the relative errors at the final time $t=1$ as well as the total time (in seconds) of the methods to integrate these systems. The results show that the two exponential integration schemes achieve the high precision of about $10^{-14}$ in all cases and obtained a higher order of convergence than we expected. An interpretation of this as the exponential integrators could be suitable for the structure of the MRDEs and capture some qualitative properties. As a comparison, we also present the results for the additive symmetric scheme of order 4 (denoted Additive4) in [36] with the same timestep and ode15s with the same accuracy. The code for the Additive4 contains parallel loops which uses 4 workers on our machine. We can see the exponential integration schemes accomplish higher computational accuracy than Additive4 and take less runtimes than ode15s.

Figures 1, 2 show plots of the F-norm of solutions obtained by ExpEuler and the reference solutions provided by ode15s for each test system. The Erow3 yields very similar behaviors with ExpEuler and is omitted here. From these figures we see that the ExpEuler follow well behaviours of the reference solutions. Although it is not very accurate in the start some time steps, the behaviours are completely corrected as the time increase. At the final time, the relative error even level out around $10^{-14}.$ This is also true in the subsequent experiments and we interpret this as exponential integration schemes being favorable for MRDEs.

Experiment 2. As a second test, we consider a finite element discretization of a heat equation arising from the optimal control of steel cooling [8]. The matrix $A\in\mathbb{R}^{n\times n}$ is symmetric and stable, $B\in\mathbb{R}^{n\times 7}$ and $C\in\mathbb{R}^{6\times n}.$ The initial value $X_{0}=L_{0}L_{0}^{\prime}$ is set as $L_{0}=0_{n\times 1}$ or random vectors $L_{0}=\textsf{rand}(n,2),$ with elements from the normal (0,1) distribution. To ensure the availability of a reference solution, we perform two sizes on the time interval [0,45] with $h=4.5$ , one with $n=371$ and the other with $n=1357.$ The reference solutions are computed by MATLAB built-in function ode45 with an absolute tolerance of $10^{-20}$ and a relative tolerance of $2.22045\cdot 10^{-14}.$ Similar the above experiment Figures 3 and 4 show the convergence behaviour of the solutions of the ExpEuler to the reference solution for the size $n=371$ and $n=1375.$ From these Figures, we observe that the ExpEuler gives a fairly good approximations of the reference solutions. We also present the numerical results of the different methods mentioned above in Table 2. We see that the proposed method performs quite well in terms of accomplished accuracy and computational time.

Experiment 3. The third experiment is the same problem as in Experiment 1 for larger-scale dimensions. We choose the same setting as in Experiment 1. We compare the low-rank approximations, LrExpEuler, Erow3, BrExpEuler to the symmetric splitting of orders 4 and 6 for systems of dimensions $n=400,900,1600,2500.$ For larger-scale stiff matrix, our tests indicate that the numerical integration formula (29) will cause low computational accuracy. In following tests, we use a relatively smaller sizes $h=0.001$ except for the BrExpEuler which is set as $h=0.01.$ Due to the systems sizes, it is infeasible to use ode15s to compute an accurate reference solutions. Instead, we use the eighth-order symmetric splitting scheme in [36] with the time step $h=0.001.$ Table 3 presents the relative errors and the corresponding computation times for each of the systems with different values of the system size $n.$ It is noted that the BrExpEuler produces the smallest errors of all methods even though the time step is larger. This is due to the poor performance of numerical integrator to $\varphi$ -function. This also shows that the exponential integrators have large computational potential and inspires us to exploit other efficiently numerical approximations to $\varphi$ -function. We plan to investigate this option in our future work. In terms of the CPU times, we observe that in some cases the exponential integration methods spend more CPU times than the symmetric splitting schemes. This might have a weakened effect when the computation is done on a single core machine as the symmetric splitting methods employ parallel loops with four cores in our laptop. Again, the LrExpEuler obtained almost the same accuracy as Erow3, which illustrates the feasibility of adaptivity. We observe that the accuracy of the resulting solutions reduce as the size of system and its stiffness increase for all methods. We attribute this mainly to the pessimistic reference solutions.

6 Conclusion

In this paper, we show how to apply exponential integrators to get approximate solutions of large stiff MRDEs. The low-rank implementation of such schemes for large-scale applications and their comparison with current state-of-the-art integrators must be addressed. Numerical experiments illustrate that the exponential integration methods can achieve convergence than the expected order. Thus the exponential integrators can provide an efficient alternative to standard integrators for large-scale stiff problems. The study of the performance and application of the higher-order exponential integration schemes and their comparative performance with implicit schemes will be presented elsewhere. We also plan to develop new more efficient algorithms to approximate the exponential and related functions of a Sylvester operator acting on an matrix.

Bibliography38

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] H. Abou-Kandil, G. Freiling, V. Ionescu and G. Jank, Matrix Riccati Equations in Control and Systems Theory, Birkhäuser, Basel, Switzerland, 2003.
2[2] A. Al-Mohy and N. Higham, A new scaling and modified squaring algorithm for matrix functions, SIAM J. Matrix Anal. Appl., 31 (2009), pp. 970-989.
3[3] A. Al-Mohy and N. Higham, Computing the action of the matrix exponential, with an application to exponential integrators, SIAM J. Sci. Comput., 33 (2011), pp. 488-511.
4[4] U.M. Ascher, R.M. Mattheij and R. G. Russell, Numerical solution of boundary value problems for ordinary differential equations, Prentice-Hall, Englewood Cliffs, NJ, 1988.
5[5] M. Behr, P. Benner and J. Heiland, Solution Formulas for Differential Sylvester and Lyapunov Equations, ar Xiv, 2018, https://arxiv.org/abs/1811.08327 .
6[6] P. Benner, M. Köhler, and J. Saak, M.E.S.S.-the matrix equations sparse solvers library, https://www.mpi-magdeburg.mpg.de/projects/mess .
7[7] P. Benner and H. Mena, Rosenbrock methods for solving Riccati differential equations, IEEE T. Automat. Contr., 58 (2013), pp. 2950-2956.
8[8] P. Benner and J. Saak. A semi-discretized heat transfer model for optimal cooling of steel profiles, In Dimension Reduction of Large-Scale Systems, volume 45 of Lect. Notes Comput. Sci. Eng.), P. Benner, V. Mehrmann, and D. Sorensen, Eds.Berlin/Heidelberg, Germany: Springer-Verlag, 2005, pp. 353-356.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Taxonomy

Exponential integrators for large-scale stiff matrix Riccati differential equations111The work was supported by the Natural Science Foundation of Jilin Province of China

Abstract

keywords:

MSC:

1 Introduction

2 Preliminaries

Theorem 1**.**

Proof.

Lemma 1**.**

Proof.

3 Exponential Rosenbrock-type integrators for MRDEs

4 Implementation issues

5 Numerical experiments

6 Conclusion

Theorem 1.

Lemma 1.