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.

enter image description here

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.
  • $\begingroup$ Are you certain you don't have any errors in your code? I just solved this example problem with an Adaptive Runge-Kutta method and it solves it fine, so I don't know why the spectral approach would have difficulties unless there's implementation error. $\endgroup$ – spektr May 31 '16 at 4:48
  • $\begingroup$ Thanks for looking into it! I used a fixed spectral order of 20, as I mentioned, so imperfect accuracy is expected. If I increase the order, the solution certainly becomes better. My issue is wanting to go for roughly uniform relative accuracy, rather than what seems to be uniform absolute accuracy. I see your point that switching to finite differences would probably help, but was looking specifically for advice on spectral methods. $\endgroup$ – Josh Burkart May 31 '16 at 6:45
  • 1
    $\begingroup$ Your problem has a known exact solution. You can get arbitrarily good accuracy by using that. What are you really interested in solving? $\endgroup$ – David Ketcheson May 31 '16 at 11:39
  • 1
    $\begingroup$ @JoshBurkart The FAQ says to please only post real problems you face, and for good reason. You're likely to get perfectly good answers (like "use the exact solution") that don't apply to your real problem. $\endgroup$ – David Ketcheson Jun 1 '16 at 15:43
  • 2
    $\begingroup$ @DavidKetcheson I did post about a real problem I am facing. I then reduced it down to a demonstrative test problem that others could reproduce, which I labeled as "example test problem". I think this is a generally useful practice when requesting assistance, since my actual problem has lots of extra complexity that distracts from the issue I'm wondering about. Since there was confusion, I'll try to edit my post to make it clearer. $\endgroup$ – Josh Burkart Jun 1 '16 at 16:50

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:

enter image description here

enter image description here

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.

  • $\begingroup$ Wow thanks so much for looking into this in detail! Very helpful to have a reference to compare to. I'll look into this more and get back to you. (I do the same thing with remapping to [-1, 1] when using Chebyshevs as a basis.) Any pointers on how to look for the kind of subtle bugs you're suspecting might be present? I have dozens of unit tests etc. but bugs are certainly possible. $\endgroup$ – Josh Burkart Jun 1 '16 at 5:19
  • $\begingroup$ @JoshBurkart No problem! Just a fun thing for me to do. But in terms of debugging, I suppose it depends on what your unit tests don't cover. So maybe first checking functionality that isn't being tested in some unit test and making sure it's doing what you expect. I would really aim at things like your basis being implemented right and your integration being correct. $\endgroup$ – spektr Jun 1 '16 at 5:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.