Research Group of Prof. Dr. M. Griebel
Institute for Numerical Simulation
Title Finite-differences on Sparse Grids for the Solution of Partial Differential Integrations
Participiants Thomas Schiekofer Michael Griebel
Keywords Finite Differences, Partial Differential Equations, Sparse Grids

Sparse Grids

We consider the solution of elliptic and parabolic partial differential equations on the domain and together with appropriate boundary conditions, respectively. Usually, the discretization of these equations on the respective domains is done by the use of full grids with equidistant meshwidths in each coordinate direction. Using the piecewise linear hierarchical basis, a full grid can be obtained by taking a rectangular scheme of the infinite tableau of subspaces into account as it can be seen on the left side of Figure 1. A discussion of the properties of the subspaces together with a comparison of accuracy vs. effort for each subspace (a detailed discussion can be found in [1]) leads to the result, that it is more reasonable to use a triangular scheme as indicated on the right side of Figure 1 which leads to a sparse grid instead of taking a rectangular scheme into account.

Figure 1: rectangular scheme (left) and triangular scheme (right) of subspaces, cut out of a theoretically infinite tableau of subspaces (l denotes the level of discretization)

We are now interested in the discretization of the above mentioned differential equations (the discretization of hyperbolic pde's is treated within the project Finite Difference Schemes on Sparse Grids for Time Dependent Problems ) using finite differences. The application of well known finite difference stencils leads to a non-consistent discretization (cf. [4]). Therefore, special finite difference operators for sparse grids have to be build. These operators are based on dimensional splitting and the (nodal to) hierarchical basis transform H and its inverse H-1. Each differential operator DjSG (which stands for first or second derivatives on sparse grids, respectively) along a coordinate direction j consists of three operations:

  • application of a hierarchical transformation in all coordinate directions but the one in which the derivative has to be applied,
  • the well-known one-dimensional finite difference stencils are applied on each single one-dimensional grid line of the sparse grid according to the local meshwidth,
  • and the last operation is the application of the back-transformation from hierarchical values to nodal values again in all coordinate directions but the one in which the finite difference stecils was just applied

Figure 2 shows this procedure in the two-dimensional case where a derivative in x-direction is considered. A more detailed description of this scheme can be found in [2] and [4].

With this is mind, DjSG simply reads . If DiSG is considered to be a discretization for the second derivative, the expression

represents the Laplacian operator. Note that, the discretized differential operators lead to non-symmetric matrices even for self-adjoint differential operators like the Laplacian. Hence, we need an iterative solver that can handle systems of linear equations with non-symmetric matrices. The solver we implemented therefore is the preconditioned BiCGstab method. Although the discretization matrices are not sparse, a matrix-vector multiplication can be implemented in O(N) where N denotes the number of grid points. In the following, some numerical examples are presented. Details about discretizations, preconditioning, properties of the discretization matrices, proofs of consistency and so on can be found in [4]. There, more numerical examples including also the Navier Stokes equations can be obtained.

elliptic PDE's

As a first example, we consider a Laplace equation on the unit square with the solution

and right hand side f=0. Additionally, Dirichlet conditions are imposed on the boundary. The numerical results can be seen in Table 1 where the error measured in the maximum norm and the L2-norm is shown. Figure 2 shows the corresponding convergence rates. Here, also the L1-norm as well as a pointwise measured error can be obtained.

Table 1: numerical results for above Laplace problem

The second example we focus at is a convection diffusion equation

with b=(2,3) and c=1. The solution of this problem is chosen to be

together with the apropriate right hand side. Again, Dirichlet boundary conditions are imposed. Table 2 shows the numerical results where again the errors measured in the maximum norm and the L2-norm can be seen. Figure 3 shows the corresponding convergence rates. Again, also the L1-norm as well as a pointwise measured error can be seen.

Table 2: numerical results for above convection diffusion problem

parabolic PDE's

We also treat parabolic equations

on the domain with boundary condition for and initioal condition for u. We discretize the time derivative by the first order explicit euler method, but other time discretizations like Crank-Nicolson or full implicit methods can also be applied. As an example we take

for the elliptic operator L with the coefficient functions

The coefficient function is chosen in a way so that

is the solution of our problem where describes the initial value for u. Hence, the solution of this problem consists of a product of sine waves in x- and y-direction which propagate with different speed in time. The following pictures where obtained on a regular sparse grid with time step 1e-4.

Figure 2: function and contour plot for t=0

Figure 3: function and contour plot for t=1.5e-1

Figure 4: function and contour plot for t=3.25e-1

Related projects