Tell me more ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.
up vote 4 down vote favorite

I have a method that is supposed to be second order accurate based on asymptotic analysis/Taylor series expansion which assumes that the solution is smooth. I am solving a PDE which has very rough coefficients (one of them is a step function) and thus the solution of the PDE is not smooth, though it does exist and is unique. As a result, I observe only linear convergence. Will further mesh refinement eventually produce quadratic convergence? How should I formulate results of the numerical method in such case? Thanks! Edit: I am solving an equation of the type: $$u_t=x^2u_{xx}+1_{\{x\geq C\}}u_y$$ for some constant $C$ which is in the range of domain of computation $[0,x_{max}]$. I am solving that with Semi-lagrangian approach.

share|improve this question
Could you provide more details about the PDE you're solving, and the finite difference scheme you're using? – Paul♦ Jul 6 '12 at 17:26
You say you're using a semi-Lagrangian approach, presumably for the advection term. How do you discretize the diffusion term? – David Ketcheson Jul 7 '12 at 6:19
You might try modifying the problem by imposing smooth coefficients, and checking that you see second order convergence in that case. – David Ketcheson Jul 7 '12 at 21:08

3 Answers

up vote 5 down vote accepted

Like Jed says, if your coefficient is a step function, then the exact solution will be continuous but have a "kink" (i.e., the gradient will be discontinuous). Now, take the Laplace equation as an example (your heat equation behaves the same way, but the analysis is a bit more involved). There, you have the following analysis ($e=u-u_h$ is the error):

$$ \|\nabla e\|^2 = (\nabla e,\nabla e) = (\nabla e, \nabla (e-\varphi_h)) $$ where the last step is simply an application of Galerkin orthogonality, for any discrete function $\varphi_h$. If you choose $\varphi_h=I_h u - u_h$ where $I_hu$ is the interpolation of the exact solution onto the finite element space, then you see the following:

$$ \|\nabla e\|^2 = (\nabla e, \nabla (u-u_h-I_hu+u_h)) = (\nabla e, \nabla (u-I_hu)) \le \|\nabla e\| \|\nabla (u-I_hu)\|. $$ In other words, $$ \|\nabla e\| \le \|\nabla (u-I_hu)\|. $$ What this means is that your finite element solution can not be better in the $H^1$ seminorm than the interpolant. Using duality arguments, you will get estimates from this that are in the $L_2$ norm and look like this: $$ \|e\| \le C h \|u-I_hu\|, $$ i.e. you now need the interpolation estimate in $L_2$.

Now, if you think about approximating a function with a "kink" with finite element shape functions using interpolation, it quickly becomes clear that if you make the mesh fine enough, then the error will ultimately be concentrated in the one element in which the kink is located (unless the kink happens to be located on a mesh line). In any case, you will find in the literature the following estimates:

  • If $u\in H^1$ (e.g., if it has a kink), then $\|u-I_hu\| \le h \|\nabla u\|$.
  • If $u\in H^2$ (e.g., if it has no kinks), then $\|u-I_hu\| \le h^2 \|\nabla^2 u\|$.

In your case, the function has a kink and so is only in $H^1$, so you can't expect to be better than ${\cal O}(h)$ because that is the order you get from the interpolation estimate, and the error estimate for the finite element solution can not be better than the interpolation estimate.

share|improve this answer
up vote 4 down vote

If your problem is elliptic with variable coefficients, a discontinuous coefficient will make the gradient discontinuous. If the solution space has a continuous gradient at that point, the method cannot be better than first order accurate. No amount of refinement will fix that unless refinement causes the solution space to allow a discontinuity at that point.

Finite element methods allow such discontinuities most naturally. If you use a standard $C^0$ conforming finite element method and align the grid with the coefficient discontinuity, you should see second (or higher, depending on the basis) order accuracy. Unfortunately, corners in coefficients cause additional trouble because the solution appears on coarse grids to have a singularity similar to that near a reentrant corner in the domain. The apparent singularity will eventually be resolved with grid refinement, but the required resolution depends on the ratio of coefficients, which may be impractically fine.

There are many ways to enrich solution spaces so the internal discontinuity is represented without needing to conform to a background mesh. These papers provide a good perspective on classical treatment of material discontinuities.

More recent methods in the finite element context allowing much more general enrichments go by names like Generalized Finite Element Method (GFEM), eXtended FEM (XFEM), and Discontinuous Enrichment Method. In case the coefficient structure is oscillatory, random, or otherwise unresolvable and not possible to conform to, there are homogenization techniques, though first order accuracy is optimal for such multiscale problems.

share|improve this answer
When you say discontinuous gradient, so you mean $u_y$ or $u_x$ has a jump? When I plot the resulting function(I wish I could post it here) it looks pretty smooth and continuous except for a small region where $u$ is continuous with respect to $x$, but not differentiable and therefore the derivative $u_x$ has a jump. In fact, it is just one point on the grid, where it is such. When I refine the mesh the impact of this point gets more negligible and also gets dumped by the magnitude of the refined mesh, thus why it would matter as the mesh goes to zero? – Kamil Jul 6 '12 at 18:06
Yes, $u_x$ or $u_y$ has a jump. The problem is that even though the effect gets smaller as you refine, it only gets better as fast as the grid is refined (first order accuracy). You have to make the spaces match (by conforming to the interface or suitable enrichment) or tolerate first order accuracy. – Jed Brown Jul 6 '12 at 18:21
:thanks, I see what's going on. I wonder if this kink only comes from the fact that the coefficient is a step function. What if I am solving $u_t=u_xx+u_y$ but the initial value is $1_{\{y\geq C\}}$. The coefficients are infinitely smooth but no diffusion in the second variable might cause propagation of the kink as well of the diffusion in the first variable takes care of it and the solution is smooth? – Kamil Jul 7 '12 at 16:00
Another quick question: how much smoothness I should require for my coefficients then if I want to achieve the quadratic convergence? Suppose they are C^1 only, or even just $C^0$. My local truncation error has terms like $u_{ttt}$ and $u_{yyx}$, so that means I have to have a lot of regularity in the coefficients to make sure those derivatives are are bounded. Looks like even if they exist in a weak sense, that is a.e. they are finite except one point I have a problem with convergence... – Kamil Jul 7 '12 at 16:29
$u_t = u_x x + u_y$ has no diffusion at all. The discontinuity $1_{\{y \ge C\}}$ remains forever (being transported in the $y$ direction) and the PDE can only be interpreted in a weak sense. If you have non-degenerate diffusion, the solution will become smooth as soon as you leave the initial condition (though perhaps only smooth at an unresolved scale). For a model elliptic problem, flux matching at the interface gives the continuity requirement for an $H^2$ solution (ignoring boundary effects), but your equation is degenerate so your solution is not even $H^1$. – Jed Brown Jul 7 '12 at 17:57
show 4 more comments
up vote 1 down vote

Unless your mesh is extremely coarse right now, refining it will not change the order of your convergence. If you are getting good linear convergence, then it is likely your mesh is fine enough for the method to be asymptotic, and it will never change, no matter how fine your mesh gets (until you bottom out your residuals). I would recommend looking at the analysis of your derivative approximation and see if the assumption is made that the derivative being approximated is continuous.

I may be mistaken about this fact, but if it is based on some polynomial approximation of the function, this may cause the quadratic convergence to fail at singularities. This will drop your convergence from second order to some lower order at the point where you state your coefficient is a step function. Dropping the order at one point, be it an internal singularity or a boundary condition reduces the overall order to that same value.

Also, another possibility I thought of, are you sure that the method is second order global, and not just second order local? It is a simple mistake to make, and would explain the one order lower convergence.

If you give more detail as to what your PDE is and what method you are using to approximate your derivatives, it may help us give a more specific answer as to what could be causing the incorrect convergence rate.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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