2
$\begingroup$

I need to solve a PDE where one of the variables is an angle, so I need to know how to deal with periodic boundary conditions.

As a warm up, I am trying to solve the Helmholtz equation in polar coordinates, where the solutions should be Bessel functions multiplied by complex exponentials.

I tried the following code:

 NDSolve[{x D[u[x, t], x] + x^2 D[u[x, t], x, x] + D[u[x, t], t, t] + x^2 u[x, t] == 0, u[b, t] == 1, (D [u[x, t], x] /. x -> b) == 0.5, u[x, 0] == u[x, 2 Pi], (D[u[x, t], t] /. t -> 0) == (D[u[x, t], t] /. t -> 2 Pi)}, u[x, t], {x, b, 100}, {t, 0, 2*Pi}]

where i have set b=0.0001 because setting b=0 gives rise to singularities for some reason. One solution to this equation is BesselJ[0,x] + BesselJ[1,x], but it is not the only one: adding any combination of higher order Bessel functions one gets infinite other solutions.

The weird thing is that I noticed, by plotting the result, that Mathematica is giving me only BesselJ[0,x] as the solution! Anyone knows what might be happening?

$\endgroup$
  • 2
    $\begingroup$ u[b, t] == 1 guarantees that there is no azimuthal dependence, which in turn allows only BesselJ[0, x] as a solution. Also, Bessel's equation is singular at x == 0, and b > 0 is necessary to avoid it. $\endgroup$ – bbgodfrey Jun 6 '18 at 1:30
  • $\begingroup$ I'm voting to close this question as off-topic because the issue it raises is not really a Mathematica issue but a matter of the OP not having grasped the mathematics involved. $\endgroup$ – m_goldberg Jun 6 '18 at 19:59
  • $\begingroup$ @bbgodfrey Would you consider turning your comment into a short answer? I agree with the voters, that this is basically not a question about Mathematica, but a misunderstanding of the PDE and numerical solving in general. On the other hand, the Helmoltz equation is not that uncommon and we answered many strongly math-related questions in the past. Therefore, I'm a bit hesitant to cast the final close vote. $\endgroup$ – halirutan Jun 8 '18 at 0:30
  • $\begingroup$ @halirutan I would be happy to do so in a few days, when I have some time. Congratulations on your election as a Moderator. $\endgroup$ – bbgodfrey Jun 8 '18 at 1:25
  • $\begingroup$ @bbgodfrey Thanks and thank you for considering an answer :) $\endgroup$ – halirutan Jun 8 '18 at 1:36
6
$\begingroup$

As I remarked in comments above, the code in the question actually represents the Laplace equation, not the Helmholtz equation. It is necessary to set b > 0, because the equation is singular at x == 0, as can be seen by solving for D[u[x, t], {x, 2}]. Why it should yield only the zero-order Bessel function can be seen as follows. Let u[x, t] be represented by f[x] Cos[n t]], with n an arbitrary integer. Inserting this into the PDE gives

Expand[#/(Cos[n t] x^2) & /@ 
FullSimplify[
    (x D[u[x, t], x] + x^2 D[u[x, t], x, x] + D[u[x, t], t, t] + x^2 u[x, t] == 0)
    /. u -> Function[{x, t}, f[x] Cos[n t]]]
]

(* f[x] - (n^2 f[x])/x^2 + f'[x]/x + f''[x] == 0 *)

which is Bessel's equation of order n. Because the x == b boundary conditions are independent of t, only the n == 0 solution is permitted. That the PDE cannot have as an answer BesselJ[0,x] + BesselJ[1,x] also can be demonstrated by.

FullSimplify[
    (x D[u[x, t], x] + x^2 D[u[x, t], x, x] + D[u[x, t], t, t] + x^2 u[x, t] == 0)
    /. u -> Function[{x, t}, BesselJ[0, x] + BesselJ[1, x]]]

(* BesselJ[1, x] == 0 *)

Incidentally, numerically solving Laplace's equation as an initial value problem generally does not work well.

$\endgroup$
  • $\begingroup$ Thank you again for taking the time. +1 $\endgroup$ – halirutan Jun 11 '18 at 8:44

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.