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

一.问题的提出

$$\begin{gathered} a(P)\frac{{\partial }^{2}u}{\partial x^2}+b(P)\frac{{\partial }^{2}u}{\partial y^2}+c(P)\frac{\partial u}{\partial x}+d(P)\frac{\partial u}{\partial y}+e(P)u(P)=q(P),p\in D \\ u(Q)=f(Q),Q\in \Gamma =\partial D \end{gathered}$$

若$a(P),b(P)>0,e(P)\le 0$求$u(P),p\in D$的数值解。

二.问题的分析

将上面的式子差分化之后得到的结果如下:
$$\begin{gathered} u=-\frac{h^2}{2}WZ+\frac{W}{4(a+b)}[(2a+hc)u_{11}+(2a-hc)u_{12}+(2b+hd)u_{13}+(2b-hd)u_{14}],P\in D \\ u(Q)=f(Q) \end{gathered}$$
其中$u=u(P),a=a(P),b=b(P),c=c(P),d=d(P),e=e(P),q=q(P),u_{1j}=u(P_{1j}),Z=Z(P)=\frac{q}{a+b},W=W(P)=[1-\frac{h^{2}e}{2(a+b)}{]}^{-1}$

且步长要满足如下条件：
$$h\le \frac{2mi{n}_{P\in D}a(P)}{ma{x}_{P\in D}|c(P)|},h\le \frac{2mi{n}_{P\in D}b(P)}{ma{x}_{P\in D}|d(P)|}$$
那么需要构造的随机游动模型如下：
$$\begin{gathered} P(u\to u_{11})=\frac{2a_{1j}+h{c}_{1j}}{4(a_{1j}+b_{1j})} \\ P(u\to u_{12})=\frac{2a_{1j}-h{c}_{1j}}{4(a_{1j}+b_{1j})} \\ P(u\to u_{13})=\frac{2b_{1j}+h{d}_{1j}}{4(a_{1j}+b_{1j})} \\ P(u\to u_{14})=\frac{2b_{1j}-h{d}_{1j}}{4(a_{1j}+b_{1j})} \end{gathered}$$
一条随机游走的路线为:
$${\rho }_P:P\to P_1\to P_2...\to P_{k-1}\to Q$$
这条路线映射的随机变量$\zeta$为：
$$\zeta =-\frac{h^2}{2}\sum _{m=0}^{k-1}(\prod _{i=0}^m{W}_i)Z_m+(\prod _{i=0}^{k-1}{W}_i)f(Q)$$
由于$E\zeta =u(P)$在此处不证明，由大数定理:
$$u(P)\sim \frac{1}{N}\sum _{j=1}^N{\zeta }_j$$

三.模型的推广

若上面的模型能化成如下的推广形式：
$$\begin{gathered} u=A{u}_{11}+B{u}_{12}+C{u}_{13}+D{u}_{14}+F,P\in D \\ u(Q)=f(Q),Q\in \partial D \end{gathered}$$
其中:
$$\begin{gathered} A=A(P)=k(\frac{a}{h^2}+\frac{c}{2h}) \\ B=B(P)=k(\frac{a}{h^2}-\frac{c}{2h}) \\ C=C(P)=k(\frac{b}{h^2}+\frac{d}{2h}) \\ D=D(P)=k(\frac{b}{h^2}-\frac{d}{2h}) \\ F=F(P)=-kq \\ k=(\frac{2a+2b}{h^2}-e{)}^{-1} \end{gathered}$$
如果我们要求$a,b\ge 0,e\le 0$,且$A,B,C,D\in [0,1],\forall P\in D,A+B+C+D\le 1$。我们就可以构造以下的随机游动模型:
$$\begin{gathered} P(u\to u_{11})=A \\ P(u\to u_{12})=B \\ P(u\to u_{13})=C \\ P(u\to u_{14})=D \\ P(u\to u)=1-(A+B+C+D) \end{gathered}$$
产生如下的游走路线:
$${\rho }_P:P=P_0\to P_1\to P_2...\to P_{k-1}\to Q$$
取的随机变量如下:
$$\zeta =\sum _{i=0}^{k-1}F(P_i)+f(Q)$$
由大数定律:
$$u(P)\sim \frac{1}{N}\sum _{j=1}^N{\zeta }_j$$

四.问题的求解

假设我们要求得方程为以下形式：
$$\frac{1}{4}\frac{{\partial }^{2}u}{\partial x^2}+\frac{1}{8}\frac{{\partial }^{2}u}{\partial y^2}+\frac{1}{2}\frac{\partial u}{\partial x}+\frac{1}{2}\frac{\partial u}{\partial y}=x+2y+1,P\in D:x^2+2y^2\le 2$$
且边界条件为:
$$u(\Gamma )=2,\Gamma =\partial D$$
这个方程的解析解是:
$$u=x^2+2y^2$$
现求该方程在$u(0.1,0.1)$处的数值计算的代码及值和第一次的随机游走过程。

现在我们求步长的要求值：$h\le 0.5$,取$h=0.05$。

function [uFinal,zeta] = possionRandomChange(xPoint,yPoint,h,N)
%UNTITLED 求a = 1/4;b = 1/8;c = 1/2;d = 1/2;e = 0;q = x+2y+1 在u(L) = 2边界条件
%   [xPoint,yPoint]表示该点坐标,h仿真步长,N仿真次数
a = @(x,y)1/4;  % 函数
b = @(x,y)(1/8);
c = @(x,y)1/2;
d = @(x,y)1/2;
e = @(x,y)0;
q = @(x,y)(x+2*y+1);
f = @(x,y)2;
u = zeros(1,N);
W = @(x,y)(1-h^2*e(x,y)/(2*(a(x,y)+b(x,y))))^(-1);
Z = @(x,y)(q(x,y)/(a(x,y)+b(x,y)));
for i = 1:N
    sumV = 0;
    xValue = xPoint;
    yValue = yPoint;
    P = [xValue ,yValue];
    Ws = 1;
    %以下是游走过程
    while(1)
        if (xValue)^2+2*(yValue)^2>=2
            P = [P;xValue,yValue];
            break;
        end
        A = (2*a(xValue,yValue)+h*c(xValue,yValue))/(4*(a(xValue,yValue)+b(xValue,yValue)));
        B = (2*a(xValue,yValue)-h*c(xValue,yValue))/(4*(a(xValue,yValue)+b(xValue,yValue)));
        C = (2*b(xValue,yValue)+h*d(xValue,yValue))/(4*(a(xValue,yValue)+b(xValue,yValue)));
        D = (2*b(xValue,yValue)-h*d(xValue,yValue))/(4*(a(xValue,yValue)+b(xValue,yValue)));
        randNumber = randsrc(1,1,[[0 1 2 3];[A C B D]]);
        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
        Ws = Ws*W(P(j,1),P(j,2));
        sumV = sumV + Ws*Z(P(j,1),P(j,2));
    end
    zeta(i) = -h^2/2*sumV + Ws*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

输入命令行指令：

[uFinal,zeta] = possionRandomChange(0.5,0.5,0.01,1000);

最后的理论结果是:

0.709698130666665;

画出来的结果是:

我只能说这个随机游走求解问题的手段效率有点太低了，收敛的时间过长了,而且和理论值的偏差也有点大，实际值为$0.71$,理论值为$0.75$。