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

[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:
<br /> \frac{T_{i,1}-T_{i,-1}}{2\delta r}=0<br />
For my flux symmetry condition? Treat it similar to how you treated the boundary condition?

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

<br /> T^{n+1}=\mathbf{A}^{-1}\mathbf{B}T^{n}+\mathbf{A}^{-1}\mathbf{c}^{n+1,n}<br />

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.

 

[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:
<br /> \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.