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

[hunt_mat]

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?

 

[mfig]

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

 

[hunt_mat]

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

 

[mfig]

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:

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]

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.

 

[hunt_mat]

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.

 

[hunt_mat]

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.