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.

  • $\begingroup$ People use finite difference for much more complex SPDEs, so I'm surprised that you think finite difference is going to fail here. I understand what you mean from nonlinear SPDE, but in practice, due to the fact that your SPDE is time-dependent you should not to be worry about those nonlinear terms at all. I solve a similar equation but for a completely different application, which is called phase-field method, see it in section 4.6, p 35: physics.mcgill.ca/~provatas/papers/Phase_Field_Methods_text.pdf. You need to be more specific if you see any problem with finite difference. $\endgroup$ Feb 17, 2020 at 3:37
  • 2
    $\begingroup$ Here's deep insight from decades of computational science: You can't find the cause of a problem if you can't clearly describe its effects. The statement "But problems arise when this is not the case." is not such a clear description and none of us will be able to point out what precisely the problem might actually be. $\endgroup$ Feb 17, 2020 at 3:38
  • $\begingroup$ @AloneProgrammer Could you elaborate a bit on what you mean by not having to worry about the nonlinear terms if the SPDE is time dependent? $\endgroup$
    – user34137
    Feb 17, 2020 at 4:39
  • $\begingroup$ This shift up or down is huge or just small? What if you take an ensemble average of your solution? Does it converge? $\endgroup$ Feb 17, 2020 at 13:32
  • $\begingroup$ Also, I don't see any sarcasm in Wolfgang's comment. I think the burden is on you to describe your problem with sufficient details to other people understand it and if they see a fit give you an answer. Otherwise, I second Wolfgang's comment and say just simply saying the problem arise without any explanation won't clear the question and would be really difficult for other people to help you here. $\endgroup$ Feb 17, 2020 at 14:12

1 Answer 1


As it turns out, the $(u_x)^2$ term is ill-defined due to the irregularity of $u$, which is caused by the space-time white noise. This causes discrete approximations to converge to solutions to different equations, in a similar phenomenon as is discussed here.

The solution is to use Martin Hairer's regularity structures, which provides a framework for renormalizing equations such as this one. In the case of coupled SPDE systems of this form, the methods outlined here can be used.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.