1

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)$.

EDIT

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

EDIT 2

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

Kyle Kanos
  • 472
  • 3
  • 11
  • 3
    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? – Bill Barth Oct 04 '13 at 02:40
  • 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. – Kyle Kanos Oct 04 '13 at 02:43
  • 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. – Godric Seer Oct 04 '13 at 14:34
  • 1
    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. – Godric Seer Oct 04 '13 at 14:43
  • @GodricSeer: I removed that section because I had formulated it incorrectly. Fixing the boundaries does indeed give me better results. – Kyle Kanos Oct 04 '13 at 14:43

0 Answers0