Crank-Nicholson scheme for transport equation

Question

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

Answer 1

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