How do you solve linear systems whose solutions decay exponentially?

Question

Consider the heat equation

$$\dot{u} = \Delta u$$

with initial conditions

$$u_0 = \delta(x)$$

for some point $x$ in the domain $\Omega$ of the problem. If $\Omega$ is $\mathbb{R}^n$, then this problem has a closed-form solution given by the Euclidean heat kernel

$$k_t(x,y) = \frac{1}{(4\pi t)^{n/2}}e^{-|x-y|^2/4t};$$

in general solutions to this problem will exhibit this same type of behavior: the magnitude of $u$ decays roughly exponentially as you travel away from the "source" $x$.

Now consider solving this problem numerically by constructing a finite-dimensional linear system whose solution approximates the solution of the original PDE (e.g., using a Galerkin method). Typical numerical linear algebra systems (which work with fixed-precision floating-point arithmetic) guarantee that the maximum error in any component of the solution will be no larger than some small fraction $\epsilon > 0$ of the largest element of the right-hand side.

For instance, in the problem outlined above the right-hand side will look like a Kronecker delta, so the absolute error will be no worse than $\epsilon$. Unfortunately, since the magnitude of the solution decays exponentially, the error $\epsilon$ may be much larger than the smallest entry of the true solution.

Question: using floating-point computations only, can one obtain a solution to a linear system that obtains a desired accuracy relative to the magnitude of each element of the solution vector, rather than the right hand side?

The fundamental problem in achieving better accuracy seems to stem from cancellation effects, i.e., if you add two numbers of very different magnitude in floating point, the smaller one is essentially ignored. I am aware of algorithms that use alternative numerical representations (e.g,. Dixon's method and so on), but this kind of answer is not particularly interesting to me due to considerations of efficiency.

Answer 1

This subject has been studied in literature --- it turns out that for some class of matrices (such as M-matrices) you can obtain componentwise accurate solutions. A good starting point is the section on "accurate floating point computation" on Jim Demmel's home page; other good search terms are "accurate linear algebra" and "componentwise error analysis". Recent literature has focused more on accurate eigenvalue computation rather than linear systems, but there are results also for them.

EDIT: And let me add that discretizing $\Delta$ with the usual finite-difference scheme results in an M-matrix.

Answer 2

I assume that you're interested in solving the boundary value problem in a situation where the boundaries are far away from the source, so that the solution is close to the analytical solution on an unbounded domain.

One option here would be to write the solution to your problem as the sum of analytical solution to the problem on an unbounded domain plus a correction term. That is,

$ u(x,y,t)=k_{t}(x,y)+v(x,y,t) $

Since the PDE is linear, your correction term $v(x,y,t)$ would also have to satisfy the heat equation. You could easily derive boundary conditions for $v$ in terms of the boundary conditions for $u$ and the values of $k_{t}(x,y)$ at the boundaries.