# Set of linear ordinary differential equations with a mass matrix

What methods are known for efficiently solving a large set of linear homogeneous ordinary differential equations of the following form? $$\mathbf{B} \frac{d\mathbf{y}}{dt} = \mathbf{A} \mathbf{y},$$ where $\mathbf{B} \in \mathbb{C}^{N\times N}$ is Hermitian positive definite, $\mathbf{A} \in \mathbb{C}^{N\times N}$ is Hermitian, $\mathbf{y}\in \mathbb{C}^N$, and $N$ is a large integer such that not all the elements of the matrix $\mathbf{A}$ or $\mathbf{B}$ can be stored on the main memory.

Is there a method that avoids computing $\mathbf{B}^{-1}$, solving a linear equation $\mathbf{B} \mathbf{x} = \mathbf{c}$ for any $\mathbf{c}$, and diagonalizing $\mathbf{B}$? There seems to be such an 'inverse free' method [1] for the related eigenvalue problem, $\mathbf{B} \lambda \mathbf{y} = \mathbf{A}\mathbf{y}$, and I would like to know if there is an extension of such a method to ODEs.

[1] e.g. Gene H. Golub, Qiang Ye, An Inverse Free Preconditioned Krylov Subspace Method for Symmetric Generalized Eigenvalue Problems, SIAM J. Sci. Comp. 24, 312 (2002).

• How large is $B$? – Wolfgang Bangerth Jun 5 '17 at 19:34
• If a more precise size matters, I am thinking of $N$ being $10^4$ to $10^6$. – norio Jun 5 '17 at 19:48
• And are $A$ and $B$ dense or sparse? You could probably solve this ODE system with the simple implicit Euler integrator and an iterative linear equation solver but it is unlikely that would be very efficient. Have you already considered that option? – Bill Greene Jun 5 '17 at 23:55
• $A$ and $B$ are both dense. I am now solving the linear equation $B dy/dt = A y$ for $dy/dt$ by an iterative method and propagating $y$ with the Runge-Kutta method. Probably this is similar to what you suggested by implicit Euler and iterative linear solver. I wanted to know if there is a better way. – norio Jun 6 '17 at 2:08
• The same inverse-free eigenvalue solver based on Krylov subspaces can also be used to solve ODEs. These are known as Krylov subspace exponential integrators. – Richard Zhang Jun 6 '17 at 2:34

A quick overview for "why it works" is because at every step these methods have to solve linear equations involving the Jacobian of the RHS. The linear problem is normally of the form $(I - \gamma J)x=b$ where $J$ is the system Jacobian, $\gamma$ and $b$ are some constants, and $x$ is the evolution of the system that must be solved for. Well, it turns out that if you just use the mass matrix $B$ by changing the system to $(B - \gamma J)x=b$, these directly solve the linear problem with the mass matrix, implicitly inverting $B$ in calculations that it would already have been doing (this is because the $I$ comes from the LHS when taking the Jacobian, but with a mass matrix instead that value is just $B$!). Doing this implementation with the non-adaptive 2nd order version of the Rosenbrock method is outlined in Shampine's paper and is not difficult, though you lose a lot of efficiency without adaptive timestepping here.