I'm trying to solve the PDE for $c(r,t)$ $$c_t=(1/r)(rJ)_r$$ using Crank-Nicolson, and I'm having difficulty with the boundary conditions. $J$ is the flux, the initial condition is $c(0,r)=c_{init}$, and the flux is 0 at both boundaries $r=r_0$ and $r=r_1$. The flux is $J=rDc_r-s\omega^2r^2c$; $D$ and $s$ may be functions of $c$, but are often constants, as they are treated in this example.

So for the Crank-Nicolson scheme I use these approximations: $$c_t=(c_{n+1,j}-c_{n,j})/\Delta t$$ $$c_{rr}=\left({1\over 2}\right)\left[{c_{n,j+1}-2c_{n,j}+c_{n,j-1}\over (\Delta r)^2}+{c_{n+1,j+1}-2c_{n+1,j}+c_{n+1,j-1}\over (\Delta r)^2}\right]$$ $$c_r=\left({1\over 2}\right)\left[{c_{n,j+1}-c_{n,j-1}\over \Delta r}+{c_{n+1,j+1}-c_{n+1,j-1}\over \Delta r}\right]$$ $$c=\left({1\over 2}\right)\left[c_{n,j}+c_{n+1,j}\right]$$ where $n=$ the time step ($t=t_0+n\Delta t$), and $j=$ the position ($r=r_0+j\Delta r$, $j=0,...,J$).

For the boundary conditions, I want $J(r_0,t)=0$ and $J(r_1,t)=0$ for all $t$.

My first attempt was to use "ghost points" at $j=-1$ and $j=J+1$. I used the centered difference at $j=0$ and $j=J$, respectively, for the derivative $c_r$, and the values $c_{n,0}$ and $c_{n,J}$: $$c_r={c_{n,1}-c_{n,-1}\over 2\Delta r}$$ $$c_r={c_{n,J+1}-c_{n,J-1}\over 2\Delta r}$$ Substituting these approximations into the expression for the flux $J$ yields $c_{n,-1}$ in terms of $c_{n,0}$ and $c_{n,1}$, and $c_{n,J+1}$ in terms of $c_{n,J-1}$ and $c_{n,J}$.

With the ghost points eliminated, the $J+1 \times J+1$ tridiagonal matrices ${\bf A}$ and ${\bf B}$ are used to solve for the vector ${\bf c}_{n+1}$ from ${\bf c}_n$ for each time step.

Yet, when the numerical solution is calculated, instead of the expected exponential curve, becoming large near $r=r_1$, $c$ becomes tiny everywhere, and the values over time make it apparent that there is flux out at both boundaries. That is the issue I am grappling with.

I thought about using a Crank-Nicolson average to approximate $c_r$ at the boundaries, but then I get four additional unknowns, and I don't know how to deal with this.

Any suggestions?

  • $\begingroup$ Please, add more information. If you discretize your origin PDE, then you need to tell how you discretize $(r J)_r$. If you modify your PDE by substituting $J$ into the PDE and even computing the derivatives with respect to $r$, please confirm it or even write the resulting form of PDE explicitly. $\endgroup$ Sep 27, 2016 at 17:55
  • $\begingroup$ Yes, I substituted for $J$ in the original PDE before discretizing anything. The PDE after the substitution is $c_t=\left({1\over r}\right)(rDc_r - s\omega^2f^2c)$. I see I made a mistake in the sign of $J$, but that doesn't affect the constraint imposed by the boundary conditions. $\endgroup$
    – Woody20
    Sep 27, 2016 at 18:04
  • $\begingroup$ The $f^2$ in the eqn for $c_t$ should be $r^2$ $\endgroup$
    – Woody20
    Sep 27, 2016 at 18:29
  • $\begingroup$ How do you discrete $(r c_r)_r$ and $(r^2 c)_r$? Are you missing $r$ in your substitution opposite to PDE? $\endgroup$ Sep 28, 2016 at 6:40
  • $\begingroup$ Sorry for the typos in my first comment. Here is the PDE after substitutions: $$c_t=Dc_{rr}+{D\over r}c_r-2s\omega^2c-s\omega^2rc_r$$ In this equation, I used the approximations for $c_t$, $c_{rr}$, $c_r$ and $c$ as shown in the OP. $\endgroup$
    – Woody20
    Sep 28, 2016 at 18:30

2 Answers 2


Let me give an answer that is a general comment on prescribed zero flux for advection-diffusion (or convection-diffusion) PDE that is an important topic and it might be (but not necessary) the problem in your situation. Zero flux boundary condition is very non-standard in this case, so you might check if there is no misinterpretation of your problem.

Your flux $J$ can be split into two parts: $$ J=J_{diff} + J_{adv} = - r D c_r + s \omega^2 r^2 c $$ Let, formally, $s \omega^2 = 0$ then $J$ is purely diffusive flux and zero flux boundary conditions at $r_0$ and $r_1$ make sense as it is parabolic PDE. The situation is different if you consider for a moment that $D=0$ and you have only advective flux with "speed" $v(r)=s \omega^2 r^2$ and your PDE would be hyperbolic. For this type of problems one prescribes boundary conditions only for "inflow" boundaries, when the flow enters the domain, in your case when $v(r_0)>0$ or $v(r_1)<0$. In the opposite case for so-called "outflow" boundaries (consider the opposite inequalities) NO boundary conditions are required as the flow simply leaves the domain.

Let me consider your constants $s$ and $\omega$ being positive. Typical boundary conditions of advection for inflow are Dirichlet boundary conditions, i.e. in your case $c(r_0,t) = c_{given}(t)$. No boundary conditions are required for the advective flux at $r=r_1$.

Now if all constants are non-zero, a standard treatment for advection-diffusion flux at boundary is that you consider them separately. For instance prescribe the zero diffusive flux everywhere and prescribe the concentration in the advective flux only at inflow boundary.

If you prescribed in your code the zero advective-diffusive flux everywhere it still make sense for the inflow, because it is equivaklent to zero diffusive flux and zero concentration at advective flux. But for the outflow it can lead to non-physical behavior, because in fact what you are prescribing is that the diffusive flux is equal minus advective flux. This is especially serious if you have advection dominated problem when so called Peclet number is large, i.e. $D<s \omega^2 r_1^2$.

If this still does not solve your problem, give us a few pictures of time development for your numerical solution, just one or two time steps from constant initial condition.

P.S. In general the zero flux B.C. in Crank-Nicolson method are considered for each time point separately, so the ghost values at $n$ and $n+1$ are eliminated separately.

  • $\begingroup$ I do not understand your final "P.S". There are ghost points for each boundary in my case, but I cannot solve for both the time $n$ and $n+1$ in terms of known values. If there is a way to do so, please tell me. $\endgroup$
    – Woody20
    Oct 3, 2016 at 20:35
  • $\begingroup$ The physical situation is a finite tube with a prescribed initial mass distribution, which redistributes over time. There is no way for mass to enter or leave the tube; thus, the no-flux boundary conditions. $\endgroup$
    – Woody20
    Oct 3, 2016 at 20:50
  • $\begingroup$ Since I am still having difficulties, can you expand on your comments regarding the "outflow" boundary condition? In the physical process, both boundaries are definitely 0-flux. The equilibrium solution , solved analytically, is a positive exponential. $\endgroup$
    – Woody20
    Oct 10, 2016 at 18:15

The main reason my solution did not work was that

$$c_r=\left({1\over 2}\right)\left[{c_{n,j+1}-c_{n,j-1}\over \Delta r}+{c_{n+1,j+1}-c_{n+1,j-1}\over \Delta r}\right]$$ is incorrect. For the centered-difference approximation for $c_r$, the denominator should be $2\Delta r$, not $\Delta r$. There is an additional factor of $1\over 2$ from the Crank-Nicolson stencil.

Making this correction fixed most, but not all, of the mass loss. There is still some problem, most likely in the application of the boundary conditions.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.