数值方法- Matlab，守恒定律的积分

Question

我想知道是否有建立健壮的库或类似fex的包来处理matlab中的标量守恒定律(例如1D)。

我目前正在处理一维非线性，非局部，守恒定律和一阶格式的扩散误差是杀了我，而且很多物理遗漏。因此，我想知道是否已经有一些强大的工具来避免自己编写一些代码(理想情况下，像boost::odeint这样的工具，用于c++中与模式无关的高阶ODE集成)。

感谢您的帮助。

EDIT:为不够清晰道歉。对于守恒定律，我指的是形式为

的一般双曲偏导数方程

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

，其中u=u(t,x)是一个密集守恒变量(如标量、一维，如质量密度、能量密度等)，F = F(t,x)是它的通量。因此，我对哈密顿系统特征(能量、电流……)的守恒性质不感兴趣(感谢@headmyshoulder的评论)。

我引用boost::odeint作为解决数学问题(ode的集成)的健壮和通用库的概念参考。因此，我正在寻找一些实现godunov类型方法等的包。

Answer 1

我目前正在研究冲击湍流模拟的新方法，并在MATLAB中进行了大量的代码测试/验证。不幸的是，我还没有找到一个通用的库，但一个基本的Godunov或musl代码相对容易实现。本文很好地概述了一些有用的方法:
[1]张建军，张建军，张建军(2000)，非线性守恒方程的高分辨率中心格式，[J]。广告样稿。物理。， 160, 214-282。PDF

下面是该论文中的几个例子，用于求解无粘Burgers方程的周期域上的一维等间距网格。这些方法很容易推广到方程系统、耗散(粘性)系统和高维系统，如[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));

Nessyahu-Tadmor方案:

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

此方法在x_j+1/2网格点处计算一个新的u值，因此它还需要在每一步进行网格移动。main函数应该是这样的:

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方案

[1]的Kurganov-Tadmor格式比NT格式有几个优点，包括较低的数值耗散和半离散形式，允许使用您选择的任何时间积分方法。使用与上面相同的空间离散化，它可以被表述为du/dt = (stuff)的ODE。该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