Achieving high relative accuracy (vs. absolute accuracy) using spectral methods

Question

Problem

I have a PDE that I'm trying to solve with spectral methods. The solution $y$ is always positive, and decays as $y \propto e^{-ax}$ for large $x$. The domain is $[0, \infty)$. (There are actually two independent variables but let's not worry about that right now.)

Unfortunately, although whatever numerical solution I obtain (using the Galerkin method) seems to have good absolute accuracy, it has poor relative accuracy towards larger $x$ due to the large dynamic range of the solution (coming from the exponential decay). I'm wondering how best to improve the relative accuracy of my solution. Any suggestions would be appreciated!

Example test problem

This is not the actual PDE I'm trying to solve (which is actually 2d), but it illustrates the difficulty I'm having:

$$(1 + e^{-x})y' = -y + 2xe^{-x},\quad y(0) = 0.$$

The analytic solution is $x^2/(1+e^x)$ (plotted as dashed red below), which is always positive. However, solving this with a Laguerre function basis set (recombined so that each basis function is $0$ at $x = 0$) of order 20 yields the blue curve below, which diverges from the analytic solution (and even becomes negative) when it becomes small.

What I've tried/considered

• I considered making the transformation $y = e^{-A(x)}$ and then solving for $A$. This would probably work, but unfortunately makes the equations nonlinear, which I was hoping not to deal with.
• I tried using Laguerre functions as a basis set, which have exponential decay built in. This didn't really work, though, as you can see above -- solutions still have low relative accuracy towards larger $x$.
• I'm planning to substitute $y = e^{-ax} z$ and then solve for $z$. I don't know $a$ a priori, but can estimate it. I think this will work all right, but I'm wondering if there's a more general method that I'm unaware of.

Answer 1

So I went ahead and implemented a code in Matlab that can solve this problem using a spectral approach, utilizing a simple polynomial basis. Using a simple polynomial basis of order 20 resulted in the following:

One thing I did do that you might not have was define the basis to exist on a domain that might have better numerical stability (especially as x increases). I made it so my basis was defined on the domain $\zeta \in [-1,1]$. Then I used the mapping between $\zeta$ and $x$ to modify the integrals and do the necessary computation needed to find the coefficients for the basis.

Now based on my solution, I question if you have some implementation error that leads to the inaccurate solution you have. I would expect your solution to be at least as good as what my simple polynomial basis produces, assuming implementation is correct.