Solving wave equation with spherical coordinates using NDSolveValue?

Questions
, , ,
1

I am trying to use NDSolveValue for a fairly simple setup (used for testing that I understand NDSolveValue which apparently I don’t). I have a wave equation (for a vector potential) where I use a plane wave propagating in the -x direction as initial conditions. It works fine in Cartesian coordinates where I only need to solve for 1 function. But if I switch to spherical coordinates where I have two functions (the r and phi components of the vector potential) then the solution is no longer correct. I wonder if it is because I am not extracting the results from NDSolveValue correctly; if I compare the initial conditions that are used by NDSolveValue then these do not correspond to the “solution” I get out at the initial time. So I am wondering if anyone might have an idea as to what is going wrong - am I for instance extracting the output from NDSolveValue incorrectly? The notebook code is shown below.

In[1]:= Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]

eqr = -r r D[Ar[t, r, phi], t, t] + 
   r r D[Ar[t, r, phi], r, r] + D[Ar[t, r, phi], phi, phi] + 
   2 r D[Ar[t, r, phi], r] - 2 r D[Aphi[t, r, phi], phi] - 
   2 Ar[t, r, phi] == 0

eqphi = -r r r D[Aphi[t, r, phi], t, t] + 
   r r r D[Aphi[t, r, phi], r, r] + D[Aphi[t, r, phi], phi, phi] r + 
   4 r r D[Aphi[t, r, phi], r] + r Aphi[t, r, phi] + 
   2 D[Ar[t, r, phi], phi] == 0

a[r_, phi_] := 1
omega := 0.2

Aric = 
 Ar[t, r, phi] == a[r, phi] Sin[phi] Cos[omega (r Cos[phi] + t)]

Aphiic = 
 Aphi[t, r, phi] ==  a[r, phi] Cos[phi] Cos[omega (r Cos[phi] + t)]/r


{Arsol, Aphisol} = 
 NDSolveValue[{eqr, eqphi, 
   DirichletCondition[{Aric, Aphiic}, t == 0]}, {Ar, Aphi}, {t, 0, 
   100}, {r, 0.001, Sqrt[2] 100}, {phi, 0, 2 Pi}, AccuracyGoal -> 24, 
  PrecisionGoal -> 10]

For any rvalue and phivalue in the intervals [0.001,Sqrt[2] 100] and [0, 2 Pi] I would expect that Arsol[0,rvalue,phivalue] and Aphisol[0,rvalue,phivalue] should be equal to the expressions indicated in Aric and Aphiic evaluated at those same rvalue, phivalue. But they are not. Why?

2

I know very little about PDEs, but your equations seem singular at r==0. Also, they do not seem compatible when t==r==0:

{eqr,eqphi,Aric}/.{t->0,r->0}

{-2 Ar[0,0,phi]+(Ar^(0,0,2))[0,0,phi] == 0,
2 (Ar^(0,0,1))[0,0,phi] == 0,
Ar[0,0,phi] == Sin[phi]}
3

Thank you for pointing out the possible issues. I don’t think it should be a problem that the equations and initial conditions naively appear to be in discordance upon evaluation at t = 0, r = 0 . As you point out, the equations+i.c.'s are singular at r = 0 so it’s not valid to simply set terms proportional to r equal to zero just because r = 0 (note also that I do not attempt to solve the equations at r = 0 - I go somewhat close to r = 0 but it was a typo in the post that r = 0 was a boundary point - the typo has been fixed now). There was a typo in one of the equations though which has now fixed but it has not made any difference.

4

I can’t really follow the notebook. Are you using a different version of Mathematica or is that just the way it displays on the web?

Nonetheless, with spherical coordinates, set up your equations to solve for the radial and angular parts separately. Don’t try to solve for all of them in the same NDSolveValue statement.

5

Thank you for your post. I use Mathematica 12.2. I haven’t included the “In[x]” parts of the lines except in the first line (accidentally) but other than that, it looks like what I see in my notebook.

I’m not sure what you mean by solving the equations for the radial and angular parts separately - do you mean using the separation of variables method often used for this type of equations or do you mean solving for A^r first and then A^phi afterwards? The latter can’t be done in general so is not very useful and the former seems overly complicated since I am asking Mathematica to give me a numerical solution and not an analytical one. I also don’t see why it should be necessary? Surely Mathematica should be able to solve this problem as is?

6

Ok, I see. I have the same problem with solving equations even in cartesian coordinates. However, I am generally looking for an analytical solution, but I sometimes simply cannot solve the equations in one DSolveValue statement. I have to get creative, and break it up somehow. I wish I would have caught this at the beginning of the weekend so I’d have time to play with it. Nonetheless, I’m curious if anyone finds a solution to your problem. I apologize I couldn’t help.

7

OK, well then that may be something I need to look at as well; for this particular problem the Cartesian version only has a single unknown so maybe I will get a similar problem in Cartesian coordinates as in spherical coordinates if I try rotating the coordinate system to have two components in Cartesian coordinates as well - I will look at this tomorrow!

8

I’ve found out that one problem with my use of NDSolveValue is that using DirichletCondition might not be appropriate; I thought I’d checked that it worked in Cartesian coordinates but I must have mixed up data files exported from Mathematica because after checking again it looks like those conditions lead to errors even in the Cartesian case. So I’ve switched (back) to using another type of initial conditions. With these the Cartesian case works and I can also add an extra component to solve for without problems (the equations aren’t coupled in Cartesian coordinates though so maybe there would still be an issue if they were coupled). The following shows what works in Cartesian coordinates:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]
eqz = D[Az[t, x], t, t]  - D[Az[t, x], x, x] == 0
eqy = D[Ay[t, x], t, t]  - D[Ay[t, x], x, x] == 0
omega := 0.2
Azic = Az[0, x] ==  - Sin[omega x]
dtAzic = Derivative[1, 0][Az][0, x] ==  -omega Cos[omega x]
Azicx = Az[t, 101] == - Sin[omega (101 + t)]
Ayic = Ay[0, x] ==  Cos[omega x]
dtAyic = Derivative[1, 0][Ay][0, x] ==  -omega Sin[omega x]
Ayicx = Ay[t, 101] == Cos[omega (101 + t)]
{Aysol, Azsol} = 
 NDSolveValue[{eqy, eqz, Ayic, dtAyic, Ayicx, Azic, dtAzic, 
   Azicx}, {Ay, Az}, {t, 0, 200}, {x, -101, 101}]

But doing the same in spherical coordinates leads to Mathematica computing until it finally just gives up, aborting the computation (seems to at least sometimes give the error that it runs out of memory). I even managed to reduce the problem to one just for A^r:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]
eqr = -r r D[Ar[t, r, phi], t, t] + r r D[Ar[t, r, phi], r, r] + 
   D[Ar[t, r, phi], phi, phi] + 4 r D[Ar[t, r, phi], r] + 
   2 Ar[t, r, phi] == 0
omega := 0.1
Aric = Ar[0, r, phi] == Sin[phi] Cos[omega (r Cos[phi])]
dtAric = Derivative[1, 0, 0][Ar][0, r, 
   phi] == -omega Sin[phi] Sin[omega (r Cos[phi])]
Arb = Ar[t, 150, phi] == 
  Sin[phi] Cos[omega (150 Cos[phi] + t)]
Arperiodic = Ar[t, r, 0] == Ar[t, r, 2 Pi]
{Arsol} = 
 NDSolveValue[{eqr, Aric, dtAric, Arb, Arperiodic}, {Ar}, {t, 0, 
   100}, {r, 0.001, 150}, {phi, 0, 2 Pi} ]

I find it quite odd that this equation should be too difficult for Mathematica? I’m currently trying to see if using a smaller interval for t and r in the solution region will help.