以比较抽象的视角来看， 大概就是用停留在状态 (i) 的概率乘上从 (i) 转移到 (j) 的概率取更新停留在状态 (j) 的概率，这个可以用矩阵表示。
设初始停留概率列向量为 (st)， 转移矩阵为 (T)， 那么答案向量就是：
[sum_{ige 0} T^i * st
\
= frac{st}{I-T}
]
其中 (I) 是单位矩阵， (T^0=I)。
求矩阵逆就算了吧， 假设答案向量是 (ans)， 上式即为：
[(I-T)*ans = st
]
这其实就是一个线性方程组， 高斯消元即可。
附：
[egin{cases}
a_{1,1}*x_1+cdots +a_{1,n}*x_1 = b_1
\
vdots
\
a_{n,1}*x_1+cdots +a_{n,n}*x_1 = b_n
end{cases}
]
用矩阵表示就是：
[left[
egin{matrix}
a_{1,1} & cdots & a_{1,n}
\
vdots&ddots&vdots
\
a_{n,1} & cdots &a_{n,n}
end{matrix}
ight]
*
left[
egin{matrix}
x_1
\
vdots
\
x_n
end{matrix}
ight]
=
left[
egin{matrix}
b_1
\
vdots
\
b_n
end{matrix}
ight]
]
具体地， 对于转移矩阵 (T)， 第 (i) 行第 (j) 列的数表示从状态 (j) 转移到状态 (i) 的概率。
参考代码：
#include<bits/stdc++.h>
using namespace std;
const int N = 25, S = 503;
const double eps = 1e-9;
int n,m,a,b,g[N][N],deg[N];
double p[N];
int s;
int id(int x,int y) { return (x-1)*n+y; }
double T[S][S];
void Ins(int x,int y,double k) {
for(int i=1;i<=s+1;++i) T[x][i] += T[y][i]*k;
}
int main()
{
scanf("%d%d%d%d",&n,&m,&a,&b);
s = n*n;
for(int i=0;i<m;++i)
{
int x,y;
scanf("%d%d",&x,&y);
g[x][y] = g[y][x] = 1;
++deg[x], ++deg[y];
}
for(int i=1;i<=n;++i)
scanf("%lf",&p[i]);
for(int i=1;i<=s;++i) T[i][i]=1;
for(int i=1;i<=n;++i)
for(int j=1;j<=n;++j) if(i!=j)
for(int i_=1;i_<=n;++i_) if(i==i_ || g[i][i_])
for(int j_=1;j_<=n;++j_) if(j==j_ || g[j][j_])
{
int now = id(i,j), nxt = id(i_,j_);
double P = 1;
if(i==i_) P *= p[i];
else P *= (1-p[i])/deg[i];
if(j==j_) P *= p[j];
else P *= (1-p[j])/deg[j];
T[nxt][now] -= P;
}
T[id(a,b)][s+1] = 1;
for(int i=1;i<=s;++i) {
int mx = i;
for(int j=i+1;j<=s;++j) if(fabs(T[j][i]) > fabs(T[mx][i])) mx=j;
swap(T[i],T[mx]);
for(int j=1;j<=s;++j) if(i^j) Ins(j,i,-T[j][i]/T[i][i]);
}
for(int i=1;i<=n;++i) {
int q = id(i,i);
printf("%.10lf ", T[q][s+1]/T[q][q]);
}
return 0;
}