I have been wondering for a while whether the correct way to design numerical integration schemes would be by tailoring them to capture the behaviour of the Green's functions of the differential equations.

I was asking myself this question when I came across strong numerical oscillations that were output from a heat equation with a finite element discretization (FE) in space and a Crank-Nicolson method (CN) in time. The oscillations occurred after the first couple of time steps in the temperature fields not far from the boundary where a very large (dirac-like) condition was imposed.

I asked myself whether these oscillations were due to the fact that the applied discretization in space and time allow only for approximate solutions that respect a separation of variables method (CN in time and FE in space); whereas the Green's functions to the heat equation are fundamentally not separable.

Therefore, to try to explain the oscillations we observed, my hypothesis is that our approximate numerical method, that seeks to estimate the fundamental solution via "special polynomials" in time and space, is not effective during the first time steps of our simulation because it is not capable of capturing the Green's function locally at those instants. Similar to the way the Gibb's phenomenon of a Fourier Series of a discontinuous function manifests.

Any clarifications and further thoughts on this are welcome, and any references to formal discussions of this are also welcome :)

Best