1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ You claim you use finite differences, but you don't have a boundary value problem, as you have only one condition at $0$. $\endgroup$ – VoB Sep 12 at 7:55
  • $\begingroup$ @VoB I think he is considering periodic BCs. Then $u(0)=0$ is sufficient. $\endgroup$ – ConvexHull Sep 12 at 9:48
  • $\begingroup$ @Pedro Secchi You mentioned a symmetric domain. Do you use periodic BCs? $\endgroup$ – ConvexHull Sep 12 at 9:49
  • $\begingroup$ 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 $\endgroup$ – Pedro Secchi Sep 12 at 10:25
  • $\begingroup$ 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 $\endgroup$ – Pedro Secchi Sep 12 at 10:31
2
$\begingroup$

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

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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