随机游走问题的神奇应用(一)

本文主要是介绍随机游走问题的神奇应用(一)，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

可以用monteCarlo方法构建一个随机游走过程来求解偏微分方程。

一.问题的提出

​ 求解二维泊松方程的第一边值问题如下:
$$\frac{{\partial }^{2}u(P)}{\partial x}+\frac{{\partial }^{2}u(P)}{\partial y^2}=q(P)P(x,y)\in D$$
​ 边界条件为:
$$u(Q)=f(Q)Q(x,y)\in \Gamma =\partial D$$

二.问题的求解

如下图所示如果我们要求在$P(x^{\ast },y^{\ast })$处的值,设求解步长为$h$。我们就将原来的方程差分化:
$$\frac{u(x+h,y)+u(x-h,y)-2u(x,y)}{h^2}+\frac{u(x,y+h)+u(x,y-h)-2u(x,y)}{h^2}=q(x,y)$$
可以化成如下差分形式:
$$u=-\frac{h^2}{4}q+\sum _{i=1}^4{u}_{1i}$$
其中$P_{1i}$的含义如下图所示:

我们假设方程的定义区域为下图，黑点为边界$\Gamma$上的点,白点为内部$D$上的点。

我们定义一个从$P(x^{\ast },y^{\ast })$开始的游走路线：
$${\rho }_P:P\to P_1\to P_2\to P_3...\to P_{k-1}\to Q\in \Gamma$$
其中的点$P_i$是由点$P_{i-1}$进过以下规则$f$得到的：
$$f:P_{i-1}\to P_i$$
$f$:我们产生一个随机数$r\in [0,1]$:

$r$	$P_{i-1}$	$P_i$
$[0,0.25]$	$P_{i-1}(x^{\ast },y^{\ast })$	$P_i(x^{\ast }+h,y^{\ast })$
$[0.25,0.5]$	$P_{i-1}(x^{\ast },y^{\ast })$	$P_i(x^{\ast },y^{\ast }+h)$
$[0.5,0.75]$	$P_{i-1}(x^{\ast },y^{\ast })$	$P_i(x^{\ast }-h,y^{\ast })$
$[0.75,1]$	$P_{i-1}(x^{\ast },y^{\ast })$	$P_i(x^{\ast },y^{\ast }-h)$

表示$P_{i-1}$分别有$\frac{1}{4}$的概率到$P_{11},P_{12},P_{13},P_{14}$，即有相同的概率前后左右随机移动,到终点为止。

那么在这样的规则下会生成一条路线${\rho }_P$。此时我们建立其以下的映射关系:
$$g:{\rho }_p\to u(P)$$
即是从$P$点出发的一条路线${\rho }_P$到该点的数值解$u(P)$的一个映射关系$g$。

现在我们直接给出这个关系：
$$g:u(P)=-\frac{h^2}{4}\sum _{i=1}^{k-1}q(P_i)+f(Q)$$
那么我们从这一次随机游走${\rho }_P$映射出了一次的$u(P)$。这个结果肯定是不精确的，我们如果设每次的结果都是一次随机变量${\zeta }_P=g({\rho }_P)$，那么可以证明的是$E{\zeta }_P=u(P)$,具体证明过程忽略。我们利用这个结论可以由大数定律:
$$u(P)\sim \frac{1}{N}\sum _{i=1}^N{\zeta }_{Pi}$$
即通过多次模拟随机游走的过程求其均值用来表示当前的解$u(P)$。

三.代码求解

假设我们要求解的是以下方程：
$$\frac{{\partial }^{2}u}{\partial x}+\frac{{\partial }^{2}u}{\partial y^2}=1$$
边界$\Gamma :x^2+y^2=2$。在此边界上：
$$u(\Gamma )=\frac{1}{2}$$
现求$u(0,1)$的值。

理论上该方程的解析解为$u(x,y)=\frac{x^2+y^2}{4}$,因此$u(0,1)=\frac{1}{4}$。在这里，我们给出求解该方程的函数并且显示当前的随机游走过程：

function [uFinal,zeta] = possionRandom(xPoint,yPoint,h,N)
%UNTITLED 求 Deltea u = 1;在u(x^2+y^2 = 2) = 0.5边界条件
%   [xPoint,yPoint]表示该点坐标,h仿真步长,N仿真次数
q = @(x,y)1;  % 函数
f = @(x,y)(0.5);
u = zeros(1,N);
for i = 1:NsumV = 0;xValue = xPoint;yValue = yPoint;P = [xValue ,yValue];%以下是游走过程while(1)if (xValue)^2+(yValue)^2>=2P = [P;xValue,yValue];break;endrandNumber = randsrc(1,1,[[0 1 2 3];[0.25 0.25 0.25 0.25]]);switch (randNumber)case 0xValue = xValue + h;yValue = yValue;case 1xValue = xValue ;yValue = yValue + h;case 2xValue = xValue - h ;yValue = yValue ;    case 3xValue = xValue;yValue = yValue - h;endP = [P;xValue,yValue];end%以下是画图过程if i == 1subplot(121);grid on;plot(P(:,1),P(:,2),'b.-');title('第一次轨迹');end%以下是计算过程[m,n] = size(P);for j = 1:m-1if isnan(q(P(j,1),P(j,2))) == 1q(P(j,1),P(j,2)) = q(h,h);endsumV = sumV + q(P(j,1),P(j,2));endzeta(i) = -h^2/4*sumV + f(P(m,1),P(m,2));u(i) = sum(zeta(1:i))/i;
end
subplot(122);
grid on ;
uFinal = u(N);
plot(1:N,u,'r');
xlabel('次数');
ylabel('进化曲线');
title('收敛过程');
end

在命令行求解$u(0,1)$的值如下,设置步长为$h=0.01$,仿真次数$N=100$：

>> [uFinal,zeta] = possionRandom(0,1,0.01,100);

得到如下图形：

最终结果为$0.2353$。可以发现还是有点误差的。

关键是这玩意如果步长设置的比较小的话就会一直游走。所以运行的时间就会比较长。

这篇关于随机游走问题的神奇应用(一)的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

---