# Mass conservation in radially symmetric parabolic PDE problems

1. Mar 29, 2010

### ndalchau

Dear all,
I'm trying to solve the 2d heat equation in a radially symmetric domain, numerically using the Crank-Nicolson method. i.e.

$$\dfrac{\partial u}{\partial t} = D\left( \dfrac{\partial^2u}{\partial r^2}+\dfrac{1}{r}\dfrac{\partial u}{\partial r}\right)$$

Applying the Crank-Nicolson method basically results in a recurrence relation:

$$(-q+\frac{z}{r})u_{i-1,j+1} + (1+2q)u_{i,j+1} - (q+\frac{z}{r})u_{i+1,j+1} = (q-\frac{z}{r})u_{i-1,j} + (1-2q)u_{i,j} + (q+\frac{z}{r})u_{i+1,j}$$

where $$q=\dfrac{\delta t}{2(\delta r)^2}$$, $$z=\dfrac{\delta t}{4\delta r}$$ and $$u_{i,j}$$ is the solution at $$r=i\delta r, t=j\delta t$$.

You can write this into a matrix equation of the form $$Au\{j+1\}=Bu\{j\}$$, which basically enables you to solve the problem. This all works very nicely for the 1d heat equation (the same differential equation but without the first spatial derivative). However, I now have the problem that the solution doesn't conserve mass. That is, given some initial condition across $$r\in (0,R)$$, the sum of the solution points decreases over time.

Clearly, if the column sums of A and B are the row vector of 1s, then I am guaranteed mass conservation. However, this is not the case for my Crank-Nicolson implementation here. This results from the $$\frac{z}{r}$$ changing values and signs between neighbouring rows, which means they don't cancel.

Anyone got any comments? Is there a way of solving these parabolic PDEs in radially symmetric domains that preserves the mass conservation law? Or is my implementation incorrect maybe?

Your assistance would be greatly appreciated,

Neil