TransWikia.com

Why can't NDSolve(Value) solve Maxwell eq. for plane wave i.c.'s when going to spherical coordinates (it works fine in Cartesian coordinates)

Mathematica Asked by Sofie on May 10, 2021

I would like to use NDSolveValue (or NDSOlve) to solve a set of PDEs in spherical coordinates. I started with a simplified example that resembles the problem I actually want to study, but where I know the solution – just to learn how NDSolve works. Specifically, I tried solving the Maxwell equations (for the vector potential) in Minkowski space (flat space/special relativity) for a plane wave traveling in the -x direction polarized in the y-direction. The script works in Cartesian coordinates but the problem I want to solve in the end has to be done in spherical coordinates so I tried transforming the simple problem of the plane wave in Minkowski space to spherical coordinates, now needing to solve for the phi and r components of the vector potential. But this doesn’t work and I cannot figure out why. I have tried somewhat randomly (out of ignorance/inexperience, based on other questions I’ve found on stackexchange) using different ways of writing the initial conditions and setting MaxStepSize etc. but nothing I’ve tried works. Different i.c.’s and MaxStepSize etc. naturally give different errors so below I have added the scripts that seem most sensible to me (but that still aren’t working).

Below is first the script using Cartesian coordinates that works (I export the results in the end because I access mathematica using a visual ssh connection which is really slow/lags a lot so I prefer studying the results in python locally on my laptop – and I use a for-loop to arrange the data to make sure I export the results in the fashion I usually read data into python.):

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]

eqy = D[Ay[t, x], t, t]  - D[Ay[t, x], x, x] == 0

omega := 0.2

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} = NDSolveValue[{eqy, Ayic, dtAyic, Ayicx}, {Ay}, {t, 0, 200}, {x, -101,101}]


xstart := -100.123;
xend := 100.123;
Nx := 400;
dx := (xend - xstart)/Nx;
tstart := 0;
tend := 199.123;
Nt := 400;
dt := (tend - tstart)/Nt;

pts = {}
For[i = 1, i <= Nx, i++, For[j = 1, j <= Nt, j++,
  tvalue = tstart + dt (j - 1); xvalue = xstart + dx (i - 1);
  pts = Join[pts, {{tvalue, xvalue, Aysol[tvalue, xvalue]}}];] ]


Export["Minkowski_tx_cart.dat", pts]

Below is the script in spherical coordinates that does not work:

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[0, r, phi] == a[r, phi] Sin[phi] Cos[omega (r Cos[phi])]

dtAric = Derivative[1, 0, 0][Ar][0, r, phi] == -omega a[r, phi] Sin[phi] Sin[omega (r Cos[phi])]
Arb = Ar[t, Sqrt[2] 100, phi] == a[Sqrt[2] 100, phi] Sin[phi] Cos[omega (Sqrt[2] 100 Cos[phi] + t)]
Arperiodic = Ar[t, r, 0] == Ar[t, r, 2 Pi]

Aphiic = Aphi[0, r, phi] ==  a[r, phi] Cos[phi] Cos[omega (r Cos[phi] )]/r
dtAphiic = Derivative[1, 0, 0][Aphi][0, r, phi] ==  -omega a[r, phi] Cos[phi] Sin[omega (r Cos[phi] )]/r
Aphib = Aphi[t, Sqrt[2] 100, phi] == a[Sqrt[2] 100, phi] Cos[phi] Cos[omega (Sqrt[2] 100 Cos[phi] + t)]/Sqrt[2]/100
Aphiperiodic = Aphi[t, r, 0] == Aphi[t, r, 2 Pi]

{Arsol, Aphisol} = NDSolveValue[{eqr, eqphi , Aric, dtAric, Arb, Arperiodic, Aphiic,  dtAphiic, Aphib, Aphiperiodic}, {Ar, Aphi}, {t, 0, 100}, {r, 0.001, Sqrt[2] 100}, {phi, 0, 2 Pi}, MaxStepSize -> 0.01 ]

start := 0.002;
rend := 100.123;
Nr := 400;
dr := (rend - rstart)/Nr;

phistart := 0;
phiend := 6.28;
Nphi := 50;
dphi := (phiend - phistart)/Nphi;

pts = {}
For[i = 1, i <= Nr, i++, For[j = 1, j <= Nphi, j++,tvalue = 0; rvalue = rstart + dr (i - 1); phivalue = phistart + dphi (j - 1);
  pts = Join[pts, {{rvalue, phivalue, Arsol[tvalue, rvalue, phivalue], Aphisol[tvalue, rvalue, phivalue],  Derivative[0, 1, 0][Arsol][tvalue, rvalue, phivalue], Derivative[0, 0, 1][Aphisol][tvalue, rvalue, phivalue],  a[rvalue, phivalue] Sin[phivalue] Cos[omega (rvalue Cos[phivalue] + tvalue)],  a[rvalue, phivalue] Cos[phivalue] Cos[omega (rvalue Cos[phivalue] + tvalue)]/rvalue}}];]] 

Export["Minkowski_t0.dat", pts]

NDSolveValue now gives the warning

During evaluation: NDSolveValue::eerri: Warning: estimated initial error on the specified spatial grid in the direction of independent variable r exceeds prescribed error tolerance.

And the script eventually just stops running without doing any further computations. If I e.g. don’t set MaxStepSize then I get the following errors:

NDSolveValue::mxsst: Using maximum number of grid points 100 allowed by the MaxPoints or MinStepSize options for independent variable r.
NDSolveValue::mxsst: Using maximum number of grid points 100 allowed by the MaxPoints or MinStepSize options for independent variable phi.
NDSolveValue::mxsst: Using maximum number of grid points 100 allowed by the MaxPoints or MinStepSize options for independent variable r.
General::stop: Further output of NDSolveValue::mxsst will be suppressed during this calculation.

NDSolveValue eventually gives up trying to solve the PDEs and aborts the computation.

I also tried specifying the initial/boundary conditions as DirichletCondition. In this case, NDSolveValue actually gives a result but when I look at the exported data I find that the result is incorrect (i.e. not a plane wave). Even Arsol and Aphisol at t = 0 are incorrect (not a plane wave). That’s why I also export the actual initial conditions Aric and Aphiic to make sure that they are correct (which they are – as far as I can tell).

To be clear, when using DirichletCondition, I modify the script according to (note that I can’t set the r-interval to go all the way to 0 because I then get complaints about division by zero and mathematica pretty quickly gives up stating that it can’t find a solution)

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, Sqrt[2] 100}, {phi, 0, 2 Pi}, AccuracyGoal -> 24, 
  PrecisionGoal -> 10]


rstart := 0.002;
rend := 100.123;
Nr := 200;
dr := (rend - rstart)/Nr;

phistart := 0;
phiend := 6.28;
Nphi := 50;
dphi := (phiend - phistart)/Nphi;

pts = {}

For[i = 1, i <= Nr, i++, For[j = 1, j <= Nphi, j++,tvalue = 0; rvalue = rstart + dr (i - 1); phivalue = phistart + dphi (j - 1);
  pts = Join[pts, {{rvalue, phivalue, Arsol[tvalue, rvalue, phivalue],  Aphisol[tvalue, rvalue, phivalue], Derivative[0, 1, 0][Arsol][tvalue, rvalue, phivalue],  Derivative[0, 0, 1][Aphisol][tvalue, rvalue, phivalue],  a[rvalue, phivalue] Sin[phivalue] Cos[ omega (rvalue Cos[phivalue] + tvalue)], a[rvalue, phivalue] Cos[phivalue] Cos[omega (rvalue Cos[phivalue] + tvalue)]/rvalue}}];]] 


Export["Minkowski_t0.dat", pts]

Note that I tried changing the accuracy and precision by adding AccuracyGoal -> 24, PrecisionGoal -> 10 in NDSolveValue to see if the incorrect result is due to numerical errors but it doesn’t seem to be the case. I also changed the values of PrecisionGoal and AccuracyGoal to as much as 30 and 20, repsectively, but it doesn’t change the outcome noticeably. I also tried increasing the interval on which NDSolveValue solves the system of PDEs so that r goes up to 400 instead of sqrt[2] 100. This does change the result slightly but not much and it’s not clear that the solution has gotten better.

I am completely stuck and I don’t understand why NDSolveValue can’t solve the equations in spherical coordinates when it works fine in Cartesian coordinate – the problem seems fairly simple. Any suggestions on what might be going wrong would be highly appreciated.

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP