Solving transport equation with pseudo-spectral and finite volume methods in MATLAB

Question

I have a linear transport equation $$ \dot{c}(x; t) + vc_x(x; t) = 0 \tag{1}$$ on an interval $[0; 2π]$ with $v = 1$, periodic boundary conditions and two different initial values $$c(x; 0) = \sin(x)$$ and $$ c(x,0) = \left\lbrace \begin{aligned} &0:x≤ \tfrac13; \\ &1: \tfrac13 < x < \tfrac23;\\ &0: x ≥ \tfrac23 ; \end{aligned}\right. $$

In both cases, we are to run until time $T = 4π$ when the solution should be identical again to the initial solution, because two full revolutions have been completed. I wish to solve this using matlab so I can get the "animated" graph as result (read further). I want to use "Heaviside" to specify the second initial value in an elegant way.

I am lost but what I know so far is that to advance this Equation $(1)$, is in essence a wave equation. So whatever initial conditions you put in, they should translate in the x-direction as time passes. From my very basic knowledge, I understand that making an array to hold the "old" values of $c(x,t)$, and another one to hold the arrays of new values $c(x,t+\Delta t)$. Resolution selected in $x$ by calling $\Delta x=2\pi/N_x$, and with some timestep $\Delta t$. Set initial conditions as $c(x,t=0)$ and as above which is then iteratively advances the given eq 1 until the time equals $4π$.

Now I specifically want to use pseudo-spectral method with implicit midpoint rule whose code I already have available to me and first order upwind Finite Volume method with forward Euler for the transport equation

Advance the equation in time by making a for-loop, and stepping the solution forward.
For example equation is df(x,t)/dt = v*f(x,t), then we can write
for k=1,nx !nx=number of points in x direction f(k,0)=whatever your initial conditions are endfor
dt=1.d-3 !dt is some arbitrary time step
for n=1,ntimes for k=1,nx f(k,n)=f(k,n-1)+vdtf(k,n-1) endfor!
now then replace the old values with the new ones for k=1,nx f(k,n+1)=f(k,n) endfor endfor

everything can be made much easier by exploiting that Matlab can update entire arrays, but I forget the exact syntax. I am a matlab beginner so please do this for me and return a single package of code that I can directly run with the resultant figures.

Answer 1

The upwind scheme for the linear advection equation $c_t + v c_x = 0$ writes $$ c_i^{n+1} = c_i^n - \frac{\Delta t}{\Delta x} \left( v^+ (c_i^n - c_{i-1}^n) + v^- (c_{i+1}^n - c_i^n)\right) , $$ using the notations $c_i^n \simeq c(i\Delta x, n\Delta t)$, and $v^+ = \max\lbrace v, 0\rbrace$, $v^- = \min\lbrace v, 0\rbrace$. It is stable under the Courant-Friedrichs-Lewy (CFL) condition $\gamma < 1$, where $\gamma = |v|\frac{\Delta t}{\Delta x}$ is the Courant number.

Here is a Matlab code for the upwind scheme with the sinusoidal initial data (the rest can be adapted as desired):

% physical parameter
v = 1;

% numerical parameters
Nx = 100;       % number of points
CFL = 0.95;     % Courant number
T = 4*pi;

% additional variables
vp = max(v,0);
vm = min(v,0);
dx = 2*pi/Nx;
dt = CFL*dx/abs(v);

% data initialization
x = linspace(0-dx,2*pi+dx,Nx+3);
t = 0;
c = sin(x);
c(1) = c(Nx+2);
c(Nx+3) = c(2);
ctemp = c;

% graph
figure(1);
clf;
h = plot(x,c);
xlim([0 2*pi]);
ylim([-1 1]);

while t<T
    % upwind integration
    for i=2:Nx+2
        ctemp(i) = c(i) - dt/dx*( vp*(c(i)-c(i-1)) + vm*(c(i+1)-c(i)) );
    end

    % boundary conditions
    ctemp(1)    = ctemp(Nx+2);
    ctemp(Nx+3) = ctemp(2);

    % actualization
    c = ctemp;
    t = t + dt;
    set(h,'YData',c);
    drawnow;
end

Here is the output:

A convenient way of replacing the sinusoidal data by the rectangular data is to write c = (1/3<x).*(x<2/3) instead of c = sin(x) in line 18 of the previous piece of code.

The pseudo-spectral method writes $$ \hat{c}^{n+1}(\xi) = R(\xi)\, \hat{c}^{n}(\xi) \, , $$ where the $\hat{\cdot}$ denotes the (time-dependent) Fourier transform in space with spatial frequency $\xi$, and $$ R(\xi) = \frac{1}{1+\text{i}v\xi\Delta t}\, ,\qquad R(\xi) = \frac{1-\text{i}v\xi\Delta t/2}{1+\text{i}v\xi\Delta t/2} $$ correspond respectively to the right rule (implicit Euler time integration) and to the trapezoidal rule. To use this method in its right rule version, the previous piece of code is adapted as follows:

• the lines of code 19-20, 36-38 which account for the periodic boundary conditions should be removed.
• here, Nx = 2^nextpow2(100) is a power of two, x = (0:Nx-1)*dx, and xi = [0:Nx/2-1 0 -Nx/2+1:-1]

• the loop in time is replaced by

  while t<T
      % Fourier transform
      ctemp = fft(c,Nx);
  
      % actualization
      R = 1./(1 + 1i*v*xi*dt);
      c = real( ifft(R.*ctemp,Nx) );
      t = t + dt;
      set(h,'YData',c);
      drawnow;
  end