3
$\begingroup$

What methods are known for efficiently solving a large set of linear homogeneous ordinary differential equations of the following form? \begin{equation} \mathbf{B} \frac{d\mathbf{y}}{dt} = \mathbf{A} \mathbf{y}, \end{equation} 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).

Improve this question
$\endgroup$
9
  • 1
    $\begingroup$ How large is $B$? $\endgroup$
    – Wolfgang Bangerth
    Jun 5, 2017 at 19:34
  • $\begingroup$ If a more precise size matters, I am thinking of $N$ being $10^4$ to $10^6$. $\endgroup$
    – norio
    Jun 5, 2017 at 19:48
  • 1
    $\begingroup$ 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? $\endgroup$
    – Bill Greene
    Jun 5, 2017 at 23:55
  • $\begingroup$ $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. $\endgroup$
    – norio
    Jun 6, 2017 at 2:08
  • 2
    $\begingroup$ 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. $\endgroup$
    – Richard Zhang
    Jun 6, 2017 at 2:34

1 Answer 1

Reset to default
5
$\begingroup$

B is known as a mass matrix. There are many methods which can directly solve with mass matrices. The main ones are implicit methods: implicit Runge-Kutta methods and the multistep BDF methods. Also Rosenbrock methods can handle mass matrices directly, but are not implicit methods.

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.

If you use Julia, DifferentialEquations.jl has access to these methods. Just change mass_matrix in the construction of the ODEProblem type. The Rosenbrock23() method will handle mass matrices like this. In addition, mass matrices can be passed to the implicit RK methods radau() and radau5() which are higher order and tend to be faster if you are looking to solve stiff problems with high accuracy. Technically Sundials' CVODE_BDF() could solve with a mass matrix, but there's an open issue for setting up the C FFI for doing so (it's high on my priority list. Full disclosure: I am a developer of DifferentialEquations.jl).

If you can handle the performance drop of going to MATLAB instead, you'll find the field is very similar. ode23s is a 2nd-order Rosenbrock method, is very similar to Rosenbrock23(), and can use constant mass matrices. ode15s is very similar to CVODE_BDF() and can take mass matrices. MATLAB does not have an implicit RK method, though it does have ode23t which is good for problems where a symplectic integrator is needed and it can handle mass matrices.

I do not see options for passing mass matrices to SciPy's ODE solvers, and have not used the SciPy wrappers. Those wrap some of the same codes as DifferentialEquations.jl so I know that in some cases the mass matrix can be solved by the underlying C/Fortran codes that they use, but it looks like that part of the API is not exposed.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you very much. Your explanation of the principle makes sense. I will look into the methods you suggested. I program now in Fortran but got interested in Julia. I will try it out for a later project. $\endgroup$
    – norio
    Jun 6, 2017 at 2:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.