# Parity for artificial dissipation term in a finite-difference solution

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?

• You claim you use finite differences, but you don't have a boundary value problem, as you have only one condition at $0$.
– VoB
Sep 12, 2020 at 7:55
• @VoB I think he is considering periodic BCs. Then $u(0)=0$ is sufficient. Sep 12, 2020 at 9:48
• @Pedro Secchi You mentioned a symmetric domain. Do you use periodic BCs? Sep 12, 2020 at 9:49
• Hi, yes, as I mentioned above I’m using u(0)=0 as a boundary condition. I’m not using a periodic BC, however. The problem I’m trying to model is the boundary layer over a symmetric airfoil, the thickness of which is supposed to be an even function. The equation is not exactly formatted as the one above (would have to add a function multiplying the derivative term which I kept out of the question for simplicity), but the mathematical considerations I made above regarding the signal are still valid Sep 12, 2020 at 10:25
• PS the boundary layer problem is not supposed to be periodic. Is the u(0)=0 BC still not enough then? I thought it would be enough to close the system Sep 12, 2020 at 10:31

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$$.