Help with approaching a hyperbolic PDE (QED)

Question

This is my first question here, so I'd like to apologize in advance if there's too little, or too much information and for my general lack of "etiquette".

I am in need of help with choosing an approach to solve a hyperbolic PDE of this form: $$ \frac{ \partial^2 F }{ \partial t^2 } + a_{zz}\frac{ \partial^2 F }{ \partial z^2 } + a_{t}\frac{ \partial F }{ \partial t } + a_{z} \frac{ \partial F }{ \partial z } + g(t,z)\,F = 0 $$ where $a_{zz}$, $a_z$, and $a_t$ are constant coefficients and $g(t,z)$ depends on both variables.

The way I got to this PDE is a bit lengthy, so in short, the unknown function $F$ is part of a solution: $$ \psi = e^{-\frac{i}{\hbar}(p \cdot x)} \cdot F ~~\text{($p$, $x$ - 4-vectors; (+,-,-,-) metric)} $$ to the Klein-Gordon equation for an electron ($q = e^-$, $m = m_e$) in two opposing laser pulses (on the spatial Oz axis), with minimal coupling. The 4-potential $A$ satisfies the gauge conditions: $\partial^\mu A_{\mu} = \partial_\mu A^{\mu} = 0$ (Lorenz) and $A \cdot n_1 = A \cdot n_2 = 0$, where $n_1 = (1,0,0,1)$, $n_2 = (1,0,0,-1)$ are the directions of propagation of the laser pulses.

The end goal is really to be able to compute $F$ for specific $p$ across some spatial domain of interest and evolve it in time, as I'll eventually use the K-G solutions to construct a wave-packet. Just in case, the actual coefficients are: $$ \begin{gather} a_{zz} = -c^2,~a_z = -2\frac{i}{\hbar}p_zc^2,~a_t=-2\frac{i}{\hbar}E \\ g(t,z) = -\frac{c^2}{\hbar^2} \left[ 2\frac{q}{c}\bigl(p \cdot A(t,z)\bigr) - \frac{q^2}{c^2}A(t,z)^2 \right] \end{gather} $$ I mention that $A$ is actually a combination of the two separate laser pulses $A = A_1 + A_2$, but I'm not too sure that the specific expression of $A$ is critical to solving the PDE.

I'm confident that the derivation is correct, but got stuck on that PDE as I don't have much experience with them and the one that I got, I also find quite intimidating. I believe there's no analytic/closed-form solution, so I would appreciate any advice and pointers towards some numerical methods and/or software solutions (apps, packages, any coding language) that are best suited for tackling this.

Answer 1

Finite differencing

In order to start numerically solving a PDE, you will want to approximate the continuous derivatives with finite differences, $$\frac{\partial\eta}{\partial \chi}\approx\frac{\eta(\chi_{i+1})-\eta(\chi_{i})}{\chi_{i+1}-\chi_{i}}\equiv\frac{\eta_{i+1}-\eta_i}{\Delta\chi}$$ which essentially follows from the definition of the derivative as a limit--in the right-most part of the equation, I've made the assumption that the spatial grid, $\chi$, is uniform so that you can replace the array of differences with a single scalar.

For your system, you would obtain, \begin{align} \frac{F_i^{n+1}+F_i^{n-1}-2F_i^n}{\Delta t^2}+a_{t}\frac{F_i^{n+1}-F_i^n}{\Delta t}\\ +a_{zz}\frac{F_{i+1}^{n}+F_{i-1}^{n}-2F_i^n}{\Delta z^2} &+a_z\frac{F_{i+1}^n-F_i^n}{\Delta z}+g(t^n,\,z_i)F_i^n = 0 \end{align} where the $n$ superscript indicates the time slice and the $i$ subscript indicates the spatial slice. You then solve the above system for $F_i^{n+1}$ on the left hand side and everything else on the right hand side, \begin{align} F_i^{n+1}\left(1+a_t\right)=f(\cdots) \end{align} which has a slight complexity in that you need the previous timestep, $F_i^{n-1}$, in order to compute the next timestep, $F_i^{n+1}$ (which really should be simple enough as you're storing it in computer memory).

Stability

You will also have to watch out for numerical stability, which can be tested using Von Neumann analysis in which you replace $F_i^n$ with the Fourier term, $$F_i^n=E^n\mathrm{e}^{ikx},\quad F_{i+1}^n=E^n\mathrm{e}^{ik(x+\Delta x)},\quad F_i^{n+1}=E^{n+1}\mathrm{e}^{ikx}$$ where the $i$ in the argument is the complex $i$, not the subscript index $i$. If you then reduce the system to generate the eigenvalue, you can find a way to ensure the growth of instabilities is damped (via constricting the eigenvalue to being less than 1).
It also might be the case that stability requires using central differences or backward differences instead of the forward differences I used. Or even using an implicit scheme (essentially all spatial derivatives are evaluated at $t^{n+1}$ instead of $t^n$) instead of the explicit scheme I've expressed.

Alternatives?

Since you have complex coefficients, you may want to consider splitting the function $F$ into a real and complex part, $F=\varphi+i\psi$, and joining like real/complex parts into separate PDEs to solve concurrently (cf. this SciComp.se post discussing the implementation for the Schrodinger equation).