Finite difference method - implementation and stability

107 Views Asked by At

I have a PDE for a function $P(r,t)$ in spherical co-ordinates of the form

$P_{t} = D(P_{rr} + (2/r)P_{r}) - a\Omega\left(\frac{P}{P + k} \right) $

where subscripts denote partial derivative with respect to the subject and $D,\Omega$ and $k$ are known constants. I also know the following conditions;

$P(r,0) = 0$ (Initial condition)

$P(0,t) = 100$ (Boundary condition at zero)

As I understand it, I can approximate the solution to such an equation by using the finite difference method and solving the algebraic expression

$P(i,j + 1) = P(i,j) + \frac{D \Delta t}{(\Delta r)^2}\left(P(i + 1,j) -2P(i,j) + P(i -1,j) \right) + \frac{2D \Delta t}{r \Delta r}\left(P(i + 1,j) - P(i,j) \right) - a\Omega \left(\frac{P(i,j)}{P(i,j) + k} \right) $

However, I am a little curious about this; surely the choice of $\Delta t$ and $\Delta r$ are not independent and have serious repercussions for the estimated value? I know for a standard diffusion equation that the Courant – Friedrichs – Lewy condition demands that $\frac{D\Delta t}{(\Delta r)^2} \leq \frac{1}{2}$ for convergence, but I imagine the situation here is more complex? I am relatively new to this area and would be grateful for any clarification before I try to code this up!