# My algorithm for the heat equation is unstable

I have implemented the 2D heat equation with what I thought was the Crank-Nicolson algorithm in the following way:

function[newDomain] = heatCN(domain, dt, Nx, Ny, hx, hy, spanX, spanY)
aux = domain;
newDomain = domain;
for i=2:Nx-1
for j=2:Ny-1
derX1 = dt * (-2.0 * domain(i,j) + domain(i+1, j) + domain(i-1,j)) / hx^2;
derY1 = dt * (-2.0 * domain(i,j) + domain(i, j+1) + domain(i, j-1)) / hy^2;
aux(i,j) = derX1 + derY1 + domain(i,j);
end;
end;

for i=2:Nx-1
for j=2:Ny-1
derX2 = dt * (-2.0 * domain(i,j) + domain(i+1, j) + domain(i-1,j)) / hx^2;
derY2 = dt * (-2.0 * domain(i,j) + domain(i, j+1) + domain(i, j-1)) / hy^2;
derXT = dt * (-2.0 * aux(i,j) + aux(i+1, j) + aux(i-1,j)) / hx^2;
derYT = dt * (-2.0 * aux(i,j) + aux(i, j+1) + aux(i, j-1)) / hy^2;
newDomain(i,j) = domain(i,j) + 0.5 * (derX2 + derY2 + derXT + derYT);
end;
end;
end;


From my knowledge, the CN method should be unconditionally stable, meaning that for any (decent) time step it should give me the right output. But the behavior I see with this program is unstable.

• My advice would be to try solving this in 1D first. You should also state your boundary conditions and initial conditions. If you have initial conditions with very sharp features then this can causes oscillations which may cause divergence. Finally, you are more likely to get a good response when you post your method mathematically rather than just showing code, for example a similar 'debugging' question where the method was emphasised and not the code, scicomp.stackexchange.com/questions/5355/… Dec 3 '13 at 0:36
• Thank you, @boyfarrell. The method I am trying to implement is just the standard 2D-diffusion Crank-Nicolson (en.wikipedia.org/wiki/…) I am on the other hand not using a matrix to solve the system, but trying to do this on the grid, thus saving space. I am though not sure if the idea of calculating the next time step in aux and then using that to calculate the new answer is the right one.. Dec 3 '13 at 8:41

The algorithm you have implemented is explicit. Crank-Nicolson is an implicit method, and thus requires a solve.

• To expand slightly on Jed's answer, the Crank Nicolson method defines the subsequent timestep as a function of both the present and future function values. Assuming you have the right boundary conditions, it breaks down into a system of equations that can be solved algebraically. Dec 3 '13 at 4:14
• Mathematically speaking, I understand that the system has to be solved in order to get to the next step. I have implemented above a grid-wise approach to the problem, so I don't have to store the matrix which might get pretty big. So I am technically calculating the next step in aux and then using that and the current step to get to the next time step. My idea is from here: en.wikipedia.org/wiki/… Dec 3 '13 at 8:38
• I think the point that others are making is that you must solve a linear system if you want the unconditional stability of the Crank-Nicolson approach. You say that you understand this point (the difference between implicit and explicit approaches) so forgive my repetition if this is unwelcome! I just wanted to clarity. Why not try with matrices? This is after all what MATLAB does best. It can even do sparse matrices which save space and are computational efficient. If the performance if poor at least you have working reference from which you can optimise to a grid-wise algorithm. Dec 3 '13 at 9:51
• I think I may have not seen your point at first. I will try it with solving the system and will see how far that will get me. The reasons for wanting to work in a grid-wise approach is the fact that I am a C/C++ programmer, so from what I see grids are more memory and access time efficient, but this might be a completely different discussion. I will try it with the matrices. Dec 3 '13 at 12:25
• @Daniel No local process can perform the solve, since the operation (the matrix inverse) is dense. You can apply this dense operation fast (in $O(n)$ time) using a hierarchical method like multigrid. If you are programming in C or C++, you may want to check out a library like PETSc, which provides many efficient solvers. Dec 3 '13 at 14:54