Problem
I have an $n\times n$ grid, and each point on the grid is assigned two values: a score, and an (inverse) speed factor. There is a "turtle" moving along the grid, and it's goal is to maximize it's score. Each time it moves to a point $x$ it gets the score $s_x$ at that point, but it takes $\epsilon_x$ turns before it can move again. The turtle can only move in the four cardinal directions (north, east, south, west) and in the "center" direction (remaining where it currently is).
I want to assign a third value $p_x$, the "potential" at $x$, which helps the turtle find an optimal spot to move to. If the grid were extended to be continuous, the potential function should be continuous with a first derivative (or gradient).
Preliminaries
A natural choice for the potential function satisfies
$$\epsilon_xp_x = s_x + \frac{\sum\limits_{y\in \text{Neighbors}(x)}p_y}{4}\qquad (1)$$
(Notice that this almost exactly defines a qtable, except for the differing $\epsilon_x$ values.)
For now, assume that $\epsilon_x$ is constant. The equation happens to then be the same as a Discrete Poisson equation. This means that physically $p_x$ is the potential field in a system where each point has a $-s_x$ charge on it. There are many ways to solve the Discrete Poisson equation, but one of the fastest is to use spectral methods. See here for using the Discrete Cosine Transform (DCT) to solve it. Our grid can be rather large ($50\times 50$), so to fit within a one-second time constraint we need an $O(n^2\log n)$ algorithm.
However, $\epsilon_x$ is not constant. $\epsilon_x$ physically represents the electric permittivity at the point $x$, so if it were a one dimensional grid we could set every $\epsilon_x$ to $1$ and insert points to mimic having a higher permittivity. However, we are on a two dimensional grid. Not only that, but $\epsilon_x$ is not necessarily an integer—it could be any positive real number. My goal is to use spectral analysis to solve equation $(1)$ with an algorithm running in $O(n^2\log n)$ time.
Progress
First, rearrange the original equation:
$$\frac{p_{i+1,j}+p_{i-1,j}+p_{i,j+1}+p_{i,j01}}{4}-p_{i,j} + \gamma_{i,j}p_{i,j} = t_{i,j}\qquad (2)$$$$\frac{p_{i+1,j}+p_{i-1,j}+p_{i,j+1}+p_{i,j01}}{4}-p_{i,j} + \gamma_{i,j}p_{i,j} = -s_{i,j}\qquad (2)$$ where $\gamma_{i,j} = 1-\epsilon_{i,j}$ and $t_{i,j} = -4s_{i,j}$.
Notice that this is the finite difference equation for
$$\nabla^2 p + \gamma p = t\qquad (3)$$$$\nabla^2 p + \gamma p = -s\qquad (3)$$
Now, $$\nabla^2(fp) = f\nabla^2p + 2\nabla f\cdot\nabla p + p\nabla^2f$$ so if I can find a function $f$ that satisfies $$2\nabla f\cdot\nabla p = 0, \frac{\nabla^2 f}{f} = \gamma$$ I would be mostly done. I just solve the Poisson equation $\nabla^2(fp) = ft$$\nabla^2(fp) = -fs$ and divide by $f$. However, I don't know how to do this in $O(n^2\log n)$ time.
Alternatively, we can go straight to spectral analysis. Let $$p = \sum \hat{p}_{j,k}e^{i(jx+ky)}$$ $$t = \sum \hat{t}_{j,k}e^{i(jx+by)}$$$$s = \sum \hat{s}_{j,k}e^{i(jx+by)}$$
Plugging this into $(3)$ we get
$$\sum -\hat{p}_{j, k}(j^2+k^2)e^{i(jx+by)} + \gamma\sum \hat{p}_{j, k}e^{i(jx+ky)} = \sum \hat{t}_{j, k}e^{i(jx+ky)}$$$$\sum -\hat{p}_{j, k}(j^2+k^2)e^{i(jx+by)} + \gamma\sum \hat{p}_{j, k}e^{i(jx+ky)} = -\sum \hat{s}_{j, k}e^{i(jx+ky)}$$
Removing the summation yields $$\hat{p}_{j,k}(j^2+k^2+\gamma) = -\hat{t}_{j, k}\implies \hat{p}_{j,k} = -\frac{\hat{t}_{j,k}}{j^2+k^2+\gamma}$$$$\hat{p}_{j,k}(j^2+k^2-\gamma) = \hat{s}_{j, k}\implies \hat{p}_{j,k} = \frac{\hat{s}_{j,k}}{j^2+k^2-\gamma}$$
But $\hat{p}_{j,k}, \hat{t}_{j,k}$$\hat{p}_{j,k}, \hat{s}_{j,k}$ are constants while $\gamma$ is not, so there is no solution.