I’m trying to come up with a time-integration procedure for the following system -
\frac{\partial c}{\partial t} = D_c \nabla^2 c -\beta c + k\delta(\mathbf{x} - \mathbf{x}_a(t) ) = f(c,\nabla^2c,\mathbf{x}_a) \tag{1}
\frac{d\phi}{dt} = -(\kappa/\gamma_R)\begin{pmatrix}-\sin \phi \\ \cos\phi\end{pmatrix}\cdot\nabla = g(\phi,\nabla c)\tag{2}
\frac{d\mathbf{x}_a(t)}{dt} = v\begin{pmatrix} \cos\phi \\ \sin\phi \end{pmatrix} = h(\phi). \tag{3}
where \mathbf{x}_a(t) is the position of some particle, c is an ambient field and D_c, \beta, k, \kappa, \gamma_R and v are all positive constants.
The time-integration is severely-complicated by the presence of the \delta and the nonlinearity in (2). One idea is to use the following pseudo-implicit Euler scheme:
\frac{c^{n+1} - c^n}{dt} = f(c^{n+1},\nabla^2c^{n+1},\mathbf{x}_a^n)
\\
\frac{\phi^{n+1}-\phi^n}{dt} = g(\phi^{n+1},\nabla c^{n+1})
\\
\frac{\mathbf{x}_a^{n+1} - \mathbf{x}_a^{n}}{dt} = h(\phi^{n+1})
however, the fact that I’m using \mathbf{x}_a^n instead of \mathbf{x}_a^{n+1} is a bit suspect and I’m not sure I could guarantee unconditional stability.
Does anyone have any ideas?