5

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

我目前正在处理一维非线性，非局部，守恒定律，一阶方案的扩散误差让我很生气，而且错过了很多物理学。因此，我想知道是否已经有一些强大的工具，以避免自己编写一些代码（理想情况下，类似boost::odeint的东西，用于 C++ 中与方案无关的高阶 ODE 集成）。

任何帮助表示赞赏。

编辑：对缺乏明确性表示歉意。对于守恒定律，我的意思是一般双曲偏导方程的形式

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

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

我引用boost::odeint了一个解决数学问题（ODE 集成）的健壮和通用库的概念参考。因此，我正在寻找一些实现 Godunov 类型方法的包，等等。

4

1 回答 1

3

我目前正在研究冲击湍流模拟的新方法，并在 MATLAB 中进行大量代码测试/验证。不幸的是，我还没有找到一个通用库来满足您的期望，但是基本的 Godunov 或 MUSCL 代码实现起来相对简单。本文对一些有用的方法进行了很好的概述：
[1] Kurganov、Alexander 和 Eitan Tadmor (2000)，New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations，J. Comp。物理。, 160, 214–282。PDF格式

以下是该论文中的一些示例，用于求解无粘性 Burgers 方程的周期域上的一维等距网格。这些方法很容易推广到方程组、耗散（粘性）系统和 [1] 中概述的更高维度。这些方法依赖于以下函数：

通量项：

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

最小模块功能：

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

方法

Nessyahu-Tadmor 方案：

二阶方案_

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 值，因此它还需要在每一步进行网格移位。主要功能应该是这样的：

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
于 2014-03-19T21:54:43.727 回答