Technical Notes > Smoothing Operators in PDE's

The Laplacian Operator as a Bandpass Filter Function

Assume that we have a function v(x) which we wish to smooth.

The Fourier expansion of this function is v(x) = ∑ Vk exp(i k x).

Let the smoothed function be u(x) = ∑ Uk exp(i k x), with k the angular velocity in radians per unit distance;

then the Laplacian of u is ∇2 u = ∑ (- k2 ) Uk exp(i k x).

We define u from the relation u - ε ∇2 u = v

then ∑ Uk exp(i k x) ( 1 + ε k2) = ∑ Vk exp(i k x).

Component by component, Uk exp(i k x) ( 1 + ε k2) = Vk exp(i k x)

Or, Uk = Vk /(1 + ε k2)

so that the kth frequency component is attenuated by a factor of 1/(1 + ε k2).

The Sampling Theorem states (McGillem and Cooper, "Continuous and Discrete Signal and System Analysis", p 164):

" A band-limited signal can be uniquely represented by a set of samples taken at intervals spaced 1/2W seconds apart, where W is the signal bandwidth in Hz."

The sampled signal is the product of the input signal and the sampling function, and the spectrum of the sampled signal is the convolution of the two transforms. The spectrum of the sampling function is a series of impulses at the harmonics of the sampling frequency (2W), and the convolution leads to a replication of the signal spectrum around each of these harmonics. If the signal bandwidth exceeds the harmonic spacing 2W, then the harmonics will overlap, and aliasing will occur.

From this we infer that if spatial data are available at a spacing of D meters, then the maximum bandwidth in the defined signal will be W = 1/(2D) cycles per meter, corresponding to k = 2πW radians per meter.

Combining these two items, we wish to infer a value of ε that will damp components of U with frequencies above W. However, the Laplacian filter does not have a sharp cutoff at any frequency, so we have some latitude in assigning ε.

Let us find ε such that the frequency component at frequency W is attenuated by a factor N, ie.

1/(1 + ε 4π2W2) = 1/N, with 1/(2W) = D.

Then ε = (N-1)/(4π2W2) = D2(N-1)/π2.

Arbitrarily choosing a frequency attenuation factor of N=2, we get ε = D2 /π2.

Smoothing Steady-state Solutions

In the solution of partial differential equation systems, it sometimes happens that auxiliary equations must be solved simultaneously with the PDE, and that these auxiliary equations have no spatial coupling, being point relations or other zero-order equations. In these cases, the finite element method works poorly, because the discretization is based on integrals over space, and oscillatory solutions can satisfy the integrals. In such systems, we are justified in adding to the equation a diffusion operator to impose a smoothing on the solution. If we have, for example,

U = F(..)

then we can replace this equation with

U - (D2/π2)∇2 U = F(..), with D the approximate spatial wavelength of acceptable oscillations.

Damping Time-dependent Systems

A similar analysis can be applied to time-dependent partial differential equations.

Suppose we have a system ∂v/∂t = f, in which the discretized equations support high frequency solutions which destabilize the numerical solution process. We wish to damp high frequency components.

Assume that v can be expanded as v(x,t) = ∑ Vk exp(i k (x-ct)), where c is a propagation velocity.

Let the smoothed function be u(x,t) = ∑ Uk exp(i k (x-ct)),

then the Laplacian of u is ∇2 u = ∑ (- k2 ) Uk exp(i k (x-ct)),

while the time derivative is ∂u/∂t = ∑ (- ikc ) Uk exp(i k (x-ct)).

We define u from the relation ∂u/∂t - ε ∇2 u = ∂v/∂t

then ∑ Uk exp(ik(x-ct)) ( -ikc + ε k2) = ∑ Vk exp(i k (x-ct))(-ikc).

Component by component, Uk = Vk (-ikc)/(ε k2 - ikc)

Or, |Uk| = |Vk |/sqrt(1 + ε2 k2/c2)

so that the kth frequency component is attenuated by a factor of

1/sqrt(1 + ε2 k2/c2).

Again defining W = 1/(2D) and seeking an attenuation factor of 2, we get

ε2 = (N2-1)c2/(4π2W2) = D2(N2-1)c2/π2 = 3D2c2/π2,

or approximately, ε = 2Dc/π.

We can now solve the equation ∂u/∂t - ε ∇2 u = f, with the expectation that u will be a frequency-filtered representation of v.

Steady-state limits of Time-dependent Equations

In some cases, a steady-state limit of a known time-dependent system is desired, but while the time-dependent equation itself is stable, the steady-state equation which results from merely setting the time derivative to zero is not. In these cases, we can replace the time derivative by -ε ∇2 u, again with the expectation that u will be a frequency-filtered representation of v.