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. I also have snapshots at different times: EDIT:

After @cos_theta comments, I obtain the following solution at $$t=2$$: • Cross-posted: mathematica.stackexchange.com/q/230429/1871 Sep 20 '20 at 14:12
• Do the oscillations reduce when the grid size is decreased? Looks like an instability to me. Sep 22 '20 at 14:09
• 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 Sep 22 '20 at 14:11

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.

• 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 Sep 20 '20 at 17:29
• 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. Sep 20 '20 at 19:04
• 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? Sep 20 '20 at 19:17
• 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$$ Sep 20 '20 at 20:51
• 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 Sep 20 '20 at 20:56