Crank-Nicholson solution to the cylindrical heat equation

• A
• hunt_mat
In summary: AT^{i+1}+AT^i)+q$$or$$\frac{dT}{dt}=-\frac{1}{2}(AT+A^TT+q)$$You can write this in your notation with an appropriate change of basis if you want, but I don't know what purpose this would serve. I thought you were just going to implement my suggestion. It should only take a couple of minutes to do that.In summary, the conversation discusses the process of solving the radially symmetric heat equation with an internal heat source to model the heating of a cylindrical battery. The code used is based on the Crank-Nicholson method and involves constructing matrices and hunt_mat Homework Helper TL;DR Summary I used a Crank-Nicholson method to solve a radially symmetric heat equation. I get some odd results. I'm not sure if it's my 1) Numerical scheme, 2) Parameters chosen or 3) Code which is wrong. Hi, I am solving the radially symmetric heat equation with an internal heat source(this is meant to model the heating of a cylindrical battery). It's meant to model heat in a cylinder with conduction to the environment, so my outer boundary condition is Newton's law of cooling. The T in the equation is excess temperature from the environment. I used a Crank-Nicholson which I think I've got right. I plotted a cross-section to see if it looked okay and it didn't basically and now I'm not sure what's going on. I've checked my code, and it matches my numerical algorithm. I think I have my physics correct but again I could be wrong. My code is: Code: %This solves the heat equation with source term in cylindrical %co-ordinates. %T_t=(k/r)*(rT_r)_r+Q(t,r) N_r=300; %Radial steps N_t=200; %Time steps r=linspace(0,1,N_r); %Radial co-ordinate t=linspace(0,20,N_t)'; %Time co-ordinate dr=r(2);dt=t(2); %Increments T=zeros(N_t,N_r); a_2=zeros(N_r-1,1); a_3=zeros(N_r-1,1); %Physical characteristics k=3; %Thermal conductivity h=2.5; %Thermal exchange constant theta=1; A=k*dt/(4*dr); B=k*dt/(2*dr^2); %Useful constants Q=10*t.^2.*exp(-t); %Construct the matrices %(i+1)th matrix a_1=-(theta+2*A)*ones(N_r,1); a_1(1)=-(4*A-theta); a_1(N_r)=-2*h*dr*(A+B/r(N_r))-2*A-theta; a_2(2:N_r-1)=A+B./r(2:N_r-1); a_2(1)=4*A; a_3(1:N_r-2)=A-B./r(2:N_r-1); a_3(N_r-1)=2*A; alpha=diag(a_1)+diag(a_2,1)+diag(a_3,-1); %ith matrix b_1=(2*A-theta)*ones(N_r,1); b_1(1)=-(theta+4*A); b_1(N_r)=2*A-theta+2*h*dr*(A+B/r(N_r)); b_2=-a_2; b_3=-a_3; beta=diag(b_1)+diag(b_2,1)+diag(b_3,-1); %Constructing the solution: b=zeros(1,N_r); for i=2:N_t b=-0.5*(Q(i,:)'+Q(i-1,:)')*dt; u=T(i-1,:)'; C=beta*T(i-1,:)'+b; x=alpha\C; T(i,:)=x'; end for i=1:N_t plot(r,T(i,:)); pause(0.1); end Can anyone see where I have made a glaring error? Attachments • cross_section.jpg 76.7 KB · Views: 357 • cylindrical_heat.pdf 112.2 KB · Views: 425 Delta2 I think the issue here is how you handle the boundary conditions, in particular, the boundary condition at r = 0. I would approach the formulation of the difference equations a little differently from the way you have done it. It would be a two step process, first using the method of lines to discretize the differential equation spatially into a coupled set of 1st order ODEs in time, and then secondly to apply the Crank Nicholson formalism to the set of ODEs. I will only show the first step here. The 2nd step is then obvious. So, applying the 1st step to the general region, I get$$\frac{dT_j}{dt}=\alpha \frac{T_{J+1}(1+1/(2j))-2T_j+T_{j-1}(1-1/(2j))}{(\Delta r)^2}+q$$where ##\alpha=\kappa/\theta## and ##q=Q/\theta## Now for the BC at r = 0. Since dT/dr is zero there, we can approximate the temperature profile in the vicinity of r = 0 by a series of quadratic accuracy$$T=A(t)+B(t)r^2\tag{1}$$, with ##T=T_0## at r = 0 and ##T=T_1## at ##r=\Delta r##. So ##A=T_0## and ##B=\frac{T_1-T_0}{(\Delta r)^2}##, From Eqn. 1, it follows that$$\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial T}{\partial r}\right)=4B(t)$$. So, for the BC at r = 0, the spatially discretized differential equation becomes:$$\frac{dT_0}{dt}=\alpha \frac{4(T_1-T_0)}{(\Delta r)^2}+q$$This is how the 2nd order accurate spatially discretized boundary condition should be applied at r = 0. Last edited: hunt_mat Hmmm, I was hoping that I could do the same thing as I did with the spherical diffusion equation. I am aware of the method of lines but that's all I am. I've seen the method for dealing with the r=0 term before and I applied it here, so why does it work for spherical and not cylindrical? I've always viewed numerics as a bit of a dark art. Delta2 hunt_mat said: Hmmm, I was hoping that I could do the same thing as I did with the spherical diffusion equation. I am aware of the method of lines but that's all I am. I've seen the method for dealing with the r=0 term before and I applied it here, so why does it work for spherical and not cylindrical? I've always viewed numerics as a bit of a dark art. It's. not clear what you are saying. Are you saying that you tried the discretization method I suggested in spherical coordinates (where the 4 becomes a 6), or are you saying that you tried something entirely different. I've used the method I've suggested a huge number of times for both spherical coordinates and cylindrical coordinates, and it's worked flawlessly in all cases. I don't really view numerics as a dark art. It seems pretty straightforward to me. If you haven't tried the method I suggested, please do. It should be easy enough to set up within your current framework. hunt_mat and Delta2 No, I computed the limit as $r\rightarrow 0$ as: $$\lim_{r\rightarrow 0}\frac{1}{r}\frac{\partial u}{\partial r}=\frac{\partial^{2}u}{\partial r^{2}}$$ So for the first point, I would simply discretise a cartesian diffusion equation and I treat everything else normally. Also for the method of lines, you basically have the equation: $$\frac{d\mathbf{T}}{dt}=\mathbf{A}\mathbf{T}$$ Where A is some matrix, could I simply write $\mathbf{A}=\mathbf{P}^{T}\mathbf{D}\mathbf{P}$, where $\mathbf{D}$ is a diagonal matrix, then if I just multiply my equation through by $\mathbf{P}$ I obtain: $$\frac{d}{dt}(\mathbf{P}\mathbf{T})=\mathbf{D}(\mathbf{P}\mathbf{T})$$ Which results in a set of decoupled equations which are simple to integrate. I have no doubt that you can get your code to work perfectly, I, on the other hand, have to fight to get something to work. hunt_mat said: No, I computed the limit as $r\rightarrow 0$ as: $$\lim_{r\rightarrow 0}\frac{1}{r}\frac{\partial u}{\partial r}=\frac{\partial^{2}u}{\partial r^{2}}$$ So for the first point, I would simply discretise a cartesian diffusion equation and I treat everything else normally. . This is not a valid discretization of the spatial derivative. The expression I gave, on the other hand actually is a valid discretization. I can write the corresponding Crank-Nicholson version of the equation for the grid point at r = 0 so that you can incorporate it into your code if that is what you need. But I was hoping you would be able to do this on your own. It should take less than 5 minutes for you to change to code. Until you do that, in my judgment, you will not be able to get a satisfactory solution. Also for the method of lines, you basically have the equation: $$\frac{d\mathbf{T}}{dt}=\mathbf{A}\mathbf{T}$$ Where A is some matrix, could I simply write $\mathbf{A}=\mathbf{P}^{T}\mathbf{D}\mathbf{P}$, where $\mathbf{D}$ is a diagonal matrix, then if I just multiply my equation through by $\mathbf{P}$ I obtain: $$\frac{d}{dt}(\mathbf{P}\mathbf{T})=\mathbf{D}(\mathbf{P}\mathbf{T})$$ Which results in a set of decoupled equations which are simple to integrate. I have no doubt that you can get your code to work perfectly, I, on the other hand, have to fight to get something to work. No. My equation using the method of lines is of the form:$$\frac{dT}{dt}=AT+q$$where q is the vector of heat source at the grid points. The Crank-Nicholson approximation to this set of ODEs is:$$\frac{T^{i+1}-T^i}{\Delta t}=\frac{1}{2}(AT^{i+1}+q^{i+1}+AT^i+q^i)$$hunt_mat Why isn't it a valid discretization? I saw this in various peer-reviewed research articles, so I naturally assumed it was perfectly fine. So basically I would only need to change a row of the $\mathbf{A}$ matrix and that should be fine? I'm battling with a finite element code at the moment which isn't working but did in the spherical case. hunt_mat said: Why isn't it a valid discretization? I saw this in various peer-reviewed research articles, so I naturally assumed it was perfectly fine. I don't want to have to go into that. Just please try the approximation that I showed you. So basically I would only need to change a row of the $\mathbf{A}$ matrix and that should be fine? I'm battling with a finite element code at the moment which isn't working but did in the spherical case. Yes, that's all you would have to do. Using the notation from previous posts, my new stencil will now become: $$-\alpha A_{j,j-1}u_{i+1,j+1}+(\theta-\alpha A_{j,j})-\alpha A_{j,j+1}u_{i+1,j+1}=\alpha A_{j,j-1}u_{i,j-1}+(\theta+\alpha A_{j,j})u_{i,j}+\alpha A_{j,j+1}+\frac{\delta t}{2}(Q_{i+1}+Q_{i})$$ Certain rows and columns will have to be adjusted for boundary conditions. This is basically the code? hunt_mat said: Using the notation from previous posts, my new stencil will now become: $$-\alpha A_{j,j-1}u_{i+1,j+1}+(\theta-\alpha A_{j,j})-\alpha A_{j,j+1}u_{i+1,j+1}=\alpha A_{j,j-1}u_{i,j-1}+(\theta+\alpha A_{j,j})u_{i,j}+\alpha A_{j,j+1}+\frac{\delta t}{2}(Q_{i+1}+Q_{i})$$ Certain rows and columns will have to be adjusted for boundary conditions. This is basically the code? I think you made so typo's in the equation you wrote. It should be exactly the same as you previously obtained, except for row j = 0 (i.e., the ##T_0## equation). Did you not understand that? No, the code will be different. The way it was discretised was different and hence the stencil will be different. Otherwise, what's the difference between the method of lines and normal finite difference? hunt_mat said: No, the code will be different. The way it was discretised was different and hence the stencil will be different. Otherwise, what's the difference between the method of lines and normal finite difference? The method of lines is just a process for deriving the finite difference discretization in a two step-process, by first discretizing with respect to the spatial parameter and then discretizing with respect to the time parameter. It should deliver the exact same Crank-Nicholson finite difference equation(s). I will derive the Crank-Nicholson finite difference equation using your parameters and demonstrate that it is exactly the same. Your starting differential equation is$$θ\frac{\partial T}{\partial t}=\frac{\kappa}{r}\frac{\partial}{\partial r}\left(r\frac{\partial T}{\partial r}\right)+Q$$​ The spatial finite difference discretization of this is:$$\theta \frac{dT_j}{dt}=\kappa\left[\frac{T_{j+1}-T_{j-1}}{2r_j\Delta r}+\frac{(T_{j-1}-2T_j+T_{j+1}}{(\Delta r)^2}\right]+Q=
\kappa \frac{T_{j+1}(1+1/(2j))-2T_j+T_{j-1}(1-1/(2j))}{(\Delta r)^2}+Qwhere, in the second term of the rhs, we have made use of the relationship that rj=jΔrrj=jΔr.

I can now see that you have an error in Eqn. 4 of your write-up. There should be a 2 in the denominator of the 1st- and the 4th terms in brackets. Do you see that? So the value of your parameter B should be ##B=\frac{\kappa \delta t}{4\delta r}##, not ##B=\frac{\kappa \delta t}{2\delta r}##.

I will continue after we get agreement on this.

Last edited:
Are you labeling the grid point at r = 0 as j = 1 and the grid point at r = R as j = N, so that there are (N - 1) grid intervals and ##\Delta r=R/(N-1)##?

I start from $j=1$ which is in-line with the MATLAB convention. So I do
Matlab:
r=linspace(0,R,N)
, this makes $r_{1}=0$.

I found some sign errors in your algebra in going from Eqn. 12 to Eqn. 14. Please check.

Using the notation before my stencil is:
$$(\theta\delta_{jk} -\alpha A_{j,k})T_{i+1,k}=(\alpha A_{j,k}+\theta\delta_{jk})u_{i,k}+\frac{\delta t}{2}(Q_{i+1}+Q_{i})$$

The top row of $\mathbf{A}$ is $(-4,4,0,\ldots,0)$, the bottom row is $(0,\dots,0,2,(h\delta r/\kappa-2))$

That's the code in a nutshell I believe.

hunt_mat said:
Using the notation before my stencil is:
$$(\theta\delta_{jk} -\alpha A_{j,k})T_{i+1,k}=(\alpha A_{j,k}+\theta\delta_{jk})u_{i,k}+\frac{\delta t}{2}(Q_{i+1}+Q_{i})$$

The top row of $\mathbf{A}$ is $(-4,4,0,\ldots,0)$, the bottom row is $(0,\dots,0,2,(h\delta r/\kappa-2))$

That's the code in a nutshell I believe.
I concur with the top row if alpha is what you previously called A in your write-up.

And what is your comment regarding my assertion about your parameter B being high by a factor of 2?

Last edited:
It looks like you have the A and B parameter definitions switched in your code.

This is the code I'm working with now:
Matlab:
%This solves the heat equation with source term in cylindrical
%co-ordinates.

%T_t=(k/r)*(rT_r)_r+Q(t,r)
N_t=1000; %Time steps

t=linspace(0,3,N_t)'; %Time co-ordinate
dr=r(2);dt=t(2); %Increments

%Solution matrix
T=zeros(N_t,N_r);

%Physical characteristics
k=1; %Thermal conductivity
h=1; %Thermal exchange constant
theta=1;
alpha=k*dt/(2*dr*2);  %Useful constants
Y_1=sech(t)*(r.^2.*(1-r).^2);
Y_2=-2*tanh(t)*(8*r.^2-9*r+2);
Q=Y_1+Y_2;
a_1=-2*ones(1,N_r);
a_2=1+dr./r(1:N_r-1);
a_3=1-dr./r(2:N_r);
A=diag(a_1)+diag(a_2,1)+diag(a_3,-1);
A(1,1)=-4;A(1,2)=4;
A(N_r,N_r-1)=2;A(N_r,N_r)=h*dr/kappa-2;
A_plus=theta*eye(N_r)-alpha*A;
A_minus=theta*eye(N_r)+alpha*A;

%Constructing the solution:
b=zeros(1,N_r);
analytical=tanh(t)*(r.^2.*(1-r).^2);
for i=2:N_t
b=0.5*(Q(i,:)'+Q(i-1,:)')*dt;
u=T(i-1,:)';
C=A_minus*T(i-1,:)'+b;
x=A_plus\C;
T(i,:)=x';
end

for i=1:N_t
plot(r,T(i,:));
hold on
plot(r,analytical(i,:));
pause(0.1);
hold off
end

hunt_mat said:
This is the code I'm working with now:
Matlab:
%This solves the heat equation with source term in cylindrical
%co-ordinates.

%T_t=(k/r)*(rT_r)_r+Q(t,r)
N_t=1000; %Time steps

t=linspace(0,3,N_t)'; %Time co-ordinate
dr=r(2);dt=t(2); %Increments

%Solution matrix
T=zeros(N_t,N_r);

%Physical characteristics
k=1; %Thermal conductivity
h=1; %Thermal exchange constant
theta=1;
alpha=k*dt/(2*dr*2);  %Useful constants
Y_1=sech(t)*(r.^2.*(1-r).^2);
Y_2=-2*tanh(t)*(8*r.^2-9*r+2);
Q=Y_1+Y_2;
a_1=-2*ones(1,N_r);
a_2=1+dr./r(1:N_r-1);
a_3=1-dr./r(2:N_r);
A=diag(a_1)+diag(a_2,1)+diag(a_3,-1);
A(1,1)=-4;A(1,2)=4;
A(N_r,N_r-1)=2;A(N_r,N_r)=h*dr/kappa-2;
A_plus=theta*eye(N_r)-alpha*A;
A_minus=theta*eye(N_r)+alpha*A;

%Constructing the solution:
b=zeros(1,N_r);
analytical=tanh(t)*(r.^2.*(1-r).^2);
for i=2:N_t
b=0.5*(Q(i,:)'+Q(i-1,:)')*dt;
u=T(i-1,:)';
C=A_minus*T(i-1,:)'+b;
x=A_plus\C;
T(i,:)=x';
end

for i=1:N_t
plot(r,T(i,:));
hold on
plot(r,analytical(i,:));
pause(0.1);
hold off
end
Y'know, it would have been very considerate of my valuable time to, from the onset, have given me a correct version of the writeup and a corresponding latest version of your code. I've wasted a lot of time looking over your flawed versions.

That was my new code which I wrote yesterday.

Thank you for your time. I will use pdepe, as it seems that doing anything in cylindrical co-ordinates is a complete nightmare which I cannot get to work.

These are my calculations for the method. I haven't made any mistakes and it simply doesn't work.

Attachments

• MOL_fd.pdf
86.8 KB · Views: 311
hunt_mat said:
These are my calculations for the method. I haven't made any mistakes and it simply doesn't work.
If you didn't make any mistakes, it would work. Try a calculation in which you use only 21 grid points and choose delta t such that it makes alpha close to 1.

I have some ideas on how you can proceed to get your CN program to work properly.

Phase 1: To test out your spatial finite difference scheme, solve the problem using forward Euler (rather than CN ) with a small enough step size that the solution remains stable.

Phase 2: Use the CN equation, but solve the matrix using predictor-corrector (successive substitution) at each time step, rather than matrix inversion. The predictor for the first time step would be forward Euler. The corrector equation would involve only the diagonal elements, while the other values at the end of the time interval would be those from the previous iteration. Iterate at each time step until the corrector has converged.

I also recommend for debugging purposes that you do all these tests with only 3 spatial grid points. The CN solution at each time step would then require the solution of only 3 simultaneous linear equations, which could be done by hand (to check that the matrix was set up and solved properly).

1. What is the Crank-Nicholson method for solving the cylindrical heat equation?

The Crank-Nicholson method is a numerical method used to solve partial differential equations, specifically the cylindrical heat equation. It is a combination of the implicit and explicit methods, making it more accurate and stable than either method alone.

2. How does the Crank-Nicholson method differ from the explicit method?

The explicit method for solving the cylindrical heat equation uses a forward difference approximation for the time derivative, while the Crank-Nicholson method uses a combination of forward and backward difference approximations. This makes the Crank-Nicholson method more accurate and less prone to instability.

3. What is the stability condition for the Crank-Nicholson method?

The Crank-Nicholson method has a stability condition of $\Delta t \leq \frac{\Delta r^2}{2\alpha}$, where $\Delta t$ is the time step, $\Delta r$ is the spatial step, and $\alpha$ is the thermal diffusivity. This means that the time step must be small enough to ensure stability and accuracy of the solution.

4. Can the Crank-Nicholson method be used for other types of equations?

Yes, the Crank-Nicholson method can be used to solve other types of partial differential equations, such as the diffusion equation and the wave equation. It is a versatile method that can be applied to a variety of problems in physics and engineering.

5. What are the advantages of using the Crank-Nicholson method over other numerical methods?

The Crank-Nicholson method is more accurate and stable than other numerical methods, such as the explicit method. It also has a higher order of convergence, meaning that the error decreases faster as the grid is refined. Additionally, it can handle complex boundary conditions and can be easily adapted to solve problems in higher dimensions.

• Differential Equations
Replies
2
Views
2K
• MATLAB, Maple, Mathematica, LaTeX
Replies
41
Views
8K
• Differential Equations
Replies
3
Views
2K
• Differential Equations
Replies
6
Views
2K
• Differential Equations
Replies
3
Views
1K
• Differential Equations
Replies
7
Views
2K
• Differential Equations
Replies
4
Views
2K
• Calculus and Beyond Homework Help
Replies
0
Views
475
• Set Theory, Logic, Probability, Statistics
Replies
62
Views
3K
• Differential Equations
Replies
8
Views
4K