2
$\begingroup$

I know very little about numerical PDE, so I apologize if this equation is ill-formed (and I appreciate any corrections).

I am currently learning about Brownian motion. A classic result is that we can use Brownian motion hitting times to solve the heat equation. More specifically, if we consider a nice domain in $\mathbb{R}^n$ with smooth boundary, and release Brownian motion from a point $x$ in the domain, then one can compute the solution to the heat equation with specified boundary conditions by averaging over the points that the Brownian motion hits (weighted according to the hitting probability). More sophisticated version of this idea are possible with more complicated PDEs (e.g. survival time, and Feymann–Kac stuff).

I would like to program a PDE solver to simulate this problem myself. That it, I would like to determine the Brownian motion hitting times by solving the corresponding PDE. (This is more accurate than just running a Monte Carlo simulation with BM because there's no Monte Carlo error.) I am most interested in dimensions $n=2,3$, although I would like to be able to go higher.

I would like to do this with very curvy shapes: ellipsoids, spheres, a U, a filled-in smiley face, etc.

My understanding is that the last method rules out finite difference methods, or at least makes them very difficult to program (due to needing some extrapolation technique I do not really understand). The books I've flipped though seem to concentrate on rectangular domains.

I was recommended FD methods as the simplest way to get started with PDE solving. What is the next-simplest method (which gives relatively good results and is easiest to program) that will allow me to solve these problems on curvy domains?

Improve this question
$\endgroup$

3 Answers 3

Reset to default
2
$\begingroup$

For complex domains like what you describe, I think that the finite element method is basically the only way to go. Spectral discretizations are nice when the spatial domain is simple, but your domains are complex. The finite volume method works well for hyperbolic conservation laws but is arguably less suited to elliptic or parabolic-type problems. I won't attempt to summarize finite element or other Galerkin-type methods here; instead I'll just refer you to my favorite books on the subject, one by Braess and the other by Gockenbach.

How simple it is to code this up can vary greatly depending on how much you want to do yourself. I've written mesh data structures, sparse matrix data structures, and assembly kernels for piecewise linear finite elements completely from scratch. This is reasonable to do and I found it very educational, but it will be time-consuming. At the other end of the spectrum is using high-level packages like FEniCS, Firedrake*, scikit-fem, sfepy, etc. To use these packages, you need to understand the weak form of the PDE and to be able to use a mesh generator. They will insulate you from needing detailed knowledge about how to turn the weak form of the PDE into a matrix-vector equation. Using these packages makes development very fast, but can leave you feeling like it's a big opaque box where you don't understand what's going on inside it. The deal.II package lies somewhere in the middle. You're partly responsible for the finite element matrix assembly, but the package provides a host of common element families, has already defined good mesh and matrix data structures, will deal with quadrature for you, etc.

Which of these choices you opt to go with is a decision you have to make. How much work do you want to invest in this and to what degree do you feel that it's important to understand the lower-level details?

No matter what, you will need a mesh generator. The two most popular choices are gmsh and Triangle.

*I contribute to Firedrake and receive grant funding to develop a downstream application built on it.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks a lot for your detailed answer. I have a basic worry. The entire point of this exercise for me is to beat the Monte Carlo error from BM simulation (like N^{-1/2} in number of samples, so pretty big). With a finite element method, there will be error from the mesh approximation. How do people convince themselves that the mesh error is negligible in practice for complicated domains? Are there theoretical (rigorous) results available? Or do they just, e.g., solve it with a sequence of finer meshes and eyeball convergence? Basically, I want to know how to do error control. $\endgroup$
    – alligator
    Feb 18, 2022 at 15:33
  • $\begingroup$ Aha so you're wondering about a posteriori error estimation and control. The blunt approach is to solve your problem at multiple resolutions and then seeing how different they are. You can be much more exact than that though. Wolfgang wrote a book on the subject which I highly recommend. It's a little harder to find but I also like this book by Ainsworth and Oden. $\endgroup$
    – Daniel Shapero
    Feb 18, 2022 at 16:29
  • $\begingroup$ Awesome, thank you! $\endgroup$
    – alligator
    Feb 18, 2022 at 17:36
4
$\begingroup$

One can certainly solve problems with curved boundaries using the finite difference method, but it is awkward.

It is simpler to use the finite element method. There, you need to subdivide your domain into a mesh, typically consisting of triangles in 2d and tetrahedra in 3d. These need not conform to the curved boundary, you just need that the vertices of the triangles/tetrahedra lie on the curved boundary -- which the commonly used mesh generators all do. Then, even though your mesh consists of straight-edged triangles/tetrahedra, you are guaranteed that the finite element solution converges to the exact solution of the PDE.

One can of course improve things (and get better accuracy) by using triangles and tetrahedra that have curved edges and faces. But that is non-trivial. You don't want to do this yourself, though you could use one of the widely used finite element libraries and see whether they support such curved cells. (One choice is deal.II, which supports curved cells -- both using polynomial approximations and using the exact geometry of your domain; I am one of the Principal Developers of this library, so consider the conflict of interest in my recommendation.)

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you. Are there results that allow one to control the error of the solution to the PDE on the mesh relative to the solution on the curvy domain? For example, theorems like: if the boundary is Lipschitz with constant C, then a certain mesh approximation will get me within 1% error. If so, this would definitely meet my needs (and I would be thankful for a pointer to a good textbook). $\endgroup$
    – alligator
    Feb 18, 2022 at 0:15
  • $\begingroup$ Only a posteriori error estimates are able to guarantee a specific percent error. These exist, certainly for the Laplace equation, for polygonal domains. I don't know whether anyone has taken the time to extend them to boundary approximation errors. $\endgroup$
    – Wolfgang Bangerth
    Feb 18, 2022 at 4:11
1
$\begingroup$

If you can define a mapping that transforms a rectangle into the domain you want to use, then you could use the simple methods described in this paper. As an example of a "very curvy" domain that can be handled, see for instance this figure.

It's focused on hyperbolic problems, but includes a description of how to discretize parabolic or elliptic problems (see Section 7.2). It's ostensibly a finite volume method, but a 2nd-order accurate finite volume method is essentially equivalent to a finite difference method.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.