The Method of Fundamental Solutions and a Taste of Python

Informal Numerical Analysis Seminar

University of Bath

Friday 11 April 2008

Introduction

• BICS talk by Jon Trevelyan last Monday (7 April)
• About method for solving the Helmholtz equation
• very inspiring
• got me thinking

Summary of talk by Jon Trevelyan

• Helmholtz equation
• boundary integral equation formulation
• standard finite elements (FE) for function on boundary
• boundary element method BEM
• enrich with plane wave expansion
• decompose solution as
• in general write solution as linear combination of basis functions
• collocation of the boundary integral equation
• gives rise to system for unknowns on the boundary
• curve (1D) for 2D, surface (2D) for 3D
• Dirichlet/Neumann
• building matrix, involves solving integrals
• his approach subdivide integration domain, Gauss along wave crests,
• rule for oscillatory integrals (steepest descent) in direction of wave
• 1D integral for 2D, 2D integral for 3D

Basis functions

• idea: can we choose beter basis functions?
• sidestep
• standard finite element (FE) method (e.g. for elliptic equations)
• Dirichlet problem
• choose basis functions that satisfy boundary conditions (BC) (but don't solve PDE)
• e.g. triangulation and hat functions for nodes in interior
• then do Galerkin or collocation inside domain
• homogeneous Helmholtz equation
• basis functions that solve PDE? (but don't satisfy boundary conditions)
• for homogeneous Helmholtz (Laplace) this is possible
• 2D Hankel function \(H_0^{(1)}(k r)\)
• spherical wave 3D \(\exp(i k r)/r\)
• point source in 3D
• cylindrical wave, line source in 3D
• essentially point source in 2D
• fundamental solutions
• approximate solution solves PDE in domain
• need to find coefficients such that boundary conditions (BC) satisfied
• collocation
• very similar in flavour to boundary element method (BEM)
• but no integration necessary

Detailed explanation of method (see notes)

Further discussion

• of course I wasn't the first to come up with this
• some searching
• memory of talk at Dundee
• method of particular solutions for PDE eigenvalue problems
• paper by Trefethen
• work by Moler, Matlab logo
• Method of Fundamental Solutions
• paper by Barnett and Betcke
• authors used Matlab
• they make the link with the boundary integral formulation (BEM)
• analysis for interior Helmholtz for disc
• convergence, conditioning, stability
• conditioning of system, size of coefficients
• solution as analytic function
• e.g. point source outside of domain
• best to choose points "inbetween" singularity and domain
• choice of curves

Advantages

• interior or exterior problems
• automatically satisfies boundary conditions at infinity (Sommerfeld radiation conditions)
• very easy to implement
• works in 2D, 3D, ...
• meshless

Disadvantages

• conditioning
• but B&B article shows that this can be alleviated by
• choosing source locations well
• also problem for other methods

Implementation

• How hard is this to implement?
• walk through python code
• show that we can very easily get very similar results to Matlab
• many libraries available
• here we use numpy, scipy and matplotlib
• module/file doc string
• importing modules
• numpy, scipy, pylab (matlab-like plotting interface to matplotlib)
• importing specific functions, submodules (, variables, classes)
• defining new functions
• scipy documentation
• tutorials, some reference documentation, api docs, wiki, source
• here api docs were useful
• scipy.special Hankel, derivative of Hankel, Bessel
• unit_circle, circle (only points and gradients needed, meshless)
• fundamental solutions
• source using Bessel function of the second kind (\(Y_0\))
• helmholtz class
• example: unit disc, bc based on Bessel point source outside
• Dirichlet/Neumann
• numerical values for parameters
• plots using matplotlib, very similar to Matlab

Demonstration

• Figures 4a, 4b
• Figures 3a, 3b, 3c
• other examples: discs with Dirichlet and Neumann conditions
• Animations

Idea for further exploration

• talk at Gene Golub memorial day by Godela Scherer
• about work with Gene Golub and Victor Pereyra on optimisation
• Gene Golub QR and SVD for LS
• generalised to separable nonlinear optimisation
• VARPRO/PORT
• paper illustrates that we want to avoid growth of coefficients
• minimise not just ||r||^2 but ||r||^2 + \nu ||c||^2
• standard regularisation
• Golub showed how to this for QR
• more stable for ill conditioned problems
• normal equations kappa^2, QR kappa
• same equations as before, but in addition n equations of the form \(\sqrt{\nu} c_j = 0\)