Matlab,保护法的整合

Question

我想知道是否有成熟的库或类似FEX的软件包来处理matlab中的标量守恒定律(比如1D).

我目前正在处理一维非线性,非局部,守恒定律,一阶方案的扩散误差正在扼杀我,而且很多物理学都被遗漏了.因此,我想知道是否已经存在一些强大的工具,以避免自己烹饪一些代码(理想情况下,类似于boost :: odeint,用于C++中与方案无关的高阶ODE集成).

任何帮助赞赏.

编辑:对缺乏清晰度的道歉.这里的守恒定律是指形式中的一般hyberbolic偏导数方程

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

其中u=u(t,x)是一个密集的守恒变量(比如标量,1D,例如质量密度,能量密度......)并且F = F(t,x)是它的通量.因此,我对汉密尔顿系统特征(能量,电流......)的保护特性并不感兴趣(感谢@headmyshoulder的评论).

我引用boost::odeint了一个强大的通用库来解决数学问题(ODE的集成)的概念性参考.因此我正在寻找一些实现Godunov类型方法的包等等.

Answer 1

我目前正在研究冲击湍流模拟的新方法，并在 MATLAB 中进行大量代码测试/验证。不幸的是，我还没有找到一个通用库来实现您所希望的功能，但是基本的 Godunov 或 MUSCL 代码相对容易实现。本文很好地概述了一些有用的方法：
\n[1] Kurganov、Alexander 和 Eitan Tadmor (2000)，非线性守恒定律和对流扩散方程的新高分辨率中心方案，J. Comp。物理。, 160, 214\xe2\x80\x93282。PDF

\n\n

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

\n\n

通量术语：

\n\n

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

\n\n

最小模函数：

\n\n

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

\n\n

方法

\n\n

Nesyahu-Tadmor 方案：

\n\n

二阶方案^​

\n\n

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

\n\n

_{该方法在 x j+1/2}网格点处计算新的 u 值，因此每一步还需要进行网格移位。主要功能应该是这样的：

\n\n

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

\n\n

库尔加诺夫-塔德莫尔方案

\n\n

[1] 的 Kurganov-Tadmor 方案比 NT 方案有几个优点，包括较低的数值耗散和允许使用您选择的任何时间积分方法的半离散形式。使用与上述相同的空间离散化，可以将其表示为 du/dt = (stuff) 的 ODE。该 ODE 的右侧可以通过以下函数计算：

\n\n

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

\n\n

该函数还依赖于了解 F(u) 的雅可比行列式的谱半径（rhodF在上面的代码中）。对于无粘性的汉堡来说，这只是

\n\n

function rho = rhodF(u)\ndFdu=abs(u);\n

\n\n

KT方案的主程序可能是这样的：

\n\n

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

\n