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

可以用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:N
    sumV = 0;
    xValue = xPoint;
    yValue = yPoint;
    P = [xValue ,yValue];
    %以下是游走过程
    while(1)
        if (xValue)^2+(yValue)^2>=2
            P = [P;xValue,yValue];
            break;
        end
        randNumber = randsrc(1,1,[[0 1 2 3];[0.25 0.25 0.25 0.25]]);
        switch (randNumber)
            case 0
                xValue = xValue + h;
                yValue = yValue;
            case 1
                xValue = xValue ;
                yValue = yValue + h;
            case 2
                xValue = xValue - h ;
                yValue = yValue ;    
            case 3
                xValue = xValue;
                yValue = yValue - h;
        end
        P = [P;xValue,yValue];
    end
    %以下是画图过程
    if i == 1
        subplot(121);
        grid on;
        plot(P(:,1),P(:,2),'b.-');
        title('第一次轨迹');
    end
    %以下是计算过程
    [m,n] = size(P);
    for j = 1:m-1
        if isnan(q(P(j,1),P(j,2))) == 1
            q(P(j,1),P(j,2)) = q(h,h);
        end
        sumV = sumV + q(P(j,1),P(j,2));
    end
    zeta(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$。可以发现还是有点误差的。

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