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.