用Fast Sweeping Method解Eikonal方程

来源：互联网发布：电脑椅推荐知乎编辑：程序博客网时间：2024/06/05 13:23

Eikonal 方程属于静态的Hamilton-Jacobi方程。静态是指方程不包含时间，或者说解不随时间变化，而Hamilton-Jacobi方程是一种双曲型的非线性偏微分方程。一般来说，非线性偏微分方程难以得到解析解，因此一般采用数值方法求解。Eikonal 方程的数值求解方法有快速行进法(Fast Marching Method)和快速扫描法(Fast Sweeping Method)，后者适用范围更广，效率更高。下面我们在数学软件Mathematica中用Fast Sweeping Method求解二维的Eikonal方程。
Mathematica代码如下：
Initialization (初始化)

gridNum = 50; (*x和y方向的离散网格数量*)center = Round[(gridNum + 1)/2.];inf = 10^6;  (*设置无穷大*)h = 1.0/(gridNum-1);  (*网格节点间的间隔，1.0 是计算范围*)Table[u[i][j] = inf; f[i][j] = 1, {i, 0, gridNum + 1}, {j, 0, gridNum + 1}, {k, 0, gridNum + 1}];(*u表示Eikonal方程的解，它是网格节点的函数，f表示网格节点处的速度*)

Boundary condition (设置边界条件，本文中都选择为中心点)

u[center][center] = 0.0;

One Iteration with Four Sweeping (一次循环包含四次扫描)

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 = Min[u[i - 1][j], u[i + 1][j]];(*u[i][j]在x方向的2个邻居*)    b = Min[u[i][j - 1], u[i][j + 1]];(*u[i][j]在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 (计算结果的可视化)

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.6}, Boxed -> False]ListContourPlot[data, Contours -> 10]

解的图像如下
这里写图片描述
前面我们对所有节点处的值进行了计算，其实还可以考虑环境中存在障碍物的情况，例如：

 ((10 <= i <= 20) && (10 <= j <= 20)) (*障碍物的区域*)

结果如下图
这里写图片描述
更复杂的障碍物需要更多次迭代，下面的形状需要4次：

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]

即：
这里写图片描述
结果如下图所示：

如果我们计算一下解的梯度向量场就会看到更有意思的现象，见下图。
400
不管从哪一点出发（障碍物内除外），一直沿着向量所指的方向前进，我们最后总能绕过障碍物到达中间的点（也就是我们初始化时设置的边界点），而且所得到的路径是最短的，如下动图。
这里写图片描述
如果我们将节点更新步骤

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

改为

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

那么Fast Sweeping Method就等价于著名的Dijkstra最短路径算法(更准确的说，是Dijkstra算法的特殊形式，因为这里我们没有显式的构造优先队列，这时可以称为波前传播算法)。这时得到的解如下图：
这里写图片描述
其对应的水平集如下图

说明
这个程序只是为了展示 Fast Sweeping Method 的原理。由于我们的离散点比较少，而且针对的是二维的Eikonal方程，所以用Mathematica实现还是很快的(50X50的网格耗时不超过0.5s)。
本文中的方法来源于论文 A Fast Sweeping Method For Eikonal Equations (Mathematics of Computation HONGKAI ZHAO)
关于Dijkstra算法和Fast Marching Method的关系可以参考http://www.numerical-tours.com/matlab/fastmarching_0_implementing/#2

1 0