comp! I'm trying to solve the advection-reaction problem

$ dg/dt = dg/dx + x\cdot g \qquad on~~x \in \Omega = (-\infty, +\infty)$

supplemented with the boundary conditions

$ \lim_{\lvert x \rvert \to\infty} g(x, t) = 0.0$

and subject to the constraint

$g(x=0, t) = 1$.

It seems to me that a non-trivial steady state solution exists for this system and is given by $g_{\infty}(x) = \exp(-0.5x^2)$. (Note that without the constraint at $x=0$ the PDE would be ill-posed as there would be infinitely many admissible solutions, each offset by a multiplicative factor.)

However, I'm unsure as to how I would go about solving this system numerically. My initial approach was to use an upwind stencil (i.e. $\delta_xg_k \approx \frac{1}{2\Delta x}(-g_{k} + g_{k+1})$) since it's an advection-like equation and then time march using implicit Euler. Moreover, I truncated the infinite domain by approximating it as $\widehat{\Omega} = [-10.0, +10.0]$ and modified the boundary conditions accordingly.

Using this scheme, I get the following results:

Numerical Solution using First Order Forward Differencing and Implicit Euler

Of course in retrospect this result isn't too surprising -- we can think of the original problem as two separate subproblems. Problem 'A' exists on the domain $x \in [-10.0, ~0.0]$ and has the outflow boundary conditions $g(-\infty, t) = 0$ and inflow boundary condition $g(0, t) = 1$ -- physically this problem makes sense, which explains why we obtain the correct response.

However, the remaining problem ('Problem B') exists on the domain $x \in [0.0, +10.0]$ and has the OUTFLOW boundary $g(0, t) = 1$ and INFLOW $g(+\infty, t) = 0.0$. Using my upwind stencil, I can never incorporate the boundary at $x = 0.0$ which explains why on the RHS we obtain the trivial solution. If I try to use a downwind or centered stencil though, the solution blows up since the corresponding semi-discrete form is unstable.

Surprisingly, even though an analytical steady state solution exists, it doesn't seem that it's possible to recover it numerically... Am I missing something in my understanding/Have you ever encountered a similar problem?

Any advice on this matter would be appreciated! Thanks in advance for your help.


2 Answers 2


The centered difference scheme should be stable for implicit Euler, so I'd try that first.

[deleted, misread question]

  • $\begingroup$ Thanks for your input, I'll definitely give centered difference another shot. However, I'm not quite sure about your remark about the wave speed changing signs -- isn't the wave speed (the coefficient in front of the spatial derivative) +1, and therefore constant? It's only the forcing function x*g that changes sign across the domain isn't it? $\endgroup$
    – user7208
    Feb 4, 2014 at 19:38
  • $\begingroup$ Oh quite right sorry, I did mis-read the question severely. I just quick sketched out the von neumann stability analysis (which you should do for this). It looks like central differencing with implicit Euler should be conditionally stable. $\endgroup$
    – Aurelius
    Feb 4, 2014 at 21:11

No, your numerical answer is correct.

The main point is that your problem is essentially advective, and particles will move to the left as time goes on, so if you start with the initial condition $g(t=0, x) = 0$, the right section (i.e. $x<0$) cannot acquire non-zero values.

Indeed your equation can be solved analytically. Consider the characteristic curves $(t(s),x(s)) = (s,-s+a)$, where $a$ is a constant. Then $dg/dt - dg/dx = (dg/dt)(dt/ds)+(dg/dx)(dx/ds) = dg(t(s),x(s))/ds$. Thus, $dg/dt-dg/dx - x g = dg(t(s),x(s))/ds - (a-s)g(t(s),x(s)) = 0$. In other words, along curves $(s,-s+a)$, the value of $g(t,x) = g(s)$ changes according to $dg(s)/ds = (a-s)g(s)$. The boundary condition at $x=0$ gives rise to "new materials", which then "travel to the left" and the amount f materials changes according to $dg(s)/ds = (a-s)g(s)$.

If at a given $s_0$, $g(s_0) = 0$, then $g(s) = 0$ for $s>s_0$. This is why. If $g(s_0) = 0$, then $dg(s_0)/ds = (a-s_0)g(s_0) = 0 $. Similarly $d^2g(s_0)/ds^2 = -g(s_0) + (a-s_0)g'(s_0) = 0$ and so on. This means the initially zero values on the right section are "transported" to the left as essentially zero values.

Thus how can the right section $x>0$ acquire non-zero values?


Your Answer

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

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