Random Walk——高斯消元法

题目

有一个 $N \times M$ 大小的格子，从(0, 0)出发，每一步朝着上下左右4个格子中可以移动的格子等概率移动。另外有些格子有石头，因此无法移至这些格子。求第一次到达 $(N-1, M-1)$ 格子的期望步数。（$2 \leq N,M\leq 10$）

分析

设 $E(x, y)$ 表示从 (x, y) 出发到终点的期望步数。

我们先考虑从 $(x, y)$  向上下左右4个方向都可以移动的情况，由于向4个方向的移动的概率是相等的，因此可以建立如下关系：

$$
\begin{aligned}
E(x, y) & =  \frac{1}{4}(E(x-1, y)+1) + \frac{1}{4}(E(x+1,y)+1) + \frac{1}{4}(E(x, y-1)+1) + \frac{1}{4}(E(x, y+1)+1)\\
&=\frac{1}{4}(E(x-1, y) + E(x+1, y) + E(x, y-1) + E(x, y+1)) + 1\\
\end{aligned}$$

如果移动不是等概率，只需要把 1/4 改成相应的数值就可以了。

如果存在不能移动的方向，我们也可以列出类似的式子。

为了使方程有唯一解，我们令有石头的格子和无法到达的格子都有 $E(x, y) = 0$。

把得到的 $N \times M$ 个方程联立，就可以解出期望步数。

#include<bits/stdc++.h>
using namespace std;

const int maxn = 10+5;
const int maxm = 10+5;
char grid[maxn][maxm+1];
int N, M;

bool vis[maxn][maxm];   //can[x][y]为true表示(x, y)能够达到终点
const int dx[4] = {-1, 0, 1, 0};
const int dy[4] = {0, 1, 0, -1};

//搜索可以达到终点的点
void dfs(int x, int y)
{
    vis[x][y] = true;
    for(int i = 0;i < 4;i++)
    {
        int nx = x + dx[i], ny = y + dy[i];
        if(nx >= 0 && nx <N && ny >= 0 && ny < M && !vis[nx][ny])
            if(grid[nx][ny] != '#')  dfs(nx, ny);
    }
}

const double eps = 1e-8;
typedef double Matrix[maxn*maxm][maxn*maxm];

//结果为A[i][n]/A[i][i]
void gauss_jordan(Matrix A, int n)
{

    int i, j, k, r;
    for(i = 0;i < n;i++)
    {
         //选绝对值一行r并与第i行交换
        r = i;
        for(j = i+1;j < n;j++)
            if(fabs(A[j][i]) > fabs(A[r][i]))  r = j;
        if(fabs(A[r][i]) < eps)  continue;      //放弃这一行，直接处理下一行
        if(r != i)  for(j = 0;j <= n;j++)  swap(A[r][j], A[i][j]);

        //与除第i行外的其他行进行消元
        for(k = 0;k < n;k++)  if(k != i)
            for(j = n;j >= i;j--)  A[k][j] -= A[k][i] / A[i][i] * A[i][j];
    }
}

void debug_print(Matrix A, int n)
{
    for(int i = 0;i < n;i++)
        for(int j = 0;j < n;j++)
            printf("%f%c", A[i][j], j == n-1 ? '\n' : ' ');
}

Matrix A;

int main()
{
    scanf("%d%d", &N, &M);
    for(int i = 0;i < N;i++)
    {
        char s[15];
        scanf("%s", s);
        for(int j = 0;j < M;j++) grid[i][j] = s[j];
    }

    dfs(0, 0);

    //构建矩阵
    for(int i = 0;i < N;i++)
        for(int j = 0;j < M;j++)
        {
            if((i == N-1 && j == M-1) || !vis[i][j])  //终点/无法到达的点，令期望步数为0
            {
                A[i*M + j][i*M + j] = 1;
                continue;
            }
            int move = 0;
            for(int k = 0;k <4;k++)
            {
                int nx = i + dx[k], ny = j + dy[k];
                if(nx >= 0 && nx < N && ny >= 0 && ny < M && grid[nx][ny] == '.')
                {
                    A[i*M + j][nx*M + ny] = -1;
                    move++;
                }
            }
            A[i*M + j][i*M + j] = A[i*M + j][N*M] = move;
        }

    //debug_print(A, N*M+1);

    gauss_jordan(A, N*M);

    printf("%.8f\n", A[0][N*M] / A[0][0]);
    return 0;
}