4
$\begingroup$

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.

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ 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/… $\endgroup$
    – boyfarrell
    Commented Dec 3, 2013 at 0:36
  • $\begingroup$ 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.. $\endgroup$
    – Daniel
    Commented Dec 3, 2013 at 8:41

1 Answer 1

Reset to default
6
$\begingroup$

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

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – Nat Wilson
    Commented Dec 3, 2013 at 4:14
  • $\begingroup$ 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/… $\endgroup$
    – Daniel
    Commented Dec 3, 2013 at 8:38
  • 1
    $\begingroup$ 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. $\endgroup$
    – boyfarrell
    Commented Dec 3, 2013 at 9:51
  • $\begingroup$ 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. $\endgroup$
    – Daniel
    Commented Dec 3, 2013 at 12:25
  • $\begingroup$ @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. $\endgroup$
    – Jed Brown
    Commented Dec 3, 2013 at 14:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.