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

hunt_mat

Homework Helper
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:
Matlab:
%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,

Attachments

• 109.9 KB Views: 32
Last edited by a moderator:
Related MATLAB, Maple, Mathematica, LaTeX, Etc News on Phys.org

Chestermiller

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

• hunt_mat

Chestermiller

Mentor
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

Homework Helper
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

Mentor
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

Homework Helper
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

Mentor
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.

• hunt_mat

Chestermiller

Mentor
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

Mentor
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

Mentor
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

Homework Helper
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 spacial 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

Mentor
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^2+(\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

Homework Helper
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

Mentor
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.

Last edited:

hunt_mat

Homework Helper
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

Mentor
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

Mentor
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

Homework Helper
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

Mentor
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

Mentor
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

Mentor
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.

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

Physics Forums Values

We Value Quality
• Topics based on mainstream science
• Proper English grammar and spelling
We Value Civility
• Positive and compassionate attitudes
• Patience while debating
We Value Productivity
• Disciplined to remain on-topic
• Recognition of own weaknesses
• Solo and co-op problem solving