Matlab, интеграция законов сохранения

Question

Мне было интересно, существуют ли надежные библиотеки или FEX-подобные пакеты для работы со скалярными законами сохранения (скажем, 1D) в matlab.

В настоящее время я имею дело с одномерными нелинейными, нелокальными законами сохранения, и диффузионная ошибка схем первого порядка убивает меня, более того, многие физики упущены. Таким образом, мне интересно, есть ли уже какой-нибудь надежный инструмент, чтобы избежать самостоятельной подготовки некоторого кода (в идеале, что-то вроде boost::odeint для независимой от схемы интеграции ODE высокого порядка в C++).

Любая помощь приветствуется.

РЕДАКТИРОВАТЬ: Извинения за отсутствие ясности. Здесь для законов сохранения я подразумеваю общие гиперболические уравнения с частными производными в виде

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

где u=u(t,x) является интенсивной консервативной переменной (скажем, скаляр, 1D, например, плотность массы, плотность энергии,...) и F = F(t,x) это его поток. Поэтому меня не интересует особенность свойств сохранения гамильтоновых систем (энергия, токи...) (спасибо @headmyshoulder за его комментарий).

Я процитировал boost::odeint для концептуальной ссылки на надежную и универсальную библиотеку для решения математической проблемы (интеграция ODE). Поэтому я ищу какой-то пакет, реализующий методы типа Годунова и так далее.

Answer 1

В настоящее время я работаю над новыми методами моделирования шоковой турбулентности и много тестирую / проверяю код в MATLAB. К сожалению, я не нашел общей библиотеки, которая бы делала то, на что вы надеетесь, но базовый код Годунова или MUSCL относительно прост в реализации. Эта статья имеет хороший обзор некоторых полезных методов:
[1] Курганов Александр и Эйтан Тадмор (2000), Новые центральные схемы высокого разрешения для нелинейных законов сохранения и уравнений конвекции-диффузии, J. Comp. Phys., 160, 214–282. PDF

Вот несколько примеров из этой статьи для одномерной равномерно распределенной сетки в периодической области для решения невязкого уравнения Бюргерса. Методы легко обобщаются на системы уравнений, диссипативные (вязкие) системы и более высокую размерность, как описано в [1]. Эти методы основаны на следующих функциях:

Поток термин:

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

Minmod функция:

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

методы

Схема Нессяху-Тадмор:

Схема 2- ^го порядка

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 );

Этот метод вычисляет новое значение u в точках сетки x _{j + 1/2,} поэтому он также требует сдвига сетки на каждом шаге. Основная функция должна выглядеть примерно так:

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

Схема Курганова-Тадмора

Схема Курганова-Тадмора [1] имеет несколько преимуществ по сравнению со схемой NT, включая меньшую числовую диссипацию и полудискретную форму, которая позволяет использовать любой метод интегрирования по времени, который вы выберете. Используя ту же пространственную дискретизацию, что и выше, ее можно сформулировать как ODE для du/dt = (stuff). Правая часть этого ODE может быть вычислена функцией:

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);

Эта функция также основана на знании спектрального радиуса якобиана F(u) (rhodF в коде выше). Для невидимых бургеров это просто

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

Основная программа схемы KT может быть что-то вроде:

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