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.
Why is my second order accurate method only converging at first order when the coefficients are rough?
pde;finite difference
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.