I have programmed a python script to do arithmetic and calculus using numpy arrays and their broadcasting capabilities and some clever Fourier Transform tricks implemented in a Field class and its methods.
The script works like this: if $u$ and $v$ are scalar fields in $\mathbb{R}^3$, it can do things like $u+v$ or $\nabla u$ and return the correct values. It also works with vectorial fields and I intend to generalize it to higher order tensors, infinite-dimensional fields and to implement tensor operations, along with indicial notation. I think something like this already exists but I wanted to try my hand at doing it anyway.
The next phase of my "field calculator" is to perform the dynamical evolution of Fields. It already has something like that implemented. For a certain field $u$, for instance, and its time derivative $\partial_t u$ I am computing the next step in the progression liek so:
$$ u(t+dt,\vec x) = u(t,\vec x)+dt\partial_t u $$
I was told that this is the Euler method. It seemed to work, when I tested it solving the Heat Equation (Equation \ref{eq:heat}).
$$ \label{eq:heat} \partial_t u = -\alpha \nabla^2 u \tag{1} $$
But it did not work so well when I tried to solve the Madelung Equations(Equations \ref{eq:madelung}) for the Hydrogen atom using the Euler method. The Madelung Equations involve two free fields, $R$ and $\phi$ as well as a derived field $\mu$(computed from $R$), and are considerably more complicated than the heat equation.
$$ \label{eq:madelung} \tag{2} \partial_t R + \nabla R \cdot \nabla \phi + \frac{1}{2} R \nabla^2 \phi = 0 \\ \partial_t \phi + \frac{1}{2}|\nabla \phi| + \mu = 0 $$
The results are incorrect by a mile, even using the finest grid my memory supports, even using the smallest timestep that my schedule allows. No matter how I tweak the parameters, the Madelung equations fail where the heat equation has succeeded.
If I try to solve the Schrödinger Equation (Equation \ref{eq:schrodinger}), then I can get it right... But not by using the Euler method. It turns out that there is a special way of solving the Schrödinger Equation using the Split-Operator method, which propagates the scalar field $\psi$ in an entirely different manner. When I use that method, I can solve the Schrödinger Equation for the Hydrogen atom and get a pretty good approximation to the ground state even with modest grid granularity.
$$ \label{eq:schrodinger} \tag{3} i\partial_t \psi = -\frac{1}{2} \nabla^2 \psi + V\psi $$
After these numerical experiments I became convinced that I had made a mistake.It appeared to me that the Heat Equation was blind luck, maybe because it was so simple. I was sure that propagating fields was indeed different from propagating particles, and that I should look for the appropriate methods instead of Eulering my way through it.
I had to find that general method to propagate fields, otherwise my dreams of creating a general-purpose field propagator were over. I liked the project and did not want to abandon it, so I went on and tried to retro-engineer the way that the Split-Operator method worked in the hopes of finding a general way to propagate the fields.
That journey led me down a rabbit hole towards Diffusion Quantum Monte Carlo. Although an interesting technique on its own, what attracted my attention was section "Stochastic Implementation" of the cited Wikipedia article which says, and I quote:
Now we have an equation that, as we propagate it forward in time and adjust $E_0$ appropriately, we find the ground state of any given Hamiltonian. This is still a harder problem than classical mechanics, though, because instead of propagating single positions of particles, we must propagate entire functions.
my emphasis on the last sentence. After that bit the article goes along and talks about Green's Functions as a way to do that. Although I should take wikipedia articles with a grain of salt, and I need to factor in my mathematical illiteracy, this sounded like that was the real deal. It seemed like I had been looking for the Green's Functions all along.
Naturally I went on to read about these functions. It turns out that each operator has its very own Green's function. This looked promising, since incorrectly using the Green's functions would possibly lead to wrong results I thought that maybe I got it to work accidentaly for the Heat equation, while not being so lucky with Schrödinger and Madelung's Equations.
The only problem now is how does that apply to my code?
Question: How can I compute the Green's functions appropriately so as to apply them to discrete, non-analytical functions defined over grids?