How do I properly implement a hyperbolic TVD flux limiter? I'm trying to implement a TVD flux limiting scheme for the hyperbolic transport equation. However, I am currently stuck at implementing the limiter. How would I properly do this?
Here's how far I got:
We start with the transport equation:
$$
\partial_t u(x,t) = -\partial_xF\left(u,x,t\right)
$$
Using a flux conservative discretization, we get:
$$
\frac{u^{k+1}_i-u^k_i}{\Delta t} = -\frac{1}{\Delta x_i}\left(F^k_{i+1/2}-F^k_{i-1/2}\right)
$$
The right and left side flux are calculated by mixing of a high-order scheme $f^H$ and a lower-order scheme $f^L$:
$$
F^k_{i\pm1/2}=\left(1-\phi_{i\pm1/2}\right)f^L_{i\pm1/2}+\phi_{i\pm1/2}f^H_{i\pm1/2}
$$
The functions $\phi_{i\pm1/2}$ are called the flux limiters, and that's where my problem is! The articles I find, all assume $F$ to be positive and use the following approach:
$$
\phi_{i+1/2}=\phi\left(r_i\right),\quad r_i=\frac{F^k_i-F^k_{i-1}}{F^k_{i+1}-F^k_i}\\
\phi_{i-1/2}=\phi\left(r_{i-1}\right),\quad r_{i-1}=\frac{F^k_{i-1}-F^k_{i-2}}{F^k_i-F^k_{i-1}}
$$
For $\phi\left(r\right)$ they use one of the following functions: https://en.wikipedia.org/wiki/Flux_limiter
I tried to implement TVD flux limiting with the low-order Osher scheme and the high-order Lax-Wendroff scheme, however my solution is somehow broken. How can I fix this?
My schemes:
$$
f^L_{i+1/2}=\min{\left(0,F^k_{i+1}\right)}-\min{\left(0,F^k_i\right)}\\
f^L_{i-1/2}=\max{\left(0,F^k_{i-1}\right)}-\max{\left(0,F^k_i\right)}\\
f^H_{i+1/2}=\frac{F^k_{i+1}+F^k_i}{2}-\frac{\Delta t}{2\Delta x_i}A^k_{i+1/2}\left(F^k_{i+1}-F^k_i\right)\\
f^H_{i-1/2}=\frac{F^k_i+F^k_{i-1}}{2}-\frac{\Delta t}{2\Delta x_i}A^k_{i-1/2}\left(F^k_i-F^k_{i-1}\right)
$$
(The $A^k_{i\pm1/2}$ are the Jacobian $\partial_u F\left(u,x,t\right)$ at $x_{i\pm1/2},u^k_{i\pm1/2}$.)
My function:
$$
F\left(u,x,t\right)=u
$$
My initial condition:
$$
u\left(x,t=0\right)=\begin{cases}
1,-\frac{1}{2}\le x\le\frac{1}{2}\\
0,\text{everywhere else}
\end{cases}$$
My result after 60 iterations with $\Delta t=0.02$:

Not what I expected...
EDIT:
If I use $\phi_{i-1/2}=\phi\left(r_i\right)$, the solution is way better:

But breaks down for negative function values:
$$
F\left(u,x,t\right)=-u
$$

 A: Let us describe the flux-limiter scheme for the linear advection equation $u_t = -u_x$, as done in the book (1) p. 176, §16.2 (see also the article (2)). 
The numerical flux is defined as a combination of the upwing and Lax-Wendroff fluxes:
$$
F_{i+1/2}^k = u_i^k + \phi(r_i) \left(1 - \frac{\Delta t}{\Delta x_i}\right) \frac{u_{i+1}^k - u_{i}^k}{2} \, ,
\qquad 
r_i =\frac{u_i^k - u_{i-1}^k}{u_{i+1}^k - u_{i}^k}\, ,
$$
where $\phi$ is deduced from a limiter function --e.g. the minmod limiter $\phi(r) = \max(0, \min(1,r))$-- and the smoothness indicator $r_i$ is the ratio of consecutive gradients.
For arbitrary wave speeds $a$ such that $u_t = -au_x$, the method has numerical flux
$$
F_{i+1/2}^k = a\frac{u_i^k + u_{i+1}^k}{2} - |a|\frac{u_{i+1}^k - u_{i}^k}{2} + \phi(r_i) \left(\operatorname{sgn} a - a\frac{\Delta t}{\Delta x_i}\right)  a \frac{u_{i+1}^k - u_{i}^k}{2} \, ,
$$
and the smoothness indicator is the ratio of slopes in the upwind direction
$$
r_i = \frac{u_{I+1}^k - u_{I}^k}{u_{i+1}^k - u_{i}^k}, \qquad I = i - \operatorname{sgn} a \, .
$$
Below is a Matlab implementation of it
% parameters
N = 400;
CFL = 0.95; % Courant number
a = 1;
tend = 1;

dx = 4/N;
x = ((-2-2*dx):dx:(2+2*dx))';
Nx = length(x);
dt = CFL*dx/abs(a);

% initialization
t = 0;
u0 = @(x) (-0.5<x).*(x<0.5);
u = u0(x);

% iterations
while (t+dt)<tend
    up = circshift(u,-1);
    upp = circshift(up,-1);
    um = circshift(u,+1);

    phi = zeros(Nx,1);
    r = ones(Nx,1);
    if a>=0
        num = u-um;
        den = up-u;
    else
        num = upp-up;
        den = up-u;
    end
    for i=1:Nx
        if abs(den(i))>=eps
            r(i) = num(i)./den(i);
        end
    end
    phi = min([ones(Nx,1), r],[],2);
    phi = max([zeros(Nx,1), phi],[],2);

    Fp05 = 0.5*(a*(u+up)-abs(a)*(up-u)) + 0.5*(1-CFL)*abs(a)*(up-u).*phi;
    Fm05 = circshift(Fp05,+1);

    u = u - dt/dx*(Fp05 - Fm05);
    u(1) = u(3);
    u(2) = u(3);
    u(end-1) = u(end-2);
    u(end) = u(end-2);

    t = t + dt;
end

% graphics
figure;
plot(x,u,'bo');
hold on
plot(x+a*t,u0(x),'k-');
xlim([x(1) x(end)]);
ylim([0 2]);

with the following output


Note that the analytical solution $u(x,t) = u(x- at,0)$ of $u_t = -au_x$ translates the initial data with speed $a$. For the prescribed rectangular initial data, we have
$$u(x,t) = \Bbb I_{at-1/2\leq x\leq at+1/2} $$
where $\Bbb I$ denotes the indicator function.
(1) R.J. LeVeque, Numerical Methods for Conservation Laws. Birkhäuser, 1992. doi:10.1007/978-3-0348-8629-1
(2) P.K. Sweby: "High resolution schemes using flux limiters for hyperbolic conservation laws", SIAM J. Numer. Anal. 21(5) (1984), 995-1011. doi:10.1137/0721062
A: is it possible to get python script for the solution in any particular scheme. From there, I would like to modulate the graph by changing different flux limiter functions. If possible please reply sooner.
