8
$\begingroup$

I am trying to solve partial differential equation in spherical coordinates $(\theta,\phi)$,but I don't know how to properly include boundary conditions of $\theta$. For $\phi$ it is periodic, but $\theta$ - I don't really know.

Here is my code:

M = 50;
z0 = 0.0;
ϕ0 = Pi;
tmax = 10;
pde = {D[Q[t, ϕ, θ], t] == -Cos[θ] D[Q[t, ϕ, θ], ϕ] + 1/M Sin[θ] D[Q[t, ϕ, θ], θ, ϕ],
       Q[t, 0, θ] == Q[t, 2 Pi, θ],
       Q[t, ϕ, 0] == Q[t, ϕ, 2 Pi],
       Q[0, ϕ, θ] == 2^-M (1 + Cos[θ] z0 + Cos[ϕ0 - ϕ] Sin[θ] Sqrt[1 - z0^2])^M};
sol = NDSolve[pde, Q, {t, 0, tmax}, {ϕ, 0, 2 Pi}, {θ, 0, 2 Pi}];

$\phi \in [0,2\pi]$ and $\theta \in [0,\pi]$. Boundary condition should reflect spherical geometry.

Analytical solution (quantum mechanics derivation):

M = 5;
Q[t_, θ_, ϕ_] := (Sum[Binomial[M, k]*(Cos[θ/2]*Cos[θ0/2])^k*(Sin[θ/2]*Sin[θ0/2])^(M - k)*Cos[(ϕ - ϕ0)*(M - k) + t/(4*M)*(2*k - M)^2], {k, 0, M}])^2 + 
(Sum[Binomial[M, k]*(Cos[θ/2]*Cos[θ0/2])^k*(Sin[θ/2]*Sin[θ0/2])^(M - k)*Sin[(ϕ - ϕ0)*(M - k) + t/(4*M)*(2*k - M)^2], {k, 0, M}])^2;

with initial conditions:

θ0 = Pi/2;
ϕ0 = Pi;

You can check for yourself:

D[Q[t, θ, ϕ], t] + Cos[θ]*D[Q[t, θ, ϕ], ϕ] - 1/M*Sin[θ]*D[Q[t, θ, ϕ], ϕ, θ] // Simplify
$\endgroup$
  • $\begingroup$ Er… so what's your boundary condition of θ, in traditional math notation? BTW I think it'll be better if you mention some background information of your PDE. $\endgroup$ – xzczd Jan 24 '15 at 6:49
17
$\begingroup$

The differential equation should reflect the fact that θ ranges between 0 and Pi. Because the spherical coordinate system becomes degenerate at these two pole, Q should be independent of ϕ there. Thus, from the differential equation, Q is time-independent at the poles. Hence, the boundary values of Q at 0 and Pi are the initial values there. Based on these considerations, the code becomes

pde = {D[Q[t, ϕ, θ], t] == -Cos[θ]*D[Q[t, ϕ, θ], ϕ] + 1/M*Sin[θ]*D[Q[t, ϕ, θ], θ, ϕ],
       Q[t, 0, θ] == Q[t, 2*Pi, θ],
       Q[t, ϕ, 0] == (1 + z0) 2^-M, Q[t, ϕ, Pi] == (1 - z0) 2^-M,
       Q[0, ϕ, θ] == 2^-M*(1 + Cos[θ]*z0 + Cos[ϕ0 - ϕ]*Sin[θ]*Sqrt[1 - z0^2])^M};
sol1 = NDSolveValue[pde, Q, {t, 0, tmax}, {ϕ, 0, 2*Pi}, {θ, 0, Pi}];

The solution can be displayed in a number of ways. For instance,

Manipulate[Plot3D[sol1[t, ϕ, θ], {ϕ, 0, 2 Pi}, {θ, 0, Pi}, PlotRange -> All, 
  AxesLabel -> {ϕ, θ, Q}, AxesStyle -> Directive[Bold, Black, FontSize -> 16]], 
  {t, 0, 10, Appearance -> "Labeled"}, ContinuousAction -> False]

A sample plot, for t = 10, is

t = 10


Kuba: couldn't resist :)

First, solved till t=40.

Manipulate[
 SphericalPlot3D[ 1 + sol1[t, ϕ, θ], {θ, 0, Pi}, {ϕ, 0, 2 Pi},
   ImageSize -> 500, PlotRange -> 2, BoxRatios -> 1, 
  ColorFunction -> (Blend["TemperatureMap", #6] &)], {t, 0, 20, .1, 
  Appearance -> "Labeled"}, ContinuousAction -> False]

enter image description here

Update

As noted by the OP in his addition to the Question, the numerical solution presented above agrees only qualitatively with his analytical solution. The accuracy of the numerical solution can be improved by reducing MaxStepSize and making the usual transformation, ` x = Cos[θ]. The revised equations

pde1 = {D[Q[t, ϕ, x], t] + D[x Q[t, ϕ, x] + (1 - x^2) D[Q[t, ϕ, x], x]/M, ϕ] == 0,
   Q[t, 0, x] == Q[t, 2*Pi, x], Q[t, ϕ, 1] == ((1 + z0)/2)^M ,
   Q[0, ϕ, x] == ((1 + x*z0 + Cos[ϕ0 - ϕ]*Sqrt[1 - x^2]*Sqrt[1 - z0^2])/2)^M};
sol1 = NDSolveValue[pde1, Q, {t, 0, tmax}, {ϕ, 0, 2*Pi}, {x, -1, 1}, 
   MaxStepSize -> {.005, .12, .054}]

now yield for t = 10 the plot

revised t=10

which agrees well with the analytical result. In fact, numerical and analytical results agree (at least for M = 50) up to about t = 15, above which they drift apart. Thus, at t = 20 the numerical result is

t=20

whereas the corresponding analytical result consists of eight identical, uniformly spaced solitons.

The numerical results are most sensitive to MaxStepSize for x. Unfortunately, reducing it below 0.054 triggers a numerical instability. I have tried numerous options for Method without improving the outcome. Of course, there may be other options that I have not tried that may produce superior results. Note, though, that I have found no sample problem in the voluminous NDSolve documentation that includes a term of the form D[Q[t, ϕ, x], x, ϕ], which appears to make this problem distinctive.

$\endgroup$
  • 3
    $\begingroup$ @Kuba, What a great addition. Thanks. $\endgroup$ – bbgodfrey Jan 24 '15 at 13:31
  • 1
    $\begingroup$ brilliant animation @Kuba! :) $\endgroup$ – WoofDoggy Jan 24 '15 at 18:28
  • $\begingroup$ Unfortunatelly, this numerical solution does not overlap with analytical form of the expression I have derived - it resembles it, but time scales and behaviour is different. I tried to increase WorkingPrecision, but an error appeared: "The precision of the differential equation (...) is less than WorkingPrecision (5.)." $\endgroup$ – WoofDoggy Jan 27 '15 at 13:05
  • $\begingroup$ What is you analytic solution? $\endgroup$ – bbgodfrey Jan 27 '15 at 13:14
  • $\begingroup$ I have updated the question $\endgroup$ – WoofDoggy Jan 27 '15 at 14:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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