My Crank-Nicolson code for my diffusion equation isn't working

[hunt_mat]

TL;DR Summary

I am solving the spherical diffusion equation in 1D using CN and for some reason it diverges.

I'm trying to solve the diffusion equation in spherical co-ordinates with spherical symmetry. I have included the PDE in question and the scheme I'm using and although it works, it diverges which I don't understand as Crank-Nicholson should be unconditionally stable for the diffusion. The code I'm using is:

%This program solves the spherical diffusion equation with constant
%diffusion coefficient, it has boundary conditions c_r(t,0)=0 and
%c_r(t,R)=J.
%The equation is T_t=(D/r^2)(r^2T_r)_r
clear;
N_r=50;N_t=1000; %N_r is the number of points in the electrodes, N_t is the
D=0.01; %This is the diffusion coefficient.
r=linspace(0,1,N_r);
t=linspace(0,5,N_t);% set time and distance scales
dr=r(2)-r(1);dt=t(2)-t(1); %Increments
s_1=D*dt/(2*dr^2);%Useful for the solution.
s_2=D*dt/(4*dr);
%The main problem is the matrix problem Au=v, the matrix A is constant, so
%we do this before the main loop.
u=zeros(N_t,N_r); %set up solution matrix
a=zeros(N_r-3,1);b=zeros(N_r-2,1);c=zeros(N_r-3,1);
for i=1:N_r-3
    a(i)=(s_1+s_2/r(i+1));
    b(i)=-(1+2*s_1);
    c(i)=(s_1-s_2/r(i+1));
end
b(end)=-(1+2*s_1);
%The beginning and end rows are different to the others due to the boundary
%conditions
A=diag(a,1)+diag(c,-1)+diag(b);
%Add in the corrections for the boundary conditions
A(1,1)=(2*s_1+4*s_2/r(2))/3;A(1,2)=-(3+2*s_1+4*s_2/r(2))/3; %Boundary conditions ar r=0;
A(end,end-1)=(2*s_1-4*s_2/r(end-1))/3;A(end,end)=-(3+2*s_1-4*s_2/r(end-1))/3; %Boundary conditions at r=R
%Add in intial condition:
u(1,:)=0;
%This is the RHS for the linear system of equations.
v=zeros(N_r-2,1);
%This is the boundary condition on r=R
J=0.02*exp(0.5*t);
%Solve the system
for i=2:N_t
    %Add in boundary conditions
    u(i,1)=(4*u(i,2)-u(i,3))/3;
    u(i,end)=(J(i)+4*u(i,end-1)-u(i,end-2))/3;
    %Construct the RHS for the solution
    for j=2:N_r-1
        A_1=-(s_1+s_2/r(j));
        A_2=(2*s_1-1);
        A_3=(s_1-s_2/r(j));
        v(j-1)=A_1*u(i-1,j+1)+A_2*u(i-1,j)+A_3*u(i-1,j-1);
    end
    v(end)=-(s_1+s_2/r(end-1))*u(i-1,end)+(2*s_1-1)*u(i-1,end-1)-(s_1-s_2/r(end-1))*u(i-1,end-2)-2*dr*(s_1+s_2/r(end-1))*J(i)/3;
    %Solve for the new time step
    w=A\v;
    u(i,2:N_r-1)=w';
end
figure
for i=1:N_t
    plot(r,u(i,:));
    pause(0.1);
end

Any suggestions are welcome,

 

[Chestermiller]

In Eqn. 5, it looks like you are missing a factor of 2 in front of the dc/dr terms.

 

[Chestermiller]

I wouldn't have discretized spatially the way you did it. I would have written: $$\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial c}{\partial r}\right)=\frac{1}{r_j^2}\left[\frac{r^2_{(j+1/2)}(c_{j+1}-c_j)-r^2_{(j-1/2)}(c_j-c_{j-1})}{(\Delta r)^2}\right]$$With this approximation, material is conserved not only by the differential equation but also automatically by the finite difference equation.

 

[hunt_mat]

Hi,

Thanks for this. The factor of 1/2 was indeed a problem, and I have changed it but that didn't seem to work. There are a few other things which I can try.

Can you explain your discretisation please? I've not seen that before.

 

[Chestermiller]

I can't go through your code because I am unfamiliar with the language you are using. But, have you put in any prints or done step-by-step debugging so that you can hand calculate the various terms and compare.

With regard to the discretization, it determines dc/dr at the half--way points between grid points, and then weights them in terms of the radii at these half-way points. It then evaluates the term r^2 dc/cr at the grid point by taking the differences between the values at the two adjacent half-way points. This is pretty straightforward; I feel like I used it a million times in practice.

The advantage of this scheme is that it automatically conserves mass. The scheme you used does not automatically conserve mass (in the finite difference scheme) because of discretization error. We stopped using that kind of scheme in the company where I worked many years ago for that reason.

 

[hunt_mat]

I don't think you've justified the use of your scheme. Can you provide a derivation of it? I don't trust the, "This scheme does X", without a rigorous justification.

 

[Chestermiller]

I don't think you've justified the use of your scheme. Can you provide a derivation of it? I don't trust the, "This scheme does X", without a rigorous justification.

Let $\xi=\frac{\partial C}{\partial r}$. Using central differences (second order accurate in $\Delta r$),
$$\frac{\partial}{\partial r}\left(r^2\xi\right)=\frac{(r+\Delta r/2)^2\xi_{[r+\Delta r/2]}-(r-\Delta r/2)^2\xi_{[r-\Delta r/2]}}{\Delta r}$$
$$\xi_{[r+\Delta r /2]}=\frac{C_{[r+\Delta r]}-C_{[r]}}{\Delta r}$$
$$\xi_{[r-\Delta r /2]}=\frac{C_{[r]}-C_{[r-\Delta r]}}{\Delta r}$$
Rigorous enough for you?

Also, what makes you such an expert to think that I (a numerical analyst with over 50 years of experience) am not justified in using such a simple (trivial) scheme? I (and may others) have used it successfully innumerable times in the past, and it has always given accurate answers. And, unlike you, I have actually been able to obtain stable solutions.

 

[bigfooted]

If you have to stick with the classical CN-method because the assignment says so:
have you tried solving the diffusion equation in cartesian coordinates?
Also note that it matters how you write the pde:
$C_t = D(\frac{1}{r^2}(r^2C_r)_r) = D(C_{rr} + \frac{2}{r}C_r)$
and you can also rewrite the pde using $u=Cr$ as
$u_t = Du_{rr}$,
so essentially you can rewrite your polar pde in the same way as the planar pde so this is the first thing to check in your code: can it solve the planar case?

Also note that Crank-Nicolson is a combination of forward Euler and backward Euler, and your first order derivative in eq. 4 is not forward Euler but central difference. This will affect convergence.

 

[Chestermiller]

If you have to stick with the classical CN-method because the assignment says so:
have you tried solving the diffusion equation in cartesian coordinates?
Also note that it matters how you write the pde:
$C_t = D(\frac{1}{r^2}(r^2C_r)_r) = D(C_{rr} + \frac{2}{r}C_r)$
and you can also rewrite the pde using $u=Cr$ as
$u_t = Du_{rr}$,
so essentially you can rewrite your polar pde in the same way as the planar pde so this is the first thing to check in your code: can it solve the planar case?

Also note that Crank-Nicolson is a combination of forward Euler and backward Euler, and your first order derivative in eq. 4 is not forward Euler but central difference. This will affect convergence.

I've never seen this transformation before. It's really quite dazzling. Wow!

Is it relatively new? Can you please provide a reference?

Chet

 

[Chestermiller]

If you have to stick with the classical CN-method because the assignment says so:
have you tried solving the diffusion equation in cartesian coordinates?
Also note that it matters how you write the pde:
$C_t = D(\frac{1}{r^2}(r^2C_r)_r) = D(C_{rr} + \frac{2}{r}C_r)$
and you can also rewrite the pde using $u=Cr$ as
$u_t = Du_{rr}$,
so essentially you can rewrite your polar pde in the same way as the planar pde so this is the first thing to check in your code: can it solve the planar case?

Also note that Crank-Nicolson is a combination of forward Euler and backward Euler, and your first order derivative in eq. 4 is not forward Euler but central difference. This will affect convergence.

Using u, how do you handle the BC at r =0?

 

[bigfooted]

I've never seen this transformation before. It's really quite dazzling. Wow!
Is it relatively new? Can you please provide a reference?

It's also relatively new to me, I actually found it in John Crank's book 'Mathematics of Diffusion', which I purchased last year (my edition is from 1967). Crank mentions as a reference Carslaw and Jaeger, but I don't have that book so I can't check if this transformation is in there.

Using u, how do you handle the BC at r =0?

For the boundary conditions, we have to use the transformation as well, so at $r=r_0$, a Dirichlet condition becomes
$C(r_0)=C_0 \rightarrow u(r_0)=r_0\cdot C_0$, so at r=0, a Dirichlet condition always transforms to the boundary condition $u=0$
A Neumann condition is more tricky because the derivative leads to:
$\frac{dC}{dr}= \frac{1}{r}\frac{du}{dr} - \frac{u}{r^2} = C_1$ so the Neumann condition becomes a Robin condition:
$r_0\cdot \frac{du}{dr} - u = r_0^2 \cdot C_1$
So when $r=r_0$, the Neumann boundary condition also transforms to the Dirichlet boundary condition $u=0$.

 

[Chestermiller]

It's also relatively new to me, I actually found it in John Crank's book 'Mathematics of Diffusion', which I purchased last year (my edition is from 1967). Crank mentions as a reference Carslaw and Jaeger, but I don't have that book so I can't check if this transformation is in there.
For the boundary conditions, we have to use the transformation as well, so at $r=r_0$, a Dirichlet condition becomes
$C(r_0)=C_0 \rightarrow u(r_0)=r_0\cdot C_0$, so at r=0, a Dirichlet condition always transforms to the boundary condition $u=0$
A Neumann condition is more tricky because the derivative leads to:
$\frac{dC}{dr}= \frac{1}{r}\frac{du}{dr} - \frac{u}{r^2} = C_1$ so the Neumann condition becomes a Robin condition:
$r_0\cdot \frac{du}{dr} - u = r_0^2 \cdot C_1$
So when $r=r_0$, the Neumann boundary condition also transforms to the Dirichlet boundary condition $u=0$.

I would do it differently. First of all, I wouldn't work in terms of u, but rather I would still work in terms of C, but I would work with the PDE in the following form:$$\frac{\partial C}{\partial t}=Dr\frac{\partial ^2}{\partial r^2}\left(\frac{C}{r}\right)$$I would finite-difference the equation based of this form of the PDE.

For the Neumann boundary condition at r = 0, I would start out by writing for small r $$C=C_0+(C_1-C_0)\frac{r^2}{(\Delta r)^2}$$So, at small r, I would have:$$\frac{2}{r}\frac{\partial C}{\partial r}+\frac{\partial ^2 C}{\partial r^2}=\frac{6(C_1-C_0)}{(\Delta r)^2}$$So, assuming I'm using the Method of Lines, at r = 0, I would have $$\frac{dC_0}{dt}=D\frac{6(C_1-C_0)}{(\Delta r)^2}$$This is 2nd order accurate in $\Delta r$. For all the other grid points, I would finite difference in terms of C/r.

 

[hunt_mat]

I'm aware of the transformation, it's well known, the problem is, that you need to know the value of $u$ on the boundary, which is some extra information you have to add in, so that's not appropriate in this case.

Bigfooted: My code was developed from the planar case, and I know that works. So the question becomes: For my first derivative term, what do I use, forward or backward difference? My understanding was that the CN method was to time average over the spatial components? The second order derivative was central difference as was the first order one? What's the matter with that?

Chestermiller: I like the derivation. Would I sort out the boundary points as I did in my previous post?

 

[Chestermiller]

Chestermiller: I like the derivation. Would I sort out the boundary points as I did in my previous post?

I recommend handling the boundary point at r = 0 the way I did in post #12. For the boundary point at r = R, I would do the following: $$\frac{\partial C}{\partial r}=f(t)$$So, temporarily using a "fictitious point" at $R+\Delta r$, we have the finite difference representation of this equation as $$\frac{C_{(R+\Delta r)}-C_{(R-\Delta r)}}{2\Delta r}=f(t)$$or, $$C_{(R+\Delta r)}=C_{(R-\Delta r)}+2\Delta r f(t)\tag{1}$$
Turning next to the spatial finite difference representation of PDE, we have
$$\frac{dC_R}{dt}=\frac{D}{R^2}\frac{(R+\Delta r/2)^2(C_{(R+\Delta r)}-C_R)-(R-\Delta r/2)^2(C_R-C_{(R+\Delta r)})}{(\Delta r)^2}\tag{2}$$Combining Eqns. 1 and 2 then yields:$$\frac{dC_R}{dt}=2\frac{D}{R^2}(R-\Delta r/2)^2\frac{(C_{(R-\Delta r)}-C_R)}{(\Delta r)^2}+2\frac{D}{R^2}\frac{(R-\Delta r/2)^2}{\Delta r}f(t)$$

 

[hunt_mat]

There are much better ways than using fictitious points, what you can do is have a one-sided derivative which I did. You don't have to do anything else. The scheme you've done is an explicit method which is only conditionally stable.

 

[Chestermiller]

There are much better ways than using fictitious points, what you can do is have a one-sided derivative which I did. You don't have to do anything else. The scheme you've done is an explicit method which is only conditionally stable.

I don't see how you can say it is explicit, since I didn't say anything about how the time integration was going to be done.

As far as the final result I presented in my previous post is concerned, all that matters is whether the scheme is 2nd order accurate, not whether or not I employed a fictitious point.

For a guy who can't even get his own code to work, you certainly are opinionated.

 

[hunt_mat]

I know enough not to play around with fictitious points...

For a simple enough equation, on a simple enough domain, why play around with fancy methods when a simple finite difference method will do. Including fictitious points is an example of this. Why include extra points when you don't need to?

I've done some numerical work before, but this particular equation seemed to be a problem when I should work.

 

[Chestermiller]

I know enough not to play around with fictitious points...

For a simple enough equation, on a simple enough domain, why play around with fancy methods when a simple finite difference method will do. Including fictitious points is an example of this. Why include extra points when you don't need to?

So the final equation I came up with doing it this way is not accurate? The fictitious point is not present in this final equation.

I've done some numerical work before, but this particular equation seemed to be a problem when I should work.

It's a trivially simple problem. You should be using an automatic stiff ODE integration package (based on the Gear method) to solve this, rather than that outdated Crank-Nicholson method.

 

[Chestermiller]

If I read your writeup correctly, you are not doing the time integration at the two end grid points, but only on the interior grid points. This is probably part of your numerical problem. You should also be doing the time integration at the two end grid points. The approach I presented in posts #12 and 14 shows how to set this up.

 

[hunt_mat]

What can I say, I applied the exact same methodology to the 1D cartesian problem and it works perfectly. I applied the method for the Crank-Nicholson method because it works. I don't use black box solvers when I need something to do it fast, which the CN method does. I think you're making a mountain out of a molehill.

With your discretisation of the second order term and my way of dealing with the boundary conditions and using the Crank-Nicholson method, would that produce a stable solution?

 

[Chestermiller]

What can I say, I applied the exact same methodology to the 1D cartesian problem and it works perfectly. I applied the method for the Crank-Nicholson method because it works. I don't use black box solvers when I need something to do it fast, which the CN method does. I think you're making a mountain out of a molehill.

With your discretisation of the second order term and my way of dealing with the boundary conditions and using the Crank-Nicholson method, would that produce a stable solution?

Sure.

 

[Chestermiller]

I wouldn't call the prestigious International Mathematics and Statistics Library of subroutines (IMSL Library) a "black box." And, if you need something to do it fast, with that solver, all you need to do is provide a subroutine which maps the vector of concentration values at the grid points into the vector of time derivatives. The ODE solver does the rest. I can't believe you are not familiar with this.

 

[dameng]

I found a small mistake, probably a typo in the pdf you attached. From Equation 1 to 5, you probably forget to times "2" to the 'dc/dr' part.
for example,
right side of Equation(5) should be D/2* [(Ci+1,j+1 - Ci+1,j-1)/rj*dr...] instead of D/2* [(Ci+1,j+1 - Ci+1,j-1)/2rj*dr...]

 

[Chestermiller]

I found a small mistake, probably a typo in the pdf you attached. From Equation 1 to 5, you probably forget to times "2" to the 'dc/dr' part.
for example,
right side of Equation(5) should be D/2* [(Ci+1,j+1 - Ci+1,j-1)/rj*dr...] instead of D/2* [(Ci+1,j+1 - Ci+1,j-1)/2rj*dr...]

See post #2.

 

[hunt_mat]

Let $\xi=\frac{\partial C}{\partial r}$. Using central differences (second order accurate in $\Delta r$),
$$\frac{\partial}{\partial r}\left(r^2\xi\right)=\frac{(r+\Delta r/2)^2\xi_{[r+\Delta r/2]}-(r-\Delta r/2)^2\xi_{[r-\Delta r/2]}}{\Delta r}$$
$$\xi_{[r+\Delta r /2]}=\frac{C_{[r+\Delta r]}-C_{[r]}}{\Delta r}$$
$$\xi_{[r-\Delta r /2]}=\frac{C_{[r]}-C_{[r-\Delta r]}}{\Delta r}$$
Rigorous enough for you?

Also, what makes you such an expert to think that I (a numerical analyst with over 50 years of experience) am not justified in using such a simple (trivial) scheme? I (and may others) have used it successfully innumerable times in the past, and it has always given accurate answers. And, unlike you, I have actually been able to obtain stable solutions.

So, with that method I have the general stencil:
$$sr_{j}^{+2}u_{i+1,j+1}-(r_{j}^{2}+s(r_{j}^{+2}+r_{j}^{-2}))u_{i+1,j}+sr_{j}^{-2}u_{i+1,j-1}=-sr_{j}^{+2}u_{i,j+1}+(s(r_{j}^{+2}+r_{j}^{-2})-r_{j}^{2})u_{i,j}-sr_{j}^{-2}u_{i,j-1}$$

where $r_{j}^{+}=r_{j}+\delta r/2$, $r_{j}^{-}=r_{j}-\delta r/2$ and $s=D\delta t/(2\delta r^{2})$

The boundary conditions can be done as the following:
$$\frac{\partial u}{\partial r}\Bigg|_{r=0}=0\Rightarrow \frac{u_{i,1}-u_{i,-1}}{2\delta r}=0$$

and

$$\frac{\partial u}{\partial r}\Bigg|_{r=1}=f(t)\Rightarrow\frac{u_{i,N+1}-u_{i,N-1}}{2\delta r}=f(t_{i})=f_{i}$$

Is this what you would suggest?

 

[Chestermiller]

So, with that method I have the general stencil:
$$sr_{j}^{+2}u_{i+1,j+1}-(r_{j}^{2}+s(r_{j}^{+2}+r_{j}^{-2}))u_{i+1,j}+sr_{j}^{-2}u_{i+1,j-1}=-sr_{j}^{+2}u_{i,j+1}+(s(r_{j}^{+2}+r_{j}^{-2})-r_{j}^{2})u_{i,j}-sr_{j}^{-2}u_{i,j-1}$$

where $r_{j}^{+}=r_{j}+\delta r/2$, $r_{j}^{-}=r_{j}-\delta r/2$ and $s=D\delta t/(2\delta r^{2})$

The boundary conditions can be done as the following:
$$\frac{\partial u}{\partial r}\Bigg|_{r=0}=0\Rightarrow \frac{u_{i,1}-u_{i,-1}}{2\delta r}=0$$

and

$$\frac{\partial u}{\partial r}\Bigg|_{r=1}=f(t)\Rightarrow\frac{u_{i,N+1}-u_{i,N-1}}{2\delta r}=f(t_{i})=f_{i}$$

Is this what you would suggest?

In post #12, I discussed what I would do at i=0. At the N'th grid point, I would combine your boundary difference equation with the regular global difference equation to eliminate the point N+1. That way, I would be solving for the time variation at i=N also.

 

[hunt_mat]

Yes, agreed. I understand this now. I

In post #12, I discussed what I would do at i=0. At the N'th grid point, I would combine your boundary difference equation with the regular global difference equation to eliminate the point N+1. That way, I would be solving for the time variation at i=N also.

Yes, I understand that now and I am getting a better feeling for what I'm doing, the stencil at r=0 is given by:
$$u_{i+1,1}-u_{i+1,0}=-(u_{i,1}-u_{i,0})$$
and the stencil at r=1 is given by the (slightly more complicated):
$$-(r_{N}^{2}+s(r_{N}^{+2}+r_{N}^{-2}))u_{i+1,N}-s(r_{N}^{+2}+r_{N}^{-2})u_{i+1,N-1}=-(-s(r_{N}^{2}-s(r_{N}^{+2}+r_{N}^{-2}))u_{i,N}+s(r_{N}^{+2}+r_{N}^{-2})u_{i,N-1})+2sr_{N}^{+2}\delta r(f_{i}-f_{i-1})$$

The problem then becomes:
$$A\mathbf{u}^{i+1}=-B\mathbf{u}^{i}+b$$

I can reduce the number of computations by multiplying $A^{-1}$, and computing things before the loop.

 

[hunt_mat]

Hi,

I have written up my calculations, and coded it up, and as you predicted, there was no blow up. However, there *is* another issue, when I set:

f=-0.375*ones(N_t,1);
c_0=0.95*ones(N_r,1);
D=1.5;

Then when I plot the boundary value for my solution as a function of time, then I should obtain a straight line but I don't. This is wrong. Where have I made my error?

 

[Chestermiller]

Hi,

I have written up my calculations, and coded it up, and as you predicted, there was no blow up. However, there *is* another issue, when I set:

f=-0.375*ones(N_t,1);
c_0=0.95*ones(N_r,1);
D=1.5;

Then when I plot the boundary value for my solution as a function of time, then I should obtain a straight line but I don't. This is wrong. Where have I made my error?

I don't understand. Please plot a graph to illustrate your issue.

 

[hunt_mat]

When I use a different code, coming from the finite element method, I get the attached solution for the boundary values. My code seems to blow up.

 

[mfig]

Here is how I would do this, FWIW. Let me know if I have made any mistakes.

 

[hunt_mat]

Whereabouts are you evaluating r on your first yellow boxed equation?

I get what you're doing and your code looks very nice. When I tried to code before, I got blow up for my code. I essentially did something similar to what you did but I treated the boundary conditions differently.

 

[mfig]

Whereabouts are you evaluating r on your first yellow boxed equation?

I get what you're doing and your code looks very nice. When I tried to code before, I got blow up for my code. I essentially did something similar to what you did but I treated the boundary conditions differently.

At the node:$r_k = k \Delta r$

I should probably include the subscript on r in future versions!

 

[hunt_mat]

Did you have a look at the code I wrote? Did you see any errors which jumped out at you? I tried your idea before as that is what made the most sense(to me at the time) but I was later informed that it violated some basic mass conservation?

Why can't I simply write:
$$\frac{T_{i,1}-T_{i,-1}}{2\delta r}=0$$
For my flux symmetry condition? Treat it similar to how you treated the boundary condition?

If I write the system of equations as:
$$\mathbf{A}T^{n+1}=\mathbf{B}T^{n}+\mathbf{c}^{n+1,n}$$
Can I then write:

$$T^{n+1}=\mathbf{A}^{-1}\mathbf{B}T^{n}+\mathbf{A}^{-1}\mathbf{c}^{n+1,n}$$

So I don't have to solve a system of linear equations at every timestep?

 

[mfig]

I am not sure what is meant by the statement that the standard Crank-Nicolson FD spatial formulation (which is the average of a simple implicit and explicit scheme) violates the conservation of mass. I have never heard this. I have used this for diffusion problems forever and had no such trouble, but maybe I have been lucky?

I think you can use what you wrote for the flux formulation at the center. The trouble is to avoid the singularity at the origin in the original PDE, which is where L'Hopital for r=0 comes in: to reformulate the PDE at that point. A subsequent application of the C-N formulation here, with your suggested symmetry consideration, works just fine, AFAIK.

I didn't see any error in your original code, but I did notice something strange at least to me. Your A matrix is not N-r-by-N_r. Compare how I constructed the M matrix in my code.