I'm currently trying to simulate the dynamics of a particle with position $\mathbf{x}_a(t)$ which moves in response to a field $c$. The field is in turn influenced by the particle. The system is governed by the following equations:-
\begin{align} \frac{\partial c}{\partial t} &= D_c \nabla^2 c -\beta c + h \delta(\mathbf{x} - \mathbf{x}_a(t) ) \tag{1} \\[1em] \frac{d\phi}{dt} &= -(\kappa/\gamma_R)\begin{pmatrix}-\sin \phi \\ \cos\phi\end{pmatrix}\cdot\nabla c \tag{2} \\[1em] \frac{d\mathbf{x}_a(t)}{dt} &= v\begin{pmatrix} \cos\phi \\ \sin\phi \end{pmatrix} \tag{3} \end{align}
where $\mathbf{x}_a(t)$ is the position of some particle and $D_c, \beta, h, \kappa, \gamma_R$ and $v$ are all positive constants. I'm solving (1) at each time step using a finite difference scheme. Time-marching is handled using an operator-splitting method similar to that used here: https://www.dealii.org/current/doxygen/deal.II/step_19.html.
The results I get from my program don't agree with published data. There are a few tricky things that complicate the problem that I suspect are responsible -
- Implementation of the Dirac delta in (1). The gradients of $c$ in (2) depend very sensitively on how the delta is implemented and on how the time-stepping is performed. To avoid a singularity in $c$, at timestep $n+1$ I evaluate gradients of $c$ at the $\mathbf{x}_a^{n}$ - that is to say, at the previous position of the particle.
When I do this, the gradients of $c$ I compute are very small ($\mathcal{O}(10^{-11}$) or smaller) and $\phi$ never changes. If I try to remedy this by artificially making the delta spike taller, gradients become too large and the observed behavior of the particle looks like a random walk (also not correct).
- Another potential problem is that I can only evaluate the delta at grid points, meaning that in practice I simply have to find the grid point closest to $\mathbf{x}_a(t)$, evaluate the delta there and gradients of $c$ as well. I suspect this could also contribute to the apparent failure of my program.
Using finite elements is not an option at the moment. Is there any way to salvage this approach?