5

matlabでスカラー保存則（たとえば1D）を処理するための堅牢なライブラリまたはFEXのようなパッケージが確立されているかどうか疑問に思いました。

私は現在、1D非線形、非局所、保存則を扱っていますが、一次スキームの拡散誤差は私を殺し、さらに多くの物理学が見落とされています。したがって、コードを自分で作成することを避けるために、堅牢なツールがすでに存在するかどうか疑問に思っています（理想的には、C++でのスキームにとらわれない高次ODE統合のためのboost:: odeintのようなもの）。

助けていただければ幸いです。

編集：明確さの欠如についてお詫びします。ここで保存則とは、次の形式の一般的な双曲線偏微分方程式を意味します。

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

ここu=u(t,x)で、は集中的に保存された変数（たとえば、スカラー、1D、たとえば質量密度、エネルギー密度など）であり、F = F(t,x)はその流束です。したがって、私はハミルトン系が特徴とする種類の保存特性（エネルギー、電流...）には興味がありません（@headmyshoulderのコメントに感謝します）。

boost::odeint数学的な問題（ODEの統合）に対処する堅牢で汎用的なライブラリの概念的なリファレンスを引用しました。したがって、Godunovタイプのメソッドなどを実装するパッケージを探しています。

4

1 に答える 1

3

私は現在、衝撃乱流シミュレーションの新しい方法に取り組んでおり、MATLABで多くのコードテスト/検証を行っています。残念ながら、私はあなたが望んでいることを実行する一般的なライブラリを見つけられませんでしたが、基本的なGodunovまたはMUSCLコードは比較的簡単に実装できます。この論文は、いくつかの有用な方法の概要を示しています。
[1] Kurganov、Alexander and Eitan Tadmor（2000）、非線形保存則と移流拡散方程式のための新しい高解像度中央スキーム、J。Comp。物理学 、160、214〜282。PDF

これは、非粘性バーガース方程式を解くための周期領域上の1D等間隔グリッドのその論文からのいくつかの例です。この方法は、[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値を計算するため、各ステップでグリッドシフトも必要です。主な機能は次のようになります。

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 に答える