BRS: Solving Killing/Geodesic Eqns with Maple's casesplit

[Chris Hillman]

For those who use GRTensorII, here is a simple and effective method for

• determining the Lie algebra of Killing vector fields
• using the result to find first integrals of the geodesic equations via Noether's theorem on symmetries

The Lie algebra of Killing vector fields "exponentiates" to give the Lie group of self-isometries on our spacetime; if this acts transitively on a family of three-dimensional submanifolds, we can obtain a complete list of first integrals.

This method requires no particular knowledge of GRTensorII.

I'll explain using a simple example, the Melvin electrovacuum, which models a cylindrically symmetric magnetostatic field, in which the EM field energy is the sole source of the gravitational field. Here is an (extensively commented) GRTensorII file you can use to define the static coframe in this spacetime:

Ndim_ := 4:
x1_ := t:
x2_ := z:
x3_ := r:
x4_ := phi:
eta11_ := -1:
eta22_ :=  1:
eta33_ :=  1:
eta44_ :=  1:
bd11_ := -(1+q^2*r^2):
bd22_ :=  (1+q^2*r^2):
bd33_ :=  (1+q^2*r^2):
bd44_ :=  r/(1+q^2*r^2):
Info_ := `Melvin non-null electrovacuum (1964); cyl chart; static obsvrs`:
# Chart -infty < t,z < infty, 0 < r < infty, -Pi < u < Pi
# (Discovered by Bonnor 1954 and Witten 1962 and then Melvin 1964)
#
# Four dimensional Lie algebra of Killing vector fields:
#	@_t			time translation
#	@_z			spatial translation in z direction
#	z @_t + t @_z		boost in z direction
#	@_phi			axial rotation
# So this is boost/rotation symmetric as well as cyl sym static.
# Geodesic equations:
# From Lagrangian we obtain quadratic invariant
#	epsilon = (1+q^2*r^2)^2*(-t1^2+z1^2+r1^2)+phi1^2*r^2/(1+q^2*r^2)^2
# From Killing vectors we obtain four linear quadratic invariants:
#	B = t1*(1+q^2*r^2)^2
#	P = z1*(1+q^2*r^2)^2
#	K = (z*t1+t*z1)*(1+q^2*r^2)^2
#	J = phi1*r^2/(1+q^2*r^2)^2
# Solving { epsilon=0,B,P,J} for {t1,z1,r1,phi1} we obtain
#   	t1 = B/(1+q^2*r^2)^2
#   	z1 = P/(1+q^2*r^2)^2
#	r1 = sqrt((B^2-P^2)*r^2-(1+q^2*r^2)^4*J^2)/(1+q^2*r^2)^2/r
#   	phi1  = J/r^2*(1+q^2*r^2)^2
# Helical null congruence
#	J/sqrt(B^2-P^2) = r/(1+q^2*r^2)^2
#	1/(1+q*r^2) (B e_1 + P e_2) + sqrt(B^2-P^2)/r e_4
#
# Physical experience of static observers:
# Acceleration vector
#   	-2*q^2*r/(1+q^2*r^2)^2 e_3
# Expansion and vorticity tensor vanish
# Vector Potential
#   	-q/2*(1+q^2*r^2) @_phi
# Axial regular magnetic field concentrated near z = 0
#   	q/(1+q^2*r^2)^2 e_2
# Einstein tensor
#   	G^(ab) = 4*q^2/(1+q^2*r^2) diag(1,-1,1,1)
# Principal gravitational field invariants
#   	RR = 64*q^4*(3*r^4*q^4-6*r^2*q^2+5)/(1+r^2*q^2)^8
#   	Rstar = 0
#   	diRiem = 256*r^2*q^8*(133-126*r^2*q^2+45*r^4*q^4)/(1+r^2*q^2)^12
#   	CC = 192*q^4*(r*q-1)^2*(r*q+1)^2/(1+r^2*q^2)^8
#   	Cstar = 0
#   	diWeyl = 768*q^8*r^2*(15*r^4*q^4-42*r^2*q^2+31)/(1+r^2*q^2)^12
# Curvature everywhere finite and concentrated near r = 0
# Tidal tensor or electro-Riemann shows less axial tension than electro-Weyl tensor
#   	E_(ab) = EW_(ab) + 4*q^2/(1+q^2*r^2)^4 diag(1,0,0)
# This describes tidal accelerations of clouds of electrically neutral test particles.
# Weyl tensor is "pure electric"
#   	EW_(ab) = 2*q^2*(1-r^2*q^2)/(1+q^2*r^2)^4 diag(-2,1,1)
# Orthogonal hyperslices have three-dimensional Riemann tensor
#   	rv^2_(323) = rv^2_(424) = 2*q^2*(q^2*r^2-1)/(q*2*r^2+1)^4
#   	rv^3_(434) = 4*q^2*(q^2*r^2-2)/(q^2*r^2+1)^4
#
# For cylindrically symmetric functions, the Laplace-Beltrami equation
# becomes the usual wave equation
#   	h_(tt) = h_(rr) + (h_r)/r 
# A simple exact solution is
#   	h = cos(a*t)*BesselJ(0,a*r)
# The corresponding vector field is
#   	sin(a*t)*BesselJ(0,a*r)/(1+q^2*r^2) e_1 
#  	-a*cos(a*t)*BesselJ(1,a*r)/(1+q^2*r^2) e_3
#
# References:
#  Bonnor, prsl A67 (1954): 225.
#  Melvin, Phys. Lett. 8 (1964): 65.

Now fire up Maple and in your session type the following commands:

qload(Melvin_nnevac_mg_stataxs_reg_cyl); 
grcalcalter(x(up),basisv(dn),basisv(up),g(dn,dn),g(up,up),1):grdisplay(_);
grcalcalter(LieD[kv0,g(dn,dn)],1):grdisplay(_);
keqs := [seq(seq(grcomponent(LieD[kv0,g(dn,dn)], [grcomponent(x(up),[j]),grcomponent(x(up),[k])]),j=1..4),k=1..4)]:
# Enter 
#	T(t,z,r,phi)
#	Z(t,z,r,phi)
#	R(t,z,r,phi)
#	Phi(t,z,r,phi)
op(%)[1];
pdsolve(%,{T,Z,R,Phi});
subs(%,[T,Z,R,Phi](t,z,r,phi));
sol := subs(_C1=c[1],_C2=c[2],_C3=c[3],_C4=c[4],%);
# Output 
#	[c[1]*z+c[4], c[1]*t+c[2], 0, c[3]]
subs(c[1]=1,c[2]=0,c[3]=0,c[4]=0,sol);
subs(c[1]=0,c[2]=1,c[3]=0,c[4]=0,sol);
subs(c[1]=0,c[2]=0,c[3]=1,c[4]=0,sol);
subs(c[1]=0,c[2]=0,c[3]=0,c[4]=1,sol);
# Output shows four dimensional Lie algebra of Killing vectors
#	z @_t + t @_z
#	@_t, @_z, @_phi
# Transitive on cylinders r=r0
#
multiply(matrix([[1,0,0,0]]), matrix( [[-(1+q^2*r^2)^2,0,0,0],[0,(1+q^2*r^2)^2,0,0],[0,0,(1+q^2*r^2)^2,0],[0,0,0,r^2/(1+q^2*r^2)^2]]), matrix( [[t1],[z1],[r1],[phi1]]))[1,1] = -EE; t1=solve(%,t1)
# Output shows t1 = EE/(1+q^2*r^2)^2
multiply(matrix([[0,1,0,0]]), matrix( [[-(1+q^2*r^2)^2,0,0,0],[0,(1+q^2*r^2)^2,0,0],[0,0,(1+q^2*r^2)^2,0],[0,0,0,r^2/(1+q^2*r^2)^2]]), matrix( [[t1],[z1],[r1],[phi1]]))[1,1] = A; z1=solve(%,z1);
# Output shows z1 = A/(1+q^2*r^2)^2
multiply(matrix([[0,0,0,1]]), matrix( [[-(1+q^2*r^2)^2,0,0,0],[0,(1+q^2*r^2)^2,0,0],[0,0,(1+q^2*r^2)^2,0],[0,0,0,r^2/(1+q^2*r^2)^2]]), matrix( [[t1],[z1],[r1],[phi1]]))[1,1] = L; phi1=solve(%,phi1); lprint(%);
# Output shows phi1 = L*(1+q^2*r^2)^2/r^2
multiply(matrix([[z,t,0,0]]), matrix( [[-(1+q^2*r^2)^2,0,0,0],[0,(1+q^2*r^2)^2,0,0],[0,0,(1+q^2*r^2)^2,0],[0,0,0,r^2/(1+q^2*r^2)^2]]), matrix( [[t1],[z1],[r1],[phi1]]))[1,1] = B; factor(%); -t*z1+z*t1=solve(%,-z*t1+t*z1): simplify(%,factor);
# Output shows	z*t1-t*z1 = B/(1+q^2*r^2)^2
[t1 = EE/(1+q^2*r^2)^2, z1 = A/(1+q^2*r^2)^2,  phi1 = L*(1+q^2*r^2)^2/r^2];
grcomponent(g(dn,dn),[t,t])*t1^2 + grcomponent(g(dn,dn),[z,z])*z1^2 + grcomponent(g(dn,dn),[r,r])*r1^2 + grcomponent(g(dn,dn),[phi,phi])*phi1^2;
subs(%%,%): simplify(%,factor): collect(%,EE): collect(%,A): collect(%,L);
lprint(%);
# Output shows
#	r1^2 = (EE^2-A^2)/(1+q^2*r^2)^3 - L^2*(1+q^2*r^2)/r^2

Here is what we did:

• We computed the Killing equations for a generic vector field
  $$T \, \partial_t + Z \, \partial_z + R \, \partial_r + \Phi \, \partial_\phi$$
  You can call the undetermined functions $T,Z,R,\Phi$ anything you like, of course.
• We put these PDEs in a list.
• The key step is feeding this list of first order PDEs to Maple's powerful "casesplit" commend, which "triangularizes" them. (The form of the Killing equations guarantees than in principle, this can be done.)
• We use Maple's "pdesolve" command solves these system.
• We extract a list of four vector fields which generate the Lie algebra,
  • the boost $z \, \partial_z + t \, \partial_z$
  • time translation $\partial_t$; this turns out to be a vorticity free timelike vector field, so the constant t slices form a family of spatial hyperslices everywhere orthogonal to the world lines of the static observers,
  • translation along axis of cylindrical symmetry $\partial_z$, a spacelike Killing vector field,
  • rotation about the axis of cylindrical symmetry $\partial_\phi$, a spacelike cyclic Killing vector field
  which confirms that the solution is static and cylindrically symmetric, and what it is more, it is a boost-rotation symmetric solution (see the review paper by Pravda and Pravdova in the "Useful Links" thread).
• Courtesy of Noether's theorem on variational symmetries of a system of PDEs determined by a Lagrangian, we immediately obtain several first integrals (more than we need, actually). You can use any notation you like for the first derivatives of a geodesic wrt an arc length parameter you like, but I find t1, z1, r1, phi1 particularly convenient for working with Lagrangians (when at various moments we need to consider derivatives as independent variables, or to work purely algebraically, as here). I used EE for energy in the Maple session so that there is no possibility of confusion with the defined constant $\exp(1)$. In particular, we obtain first integrals for $\dot{t}, \; \dot{z}, \; \dot{\phi}$, but we still need a first integral for $\dot{r}$.
• We plug these into the line element $0 = ds^2$ to obtain the first integral for $\dot{r}$, in the case of null geodesics.

The grand result is that we have the following list of first integrals for null geodesics:
$$\begin{array}{rcl} \dot{t} & = & \frac{E}{(1+q^2 \, r^2)^2} \\ \dot{z} & = & \frac{A}{(1+q^2 \, r^2)^2} \\ \dot{r} & = & \pm \sqrt{\frac{E^2-A^2}{(1+q^2 \, r^2)^3} - \frac{L^2 \, (1+q^2 \, r^2)}{r^2}} \\ \dot{\phi} & = & \frac{L \, (1+q^2 \,r^2)^2}{r^2} \end{array}$$
where $E, \, A, \, L$ are the invariants of null geodesic motion, and where the choice of sign determines an ingoing or outgoing motion.

We can immediately write down a null geodesic vector field whose integral curves form a null geodesic congruence:
$$\vec{k} = \frac{1}{(1+q^2 \, r^2)^2} \; \left( E \, \partial_t + A \, \partial_z \right) \pm \sqrt{ \left( \frac{E^2-A^2}{(1+q^2 \, r^2)^3} - \frac{L^2 \, (1+q^2 \, r^2)}{r^2} \right)} \; \partial_r + \frac{L \, (1+q^2 \,r^2)^2}{r^2} \; \partial_\phi$$
Far from the axis of cylindrical symmetry, this reduces to
$$\vec{k}_{\rm far} = \left( E \, \partial_t + A \, \partial_z \right) \pm \sqrt{E^2-A^2} \; \partial_r + \frac{L}{r^2} \; \partial_\phi$$
which corresponds to a plane wave in cylindrical coordinates in flat spacetime, with a constant wave vector pointing in an arbitrary direction (determined by E,A,L). This also shows that for our null geodesics, we can loosely interpret E as "the energy of a photon measured at spatial infinity", and likewise L as "linear momentum parallel to the axis of cylindrical symmetry" and "angular momentum about the axis of cylindrical symmetry". This "photon" doesn't remain "at spatial infinity", but because E, A, L are constants of geodesic motion, these parameters still have meaning as it approaches its point of closest approach to the axis of cylindrical symmetry, even though we interpret them in terms of measurements made when the "photon" was "at spatial infinity".

By using $-1 = ds^2$ instead, you can find first integrals for an arbitrary timelike geodesic in the Melvin electrovacuum. Using this you can study such interesting phenomena as the light bending and geodesic precession due to the gravitational attraction of the magnetic field, which is strongly concentrated near the symmetry axis.

In spacetimes with smaller Lie groups of self-isometries, it may not be possibly to obtain in such trivial fashion a complete list of first integrals for null, timelike, or spacelike geodesics, but one almost always obtains some valuable information. The case of the Kerr vacuum is remarkable because, although we have only a two dimensional Lie group of self-isometries, and thus we can expect only three first integrals for (null, timelike or spacelike) geodesicsl. In this case, it turns out that there also exists a more subtle symmetry, given by a nontrivial Killing tensor field. Brandon Carter was able to exploit this to solve the geodesic equations for the Kerr vacuum. Trivial Killing tensor fields are common and useless; nontrivial Killing tensor fields seem to be quite rare, in fact, IIRC no example is known other than the Kerr vacuum!

At this point, if you want to examine the Killing equations with which we started, you can try this

declare(T(t,z,r,phi), Z(t,z,r,phi), R(t,z,r,phi), Phi(t,z,r,phi);
for guy in keqs do guy od;

The first command instructs Maple to abbreviate partials by subscripts. The second prints out each equation, one line per equation. From this you can see that as I said, this is a first order system of PDEs. Compare with the "triangularized" version output by casesplit!

I should say that casesplit uses nontrivial Groebner-basis like methods in differential rings, which we can think of as a sophisticated analog of using Gaussian reduction to "triangularize" a system of linear equations. In both cases, the point is that we can systematically solve a triangularized system by solving for the variable for which we have obtained equations involving only that variable, then plugging the result into the equations involving that variable and one other, solving, and so forth until we are done.

There is a general principle in computational algebra which says that the fastest methods tend to be non-deterministic. Thus, casesplit may sometimes present results in a different order, so the directions I gave above must be interpreted with this in mind. Also, casesplit is not perfect and sometimes the above method will fail to find all the solutions it should, so be on the lookout for such problems.

 

[Chris Hillman]

BRS: Solving the EFE with Maple's casesplit

An even more obvious application of casesplit is solving the EFE!

For example, if you have written out a frame field expressed in terms of two undetermined functions p,q of some of the coordinates in the chart you are using, say t,r, where you suspect p plays the role of a master function and q a secondary metric function, you can try to find vacuum solutions

grcalcalter(G(bup,bup),1):grdisplay(_);
[seq(seq(grarray(G(bup,bup))[j,k],j=1..4),k=1..4)];
casesplit(%,[[q],[p]]);

This method is easily extended to non-null electrovacuums; for example, here are some simple coframe Ansatze we can try for a static cylindrically symmetric spacetime written wrt various charts:

• a cylindrical chart:
  $$\begin{array}{rcl} \sigma^1 & = & -u(r) \, dt \\ \sigma^2 & = & v(r) \, dz \\ \sigma^3 & = & v(r) \, dr \\ \sigma^4 & = & r \, d\phi \end{array}$$
• a cylindrical, spatially isotropic type chart:
  $$\begin{array}{rcl} \sigma^1 & = & -u(r) \, dt \\ \sigma^2 & = & v(r) \, dz \\ \sigma^3 & = & v(r) \, dr \\ \sigma^4 & = & v(r) \, r \, d\phi \end{array}$$
• yet another possibility
  $$\begin{array}{rcl} \sigma^1 & = & -u(r) \, dt \\ \sigma^2 & = & u(r) \, dz \\ \sigma^3 & = & u(r) \, dr \\ \sigma^4 & = & v(r) \, r \, d\phi \end{array}$$

If we seek a non-null electrovacuum solution featuring a magnetostatic field parallel to the axis of cylindrical symmetry, we know that wrt the static frame the Einstein tensor should have the simple form
$$G^{\hat{a}\hat{b}} = h(r) \; \operatorname{diag} \, (1,-1,1,1)$$
so we can compute the Einstein tensor from our coframe and write down a simple system of equations,

[ grcomponent(G(bup,bup),[1,1])+grcomponent(G(bup,bup),[2,2]), grcomponent(G(bup,bup),[1,1])-grcomponent(G(bup,bup),[3,3]), grcomponent(G(bup,bup),[1,1])-grcomponent(G(bup,bup),[4,4]) ];

which we can then simplify using casesplit and then we can solve the simplified system. (In this example, using casesplit is overkill, but in more complicated problems it can be very useful.) To be sure, this only produces a metric tensor whose Einstein tensor has the correct form for a magnetostatic field parallel to the axis of cylindrical symmetry; we still have to use an EM potential one-form Ansatz to solve the curved-spacetime Maxwell equations on this spacetime, and finally to match the EM contribution to the stress-energy tensor to the Einstein tensor, in order to produce the desired magnetic field needed to complete the solution:
$$\begin{array}{rcl} ds^2 & = & (1+q^2 \, r^2)^2 \; \left( -dt^2 + dz^2 + dr^2 \right) + \frac{r^2 \, d\phi^2}{(1+q^2 \,r^2)^2}, \\ \vec{A} & = & -q/2 \; (1+q^2 \,r^2) \; \partial_\phi, \\ \vec{B} & = & \frac{q}{(1+q^2 \, r^2)^2)} \; \vec{e}_2, \\ && -\infty < t, \, z < \infty, \; 0 < r < \infty, \; -\pi < \phi < \pi \end{array}$$

But where "casesplit" really shines is in more powerful examples, where our goal is to reduce the EFE to a system of PDEs whose solutions are readily studied by analytic means. A superb example is provided by the Weyl vacuums (see the BRS thread on these), which are defined by the system
$$u_{zz} + u_{rr} + \frac{u_r}{r} = 0, \; \; v_z = 2 \, r \, u_z \, u_r, \; \; v_r = r \, ( u_r^2-u_z^2)$$
which we can find by plugging in Weyl's metric Ansatz and reducing the vacuum EFE using casesplit (make sure to ask casesplit to try to express v in terms of u!).

As I remarked in the Weyl vacuums BRS thread, it is perhaps not clear at a glance that the second and third equations are mutually consistent. To find out

[ diff(v(z,r),z) = 2*r*diff(u(z,r),z)*diff(u(z,r),r), diff(v(z,r),r) = r*(diff(u(z,r),r)^2-diff(u(z,r),z)^2) ]
casesplit(%, [[v],[u]]);
op(%[1])[1];
expand(%[3],`+`);

This shows that the condition on u in order for these two PDEs to be consistent is precisely the first equation in the Weyl system. So the Weyl system is self consistent and it says: choose any harmonic axisymmetric function u and then find v by quadrature from the second and third equations.

 

[Chris Hillman]

BRS: Solving Killing/Geodesic Eqns... A Killing tensor reference

For those interested in Killing tensors, a new eprint
http://arxiv.org/abs/1003.0019
David Garfinkle and E.N. Glass,
Killing Tensors and Symmetries

A new method is presented for finding Killing tensors in spacetimes with symmetries. The method is used to find all the Killing tensors of Melvin's magnetic universe and the Schwarzschild vacuum. We show that they are all trivial. The method requires less computation than solving the full Killing tensor equations directly, and it can be used even when the spacetime is not algebraically special.

Both authors are well known and very experienced in this field.

The background you need to understand this eprint may be found in Stephani and MacCallum, Differential Equations: their solution using symmetries, Cambridge University Press, 1989.

It must seem like a curious coincidence that the authors choose the same example (Melvin electrovacuum), but it's not so strange if you know that--- IIRC--- Garfinkle is a former student of Melvin and has generalized the Melvin electrovacuum to include an (exact!) gravitational wave traveling along the axis, with plus polarization (wrt a suitable frame), with arbitrary time variation in the amplitude, but decaying in a specific way with distance from r=0, and blowing up near r=0.

Here is a ctensor (Maxima) file you can run in batch mode under wxmaxima:

/* 
Garfinkle-Melvin nnevac; double null chart

Einstein tensor
	G^(ab) = -f q^2/(1+q^2 r^2)^4 diag(1,-1,1,1)
shows this is a nonnull electrovacuum.

Weyl spinor has components (wrt NP tetrad corresponding to given frame)
        Psi_0 = f (1-q^2*r^2)/r^2/(1+q^2*r^2)^3         plus pol gw
        Psi_2 = -2*q*(1-q^2*r^2)/(1+q^2*r^2)^4          Coulomb
Notice the wave has arbitrary amplitude profile (time variation) and is strongest near r=0.
*/
load(ctensor);
cframe_flag: true;
ratchristof: true;
ctrgsimp: true;
/* define the dimension */
dim: 4;
/* list the coordinates */
ct_coords: [u,v,r,phi];
/* dependent and independent variables */
depends(f,v);
constant(q);
/* define background metric */
lfg: ident(4);
lfg[1,1]: -1;
/* define the coframe */
fri: zeromatrix(4,4);
fri[1,1]: -(1+q^2*r^2)/sqrt(2);
fri[1,2]:  (1+q^2*r^2)/sqrt(2)-f*log(r)*(1+q^2*r^2)/sqrt(2);
fri[2,1]:  (1+q^2*r^2)/sqrt(2);
fri[2,2]:  (1+q^2*r^2)/sqrt(2)+f*log(r)*(1+q^2*r^2)/sqrt(2);
fri[3,3]:  (1+q^2*r^2);
fri[4,4]:  r/(1+q^2*r^2);
/* setup the spacetime definition */
cmetric();
/* display matrix whose rows give coframe covectors */
fri;
/* compute a matrix whose rows give frame vectors */
fr;
/* metric tensor g_(ab) */
lg;
/* compute g^(ab) */
ug: invert(lg);
christof(false);
/* Compute fully covariant Riemann components R_(mijk) = riem[i,k,j,m] */
lriemann(true);
/* Compute R^(mijk) */
uriemann(false);
/* Compute Ricci componets R_(jk) */
ricci(true);
/* Compute trace of Ricci tensor */
tracer;
/* Compute R^(jk) */
uricci(false);
/* Compute and display MIXED Einstein tensor G^a_b */
/* For (-1,1,1,1) sig Flip sign of top row to get G^(ab) */
einstein(false);
cdisplay(ein);
/* WARNING! leinstein(false) only works for metric basis! */
/* Compute Kretschmann scalar */
rinvariant();
expand(factor(%));
/* Compute NP tetrad and Weyl scalars */
nptetrad(true);
weyl(false);
psi(true);
factor(psi[0]);
factor(psi[2]);
petrov();

Note that this gives the solution in terms of a double null chart in which the line element becomes
$$\begin{array}{rcl} ds^2 & = & 2 (1+ q^2 \, r^2) \; du \, dv + (1+q^2 \, r^2)^2 \; dr^2 + \frac{r^2 \; d\phi^2}{(1+q^2 \, r^2)^2} + 2 \; f \; (1+q^2 \, r^2)^2 \; \log(r) \; dv^2 \\ && -\infty < u, \; v < \infty, \; 0 < r < \infty, \; -\pi < \phi < \pi \end{array}$$
where f is an arbitrary smooth function of v only. Here, the real constant q controls the amplitude of the magnetostatic field, while f controls the amplitude profile of the wave. When f is identically zero, this reduces to the Melvin electrovacuum, and when q is zero, it reduces to a certain EK 4 plus polarized axisymmetric vacuum pp-wave. (EK4 refers to the classification by Ehlers and Kundt of all vacuum pp-waves; these are exact Petrov type N gravitational wave solutions generalizing the well known plane waves of Baldwin and Jeffery.)

If you have installed Maxima, you can at least use this to verify that the Einstein tensor has the correct form for a nonnull electrovacuum with axial EM field
$$G^{ab} = \frac{4 \, q^2}{(1+q^2 \, r^2)^4} \; \operatorname{diag} (1,-1,1,1)$$
(which is independent of f) and that the electroriemann and magnetoriemann tensors do exhibit plus polarized radiation; the trick is separate the terms containing f, like so:
$$\begin{array}{rcl} E\left[\vec{e}_1\right]_{22} & = & \frac{4 \, q^4 \, r^2}{(1+q^2 \, r^2)^4} \\ E\left[\vec{e}_1\right]_{33} & = & \frac{2 \, q^2 \; (1- q^2 \, r^2)}{(1+q^2 \, r^2)^4} + \frac{f \; (1-q^2 \, r^2)/2/r^2}{ (1+q^2 \, r^2)^3} \\ E\left[\vec{e}_1\right]_{44} & = & \frac{2 \, q^2 \; (1- q^2 \, r^2)}{(1+q^2 \, r^2)^4} - \frac{f \; (1-q^2 \, r^2)/2/r^2}{ (1+q^2 \, r^2)^3} \\ B\left[\vec{e}_1\right]_{ab} & = & \frac{-f \; (1-q^2 \, r^2)/2/r^2}{(1+q^2 \, r^2)^3} \; \left[ \begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{array} \right] \end{array}$$
Here
$${E\left[\vec{U}\right]}_{ab} = R_{ambn} \, U^m \, U^n, \; \; {B\left[\vec{U}\right]}_{ab} = {^{\ast}R}_{ambn} \, U^m \, U^n$$
The same information is contained more compactly in the Weyl spinor computed using the NP tetrad corresponding to the given coframe:
$$\Psi_0 = \frac{f \; (1-q^2 \, r^2)/r^2 } {(1+q^2 \, r^2)^3}, \; \Psi_2 = \frac{-2 \, q^2 \; (1-q^2 \, r^2)}{(1+ q^2 \, r^2)^4}$$
(Warning! Ctensor will use a -+++ signature here, consistent with our coframe definition, but the literature and GRTensorII insist upon using a +--- signature for any computations in NP formalism, so there is a sign discrepancy between the spinor components reported by Ctensor viz GRTensorII.) Here $\Psi_2$ gives the Coulomb term, the Petrov type D component, the static part of the gravitational field, which is generated by the mass-energy of the magnetostatic field, which is strongly concentrated near r=0. And $\Psi_0$ gives the Petrov type N component (tranverse gravitational radiation), the amplitude of the dynamic part of the field due to the gravitational wave, which has wave vector $\partial_v = 1/\sqrt{2} \; \left( \partial_t + \partial_z \right)$ and is traveling along the axis r=0 (as shown by the fact that the wavelike components in the electroriemann and magnetoriemann tensors are orthogonal to $\vec{e}_2$). Since this spinor component is real, the radiation is plus polarized (also seen in the electroriemann and magnetoriemann tensors) wrt our frame/tetrad. Since both a type N and type D component are present, the Weyl tensor is Petrov type II.

Compare the electroriemann tensor (aka electrogravitic tensor, aka tidal tensor) with the electroweyl tensor (confusingly, also sometimes called "tidal tensor", an inappropriate name in a nonvacuum)
$${\tilde{E}\left[\vec{U}\right]}_{ab} = C_{ambn} \, U^m \, U^n$$
For clarity, momentarily put f = 0, i.e. restrict to the case of the Melvin electrovacuum (no gravitational wave). Then the axial tidal acceleration on a small cloud of electrically neutral test particles is determined by
$${E\left[\vec{U}\right]}_{22} = \frac{ 4 q^4 \, r^2}{(1+q^2 \, r^2)^4}$$
not
$${\tilde{E}\left[\vec{U}\right]}_{22} = \frac{ -4 q^2 \; (1- q^2 \, r^2)}{(1+q^2 \, r^2)^4}$$
The discrepancy is caused by tidal effects due to the "immediate presence" of EM field energy, which the Weyl tensor doesn't take into account. This discrepancy is a bit awkward, since otherwise considerations from representation theory would impel us to first decompose the Riemann tensor into the Ricci scalar, the traceless Ricci tensor, and (completely traceless part) the Weyl tensor, and only then to decompose the Weyl tensor wrt a timelike (possibly nongeodesic) unit vector field $\vec{U}$ into electroweyl and magnetoweyl (second rank tensors living in a spatial hyperplane element orthogonal to $\vec{U}$; if $\vec{U}$ has nonzero vorticity, these elements cannot be "knit" into orthogonal hyperslices).

Here is a GRTensorII (Maple) definition of the Garfinkle-Melvin electrovacuum:

Ndim_ :=  4:
x1_ :=  u:
x2_ :=  v:
x3_ :=  r:
x4_ :=  phi:
eta11_ := -1:
eta22_ :=  1:
eta33_ :=  1:
eta44_ :=  1:
bd11_ := -(1+q^2*r^2)/sqrt(2):
bd12_ :=  (1+q^2*r^2)/sqrt(2)-f(v)*log(r)*(1+q^2*r^2)/sqrt(2):
bd21_ :=  (1+q^2*r^2)/sqrt(2):
bd22_ :=  (1+q^2*r^2)/sqrt(2)+f(v)*log(r)*(1+q^2*r^2)/sqrt(2):
bd33_ :=  (1+q^2*r^2):
bd44_ :=  r/(1+q^2*r^2):
Info_ := `Garfinkle-Melvin electrovacuum (1992); double null chart`:
# Chart covers -infinity < u,v,z < infinity, -Pi < phi < Pi
# Transformation from cylindrical chart is
#   u = (t+z)/sqrt(2),  v = (t-z)/sqrt(2)
# N.B.: here f(v) is arbitrary bounded smooth function (profile of
# gravitational traveling wave).
#
# Vector potential q*(1+q^2*r^2) @_phi
# Purely axial magnetic field 
#                      2*q/(1+q^2*r^2)^2 e_2
# EM stress energy tensor matches Einstein tensor
#   4*q^2/(1+q^2*r^2)^4 diag(1,-1,1,1)
# Principle Lorentz invariants
#   FF = 8*q^2/(1+q^2*r^2)^4
#   FstarF = 0
# In addition to this Melvin background we have a polarized 
# gravitational traveling wave moving in the direction e_2.
# This wave does not disturb the axial magnetic field.
# The magnetoriemann tensor is
#   B_(34) = -f(v)*(1-q^2*r^2)/2/r^2/(1+q^2*r^2)^3
# The electroriemann tensor is 
#   E_(22) =  4*q^2*r^2/(1+q^2*r^2)^4
#   E_(33) = -2*q^2*(1-q^2*r^2)/(1+q^2*r^2)^4 
#            + f(v)*q^2*(1-q^2*r^2)/2/r^2/(1+q^2*r^2)^3
#   E_(44) = -2*q^2*(1-q^2*r^2)/(1+q^2*r^2)^3
#            - f(v)*q^2*(1-q^2*r^2)/2/r^2/(1+q^2*r^2)^3
# where new term f(v) shows plus polarized gravitational radiation.  
# The wave component vanishes as r -> infinity but blows up as r -> 0.
#
# Reference:
#  David Garfinkle & M. A. Melvin, Phys. Rev. D 45 (1992): 1188

If you have Maple and GRTensorII installed, you can use this to verify that the magnetic field is
$$\vec{B} = \frac{2q}{(1+q^2 \, r^2)^2} \; \vec{e}_2$$
The energy density of the magnetostatic field is
$$\varepsilon = \frac{q^2}{2 \, \pi \; (1+q^2 \, r^2)^4}$$

A curious feature of this solution is that the wave does not disturb the magnetostatic field. Physically speaking, what could cause such a field? The source of magnetostatic field is presumably some steady electrical currents "at $r=0, \; z = \pm \infty$", and the wave must result from the motion of uncharged matter "at $r=0, \; z = -\infty$", in just the right way to cause a wave with just the right characteristics to leave the magnetostatic field undisturbed. So its a bit artificial, while certainly very interesting as a simple explicit example of an exact solution in which "something actually happens".

It's quite hard to find exact solutions in which "an actual physical interaction occurs", and this one probably doesn't count. Someday I hope to get around to trying to explain enough about colliding plane wave (CPW) solutions that interested SA/Ms will be able to model interesting questions like possible "generic black hole interiors" (as treated in gtr).