Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It only takes a minute to sign up.
Sign up to join this community
When solving a system of ODEs, as part of a preconditioner, I get the system
$(A + D(t))x = b(t)$
where $A$ is a sparse matrix and $D(t)$ is diagonal. I'm currently solving this by taking the LU-decomposition of the system for each $t$. $A$ is constant throughout the solution.
Is there an algorithm that would allow me to solve this efficiently without recomputing a decomposition at every step?
I do not believe that a direct solve is possible in your situation, there are simply too many applications with too many users which would benefit too much if it was possible. In particular, the problem of computing, say, the transfer function $T(s) = C(A - sI)^{-1}B+D$ of a linear time invariant system $$ \begin{align} x'(t) &= Ax(t) + Bu(t) \\ y (t) &= Cx(t) + Du(t) \end{align}$$ or solving the Lyapunov matrix equation $$AX + XA^T = BB^T$$ using the alternating direction implicit (ADI) method would be greatly facilitated by any algorithm which would allows to pass cheaply from a factorization of $A$ to a factorization of $A - sI$.
On the off chance that $D(t)$ has low rank, then the Sherman-Morrison-Woodbury formula would allow you to quickly solve systems with coefficient matrix $A + D(t)$ provided you can quickly solve system with coefficient matrix $A$. However, I gather that this is not likely.
I suggest that you view the solution $x$ of your linear system $$(A + D(t)) x = b(t)$$ as a function of $t$. It is likely that $x$ can be differentiate a few times with respect to $t$ which suggests that future values of $x$ can be extrapolated from past values. In particular, if $h$ is small enough, then $$ x_0 = x(t), \quad \text{or}\quad x_0 = 2x(t) - x(t-h)$$ represent simple and easy to implement initial guesses for the value of $x(t+h)$. More sophisticated approximation are possible, but they are not necessarily worth the effort.
I suggest further that your solve your linear system using a preconditioned Krylov subspace method. Instead of computing a new preconditioner for every timestep, I suggest that you recycle the preconditioner a few times and compute a new preconditioner only when the current one appears ineffective, i.e., the iteration count becomes too large.
I cannot promise that these options will give you the speed you are looking for, but at least they provide some parameters which can be tuned: the choice of the initial guess, the preconditioner and the number of time-steps before you compute a new preconditioner.
No (unless $D$ has few nonzeros). Variants of this question get asked a lot here, and the answer is invariably the same.
