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}$.

  • 1
    $\begingroup$ 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 $\endgroup$
    – whpowell96
    Jul 4, 2023 at 14:34
  • $\begingroup$ @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? $\endgroup$
    – lightxbulb
    Jul 4, 2023 at 15:00
  • 1
    $\begingroup$ 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 $\endgroup$
    – whpowell96
    Jul 4, 2023 at 15:22
  • 1
    $\begingroup$ The paper I linked is from the late 90s. SotA = state of the art $\endgroup$
    – whpowell96
    Jul 4, 2023 at 21:08
  • 1
    $\begingroup$ 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. $\endgroup$
    – davidhigh
    Jul 6, 2023 at 20:06


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.