I am applying the Crank-Nicolson-Method to solve the diffusion equation in 1D.
I need to solve the implicit equation $$Au^{j+1} = Bu^{j} \equiv b$$
Note that $$A= \begin{bmatrix} \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \mathstrut^{.^{.^{.}}} \\ \dotsc & -\alpha & 2(1+\alpha) & -\alpha & 0 & 0 & 0 & 0 & \dotsc \\ \dotsc & 0 & -\alpha & 2(1+\alpha) & -\alpha & 0 & 0 & 0 & \dotsc \\ \dotsc & 0 & 0 & -\alpha & 2(1+\alpha) & -\alpha & 0 & 0 & \dotsc \\ \dotsc & 0 & 0 & 0 & -\alpha & 2(1+\alpha) & -\alpha & 0 & \dotsc \\ \dotsc & 0 & 0 & 0 & 0 & -\alpha & 2(1+\alpha) & -\alpha & \dotsc \\ \mathstrut^{.^{.^{.}}} & \vdots &\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \\ \end{bmatrix}$$
The strategy is as follows: compute the LU factorization of $A,$ then store it, using the function dgttrf from LAPACK. Then, I can solve the implicit equation with the function dgttrs.
However, I have different species of particles, and so the alpha is different for all of them. I would therefore need to calculate the LU factorization for every alpha. Note though that $$A=2\mathbf{I} + \alpha D.$$
This leads me to wonder: is it possible to just compute the LU factorization for D, store this, and then use this to determine the LU factorization for $A=2\mathbf{I} + \alpha D?$ In other words, are there any nice ways to convert: $$2\mathbf{I} + \alpha PLU \to P'L'U'$$ Where $P', L',$ and $U'$ are new factorizations?
The answer is almost certainly no.
The key fact to notice is that your matrix $A = I + \alpha D$ is a rank $n$ perturbation of the matrix $\alpha D$. You cannot profit from applying the Sherman-Morrison-Woodbury formula in this situation. It is advantageous to apply this formula when $A = B + \alpha D$ and $B = UV^T$ and $U$, $V$ have $k \ll n$ columns.
Naturally, this does not prove that what you seek is impossible. However, there are simply so many people that would benefit so much from this that it seems unlikely that such a scheme hasn't been discovered at this point.
There is a very small chance that you could profit from the fact that $LU$ factorization of $A$ depends continuously on the parameter $\alpha$. If $\alpha_1 \approx \alpha_2$ and you have an LU factorization of $A(\alpha_1)$, then you may be solve $A(\alpha_2)x=b$ using the $LU$ factorization of $A(\alpha_1)$ and iterative refinement. Realistically, this will not yield a speedup on a modern computer as flops are cheap, but memory operations are expensive.
If speed is a critical issue, then your simplest course of action is to vectorize the solution of your differential equations using SIMD operations (search for AVX(2) instructions set, vector intrinsics, automatic vectorization + the name of your favorite compiler). Instead of solving 1 differential equation at a time you will solve, say, 8 equations at the same time. You will need control over the memory alignment of your data so that all arrays start at a cache line.
If you only have a single $\alpha$ then SIMD operations can still be applied if you adapt one of the standard algorithms for solving narrow banded linear systems in parallel. Examples include: cyclic reduction, the Spike algorithm.