Hi. I am willing to apply the Alternating Direction Method in a 3 dimensional domain, for the diffusion equation with non constant coefficients:

$\displaystyle \displaystyle \frac{1}{c}\frac{\partial \phi(\mathbf{r},t)}{\partial t}- \nabla \cdot \left ( \kappa (\mathbf{r})\nabla \phi(\mathbf{r},t) \right) + \mu(\mathbf{r})\phi(\mathbf{r},t)=q(\mathbf{r},t)$

However, I'm not sure if ADI methods are stable for this equation in 3 dimensions. Some schemes, like the Peaceman-Rachford scheme are not stable in 3D even for the constant coefficient case. Most theorem I've found in the literature, requiere the commutativity of the operator to be split, which could be written as:

$\displaystyle \displaystyle \mathcal{\hat L}=\mathcal{\hat L}_x+\mathcal{\hat L}_y+\mathcal{\hat L}_z=\sum_i \mathcal{\hat L}_i$

such that:

$\displaystyle \mathcal{\hat L}_i=-\frac{\partial}{\partial x_i}\left(\kappa(\mathbf{r}) \frac{\partial \phi(\mathbf{r},t)}{\partial x_i}\right)+\frac{\mu(\mathbf{r})}{3}$

So, it is clear that the operators don't commute, i.e.:

$\displaystyle \left[ \mathcal{\hat L}_i,\mathcal{\hat L}_j \right]\neq 0$ for $\displaystyle i \neq j$

So, I have built a scheme to solve this equation in 3D, I have already written a program in 2D that seems to work (however, I've found that changing the value of c the scheme gets unstable, but I have to check for coding mistakes). In 3D I have splitted the operator, and I am ready to write the code, but before I do that I would like to know if the scheme would be stable. I have read the original paper by Douglas, and others, but it is not clear to me if this will be stable or not in 3D (the theorems that guarantees the stability and convergence of the scheme can't be applied due to the non commutative of the operators).

I have used a Crank-Nicolson scheme, which I think is stable, but I'm not sure. The theorems of Douglas says that if the original scheme which will be splitted is stable, under some assumptions (one is the commutativity I've mentioned), the splitted scheme will also be stable. However, this is not a necessary condition, so I don't know what should happen with my particular problem. I could just write the code and test it experimentally, but if then I make some other silly mistake in coding, I could think that the scheme isn't stable when it actually is, and if it isn't stable I would like not to invest time in coding it.

For the discretized differential operator I have used:

$\displaystyle -\frac{\partial}{\partial x_j}\left(\kappa(\mathbf{r}) \frac{\partial \phi(\mathbf{r},t)}{\partial x_j}\right) \sim -\frac{1}{\Delta x_j^2} \left( \kappa_{j+1/2,k,l} ( \phi_{j+1,k,l}-\phi_{j,k,l} )-\kappa_{j-1/2,k,l} (\phi_{j,k,l}-\phi_{j-1,k,l}) \right)$

Can anyone help me with this before I proceed?