I'm looking at Murray's book: Mathematical biology: an introduction , first volume, pag. 404

In particular, I'm interested to solve the following PDE: $$\partial_t u = \partial_x (\text{sign}(x) u) + \partial_x (u^2\partial_x u)$$

and I chose as initial condition $u_0(x)=e^{-x^2}$ and boundary conditions $u(-L,t)=u(L,t)=0$

I want to use linear finite elements to solve it, so I followed the approach in this question that I asked yesterday.

Using this, I obtain the following plot at $t=2$, and I'd like a check.

enter image description here

I also have snapshots at different times:

enter image description here


After @cos_theta comments, I obtain the following solution at $t=2$: enter image description here

  • 1
    $\begingroup$ Cross-posted: mathematica.stackexchange.com/q/230429/1871 $\endgroup$
    – xzczd
    Sep 20, 2020 at 14:12
  • $\begingroup$ Do the oscillations reduce when the grid size is decreased? Looks like an instability to me. $\endgroup$
    – cos_theta
    Sep 22, 2020 at 14:09
  • $\begingroup$ Yes, they descrease always more. Also with Mathematica it's the same. I asked a new question (scicomp.stackexchange.com/questions/35967/…) using your comments. Could you give a check please? @cos_theta $\endgroup$
    – Vefhug
    Sep 22, 2020 at 14:11

1 Answer 1


First off, I'd note that your initial condition doesn't satisfy the boundary conditions, so you might want to instead use $u_0(x) = e^{-x^2} - e^{-L^2}$.

A great sanity check for problems like yours is the conservation property -- the total mass of $u$ should stay the same.

$$\begin{align} \frac{d}{dt}\int_{-L}^Lu\, dx & = \int_{-L}^L\frac{\partial u}{\partial t}dx \\ & = \int_{-L}^L\frac{\partial}{\partial x}\left(vu + u^2\frac{\partial u}{\partial x}\right)dx \\ & = \left(vu + u^2\frac{\partial u}{\partial x}\right)\Big|_{x=-L}^{x=L} \\ & = 0, \end{align}$$

because you've assumed that $u = 0$ at both endpoints. Here I've written $v = \text{sign}(x)$ for the advection field, but this relation would hold true regardless of what $v$ was as long as you had the same nonlinear diffusion coefficient. From the plot you've shown, it looks as if the numerical solutions you're obtaining are monotonically decreasing, which would violate the conservation property. That suggests that there's an error in your numerical implementation somewhere.

When I run into problems like these, I usually try and come up with a simpler system and see if I can solve that first. For example, what happens if you take out the advection term? The PDE

$$\partial_tu = \partial_x(u^2\partial_xu)$$

is challenging enough by itself -- it's a free boundary problem. Similarly, what happens if you take out the diffusion term and then smooth over the advection field? Can you get a good approximation to the solutions of

$$\partial_tu = \partial_x(\tanh(x/\epsilon)u)$$

for different values of $\epsilon$? Start with $\epsilon = L / 2$ and then see how things go as you decrease it to be equal to the mesh spacing $\delta x$. You might even be able to write down an analytical solution using the method of characteristics. Both of these simplified problems have conservation principles and other intrinsic mathematical properties that you can use as sanity checks.

  • $\begingroup$ Thanks. In the question I asked yesterday, I solved $\partial_t u = \partial_x (u^2 \partial_x u)$ (scicomp.stackexchange.com/questions/35940/…), so I just want to extend that FEM method to my case. What changes is the weak formulation, because I have the integral $$\int_{-L}^{L} \partial_x (\phi_j \text{sgn}(x)) \phi_i dx$$ and I wrote it as $$\int_{-L}^{L} \text{sgn}(x) \phi_j^{'} \phi_i dx$$ so in some sense I didn't consider the derivative of $\text{sgn(x)}$. I think this is the source of the problem, but I don't know how to handle it $\endgroup$
    – Vefhug
    Sep 20, 2020 at 17:29
  • 1
    $\begingroup$ First, split the integral in order to correctly treat the discontinuity at $x = 0$: $$\int_{-L}^L \partial_x\left(\mathrm{sgn}(x)\phi_j(x)\right)\phi_i(x)\,\mathrm{d}x=\lim_{a\to 0-}\int_{-L}^a \partial_x \left(\mathrm{sgn}(x)\phi_j(x)\right)\phi_i(x)\,\mathrm{d}x+\lim_{b\to 0+}\int_{b}^L\partial_x\left(\mathrm{sgn}(x)\phi_j(x)\right)\phi_i\,\mathrm{d}x$$ Integrate by parts: $$\int_{-L}^a \partial_x\left(\mathrm{sgn}(x)\phi_j(x)\right)\phi_i(x)\,\mathrm{d}x=\left[\mathrm{sgn}(x)\phi_j(x)\phi_i(x)\right]_{-L}^a-\int_{-L}^a\mathrm{sgn}(x)\phi_j(x) \partial_x\phi_i(x)\,\mathrm{d}x$$ Apply BCs. $\endgroup$
    – cos_theta
    Sep 20, 2020 at 19:04
  • $\begingroup$ Well, so the first term of your r.h.s. is $- \text{sgn}(a) \phi_{J}(a)$ where $\phi_J(x)$ is the hat function with value $1$ at $x=0$, right? @cos_theta But I have some issues in understanding what is the true value: I'd say that as $a \rightarrow 0^{-}$, $- \text{sgn}(a) \phi_{J}(a) \rightarrow +1$ correct? $\endgroup$
    – Vefhug
    Sep 20, 2020 at 19:17
  • $\begingroup$ Let's consider $\phi_i=\phi_j$, where the hat function is centered at $x=0$ with value $1$ and slope $\pm 1/h$ (otherwise, the non-integral term vanishes): $$\mathrm{sgn}(x)\phi_j(x)\phi_i(x) \to -1\quad(x\to 0^+)\\ \mathrm{sgn}(x)\phi_j(x)\phi_i(x) \to 1\quad(x\to 0^-)\\ -\int_{-h}^0\mathrm{sgn}(x)\phi_j(x)\partial_x\phi_i(x)\,\mathrm{d}x=\frac{1}{2}\\-\int_{0}^h\mathrm{sgn}(x)\phi_j(x)\partial_x\phi_i(x)\,\mathrm{d}x=\frac{1}{2}$$ So, we get $$\int_{-h}^h \partial_x\left(\mathrm{sgn}(x)\phi_j(x)\right)\phi_i(x)\,\mathrm{d}x=-1$$ $\endgroup$
    – cos_theta
    Sep 20, 2020 at 20:51
  • $\begingroup$ Yes, I worked out the computations and I found the same thing you wrote. So this is the diagonal entry of the matrix (for the point corresponding to $x=0$, which should be in the middle of my discretization points). The other entries of the matrix should be the usual ones, right? @cos_theta $\endgroup$
    – Vefhug
    Sep 20, 2020 at 20:56

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.