用Fast Sweeping Method解Eikonal Equation方程

Eikonal 方程属于静态(即不含时间)Hamilton-Jacobi方程(是一种非线性偏微分方程)，最高效的求解方法有快速行进法(Fast Marching Method)和快速扫描法(Fast Sweeping Method)。后者的效率更高。下面在Mathematica中用Fast Sweeping Method求解二维的Eikonal方程。

这里以2种边界条件为例：

1. 边界条件是中心点





2. 边界条件是四周的边界




代码如下：

(Initialization)

[code]gridNum = 50;
center = Round[(gridNum + 1)/2.0];
inf = 10^6;
hx = hy = h = 1.0/(gridNum - 1);
Table[u[i][j] = inf; f[i][j] = 1, {i, 0, gridNum + 1}, {j, 0, gridNum + 1}, {k, 0, gridNum + 1}];

Boundary condition 1

[code]u[center][center] = 0;

Boundary condition 2

[code]Table[u[1][j] = 0, {j, gridNum}];
Table[u[gridNum][j] = 0, {j, gridNum}];
Table[u[i][1] = 0, {i, gridNum}];
Table[u[i][gridNum] = 0, {i, gridNum}];

One Iteration with Four Sweeping

[code]For[k = 1, k <= 4, k++,
 Which[k == 1,
  isweepStart = gridNum;
  istep = -1;
  jsweepStart = gridNum;
  jstep = -1,
  k == 2,
  isweepStart = gridNum;
  istep = -1;
  jsweepStart = 1;
  jstep = 1,
  k == 3,
  isweepStart = 1;
  istep = 1;
  jsweepStart = gridNum;
  jstep = -1,
  k == 4,
  isweepStart = 1;
  istep = 1;
  jsweepStart = 1;
  jstep = 1
  ];
 For[i = isweepStart, 1 <= i <= gridNum, i = i + istep,
  For[j = jsweepStart, 1 <= j <= gridNum, j = j + jstep,
   a = uxmin = Min[u[i - 1][j], u[i + 1][j]];(*x方向的2个邻居*)
   b = uymin = Min[u[i][j - 1], u[i][j + 1]];(*y方向的2个邻居*)
   ubar = inf;
   If[Abs[a - b] >= f[i][j]*h,(*判断upwind方向,根据方向选择用1个还是2个邻居*)
    ubar = Min[a, b] + f[i][j]*h;
    ,
    ubar = N[(a + b + Sqrt[2*(f[i][j]*h)^2 - (a - b)^2])/2];
    ];
   u[i][j] = Min[ubar, u[i][j]];
   ]
  ]
 ]

Visualization

[code]data = Table[u[i][j], {i, 1, gridNum, 1}, {j, 1, gridNum, 1}];
ListPlot3D[data, PlotRange -> {{0, gridNum}, {0, gridNum}, {0, gridNum*h}}, Mesh -> None, Axes -> False, ColorFunction -> "Rainbow", BoxRatios -> {1, 1, 0.8}, Boxed -> False]
ListContourPlot[data, Contours -> 10]

还可以试试其他边界条件，例如：

[code]u[15][15] = 0; u[15][35] = 0; u[35][15] = 0; u[35][35] = 0;

解如下图所示：




或者

[code]Table[If[i + j == 50 || Abs[i - j] == 0, u[i][j] = 0], {i, gridNum}, {j, gridNum}];

又或者

[code]Table[If[15 < Sqrt[(i - 25)^2 + (j - 25)^2] < 16, u[i][j] = 0], {i, gridNum}, {j, gridNum}];

我们还可以考虑环境中存在障碍物的情况，例如：

[code] ((10 <= i <= 20) && (10 <= j <= 20))(*障碍物区域*)

结果如下图




更复杂的障碍物(需要更多次迭代，这里需要4次)：

[code]Clear[obstacle];
center = center - 1;
R1 = RegionDifference[Disk[{center, center}, 15], Disk[{center, center}, 13]];
R2 = RegionDifference[R1, Polygon[{{center, center - 2}, {center + 16, center - 2}, {center + 16, center + 2}, {center, center + 2}}]];
obstacle[1] = DiscretizeRegion[R2];
R1 = RegionDifference[Disk[{center, center}, 8], Disk[{center, center}, 6]];
R2 = RegionDifference[R1, Polygon[{{center, center - 2}, {center - 15, center - 2}, {center - 15, center + 2}, {center, center + 2}}]];
obstacle[2] = DiscretizeRegion[R2];
obstacles = Table[obstacle[i], {i, 2}];
obstacleRegionAll = RegionUnion[obstacles]

即：




结果如下2幅图：




另一视图




如果将更新步骤

[code]If[Abs[a - b] >= f[i][j]*h
  ,
  ubar = Min[a, b] + f[i][j]*h;
  ,
  ubar = N[(a + b + Sqrt[2*(f[i][j]*h)^2 - (a - b)^2])/2];
  ];

改为

[code]ubar = Min[a, b] + f[i][j]*h;

那Fast Sweeping Method就变成了著名的Dijkstra algorithm(更准确的说，是Dijkstra 算法的特殊形式，因为这里没有显式的构造优先队列，这时可以称为波前传播算法)。

以第一种边界条件得到的解如下图：




其对应的水平集如下图




说明

这个程序只是为了展示Fast Sweeping Method 的原理。由于我们的离散点比较少，而且考虑二维的Eikonal方程，所以用Mathematica实现还是很快的(耗时不超过0.5s)。关于Dijkstra算法和Fast Marching Method的关系可以参考http://www.numerical-tours.com/matlab/fastmarching_0_implementing/#2