For a project I'm doing, I have to numerically solve a nonlinear parabolic stochastic partial differential equation, of the form
$$ u_t = u_{xx} + f(u)(u_x)^2 + a(u) + b(u)W(t, x), $$
where primes are derivatives with respect to $x$, $W(t, x)$ is space-time white noise, $f$, $a$ and $b$ are smooth and in general nonlinear. The equation is usually solved with a Dirichlet boundary condition at $x = 0$ and a Robin-type boundary at $x = 1$, of the form $u'(t, 1) = g(u(t, 1))$.
Now, if $f(u) = 0$, this is easy enough to solve; I've used finite differences as well as finite elements to do so. But problems arise when this is not the case. The software I'm using approximates the derivatives with finite differences, which I'm surprised works at all.
As I understand (hand waviness incoming), the noise introduced by using finite differences is sort of 'canceled out' when averaging over many ensembles when the derivatives are linear. The quadratic term now amplifies the finite difference error even more, and it no longer cancels when taking averages.
Are there any methods for dealing with nonlinearities like this in SPDEs? I've been scouring the internet for the last couple of days, but can't seem to find anything that is directly relevant.
I posted this on /r/math before, and got a couple of suggestions:
- TVD schemes as used in CFD for hyperbolic equations,
- Nonlinear finite elements,
- Invariant embedding,
- Nonlinear Feynman-Kac formulas,
- Using an adaptive grid.
Thanks in advance!
EDIT: As per request, some details on what I suspect is going wrong:
When solving the above equation for a given grid spacing $h$ in $x$, I can obtain a solution just fine. However, decreasing $h$ and solving again, the solution is different; it has the same general shape each time but will be shifted up or own depending on the function $f$.
Now, this only happens $f(u) \neq 0$. I have also tested equations where with a term linear in $u_x$ rather than quadratic, and they do not show this behavior.
The reason I suspect the finite difference approximation to $u_x$ is simply that $u$, being the solution to an SPDE, is not differentiable. Taking the FD derivative anyways results in very noisy data, which is then amplified by taking the square.
I do not have much experience with SPDEs, so I may be completely off base here. If anyone has any other suggestions on what might be at play here, I would certainly appreciate it.
EDIT 2: As a concrete example, the simplest SPDE of this type we're trying to solve is
$$ u_t = u_{xx} - (u_x)^2/u + \sigma u W(t, x), $$
with $u(t, 0) = 0$ and $u_x(t, 1) = \left(a - \frac{\sigma^2}{2}\right)u(t, 1)$. We know that the analytical mean solution for $t \rightarrow \infty$ to be
$$ \langle u(\infty, x) \rangle = \exp(ax). $$
The following image shows the result we obtain directly solving this equation, with $\sigma = \frac{1}{2}\sqrt{2}$ and $a = -1$. The grid spacing uses $N = 30$ points in $x$ (in the application we're looking at, $t \mapsto \tau$, $x \mapsto t$), and $10^4$ steps in time.
The dotted line represents the analytical steady-state solution, while the solid lines are the mean $\pm$ the sampling error. Error bars are an estimate of the time step error.
We have also taken an alternative approach (we just got this result yesterday): since $u\partial_{xx} \ln u = u_{xx} - (u_x)^2/u$, we solved
$$ u_t = u(\ln u)_{xx} + \sigma u W(t, x). $$
Solving this directly gives us the solution shown below, which agrees with the analytical result within the sampling error.
