2
$\begingroup$

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

enter image description here

$\endgroup$
11
  • 1
    $\begingroup$ Can you add the plots of the numerical solution and the exact solution? $\endgroup$
    – Anton Menshov
    Apr 8 '20 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 '20 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 '20 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 '20 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 '20 at 9:57
2
$\begingroup$

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

numerical solution over reference solution

$\endgroup$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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