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

• MATLAB
• hunt_mat
In summary, the code provided solves the spherical diffusion equation with constant diffusion coefficient using a Crank-Nicholson scheme. However, the code diverges and the author is unsure why, as Crank-Nicholson should be unconditionally stable for diffusion. The author has included the PDE in question, as well as the scheme and code used. In addition, the author is open to suggestions and has asked for a justification of the scheme used.
Why wouldn't you think my A matrix wasn't an $N_{r}\times N_{r}$ matrix? The length of a is clearly $N_{r}$ and the command diag(a) makes an $N_{r}\times N_{r}$ matrix?

Can you do the inverse out of the loop as I suggested?

hunt_mat said:
Why wouldn't you think my A matrix wasn't an $N_{r}\times N_{r}$ matrix? The length of a is clearly $N_{r}$ and the command diag(a) makes an $N_{r}\times N_{r}$ matrix?

Because it isn't. Here is your original code which I only modified by changing N_r and N_t to 5 so I could inspect them. I also added a line that shows the size of A and N_r at the command line.

As you can see, A is (N_r-2)-by-(N_r-2).

Is this the first version of the code? I posted a second version a few posts up.

hunt_mat said:
Is this the first version of the code? I posted a second version a few posts up.

Yes, as I said... I was looking at your original code:

mfig said:
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.

The reason I was looking at your original code is in your reply #34 you mentioned that you had previously tried what I did, but you were dissuaded from using it over some conservation problem. I have not looked into your latest code because the derivation you provided in a .pdf is more complex than the classical C-N method I used, and needlessly so, IMO.

I may get a chance later to take a look, but I have meetings and office hours for the rest of today (Numerical Analysis, BTW )

hunt_mat
mfig said:
Yes, as I said... I was looking at your original code:
The reason I was looking at your original code is in your reply #34 you mentioned that you had previously tried what I did, but you were dissuaded from using it over some conservation problem. I have not looked into your latest code because the derivation you provided in a .pdf is more complex than the classical C-N method I used, and needlessly so, IMO.

I may get a chance later to take a look, but I have meetings and office hours for the rest of today (Numerical Analysis, BTW )
I think it was my boundary conditions that were causing the problems and I wanted to make sure that I could code something up that worked, and I think I managed to. I did something reasonable and the code worked nicely and there was nothing blowing up, so I am going to compare the results from my code and yours to see if they get the right thing and if they are then that's good. The next thing to consider is the speed of execution. I'll be calling this routine potentially thousands of times and I want it to be as fast as possible.

I'm also interested as to why my other code isn't working.

I coded up my original method with the new way of doing the boundary conditions and I got something which both ran and gave reasonable results and then I wrote it up as a function and tested it against a fancy adaptive timestep code and I give the comparison on my code (in blue) and the fancy one in red. I'm not sure what is gong on, it could be a question of units.

#### Attachments

• Comparison.jpg
22.7 KB · Views: 236
mfig
I'm harping back on this again. I take the equation:
$$\frac{\partial u}{\partial t}=\frac{D}{r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial u}{\partial r}\right)$$
With boundary conditions:
$$\frac{\partial u}{\partial r}\Bigg|_{r=0}=0\qquad D\frac{\partial u}{\partial r}\Bigg|_{r=1}=q(t)$$

I take the original equation, multiply it through by $r^{2}$ and integrate from 0 to 1 to obtain the following:
$$\frac{d}{dt}\int_{0}^{1}r^{2}u(t,r)dr=q(t)$$

Now my code should be able to reproduce this little fact. If I take $q(t)$ constant for example I can get an exact equation for the average value as being linearly increasing (if $q(t)$ is positive). Sadly I don't seem to get this. If my algorithm is of the form:
$$A[\mathbf{u}]^{i+1}=B[\mathbf{u}]^{i}+\mathbf{c}$$ then I can directly check this property by doing:
$$\mathbf{r}^{2}\star A[\mathbf{u}]^{i+1}=\mathbf{r}^{2}\star B[\mathbf{u}]^{i}+\mathbf{r}^{2}\star\mathbf{c}$$
where $\star$ denotes pointwise multiplication and boldface denotes a vector, I get something attaached.

#### Attachments

• av_graph.jpg
11.4 KB · Views: 101
Last edited:

• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
406
• Differential Equations
Replies
23
Views
4K
• Differential Equations
Replies
2
Views
2K
• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
1K
• Biology and Chemistry Homework Help
Replies
3
Views
303
• MATLAB, Maple, Mathematica, LaTeX
Replies
8
Views
2K
• MATLAB, Maple, Mathematica, LaTeX
Replies
3
Views
3K
• MATLAB, Maple, Mathematica, LaTeX
Replies
4
Views
881
• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
1K
• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
3K