# Crank-Nicholson scheme for transport equation

This is my attempt to find the approximate solution of the folowing transport equation $$\left\{\begin{array}{ll} \partial_{t} u+\partial_{x} u= (x^2-x)t+x^3/3-x^2/2, & t \in(0,0.4), x \in(0,1) \\ u(0, x)=2, & x \in(0, 1) \\ \frac{\partial u}{\partial n}=0 & t \in(0,0.4), x=0\quad or \quad1 \end{array}\right.$$ by using the following Crank-Nicholson scheme $$U_{i}^{j+1}+ \frac{k}{4h}({U_{i+1}^{j+1}-U_{i-1}^{j+1})=U_{i}^{j}- \frac{k}{4h}(U_{i+1}^{j}-U_{i-1}^{j}})+\frac{k}{2}\left[f_{i}^{j+1}+f_{i}^{j}\right]$$ it can work but it's too far from the exact solution $$u(t,x)=(x^3/3 -x^2/2)t + 2$$ I tried to found why but it turned out to be very difficult.

%% Left-hand side A*U

A=diag(-cte/4*ones(N-2,1),0)+diag(ones(N-3,1),1)+diag(cte/4*ones(N-4,1),2);
A = A([1:N-4],:);
A = sparse(A);

%% Fill in the initial condition

for i=1:N-2
uCN1(i)= g(a+(i+1)*h);
end

uCN11= [uCN1(1); uCN1 ; uCN1(N-2)];  % Neumann boundary

uCN2=uCN1;
uCN22=uCN11;

%% uCN2
F=zeros(N-4,1);

%% Left-hand side C*U+F
C=diag(cte/4*ones(N-2,1),0)+diag(ones(N-3,1),1)+diag(-cte/4*ones(N-4,1),2);
C = C([1:N-4],:);
C = sparse(C);

for j=1:M-1

for i=1:N-4
F(i)=k*(f(S(i+2),T(j))+f(S(i+2),T(j+1)))/2;
end

F1=F;

B1=C*uCN2+F1;  % Rght-hand side

% Solve for linear equation
uCN2=A\B1;

uCN22= [uCN2(1); uCN2 ; uCN2(N-2)];

uCN=uCN2;

uCN=uCN22;   % approximate solution
plot(S,uCN,'c');
R=j*k;
title([' U at at T=' num2str(R)])
pause(k)
end


EDIT

for the exact solution, I wrote this code:

for i=1:N
uEX(i)=(S(i)^3/3 -S(i)^2/2)*(0.4) +2; % I forgot +2
end


• Can you add the plots of the numerical solution and the exact solution? Apr 8 '20 at 10:39
• This now meets all the points I mentioned on the SO post, even without pictures (but pictures, moderately sized, are always nice). Apr 8 '20 at 10:43
• I see a discrepancy between the initial condition $u(0,x)=2$ and the solution formula giving $u(0,x)=0$. If you add the constant to get $u(t,x)=2+t(x^3/3-x^2/2)$, this problem would vanish. Is this the difference you observe? Apr 8 '20 at 10:47
• Thank you again, things always become nicer with your directions :))), sorry it was my mistake, I didn't add two in the expression of the exact solution Apr 8 '20 at 12:18
• It defeats somewhat the purpose if you have an order 2 method but only implement the boundary conditions with order 1. This reduces the order of the method to 1 in space direction, to $O(h+k^2)$, but this is no reason for the observed divergence. With those explanations I'll look again at your code if I find some error there. Apr 9 '20 at 9:57

Your problem seems to be the implementation of the boundary conditions. Apart from that you seem to not compute existing variables efficiently.

If I understand your initialization right, then $$x_2=a+2h$$ to $$x_{N-1}=a+(N-1)h$$ is the middle segment of the space discretization, so that $$x_0=a=0$$ implies $$x_{N+1}=b=1$$ and $$h=\frac1{N+1}$$. Shifting to index-one based arrays means that the state vector in space direction is U(1:N+2).

## Analysis of your first order implementation

Your setup is that S=a:h:b and N=lenght(S). The Neumann boundary conditions are implemented as $$u(t,x_0)=u(t,x_1)$$ and $$u(t,x_{N-1})=u(t,x_{N-2})$$. Then the state vector at time $$t_j$$ is $$u^j_i=u(t_j,x_i)$$, $$i=1,...,N-2$$, the initialization would have to be

U(1,:) = g(S(2:N-1));


(or via loop, use S(i+1) for the x value). In the C-N equation, the difference part on the right side acts on

[ U(j,1) U(j,:) U(j,N-2) ]


either via matrix multiplication or more simply due to the Toeplitz structure, via convolution. The function part can likewise constructed from the vectorized evaluations f(T(j),S(2:N-1)). (One might think about introducing Si=S(2:N-1) for the inner part of this array). Then the system for the left side is $$\begin{bmatrix} -1-α&α\\ -α&1&α\\ &\ddots&&\ddots\\ &&-α&1&α\\ &&&-α&1+α\\ \end{bmatrix} \begin{bmatrix} u^{j+1}_1\\u^{j+1}_2\\\vdots\\u^{j+1}_{N-3}\\u^{j+1}_{N-2} \end{bmatrix} = \begin{bmatrix} r_1\\r_2\\\vdots\\r_{N-3}\\r_{N-2} \end{bmatrix}$$ This matrix has dimension $$N-2$$, I do not understand how you can reduce it to size $$N-4$$ without accounting for the first and last row and their corresponding equations. You compute something for the inside of the inside, and have it only loosely connected to the boundary.

## Second order implementation

I prefer the convention where a subdivision has $$N$$ segments and thus $$N+1$$ node/sampling points $$x_i=S(i+1)$$.

S = linspace(a,b,N+1); h=S(2)-S(1);


The easiest way to implement the boundary conditions to second order is to use ghost cells. This means that one observes $$u(t,-h)\simeq u(t,+h)$$ and the same at the other boundary. In the implementation this means that the extended array standing for the extended node sequence $$x_{-1},x_0,...,x_{N},x_{N+1}$$ or (virtually) S(0:N+2) should be

[ U(2)  U  U(N) ]


The difference operator acting on this array could be realized as

conv ([ U(2)  U  U(N) ], [ -alf, 1, alf ], shape="valid")


as convolution applies the second factor in reversed order. shape="valid" removes the outer two elements on each side and thus restores the shape of U in the result.

To solve for the next step the system to be solved, using $$u^j_1-u^j_{-1}=0$$ and $$u^j_{N+1}-u^j_{N-1}=0$$, is $$\begin{bmatrix} 1&0\\ -α&1&α\\ &\ddots&&\ddots\\ &&-α&1&α\\ &&&0&1 \end{bmatrix} \begin{bmatrix} u_0\\u_1\\\vdots\\u_{N-1}\\u_{N} \end{bmatrix} = \begin{bmatrix} r_0\\r_1\\\vdots\\r_{N-1}\\r_{N} \end{bmatrix}$$ where $$r_i$$ is the right side of the C-N equation, that is, everything not depending on the $$u^{j+1}_i$$.

This can be shortened by reading off $$u_0=r_0$$ and $$u_{N}=r_{N}$$, $$\begin{bmatrix} 1&α\\ -α&1&α\\ &\ddots&&\ddots\\ &&-α&1&α\\ &&&-α&1\\ \end{bmatrix} \begin{bmatrix} u_1\\u_2\\\vdots\\u_{N-2}\\u_{N-1} \end{bmatrix} = \begin{bmatrix} r_1\\r_2\\\vdots\\r_{N-2}\\r_{N-1} \end{bmatrix} +α \begin{bmatrix} r_0\\ \\\vdots\\ \\-r_{N} \end{bmatrix}$$

In GNU octave the following works

function scicomp34793_crank_nicholson_4_transport
clear all; clf;
N = 10;
M = 20;
S = linspace(0,1,N+1)
h = S(2)-S(1);
T = linspace(0, 0.4, M+1)
k = T(2)-T(1);

v = 1;
f=@(t,x) (x.^2-x).*t+x.^3/3-x.^2/2;
g=@(x) 2+0*x;

ref=@(t,x) 2+t.*(x.^3/3-x.^2/2);

alf = v*k/(4*h);

U(1,:) = g(S);

C = spdiags([ -alf*ones(N-1,1), ones(N-1,1), alf*ones(N-1,1)], [-1,0,1], N-1,N-1);

for j=1:M
Ujdiff = conv([U(j,2)  U(j,:)  U(j,N)],[-alf,1,alf], shape="valid");
R = 0.5*k*(f(T(j),S)+f(T(j+1),S)) + Ujdiff;
R(2)  += alf*R(1);
R(N) -= alf*R(N+1);
U(j+1,:) = [ R(1) R(2:N)/C' R(N+1) ];
end%for
X = linspace(0,1,301);
for j=1:M+1
plot(X,ref(T(j),X),'b', 'LineWidth',6);
ylim([1.9,2.03]);
hold on;
plot(S,U(j,:),'-oy', 'LineWidth',2);

title([' U at at T=' num2str(T(j))])
hold off;
pause(k*20)
end%for

end%function


numerical solution over reference solution