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 = A([1:N-4],:);
A = sparse(A);

%% Fill in the initial condition 

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

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


%% uCN2

%% Left-hand side C*U+F
C = C([1:N-4],:);
C = sparse(C);

for j=1:M-1

    for i=1:N-4


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

    % Solve for linear equation

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


    uCN=uCN22;   % approximate solution
    title([' U at at T=' num2str(R)])


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

enter image description here

  • 1
    $\begingroup$ Can you add the plots of the numerical solution and the exact solution? $\endgroup$
    – Anton Menshov
    Apr 8, 2020 at 10:39
  • 1
    $\begingroup$ This now meets all the points I mentioned on the SO post, even without pictures (but pictures, moderately sized, are always nice). $\endgroup$ Apr 8, 2020 at 10:43
  • 1
    $\begingroup$ 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? $\endgroup$ Apr 8, 2020 at 10:47
  • $\begingroup$ 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 $\endgroup$ Apr 8, 2020 at 12:18
  • 1
    $\begingroup$ 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. $\endgroup$ Apr 9, 2020 at 9:57

1 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) ];
  X = linspace(0,1,301);
  for j=1:M+1
    plot(X,ref(T(j),X),'b', 'LineWidth',6);
    hold on;
    plot(S,U(j,:),'-oy', 'LineWidth',2);

    title([' U at at T=' num2str(T(j))])
    hold off;


numerical solution over reference solution

numerical solution over reference solution


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.