Re: Heat Equation on a surface sphere using NDSolve?
- To: mathgroup at smc.vnet.net
- Subject: [mg127528] Re: Heat Equation on a surface sphere using NDSolve?
- From: Roland Franzius <roland.franzius at uos.de>
- Date: Thu, 2 Aug 2012 04:49:41 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- Delivered-to: l-mathgroup@wolfram.com
- Delivered-to: mathgroup-newout@smc.vnet.net
- Delivered-to: mathgroup-newsend@smc.vnet.net
- References: <jutlae$ciq$1@smc.vnet.net>
Am 27.07.2012 10:59, schrieb georgesabitbol4 at gmail.com:
> 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
with
Laplacian = d_u (1-u^2) d_u = d_u d_u - 2 u d_u + 1/(1-u^2) d_phi d_phi
equns=Simplify@
{
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],
Derivative[F][t,-1.0,p]==0,
Derivative[F][t,1.0,p]==0,
F[0,u,p]==F0[u,p]}
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
completely.
--
Roland Franzius