matlabnumerical-methodsconservation-laws

Matlab, integration of conservation laws


I was wondering if there are established robust libraries or FEX-like packages to deal with scalar conservation laws (say 1D) in matlab.

I am currently dealing with 1D non-linear, non-local, conservation laws and the diffusive error of first order schemes is killing me, moreover a lot of physics is missed. Thus, I am wondering if there is some robust tool already there so to avoid cooking some code myself (ideally, something like boost::odeint for scheme agnostic high order ODE integration in C++).

Any help appreciated.

EDIT: Apologies for the lack of clarity. Here for conservation laws I mean general hyberbolic partial derivative equations in the form

 u_t(t,x) + F_x(t,x) = 0

where u=u(t,x) is an intensive conserved variable (say scalar, 1D, e.g. mass density, energy density,...) and F = F(t,x) is its flux. Therefore, I am not interested in the kind of conservation properties Hamiltonian systems feature (energy, currents...) (thanks to @headmyshoulder for his comment).

I cited boost::odeint for a conceptual reference of a robust and generic library addressing a mathematical issue (integration of ODEs). Therefore I am looking for some package implementing Godunov-type methods and so on.


Solution

  • I am currently working on new methods for shock-turbulence simulations and doing lots of code testing/validation in MATLAB. Unfortunately, I haven't found a general library that does what you're hoping, but a basic Godunov or MUSCL code is relatively straightforward to implement. This paper has a good overview of some useful methods:
    [1] Kurganov, Alexander and Eitan Tadmor (2000), New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations, J. Comp. Phys., 160, 214–282. PDF

    Here are a few examples from that paper for a 1D equally spaced grid on a periodic domain for solving inviscid Burgers equation. The methods easily generalize to systems of equations, dissipative (viscous) systems, and higher dimension as outlined in [1]. These methods rely on the following functions:

    Flux term:

    function f = flux(u)
    %flux term for Burgers equation: F(u) = u^2/2;
    f = u.^2/2;
    

    Minmod function:

    function m = minmod(a,b)
    %minmod function:
    m = (sign(a)+sign(b))/2.*min(abs(a),abs(b));
    

    Methods

    Nessyahu-Tadmor scheme:

    A 2nd order scheme

    function unew = step_u(dx,dt,u)
    %%%   Nessyahu-Tadmor scheme
    
    ux = minmod((u-circshift(u,[0 1]))/dx,(circshift(u,[0 -1])-u)/dx);
    
    f = flux(u);
    fx = minmod((f-circshift(f,[0 1]))/dx,(circshift(f,[0 -1])-f)/dx);
    
    umid = u-dt/2*fx;
    fmid = flux(umid);
    
    unew = (u + circshift(u,[0 -1]))/2 + (dx)/8*(ux-circshift(ux,[0 -1])) ...
          -dt/dx*( circshift(fmid,[0 -1])-fmid );
    

    This method computes a new u value at xj+1/2 grid points so it also requires a grid shift at each step. The main function should be something like:

    clear all
    
    % Set up grid
    nx = 256;
    xmin=0; xmax=2*pi;
    x=linspace(xmin,xmax,nx);
    dx = x(2)-x(1);
    
    %initialize
    u = exp(-4*(x-pi*1/2).^2)-exp(-4*(x-pi*3/2).^2);
    
    %CFL number:
    CFL = 0.25;
    
    t = 0;
    dt = CFL*dx/max(abs(u(:)));
    while (t<2)
    
        u = step_u(dx,dt,u);
        x=x+dx/2;
    
        % handle grid shifts
        if x(end)>=xmax+dx
            x(end)=0;
            x=circshift(x,[0 1]);
            u=circshift(u,[0 1]);
        end
    
        t = t+dt;
    
        %plot
        figure(1)
        clf
        plot(x,u,'k')
        title(sprintf('time, t = %1.2f',t))
        if ~exist('YY','var')
            YY=ylim;
        end
        axis([xmin xmax YY])
        drawnow
    end
    

    Kurganov-Tadmor scheme

    The Kurganov-Tadmor scheme of [1] has several advantages over the NT scheme including lower numerical dissipation and a semi-discrete form that allows the use of any time integration method you choose. Using the same spatial discretization as above, it can be formulated as an ODE for du/dt = (stuff). The right hand side of this ODE can be computed by the function:

    function RHS = KTrhs(dx,u)
    %%% Kurganov-Tadmor scheme
    
    ux = minmod((u-circshift(u,[0 1]))/dx,(circshift(u,[0 -1])-u)/dx);
    uplus = u-dx/2*ux;
    uminus = circshift(u+dx/2*ux,[0 1]);
    a = max(abs(rhodF(uminus)),abs(rhodF(uplus)));
    RHS = -( flux(circshift(uplus,[0 -1]))+flux(circshift(uminus,[0 -1])) ...
             -flux(uplus)-flux(uminus) )/(2*dx) ...
          +( circshift(a,[0 -1]).*(circshift(uplus,[0 -1])-circshift(uminus,[0 -1])) ...
             -a.*(uplus-uminus) )/(2*dx);
    

    This function also relies on knowing the spectral radius of the Jacobian of F(u) (rhodF in the code above). For inviscid Burgers this is just

    function rho = rhodF(u)
    dFdu=abs(u);
    

    The main program of the KT scheme could be something like:

    clear all
    
    nx = 256;
    xmin=0; xmax=2*pi;
    x=linspace(xmin,xmax,nx);
    dx = x(2)-x(1);
    
    %initialize
    u = exp(-4*(x-pi*1/2).^2)-exp(-4*(x-pi*3/2).^2);
    
    %CFL number:
    CFL = 0.25;
    
    t = 0;
    dt = CFL*dx/max(abs(u(:)));
    while (t<3)
    
    
        % 4th order Runge-Kutta time stepping
        k1 = KTrhs(dx,u);
        k2 = KTrhs(dx,u+dt/2*k1);
        k3 = KTrhs(dx,u+dt/2*k2);
        k4 = KTrhs(dx,u+dt*k3);
        u = u+dt/6*(k1+2*k2+2*k3+k4);
    
        t = t+dt;
    
        %plot
        figure(1)
        clf
        plot(x,u,'k')
        title(sprintf('time, t = %1.2f',t))
        if ~exist('YY','var')
            YY=ylim;
        end
        axis([xmin xmax YY])
        drawnow
    end