# Solving Maxwell's Equations in Mathematica

For the sake of demonstration (for my students) and practice, I wanted to numerically calculate the electric and magnetic field (3D) for a sphere of uniform density moving slightly to the right and back. Just to show/see that the changes in the E-M fields propagate at a fixed speed, as predicted by Maxwell's theory, but because of the Gauge Invariance of the potentials, they don't necessarily always follow that restriction.
I'm trying to have Mathematica solve Maxwell's equations numerically for the scalar and vector potentials for the given situation (with no success). Could you guys help me point out where I'm messing up? Thanks!
PS: I've elected to have the change distribution and it's motion being defined as an interpolating function, which gives me room to improve accuracy as necessary, for testing sake.
Here is the code as neatly as I could put it :)

    (*Max time of the simulation*)
tMax=10;
(*x coordinate shift of the charged sphere, you can run Plot[\[ScriptX][t],{t,0,tMax}] to see what it looks like*)
\[ScriptX]=
Interpolation[{{0, 0},{3, 0},{5, 2},{7, 0},{tMax, 0}},InterpolationOrder->1];
(*Charge distribution as a interpolating function in 4D (x,y,z,t), all points are zero except the inside of the sphere which is 1, you can run Manipulate[Plot[ρ[x,0,0,t],{x,-5,5}],{t,0,tMax}] to see what it looks like*)
ρ=Interpolation[Flatten[Table[{{x,y,z,t},If[(x-\[ScriptX][t])^2+y^2+z^2<=1,1,0]},{x,-5,5,1/2},{y,-5,5,1/2},{z,-5,5,1/2},{t,0,tMax,1/2}],3],InterpolationOrder->1];
(*Current Density, defined as \[ScriptJ]=ρv*)
\[ScriptJ]=ρ[x, y, z, t]{\[ScriptX]'[t],0,0};
(*scalar potential*)
ϕ=φ[x,y,z,t];
(*Vector potential*)
\[ScriptCapitalA]={\[ScriptCapitalA]x[x,y,z,t],\[ScriptCapitalA]y[x,y,z,t],\[ScriptCapitalA]z[x,y,z,t]};
(*Electric field in terms of scalar and vector potentials*)
ℰ=-\!$$\*SubscriptBox[\(∇$$,$${x,y,z}$$]ϕ\)-\!$$\*SubscriptBox[\(∂$$,$$t$$]\[ScriptCapitalA]\);
(*Magnetic field in terms of vector potential*)
ℬ=\!$$\*SubscriptBox[\(∇$$,$${x,y,z}$$]\[Cross]\[ScriptCapitalA]\);
(*Relevant maxwell's equations
The equation ∇.ℬ=0 is trivial form the definition of ℬ in terms of \[ScriptCapitalA], if you try to use it, it simply returns "True", which it is.
Same for ∇\[Cross]ℰ=-\!$$\*SubscriptBox[\(∂$$, $$t$$]ℬ\), so those can be ignored if we're dealing with potentials*)
Eqs={
\!$$\*SubscriptBox[\(∇$$,$${x,y,z}$$].ℰ\)==ρ,
\!$$\*SubscriptBox[\(∇$$,$${x,y,z}$$]\[Cross]ℬ\)==\[ScriptJ]+\!$$\*SubscriptBox[\(∂$$,$$t$$]ℰ\)
};
NDSolve[Eqs, {φ, \[ScriptCapitalA]x, \[ScriptCapitalA]y, \[ScriptCapitalA]z},{x,-5,5},{y,-5,5},{z,-5,5},{t,0,tMax}]


(it's neater on the original format, but copying and pasting it into a mathematica notebook makes it not-so-ugly)
If you run it, you'll see that we get the error: "Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable". Maybe I'm missing on some boundary conditions or initial conditions? But I cannot think of what those would be. I know that Gauge cannot be the problem, since it would not show in the solution anyway. Any help you could provide would be great!
Thanks!
PS: i know there are other tools I can use to do a numerical simulation, and I'm willing to learn them if that's what I have to do, but I'd be great if I could solve this and future problems with Mathematica :)

• Such a straightforward approach will not do anything. It is necessary to formulate a boundary value problem and use FEM or BEM. I solved a similar problem for a flat capacitor moving at a constant speed. Sep 21, 2019 at 5:35
• If you're interested, I have some notebooks where I implemented basically everything is explained in this series of lectures on FDTD. Sep 21, 2019 at 8:09
• The missing of i.c. and b.c. is certainly a problem, but it's not too hard to resolve (if you don't need things like PML). What's more serious is, NDSolve cannot handle terms like Derivative[0, 1, 0, 1][\[ScriptCapitalA]y][x, y, z, t] (mixture of temperal and spatial derivative) well, at least now. Related: mathematica.stackexchange.com/q/184281/1871 I believe one should first solve for electric field and magnetic field based on Faraday's law and Maxwell–Ampère law first using NDSolve and calculate for φ etc. afterwards. Sep 21, 2019 at 8:18
• @ Fortsaint Could you share the notebooks on FDTD with me? Sep 21, 2019 at 13:33
• @fortsaint I would like to see these as well, perhaps a question about implementing efficient fdtd methods? Sep 22, 2019 at 0:35

We can solve this problem numerically using the Lienard - Wiechert potentials discussed here and here. First we should note that 4 dimensional region $$0\le t\le 10, -5\le x\le 5, -5\le y\le 5, -5\le z\le 5$$ in a case of FDM scheme needs to be map on over 10^8 points at step $$dt=dx=dy=dz=0.1$$. Also this problem has moving boundary or/and moving source (ball). As it well known the point like moving charge generates special kind of retarded potentials dependent on retarded time. To compute retarded time we need to know coordinate and speed of point charge as, for example, in a form of interpolation function

tMax = 10;
(*x coordinate shift of the charged sphere*)
\[ScriptX] =
Interpolation[{{0, 0}, {3, 0}, {5, 1.9}, {7, 0}, {tMax, 0}},
InterpolationOrder -> 1];
vx = Interpolation[Table[{t, \[ScriptX]'[t]}, {t, 0, 10, .01}],
InterpolationOrder -> 1];
xb[t_] := Piecewise[{{\[ScriptX][t], 0 <= t <= 10}, {0, True}}];
yb[t_] := 0; zb[t_] := 0;
vb[t_] := {Piecewise[{{vx[t], 0 <= t <= 10}, {0, True}}], 0, 0};


Now we can define retarded time as follows

R1 = {x, y, z};
R2 = {xb[t], yb[t], zb[t]};
R = R1 - R2;
tr[t_?NumericQ, x_?NumericQ, y_?NumericQ, z_?NumericQ] :=
r /. FindRoot[
r == t - Sqrt[(x - xb[r])^2 + (y - yb[r])^2 + (z - zb[r])^2], {r,
t - Sqrt[(x - xb[t])^2 + (y - yb[t])^2 + (z - zb[t])^2]},
Method -> {"Newton", "StepControl" -> "TrustRegion"}];


Equations for the Lienard-Wiechert scalar and vector potential

LWS = 1/(Sqrt[R . R] - Dot[R, vb[t]]) /. t -> tr[t, x, y, z];
LWV = vb[t]/(Sqrt[R . R] - Dot[R, vb[t]]) /. t -> tr[t, x, y, z];


Using these potentials we can define interpolation function for LWS,LWV in a form

lst1 = Table[{{t, x, y, z}, LWS}, {t, 0, 10, 1}, {x, -5, 6.01,
1.01}, {y, -5, 6, 1.01}, {z, -5, 6, 1.01}]; // AbsoluteTiming
lst2 = Table[{{t, x, y, z}, LWV[[1]]}, {t, 0, 10, 1}, {x, -5, 6,
1.01}, {y, -5, 6, 1.01}, {z, -5, 6, 1.01}]; // AbsoluteTiming
phi = Interpolation[Flatten[lst1, 3]];
Ax = Interpolation[Flatten[lst2, 3]];


As we can see the step of discretization is about 1.01 only instead of 0.1. Visualization

{ContourPlot[phi[t, x, 2, 2], {t, 0, 10}, {x, -5, 5}, Contours -> 20,
ColorFunction -> "BrightBands", PlotLegends -> Automatic,
PlotRange -> All],
ContourPlot[Ax[t, x, 2, 2], {t, 0, 10}, {x, -5, 5}, Contours -> 20,
ColorFunction -> "BrightBands", PlotLegends -> Automatic,
PlotRange -> All]}


We can compute electric and magnetic field directly as

Ex = -{Derivative[1, 0, 0, 0][Ax][t, x, y, z], 0, 0} -
Derivative[0, 1, 0, 0][phi][t, x, y, z];
Ey= -
Derivative[0, 0, 1, 0][phi][t, x, y, z];
Ez= -
Derivative[0, 0, 0, 1][phi][t, x, y, z];
B = Curl[{Ax[t, x, y, z], 0, 0}, {x, y, z}];


Or we can use retarded field expressions in a form

Vt = D[R2, t, t];
Ev = ((1 - V0 . V0)*(R - V0*Norm[R]) +
Cross[R, Cross[(R - V0*Norm[R]), Vt]])/(Norm[R] - R . V0)^3 /.
t -> tr[t, x, y, z]; R0 = R /. t -> tr[t, x, y, z];
B = Cross[R0,Ev]/Norm[R0];


Expressions above computed for a point charge. To extend solution to the moving charged ball we use Monte-Carlo method with about 50000 points in a ball

Phi[t_?NumericQ, x_?NumericQ, y_?NumericQ, z_?NumericQ] :=
NIntegrate[
phi[t, x - x1, y - y1, z - z1] Boole[
Element[{x1, y1, z1}, Ball[]]], {x1, -1, 1}, {y1, -1, 1}, {z1, -1,
1}, Method -> "MonteCarlo"];
A[t_?NumericQ, x_?NumericQ, y_?NumericQ, z_?NumericQ] :=
NIntegrate[
Ax[t, x - x1, y - y1, z - z1] Boole[
Element[{x1, y1, z1}, Ball[]]], {x1, -1, 1}, {y1, -1, 1}, {z1, -1,
1}, Method -> "MonteCarlo"];


Finally we compute Phi[t,2,2,2] and compare it with 4 Pi/3 phi[t,2,2,2]

p0 = Table[{t, Phi[t, 2, 2, 2]}, {t, 0, 10, .5}];

Show[Plot[{4 Pi/3 phi[t, 2, 2, 2], Interpolation[p0]@t}, {t, 0, 10},
PlotRange -> All, PlotLegends -> {"Point", "Ball"}],
ListPlot[p0, PlotStyle -> Red]]