MathGroup Archive 2012

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Heat Equation on a surface sphere using NDSolve?

  • To: mathgroup at
  • Subject: [mg127528] Re: Heat Equation on a surface sphere using NDSolve?
  • From: Roland Franzius <roland.franzius at>
  • Date: Thu, 2 Aug 2012 04:49:41 -0400 (EDT)
  • Delivered-to:
  • Delivered-to:
  • Delivered-to:
  • Delivered-to:
  • References: <jutlae$ciq$>

Am 27.07.2012 10:59, schrieb georgesabitbol4 at
> Hello Folks,
> I am new to Mathematica, but it seems to be the most suitable tool for my issue.
> I have to solve a equation which is similar to heat equation on a sphere surface. The exact expression in a spherical coordinates system is (in LaTeX):
> \frac{\partial{h}}{\partial t} = \frac{1}{n} \left[ \frac{1}{r^2 \sin^2{\phi}} \frac{\partial}{\partial \theta} \left( K \frac{\partial h}{\partial \theta} \right) + \frac{1}{r^2 \sin \phi} \frac{\partial}{\partial \phi} \left( K \sin \phi \frac{\partial h}{\partial \phi}\right) + s(\theta,\phi,t)\right]
> K is kind of a diffusion coeff.
> s a source term (varying in time and space)
> and n a constant (although it might depends on \phi and \theta).
> \theta and \phi are defined here in a mathematical way.
> r the radius (constant as we deal with a perfect sphere).
> (no variation regarding the radius, we only look at the surface of the sphere).
> I found that in Mathematica, and in a Cartesian frame, such equation (somewhat not exactly the same but similar) would be:
> N = 1
> K = 34
> myfunc = NDSolve[{D[h[x, y, t], t] ==
>   1/N*(D[h[x, y, t], x, x] + D[h[x, y, t], y, y]), h[x, y, 0] == 0,
>   h[0, y, t] == K*t, h[9, y, t] == K*t, h[x, 0, t] == K*t,
>   h[x, 9, t] == K*t}, h, {x, 0, 9}, {y, 0, 9}, {t, 0, 6}];
> I'm kind of stuck right now, because I do not know how can I switch to the spherical frame. I would appreciate any ideas/suggestions.
> Thanks a lot and sorry for bothering you guys with this.

Looks a bit confusing.

The Laplacian on the unit sphere with metrics

ds^2 = dtheta^2 + sin^2 theta dphi^2

is  with s -> sin theta

Laplacian = s^-1 d_theta s d_theta + s^-2 d_phi d_phi

and the general heat equation with production and destruction terms is

d_t F(t,theta,phi) =  diffconst Laplacian( F(t,theta,phi) )
                       -  F(t,theta,phi)/lifetime + productionrate

to avoid the sine switch to cylinder coordinates

u = Cos theta,  du = -Sqrt(1-u^2) dtheta,  d_theta = - Sqrt(-1-u^2) d_u


Laplacian = d_u (1-u^2) d_u = d_u d_u - 2 u d_u + 1/(1-u^2) d_phi d_phi

D[F[t,u,p],t] == K D[(1-u^2)D[F[t,u,p],u]],u] +
   K/(1-u^2) Derivative[0,0,2][F][t,u,phi]
   - a F[t,u,p] + b,
   F[t,u,2 Pi]==F[t,u,0],

NDSolve[equns,{t, 0.0, tmax}, {u,-1.0, 1.0}, { p, 0.0, 2.0 N[Pi]}]

The polar points u=+-1 are replaced by reflecting boundary conditions. 
This guarantees smothness of the solution at the poles and conservation 
of probability.

In case the start distribution F0 is rotational symmetric around some 
axis you can choose this direction as north and discard the variable phi 


Roland Franzius

  • Prev by Date: Re: Data Format
  • Next by Date: Re: Mathematica as a New Approach to Teaching Maths
  • Previous by thread: Re: mac new OS "mountain lion", compile function,
  • Next by thread: Strange Compile runtime error