# Efficient heat diffusion implementation with varying coefficients

I have the following heat diffusion equation: \begin{alignat}{3} \partial_t u(t, \vec{x}) &= g(\vec{x})\Delta u(t,\vec{x}), &\quad& \vec{x} \in\Omega, \, t\in(0,\infty],\\ \partial_n u(t,\vec{x}) &= 0, &\quad& \vec{x}\in\partial\Omega,\, t\in(0,\infty], \\ u(0,\vec{x})&=f(\vec{x}), &\quad& \vec{x}\in\Omega. \end{alignat}

I want to solve this on a regular grid, thus I use the finite difference method to discretise $$\Delta$$ with the Neumann boundary conditions with a matrix $$-W$$, and I discretise the functions $$f, g, u(t)$$ as the vectors $$f,g\in\mathbb{R}^n$$ and $$u:\mathbb{R}^+_0\to \mathbb{R}^n$$, defined on the grid with $$n$$ points. Setting $$G=\operatorname{diag}(g)$$ I can now rewrite the PDE as a system of ODEs:

$$\partial_t u(t) = -GWu(t), \quad u(0) = f.$$

The solution of this ODE at time $$T$$ is given as $$u(T) = \exp(-TGW)f$$. How would I go about implementing this efficiently in practice?

One option is of course to use the series expansion and truncate for a large enough $$m$$: $$\exp(-TGW)f = \sum_{k=0}^{\infty} \frac{(-TGW)^k}{k!}f \approx \sum_{k=0}^{m}\frac{(-TGW)^k}{k!}f.$$

Another option is to use the eigendecomposition of $$-GW = P\Lambda P^{-1}$$ to compute the matrix directly: $$\exp(-TGW) = P\begin{bmatrix} \exp(-T\lambda_1) & & \\ & \ddots & \\ & & \exp(-T\lambda_n) \end{bmatrix} P^{-1}.$$

However the latter sounds like a pretty bad deal compared to the equation with constant coefficients, since for the latter I have the eigendecomposition $$-W = V \Lambda' V^T$$ where I know the form of $$V$$ and $$\Lambda'$$ explicitly and there's even an $$O(n\log n)$$ algorithm for computing the application of $$V$$ and $$V^T$$ to a vector (it's a fast DCT-II and DCT-III transform). Is there any way the knowledge about $$V$$ and $$\Lambda'$$ can be used in the variable coefficients setting?

Edit: I guess I could set $$u(t) = G^{1/2}v(t)$$, and then I get (for $$G$$ invertible): $$\partial_t G^{1/2} v(t) = -GWG^{1/2}v(t) \implies \partial_t v(t) = -G^{1/2}WG^{1/2} v(t).$$ Then at least I get a symmetric matrix, so the eigenvalue and eigenvectors may be computed more efficiently. On the other hand I am not sure whether there is an obvious way to relate $$V$$ and $$\Lambda'$$ to the eigendecomposition of $$-G^{1/2}WG^{1/2}$$.

• How big is $n$ and how many $f$'s would you like to apply this to? If $n$ is large enough and you don't have many $f$'s then there are Krylov-based methods for computing matrix exponentials that would probably be pretty fast for this matrix Jul 4, 2023 at 14:34
• @whpowell96 Sometimes $n$ is small enough for to do the eigendecomposition directly, sometimes it is so large that there is no way I would be able to even store the eigendecomposition (maybe only the largest eigenvalues/vectors of it). The right hand side $f$ is only one. Do you have any suggestion for a Krylov-based algorithm for matrix exponentials? I wasn't even aware this existed, although it makes sense since generally I would use Arnoldi for the eigenvalues. Do you have any idea whether the knowledge of $V$ and $\Lambda'$ can somehow be leveraged here? Jul 4, 2023 at 15:00
• I'm not aware of state of the art methods but a description of approximating a matrix exponential matvec with an Arnoldi iteration is given in On Krylov Subspace Approximations to the Matrix Exponential by Hochbruck and Lubich. Also, I would double check your finite difference scheme, as spatially-varying diffusion coefficients are typically modeled as $\nabla\cdot(g\nabla u)$, not $g\nabla^2 u$. This is necessary to keep the operator Hermitian Jul 4, 2023 at 15:22
• The paper I linked is from the late 90s. SotA = state of the art Jul 4, 2023 at 21:08
• What is your intention when you avoid the time discretization? In fact, solving an ODE is a standard technique in approximating the (action of the) matrix exponential. See e.g. the famous Nineteen Dubious Ways to Compute the Exponential of a Matrix paper. Jul 6, 2023 at 20:06