I modeled the heat equation, $$ u_t = au_{xx} $$ using the common 2nd-order Crank-Nicolson scheme, $$ \frac{u^{n+1}_i-u^{n}_i}{dt} = \frac{a}{2\,dx}\left(u_{i-1}^{n+1}+u_{i+1}^{n+1}-2u_i^{n+1} + u_{i-1}^n+u_{i+1}^n-2u_i^n\right) $$ and the standard 4th-order Crank-Nicolson scheme, $$ \frac{u^{n+1}_i-u^{n}_i}{dt} = \frac{a}{24\,dx}\left(-u_{i-2}^{n+1}+16u_{i-1}^{n+1}-30u_{i}^{n+1}+16u_{i+1}^{n+1}-u_{i+2}^{n+1}\right) \\ \quad+\frac{a}{24\,dx}\left(-u_{i-2}^{n}+16u_{i-1}^{n}-30u_{i}^{n}+16u_{i+1}^{n}-u_{i+2}^{n}\right) $$ Outside of figuring out how to solve the pentadiagonal matrix (without resorting to libraries like LAPACK), this was fairly straight-forward.

But my results were strange: the 2nd-order scheme was more accurate to the exact solution (in 3 different test cases!) than the 4th-order scheme. Is there any reason we should expect a 2nd-order accurate scheme to perform better than the 4th-order scheme?

If it matters, $a=1$, $x\in(0,1)$ and $u(0,t)=u(1,t)=0$ for all examples:

  1. $u(x,0)=6\sin(\pi x)$ leading to the analytic solution $u(x,t)=6\sin(\pi x)\exp(-\pi^2t)$
  2. $u(x,0)=12\sin(9\pi x) - 7\sin(4\pi x)$ leading to $u(x,t)=12\sin(9\pi x)\exp(-9\pi^2t) - 7\sin(4\pi x)\exp(-9\pi^2t)$
  3. $u(x,0)=10$ leading to $u(x,t)=\sum_{k=1}^{100}\frac{40}{(2k-1)\pi}\exp\left(-(2k-1)^2\pi^2t\right)\sin\left((2k-1)\pi x\right)$

With this last one, though, both did pretty poorly at the boundaries, but the 2nd-order method was closer to the exact in the range $x\in(0.1,0.9)$.


You can see that it's an order of magnitude difference between the two schemes.

two sine term enter image description here

single sine term enter image description here

constant 10 enter image description here


I ran this with different $N$'s and still found the 2nd order version to be more accurate than the 4th order.

  • 3
    $\begingroup$ The errors at the boundaries can pollute your solution pretty badly. Have you done a mesh refinement study? What does the error look like if you plot it? $\endgroup$
    – Bill Barth
    Oct 4, 2013 at 2:40
  • $\begingroup$ I'm not too worried about the errors in the last example because the function I am using does not have such a discontinuity. I have not looked at different mesh sizes, only the different $u(x,0)$'s. I'll post the error I get in a moment. $\endgroup$
    – Kyle Kanos
    Oct 4, 2013 at 2:43
  • $\begingroup$ Not just the last one. Even the 2nd one looks like you have pretty serious instabilities near the edges. If your mesh isn't fine enough, then the 4th order scheme may be unstable. As Bill Barth says, those errors can propogate into the body of your solution. $\endgroup$ Oct 4, 2013 at 14:34
  • 1
    $\begingroup$ In regards to your second edit, can you add the formulation at the boundaries? If you accidentally used a 1st order method for the boundary in the 4th order scheme, the whole thing becomes first order. $\endgroup$ Oct 4, 2013 at 14:43
  • $\begingroup$ @GodricSeer: I removed that section because I had formulated it incorrectly. Fixing the boundaries does indeed give me better results. $\endgroup$
    – Kyle Kanos
    Oct 4, 2013 at 14:43


Your Answer

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

Browse other questions tagged or ask your own question.