1

I am implementing a program which solves differential equations - 1d diffusion.

I am using the Crank - Nicolson method whose accuracy should be second-order for time and second-order for space.

Unfortunately my results are second-order for time and first-order for space. How is that possible? Or did I mess up something and is it not possible?

 I know f.e. that if it should be second-order and becomes third-order that means that T = Ahp + Bhp+1 and etc.
if Ahp equals 0 then Bhp+1 will become dominant and we have p+1 order of accuracy.

I found one small issue in my code so I am closing the question. Thanks @PierreCarre and @Uranix for your help

  • Could you please elaborate on how you study the convergence? Do you change the spatial step while the time step being fixed, or you change them both at once? – uranix May 24 '22 at 07:05
  • Do you approximate the boundary conditions with $O(h^2)$? It is tricky when they contain spatial derivatives (Neumann) – uranix May 24 '22 at 07:06
  • spatial step is divided by 2 and time is always equal h*h. – GenLancelot May 24 '22 at 09:56
  • In my case i am elaborating convergence using log10(h) vs log10(maxerror in max time) – GenLancelot May 24 '22 at 09:58
  • In my case boundary conditions are dirichlet boundary conditions – GenLancelot May 24 '22 at 10:01
  • @GenLancelot Can you post the exact problem you are solving, initial and boundary conditions included? – PierreCarre May 24 '22 at 12:56
  • Problem :\frac{\partial U(x,t)}{\partial t} = D[\frac{\partial^2 U(x,t)}{\partial x^2} + \pi^2sin(\pi x)))] Edit: i edited post with img cuz latex is not working idk why – GenLancelot May 24 '22 at 21:44

2 Answers2

1

Crank-Nicolson method is a second order in time and space, that is its error is $\varepsilon = A \Delta t^2 + B \Delta x^2 + \text{higer order terms}$. The leading term is $A \Delta t^2 + B \Delta x^2$.

If you take $\Delta t = \Delta x^2$ then the leading term becomes $A\Delta x^4 + B\Delta x^2$ and $A\Delta x^4$ becomes one of the higher order terms. The main error is given simply by $B\Delta x^2$. I also cannot understand why do you think your results are second-order in time. As I said, when $\Delta t = \Delta x^2$ the temporal error is very small and can be neglected compared to the spatial error.

The proper convergence study for the Crank-Nicolson method requires fixing $\Delta t = K \Delta x$ relation for some constant $K$. Then the leading term of the error decreases exactly in 4 times when you divide the step size by 2. Note that you need to divide both of the steps (the temporal and the spatial) by the factor of 2.

uranix
  • 7,503
1

The CN scheme results in solving at each time step the linear system $$ \left(I-\frac{D \Delta t}{2\Delta x^2} A\right)U^{n+1}=\left(I+\frac{D \Delta t}{2\Delta x^2} A\right)U^n + \Delta t F, $$

where $A$ is the usual central difference matrix and $F$ the source term $ f_i = D \pi^2 \sin(\pi x_i)$. Taking $\Delta t = \Delta x$, which is fine because the scheme is unconditionally stable, you get the table ($D=1$): $$ \begin{array}{|c|c|}\hline \Delta x, \Delta t & Max. Error\\ \hline 0.1 & 0.286 \times 10^{-1}\\ \hline 0.01 & 0.316 \times 10^{-3}\\ \hline 0.001 & 0.321 \times 10^{-5}\\ \hline \end{array} $$

where the error is measured as $\displaystyle \max_{i,j}|U^i_j - u(x_j, t_i)|$. Using the last two lines to estimate the order of convergence, you get $$ \alpha \approx \dfrac{\log \left(\dfrac{0.321\times 10^{-5}}{0.32\times 10^{-3}}\right)}{\log\left(\dfrac{0.001}{0.01}\right)}=2.011 $$

So, it would be useful if you post your exact discretization and code, because there must be some mistake.

PierreCarre
  • 20,974
  • 1
  • 18
  • 34