Parity for artificial dissipation term in a finite-difference solution

Question

I have a doubt regarding the signal of the dissipative term in a finite difference solution for an equation of the form

$$ \frac{\partial u}{\partial x}+f(x)=0, u(0)=0 $$

In which $f$ is an odd function and the solution (judging by the modelled phoenomena) should be an even function in a symmetric domain (say, $D=[-1; 1]$).

If I'm using the finite difference approximation

$$ \frac{\partial u}{\partial x}=\Delta_h u_i=\frac{u_i-u_{i-1}}{x_i-x_{i-1}} $$

The second order Taylor series decomposition of $u_i$ and $u_{i-1}$ results in a dissipative term proportional to $\Delta x \frac{\partial^2 u}{\partial x^2}$. Thus I'm trying to solve the equation by defining the residual as

$$ r_i=\Delta_hu_i+f(x_i)+\alpha \Delta^2_hu_i \times \Delta x $$

With $\alpha$ being a hyperparameter by which I can regulate how much artificial dissipation I want to apply.

I'm confused about the parity of the artificial dissipation term, however. Without it, the residual is an odd function in $D$; however, if $u$ is an even function as I said above, $\Delta x \frac{\partial^2 u}{\partial x^2}$ should also be one: therefore after it is added to the residual function to account for dissipation, the residual will not be odd when provided with an even $u$. In my specific case, this leads to $u$ converging to a physically impossible, uneven funcion after some Newton iterations when $\alpha\neq 0$.

Am I getting the idea of the dissipation term wrong?

Answer 1

I think you can understand this using the concept of a modified equation. As you have shown, your discretization

$$ \frac{u_i - u_{i-1}}{h} + f(x) = 0 $$

is an approximation of your differential equation $u_x + f(x) = 0$ with accuracy $\mathcal{O}(h)$. Now consider the Taylor expansion of your exact solution

$$ u(x_{i-1}) = u(x_i) - h u_x(x_i) + \frac{1}{2} h^2 u_{xx}(x_i) + \mathcal{O}(h^3) $$

and plug it into your finite difference approximation. This results in

$$ u_x(x_i) + f(x) + \frac{h}{2} u_{xx}(x_i) + \mathcal{O}(h^2) = 0 $$

Therefore, while your approximation is an $\mathcal{O}(h)$ approximation of your original differential equation, it is an $\mathcal{O}(h^2)$ (that is: more accurate!) approximation to the modified differential equation

$$ u_x + \frac{h}{2} u_{xx} + f(x) = 0 $$

Note that in the limit $h \to 0$, both become identical but for every finite $h$ they will be different. Using https://www.wolframalpha.com/, we can solve the modified equation for the case $f(x) = 1$ to illustrate the effect this has. The solution is

$$ u(x) = c_1 \frac{h}{2} \exp(-\frac{2 x}{h}) - x + c_2 = 0 $$

Note that if $h = 0$, you get $u(x) = -x + c_2$ which is the correct solution to $u_x + 1 = 0$. Your boundary condition $u(0) = 0$ alone is not enough to fix both $c_1$ and $c_2$ - thus, your scheme will somehow implicitly fix a second BC, I am not sure how. But unless it is precisely in a way that $c_1 = 0$, your actual solution will look like what is shown in the figure.

Note how the exponential term makes it a non-even function (neither is it odd). So while your exact solution in the limit $h\to0$ is even, your numerical scheme does not replicate this property for finite mesh sizes $h > 0$.