Quick von Neumann stability analysis

von Neumann stability analysis is a great tool to assess the stabilty of discretizing schemes for PDEs. But too often, imho, the discussion is too convoluted. Here, I try to provide a shortcut.

In a nutshell:

• The standard approach involves Fourier techniques, involving (of course) complex numbers
• The real part of these numbers is analysed, with some trigonometric expression resulting, identifying the troublesome modes
• I claim this mode can be identified in advance, which makes the whole Fourier procedure unnecessary

BTW: it’s pronounced “fon no ee man”

Discrete Fourier transform and Fourier series

This is quite silly, but the relationship between the discrete Fourier transform (DFT) and the Fourier series (FS) is clouded by annoying factors. I will try to connect both in this article. The motivation is to employ DFT techniques in a computer simulation. In the latter, one usually has a finite simulation box, which makes Fourier series the most interesting (a connection to the Fourier transform may also be made, see below).

Sound waves with attenuation

Just a simple derivation of the role of attenuation in the standard sound wave equation. Original work: Stokes, 1845.

Starting with the Navier-Stokes momentum equation

where  is a Lamé viscosity coefficient. The bulk viscosity coeficient  is defined as . The last term  is often neglected, even in compressible flow, but sound attenuation is one of the few cases where it may have some influence. All viscosities are assumed to be constant, but in this case this is a safe assumption, since we are going to assume small departures about equilibrium values.

Sparse Poisson problem in eigen

Back to scientific computing. Lately, I have been using the Eigen libraries for linear algebra. They were recommended by the CGAL project, and they indeed share a common philosophy. Thanks to the rather high C++ level they can accomplish this sort of magic:

  int n = 100;
  VectorXd x(n), b(n);
  SpMat A(n,n);

  fill_A(A);

  fill_b(b);

  // solve Ax = b
  Eigen::BiCGSTAB<SpMat> solver;
  //ConjugateGradient<SpMat> solver;

  solver.compute(A);

  x = solver.solveWithGuess(b,x0);

Notice that A is a sparse matrix! I am next describing how to use this in order to solve the 1D Poisson equation.

Approximation for weighted residuals

Intro

We ended up the article about weighted residuals pointing out three possible approaches to the problem. We expand them here.

The Rayleigh-Ritz method

Intro

This method starts straight from the weak form of a differential equation. It follows by taking the same functional space for the solution functions and the weight functions. This makes it similar to the straight Galerkin method – indeed it is often equivalent, as we will discuss in a future article.

The weak form

Intro

One can write a general differential equation as:

An example, which we will often consider is

ex.

to model e.g. heat conduction with a space-dependent heat conduction coefficient and a heat production term . (Notice the “ex” when we are talking about a particular example.)

This is pretty technical, so please stop reading here if you are already lost!

The method of weighted residuals

Intro

One can write a general differential equation as:

An example, which we will often consider is

ex.

to model e.g. heat conduction with a space-dependent heat conduction coefficient and a heat production term . (Notice the “ex” when we are talking about a particular example.)

This is pretty technical, so please stop reading here if you are already lost!

Math fonts for (La)TeX

How could I have overlooked this Survey of Free Math Fonts for TeX and LaTeX! An excellent read, it shows the author has himself been involved in font design. As I mentioned in my entry on XeTeX,  it is great to have many fonts for text (and XeTeX really takes a leap forward in this direction). But, fonts in mathematics should match.

Natural coordinates

Here is an interesting application of Voronoi tesselations / Delaunay triangulations (see previous post The alpha shapes for another one.) Suppose you have a set of points carrying some information, let’s call them particles; the simplest case is just a scalar number representing some field whose value is only known for the particles. Then, you want to compute the value of the field in some other, arbitrary point. (In the simplest instance, you may want to plot the field, and you need to know the values at nodes of a grid.)