From: jonathan@geop.ubc.ca (Jonathan Thornburg)
Newsgroups: sci.math.num-analysis
Subject: Re: nonlinear wave equation on a sphere
Date: 11 Nov 1992 19:14:12 GMT
Organization: U of BC Astronomy  +  U of Texas at Austin Physics/Relativity
Distribution: inet
NNTP-Posting-Host: rubis.astro.ubc.ca
Summary: it's tricky
Keywords: partial differential equation PDE sphere finite difference pseudo spectral

In article <0096373F.07D40080@MV3600.BMEN.TULANE.EDU>
ramesh@MV3600.BMEN.TULANE.EDU writes:
>	I was wondering about how one numerically solves a nonlinear pde 
>(like the wave equation, or a nonlinear Helmholz equation) on a closed
>object like a sphere.  
>
>	My previous experience has been that a spectral method (like
>ORSZAG) is required.  This is eeasy to do using Fourier series and the
>appropriate trignoetric identities.  However, on a sphere, sphereical harmonics
>have to be used, for which there aren't the same type of identities.

I know this isn't what you wanted, but permit me to suggest
	\mybibitem{Boyd, J. P.}{1989}
		``Chebyshev \& Fourier Spectral Methods'',
		Springer-Verlag Lecture Notes in Engineering {\bf 49},
		Springer-Verlag,
		ISBN 3-540-51487-2 (Berlin),0-387-51487-2 (New York).
as being a truly superb book on this topic, with a very nice chapter
on "sphere" problems.  It also has quite a bit of useful stuff on
finite differencing methods.



>	Does anyone know how to set up a finite difference scheme for a 
>closed object like a sphere?

The basic answer is "very carefully" :-).

Seriously, it's tricky.

There are several ways to do it, depending on how many dimensions
your "sphere" is in, and how many of those you're treating numerically
(eg axisymmetry in 3D means you only treat r and theta numerically;
phi can be done analytically).  There are several problems which arise
if you just plunge ahead with a straightforward FD scheme in (r,theta,phi).
In the following laundary list of troubles I have personally struggled
with, I'll assume 3D, and I'll also assume that the underlying PDEs are
nonsingular in (x,y,z) coordinates:
- The continuum PDEs in polar coordinates probably have 0/0 terms
  on the $z$ axis.  You can generally fix these via L'Hopital's rule.
- If nature is unkind to you, the continuum PDEs in polar coordinates
  may have $\infty - \infty$ cancellations on the $z$ axis!  These
  can in theory be fixed, but it may take a lot of algebraic manipulation
  to do so.  (Perhaps try a symbolic algebra system?)
- If you have *tensor* PDEs, you may well have to worry about regularity
  conditions on the $z$ axis.  Enforcing these numerically is very
  tricky.  For an example of this (in numerical general relativity), see
	\mybibitem{Evans, C. R.}{1984}
		``A Method for Numerical Relativity:
		  Simulation of Axisymmetric Gravitational Collapse
		  and Gravitational Radiation Generation'',
		University of Texas at Austin Ph.D Thesis,
		University Microfilms Order \#DA 85-08264.
	\mybibitem{Evans, C. R.}{1989}
		``Enforcing the Momentum Constraints
		  During Axisymmetric Spacelike Simulations'',
		pp.~194--205 in \mycite{Evans, Finn, and Hobill}{1989}.
	\mybibitem{Evans, C. R., Finn, L. S., and Hobill, D. W., Eds.}{1989}
		``Frontiers in Numerical Relativity'',
		Cambridge University Press, Cambridge (UK),
		ISBN 0-521-36666-6.
- In axisymmetry there is another way to regularize polar-coordinate
  PDEs, involving a local coordinate transformation to a hybrid
  rectangular/polar coordinate system, but it's too (ugly and)
  complicated to describe here.  Contact me for the relevant appendix
  to my thesis if you're interested.
- For initial value (hyperbolic/parabolic) problems,
  assuming you get a regular system of FD equations on an (r,theta,phi)
  type grid, any explicit time integration scheme is going to suffer
  from a severe CFL time step limit caused by the close spacing of
  the grid points at near-polar latitudes.  Boyd talks about ways
  to fix this, and gives references to some of the tricks people
  use in numerical weather forecasting codes.
- Under some circumstances you can just ignore the spherical symmetry
  and use xyz coordinates and grids.  This is easy to program and
  mostly avoids the above hassles, but may be expensive.  You also
  have to worry much more about stability in hyperbolic/parabolic
  problems.

All in all, it's a messy subject.  Alas, nature is full of
topologically-2-spherical objects...

- Jonathan Thornburg
  <jonathan@geop.ubc.ca> through mid-November 92, then
  <jonathan@einstein.ph.utexas.edu> or <jonathan@hermes.chpc.utexas.edu>
  [for a few more months] UBC / {Geophysics & Astronomy, Physics}
  [then through Aug/92] U of Texas at Austin / Physics / Center for Relativity