Numerical procedures for PDEs using Feynman-Kac Formula

Numerical procedures for PDEs
using the Feynman-Kac Formula


Project Leader: Prof. R. Jeltsch

Researchers: Prof. R. Jeltsch, Dr. W.P. Petersen, F. Buchmann

Date: 7.3.2002



Description:

We are developing algorithms for solving Dirichlet problems in bounded domains in n-space. The approach is probabilistic, which has the advantage that generalizations to high dimensions are in principle quite simple in concept. What is involved is to simulate stochastic differential equations in n dimensions and compute the partial differential equation solution using variants of the Feynman-Kac formula. Potential terms are also computed as numerical solutions to differential equations. Two difficulties are evident in this approach to the Dirichlet problem: finding the boundary in a reasonable number of steps, and determining accurately the process value and its exit time. A two layer algorithm is currently being used - a boundary finding layer, and an exit layer whose width is smaller than the weak order of approximation to the stochastic differential equation. Upon entering the boundary finding layer, a step size is chosen to prevent excursions. This means that the bounded approximation to the Brownian increment will keep the process within the domain or reach the boundary. Two approaches have been tested. A walk on spheres method requires modification of the stochastic calculus to confine the random walks to n-spherical surfaces. Our current favorite approach is more amenable to Runge-Kutta methods and is called walk on cubes . Using walk on cubes, the increment of each component of the Brownian movement is a three choice random variable. The step size is determined by the limitation that the process can at most reach the boundary of the domain.

Our results so far show the following features. First, the number of steps needed to find the boundary are O(1/h), where h is the starting step-size. Next, statistical errors due to a finite sample size N are the canonical O(N-1/2 ), independent of the dimension n of the space. In fact, it seems that the coefficient of N-1/2 seems to actually decrease with n. The number of steps required to find the boundary appears to be less sensitive to the initial step size as the dimension increases. All of these results are quite encouraging. Tests on Poisson problems with Dirichlet boundary conditions in hyperspherical domains have been run on problems as large as n=128, and have been carefully studied for n <= 64. Problems with hypercubic domains are currently under way. The algorithms are easy to implement, grid-free, and ideally suited for cluster architectures as they are embarrasingly parallel. Each sample path is completely independent of every other path, so any ordering of subsample sizes can be run on independent processors essentially without communication. Only the partial differential equation solution requires communication when all the paths have exited in O(1/h) steps.

In the following, results are shown for (1/2) Laplace(u) + 1 = 0 in various dimensions in a hyperspherical domain. All simulations were run on the Beowulf cluster Asgard.

Absolute error |u(0)-uN(0)| versus initial stepsize h. The boundary value is 0 on the hypersphere. The sample size is N=104
in the following dimensions: n=1 (circle), 2 (asterix), 4 (diamond), 8 (cross), 16 (triangle), 32 (+),64 (inverted triangle).






Early results are available in a Research Report No. 2002-01, Feb. 2002.



Contacts:
Dr. Wesley P. Petersen Fabian Buchmann
Seminar für Seminar für
Angewandte Mathematik Angewandte Mathematik
ETH Zürich ETH Zürich
CH-8092 Zürich CH-8092 Zürich
Switzerland Switzerland

Tel.: 0041 1 632 5575 Tel.: 0041 1 632 3533
Fax.: 0041 1 632 11 04 Fax.: 0041 1 632 11 04
Email: wpp@sam.math.ethz.ch Email: fab@sam.math.ethz.ch