Error in Crank-Nicolson scheme for diffusion equation

finite difference methodsfinite differencesnumerical methodsparabolic pdepartial differential equations

I'm solving the diffusion equation

$$\frac{\partial u}{\partial t}=D\frac{\partial^2u}{\partial^2x}$$

subject to the BCs $\partial_xu(x=0)=0$ and $u(x)=1$, using the Crank Nicolson scheme. For the middle points I'm using

$$\frac{u_i^{n+1}-u_i^n}{\Delta t}=\frac{1}{2}D\frac{u_{i+1}^n-2u_i^n+u_{i+1}^n}{\Delta x^2}+\frac{1}{2}D\frac{u_{i+1}^{n+1}-2u_i^{n+1}+u_{i+1}^{n+1}}{\Delta x^2},$$

whereas in the left point, I make use of the fact that $\frac{u_1-u_{-1}}{2\Delta x}=0\implies u_1=u_{-1}$ and thus write

$$\frac{u_0^{n+1}-u_0^n}{\Delta t}=\frac{1}{2}D\frac{2(u_{1}^n-u_0^n)}{\Delta x^2}+\frac{1}{2}D\frac{2(u_{1}^{n+1}-u_0^{n+1})}{\Delta x^2}.$$

When an initial condition $u(0,x)$ is given, I can also find the analytical solution of this equation and thus compare the simulation with theory.

I am evaluating the absolute error $||\mathbf{T_{simul}}-\mathbf{T_{anal}}||$ at each time $t$.

When $t\rightarrow\infty$, the solution goes to zero with machine precision, which makes sense.

Now let's say I pick up a time $t$, when heat is still diffusing, and evaluate the error, and then I make the same thing for a different time-step $\Delta t$.

The resultant log-plot of the error with $\Delta t$ is the following:

The points fit very well with a linear function, but because this is a second-order scheme in $\Delta t$, I was expecting that the error would go with $(\Delta t)^2$ and not $\Delta t$!

EDIT: As suggested in the comments, the log-log plot is:

(This fit is made to a function of the type $a+b(\Delta t)^2$. I have been told however that the fit should not contain the constant term.)

FINAL EDIT:

With the help of @LutzL I have realized I was making one mistake in my code, namely I was comparing theory and simulation at times separated precisely by $\Delta t$.

Besides that, I have found the following very interesting results:

1)
If I use as initial condition

$$u(x,0)=1.0+\cos\left(\frac{\pi}{2}x\right),$$

which gives the exact solution (with $D=1$)

$$u(x,t)=1.0+e^{-\pi^2t/4}\cos\left(\frac{\pi}{2}x\right)$$

in the domain $x\in[0,1]$, I get that the error indeed goes with $(\Delta t)^2$:

The blue line gives the fit $4.106(\Delta t)^2$

2)
Now if I use as initial condition

$$u(x,0)=1.0+\cos\left(\frac{x}{2}\right),$$

which gives the exact solution (with $D=1$)

$$u(x,t)=1.0+e^{-t/4}\cos\left(\frac{x}{2}\right)$$

in the domain $x\in[0,\pi]$, I get the following trend for the error:

where the blue line is the fit $0.00135462 + 0.00283811x + 0.0239431x^2$

It is strange that the simple change of variable $\pi x\rightarrow x$ makes this behaviour. Note also that in both cases 1) and 2) I have kept $\Delta r=0.001$ fixed.

Best Answer

If you did everything right the matrix should have been constructed as follows or similar. The dynamic is in the inner points. Thus the approximation of $F$ gets the vector $(u_1,...,u_N)$ as input. The right boundary condition gives $u_{N+1}=1$. To reconstruct $u_0$ use the forward differentiation formula $$u_x(x=0)=\frac{3u_0-4u_1+u_2}{Δx}+O(Δx^2)$$ so that $u_0=\frac13(4u_1-u_2)$. With that the second derivative at the first position is approximated by $$ (u_{xx})_1=\frac{\frac13(4u_1-u_2)-2u_1+u_2}{Δx^2}=\frac{-2u_1+2u_2}{3Δx^2} $$ This gives the matrix for $F(u)=Au+b$ as $$ A=D\frac{Δt}{Δx^2} \pmatrix{ -\frac23&\frac23&0&\dots\\ 1&-2&1&0&\dots\\ &\ddots&\ddots&\ddots\\ \dots&0&1&-2&1\\ &\dots&0&1&-2 }, ~~~ b=D\frac{Δt}{Δx^2}\pmatrix{0\\0\\\vdots\\0\\1}. $$

Building a test case

Continuous problem with exact solution

$u_t=u_{xx}$, that is, $D=1$, with $u_x(0,t)=0$, $u(\pi,t)=1$ over $t\in[0,1]$.

Separation of variables tells us that terms $e^{-\omega^2t}\cos(\omega t)$ are simple functions that satisfy homogeneous condition automatically, using a constant for the shift gives as solution example $$u(x,t)=1+a_1e^{-t/4}\cos(x/2)+a_3e^{-9t/4}\cos(3x/2).$$

Code snippets (python, using scipy.sparse)

the test solution

a1, a3 = 2,1;
def u_exact(x,t): return 1+a1*np.exp(-0.25*t)*np.cos(0.5*x)+a3*np.exp(-2.25*t)*np.cos(1.5*x)

construction of the system matrices

I = sparse.eye( M, format="csc" )
A = sparse.diags([1, -2, 1], [-1, 0, 1], shape=(M, M), format="csc");
A[0,0]=-2.; A[0,1]=2.;
b = np.zeros(M); b[-1]=1

the time iteration

x = np.linspace(0,np.pi,M+1); dx = x[1]-x[0];
u1 = u_exact(x,0)[:-1];
t = np.linspace(0,T,N+1); dt = t[1]-t[0];
p = dt/dx**2;
for n in range(N):
    u1 = sparse.linalg.spsolve(I-0.5*p*A, (I+0.5*p*A).dot(u1)+p*b)

Numerical results

Using the left side BC treatment like in the question, I get the following table of results. Horizontally $M$ varies over the size of the space subdivision, while vertically $N$ gives the number of time steps.

In each cell the upper value is the error against a run with half the time step, which is an estimate of the error against exact time integration. The lower value is the error against the exact solution. The errors displayed is an approximation of the $L^2$ norm, that is, the Euclidean norm divided by $\sqrt M$. One can see that where $\frac{Δt}{Δx^2}$ gets small, the error against the exact solution becomes constant.

\begin{array}{|c|l|l|l|l|l|}\hline & ~M=10& ~M=50& ~M=200& ~M=1000& ~M=5000 \\ \hline N=10 & \begin{split} 5.5166e-04 \\ 2.6801e-03 \end{split} & \begin{split} 5.3780e-04 \\ 5.9059e-04 \end{split} & \begin{split} 5.3395e-04 \\ 7.0361e-04 \end{split} & \begin{split} 5.3286e-04 \\ 7.0966e-04 \end{split} & \begin{split} 5.3264e-04 \\ 7.0967e-04 \end{split} \\ \hline N=20 & \begin{split} 1.3755e-04 \\ 3.2229e-03 \end{split} & \begin{split} 1.3409e-04 \\ 5.6237e-05 \end{split} & \begin{split} 1.3313e-04 \\ 1.6967e-04 \end{split} & \begin{split} 1.3286e-04 \\ 1.7680e-04 \end{split} & \begin{split} 1.3280e-04 \\ 1.7703e-04 \end{split} \\ \hline N=40 & \begin{split} 3.4365e-05 \\ 3.3587e-03 \end{split} & \begin{split} 3.3500e-05 \\ 8.4139e-05 \end{split} & \begin{split} 3.3260e-05 \\ 3.6555e-05 \end{split} & \begin{split} 3.3193e-05 \\ 4.3943e-05 \end{split} & \begin{split} 3.3179e-05 \\ 4.4224e-05 \end{split} \\ \hline N=80 & \begin{split} 8.5898e-06 \\ 3.3927e-03 \end{split} & \begin{split} 8.3736e-06 \\ 1.1693e-04 \end{split} & \begin{split} 8.3137e-06 \\ 3.5075e-06 \end{split} & \begin{split} 8.2968e-06 \\ 1.0751e-05 \end{split} & \begin{split} 8.2933e-06 \\ 1.1045e-05 \end{split} \\ \hline N=160 & \begin{split} 2.1474e-06 \\ 3.4011e-03 \end{split} & \begin{split} 2.0933e-06 \\ 1.2518e-04 \end{split} & \begin{split} 2.0783e-06 \\ 5.1960e-06 \end{split} & \begin{split} 2.0741e-06 \\ 2.4544e-06 \end{split} & \begin{split} 2.0733e-06 \\ 2.7519e-06 \end{split} \\ \hline N=320 & \begin{split} 5.3683e-07 \\ 3.4033e-03 \end{split} & \begin{split} 5.2333e-07 \\ 1.2725e-04 \end{split} & \begin{split} 5.1958e-07 \\ 7.2298e-06 \end{split} & \begin{split} 5.1852e-07 \\ 3.8307e-07 \end{split} & \begin{split} 5.1832e-07 \\ 6.7867e-07 \end{split} \\ \hline \end{array}

Speculation

Without knowing details of the code used, there are two typical places that come to mind that might introduce a first order error:

a wrong coefficient somewhere, reducing the method to first order. This is unlikely as the matrix and function in general has a simple structure. If there were an error one would expect more grave consequences.
an error in the implementation of the time iteration. Such as
- iterating N times but setting dt=T/(N-1), or
- having T and dt as input and iterating while t<T without caring that the final t can be anywhere between T and T+dt.
Such error results in a time difference between exact and numerical solution of order O(dt) which gives a difference in the function values of also O(dt).

Related Solutions

LeapFrog scheme for the Advection equation

You need to cycle all of the arrays in the loop, at the moment you leave unp0 with the original value of u0. Also, as usual, avoid loops that you can express as vector operations, the execution will be faster.

 for n=1:tsteps
     u2 = u0;
     u2(2:N+2) -=  (dt/dx)*(u1(3:N+3)-u1(1:N+1))
     u2(1)=u2(3);
     u2(N+3)=u2(N+1);

    t=t+dt;
    u0=u1; u1=u2;
    exact=F(x-t);
    plot(x,u1,'bo-',x,exact,'r-')

 end

Crank-Nicolson Scheme Not giving Accuracy of Order 2

You are computing the Euclidean norm of the vector that contains the pointwise errors. Each component is essentially $O(h^2)$ where $h$ is the spatial stepsize. In your minimal working example the time step has been reduced to the point where it makes a neglible contribution to the error. Therefore your error behaves as $$E(m,h)=\sqrt{m}h^2$$ where $m$ is the length of the vector and $h$ and $mh$ is a constant, the length of your interval. Your order is computed as the ratio $$E(2m,h/2)/E(m,h) = \frac{\sqrt{2m} \frac{1}{4} h^2}{\sqrt{m} h^2} = \sqrt{2} \frac{1}{4} \approx 0.35.$$ You are getting exactly the results that you should.

Your objective is to compute an approximation of a function $$x \rightarrow f(T,x)$$ for some fixed value $t=T$ of the time. You have approximations $a_j$ of $f(T,x_j)$ for equidistant $x_j$, i.e., $x_j = x_0 + jh$. The global error is $$e_j = f(T,x_j) - a_j = O(k^2 + h^2).$$ You have reduced the temporal stepsize $k$ to the point where it is irrelevant. I would like you to compute either $$E_1 = h \sum_{j=0}^{m-1} |e_j|$$ or $$E_2 = \sqrt{h \sum_{j=0}^{m-1} e_j^2}.$$ Both expressions are of class $O(h^2)$ because $mh = x_R - x_L$ is a constant. [Note: You should be using $m$ to compute $h$, rather than the other way around using the ceiling function. It is not a significant issue here, it will show up in a test based on Richardson's techniques.] Your original/new measure is $$\tilde{E} = \sqrt{\sum_{j=0}^{m-1} e_j^2} = O(\sqrt{mh^4}) = O(\sqrt{m}h^2) = O(h^{\frac{3}{2}}).$$