Smoothing the diffusion coefficient to improve convergence

I have been reading a book by Thomee and he considers the case of $u_t=(au_x)_x$, for the case of $a$ possibly being discontinuous. Then he says that the problems with convergence might occur, and thus to fix that one can smooth the coefficient $a$ and restore convergence rate.

Please correct me if I am wrong but discontinuous coefficients affect the smoothness of the solution of a pde at the first place, and therefore affects the local truncation error of a numerical method. How does smoothing of the coefficient affects convergence if the convergence depends on local truncation error, which depends on the properties of the function itself and its derivatives by the original stated problem which we can't control, because it is an input to the problem?

• Please include a link or reference to Thomee's work so we know exactly what your question is based on. – David Ketcheson Jul 14 '12 at 18:28
• this is a book "Handbook of Numerical analysis, part 1" by Ciarlet and Lions. Thomee has a first chapter there. My question is related to page 172, where he considers discontinuous coefficient $a$ with finite number of discontinuities. He actually refers to the work of Samarskij and Fryazinov that the convergence is improved if $a$ is smoothed, but I would like to know where my argument is wrong. – Kamil Jul 16 '12 at 1:27

I don't know what exactly Thomee had in mind, but if you can't get solution at all for a discontinuous coefficient and a particular discretization scheme, you can smooth the coefficient on a lengthscale proportional to the mesh size $h$ and obtain a problem that's solvable. Then you make the smoothing distance go to zero as you refine the mesh and you obtain a sequence of numerical solutions that approximate a sequence of exact solutions of perturbed problems that converge to the solution of the original problem. That's a common trick.

• So, am I understanding you right that the solution might be discontinuous due to the bad coefficients and the scheme has a hard time to converge at points where discontinuities occur? And what you propose is the trick to make the scheme to converge in a limit to the discontinuous solution of a pde? – Kamil Jul 16 '12 at 2:37
• Yes, something like this. – Wolfgang Bangerth Jul 16 '12 at 17:03

Strong form discretizations are extremely cumbersome for the case of discontinuous coefficients. It is better to use weak form methods that do not require differentiating the solution. Furthermore while the product $a u_x$ is continuous even for discontinuous coefficients, $u_x$ alone is not. Therefore even evaluating one derivative of the solution is something for a numerical method to be very careful with. The standard (robust) approaches are $H(\textrm{div})$ mixed finite elements, mimetic finite differences, and primal finite volume/difference methods with a choice of quadrature making them equivalent to certain mixed methods.

I started my reply with the note above because a numerical discretization must be stable and conserve the correct properties for it to be worth trying to solve the resulting algebraic system. Make sure the spatial discretization is robust for the chosen purpose before proceeding further. Local conservation, for example, is critical for any problem exhibiting non-smoothness. Also confirm that your quadrature respects the non-smoothness.

When the coefficient has discontinuities, but when the jumps conform to the grid, you can use discrete spaces that have a discontinuous gradient at the interface, thus allowing accurate approximation without elaborate tricks. This is most natural in the finite element context, which partly explains its popularity for structural mechanics.

For the case of the coefficient $a$ being rough at a subgrid scale, there are numerous homogenization techniques. For problems with scale separation in 1D, the simple approach of replacing the rough coefficient with its harmonic mean on each element is optimal. For multiple dimensions or for problems that do not have separated scales, more sophisticated techniques are needed.

In all cases, simple smoothing of the coefficients degrades accuracy significantly, and you will have no hope of accurate large-scale solutions until the mesh resolves the fine scale structure and the scheme will generally not be convergent in the energy norm.

• Thanks, but that's quite sophisticated answer for me as I know about finite elements very little. However, I understand that bad coefficients affect the convergence results simply because the solution itself is not smooth and it is hard to make a numerical method converge to that solution at those points. But do I have a hope to have a convergence in any norm because it does depends on the regularity solution and if the derivatives are not bounded the smoothing of the coefficient for the numerical scheme is not going to help to make the solution smooth. Where is my logic wrong? – Kamil Jul 16 '12 at 2:45