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

变概率椭圆型方程的第一边值问题的随机游走求解

一.问题的提出
二.问题的分析
三.模型的推广
四.问题的求解

一.问题的提出

a(P)∂2u∂x2+b(P)∂2u∂y2+c(P)∂u∂x+d(P)∂u∂y+e(P)u(P)=q(P),p∈Du(Q)=f(Q),Q∈Γ=∂Da(P)\frac{\partial ^2u}{\partial x^2}+b(P)\frac{\partial ^2u}{\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 a(P)∂x2∂2u+b(P)∂y2∂2u+c(P)∂x∂u+d(P)∂y∂u+e(P)u(P)=q(P),p∈Du(Q)=f(Q),Q∈Γ=∂D

若a(P),b(P)>0,e(P)≤0a(P),b(P)>0,e(P)\leq 0a(P),b(P)>0,e(P)≤0求u(P),p∈Du(P),p\in Du(P),p∈D的数值解。

二.问题的分析

将上面的式子差分化之后得到的结果如下:
u=−h22WZ+W4(a+b)[(2a+hc)u11+(2a−hc)u12+(2b+hd)u13+(2b−hd)u14],P∈Du(Q)=f(Q)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) u=−2h2WZ+4(a+b)W[(2a+hc)u11+(2a−hc)u12+(2b+hd)u13+(2b−hd)u14],P∈Du(Q)=f(Q)
其中u=u(P),a=a(P),b=b(P),c=c(P),d=d(P),e=e(P),q=q(P),u1j=u(P1j),Z=Z(P)=qa+b,W=W(P)=[1−h2e2(a+b)]−1u = 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^2e}{2(a+b)}]^{-1}u=u(P),a=a(P),b=b(P),c=c(P),d=d(P),e=e(P),q=q(P),u1j=u(P1j),Z=Z(P)=a+bq,W=W(P)=[1−2(a+b)h2e]−1

且步长要满足如下条件：
h≤2minP∈Da(P)maxP∈D∣c(P)∣,h≤2minP∈Db(P)maxP∈D∣d(P)∣h\leq\frac{2min_{P\in D }a(P)}{max_{P \in D}|c(P)|},h\leq\frac{2min_{P\in D }b(P)}{max_{P \in D}|d(P)|} h≤maxP∈D∣c(P)∣2minP∈Da(P),h≤maxP∈D∣d(P)∣2minP∈Db(P)
那么需要构造的随机游动模型如下：
P(u→u11)=2a1j+hc1j4(a1j+b1j)P(u→u12)=2a1j−hc1j4(a1j+b1j)P(u→u13)=2b1j+hd1j4(a1j+b1j)P(u→u14)=2b1j−hd1j4(a1j+b1j)P(u \rightarrow u_{11}) = \frac{2a_{1j}+hc_{1j}}{4(a_{1j}+b_{1j})}\\ P(u \rightarrow u_{12}) = \frac{2a_{1j}-hc_{1j}}{4(a_{1j}+b_{1j})}\\ P(u \rightarrow u_{13}) = \frac{2b_{1j}+hd_{1j}}{4(a_{1j}+b_{1j})}\\ P(u \rightarrow u_{14}) = \frac{2b_{1j}-hd_{1j}}{4(a_{1j}+b_{1j})}\\ P(u→u11)=4(a1j+b1j)2a1j+hc1jP(u→u12)=4(a1j+b1j)2a1j−hc1jP(u→u13)=4(a1j+b1j)2b1j+hd1jP(u→u14)=4(a1j+b1j)2b1j−hd1j
一条随机游走的路线为:
ρP:P→P1→P2...→Pk−1→Q\rho_P :P\rightarrow P_1 \rightarrow P_2...\rightarrow P_{k - 1}\rightarrow Q ρP:P→P1→P2...→Pk−1→Q
这条路线映射的随机变量ζ\zetaζ为：
ζ=−h22∑m=0k−1(∏i=0mWi)Zm+(∏i=0k−1Wi)f(Q)\zeta = -\frac{h^2}{2}\sum_{m = 0}^{k- 1}(\prod_{i = 0}^mW_i)Z_m + (\prod_{i = 0}^{k - 1}W_i)f(Q) ζ=−2h2m=0∑k−1(i=0∏mWi)Zm+(i=0∏k−1Wi)f(Q)
由于Eζ=u(P)E\zeta = u(P)Eζ=u(P)在此处不证明，由大数定理:
u(P)∼1N∑j=1Nζju(P) \sim \frac{1}{N}\sum_{j = 1}^N\zeta_j u(P)∼N1j=1∑Nζj

三.模型的推广

若上面的模型能化成如下的推广形式：
u=Au11+Bu12+Cu13+Du14+F,P∈Du(Q)=f(Q),Q∈∂Du = Au_{11}+Bu_{12}+Cu_{13}+Du_{14}+F,P \in D\\ u(Q) = f(Q),Q \in \partial D u=Au11+Bu12+Cu13+Du14+F,P∈Du(Q)=f(Q),Q∈∂D
其中:
A=A(P)=k(ah2+c2h)B=B(P)=k(ah2−c2h)C=C(P)=k(bh2+d2h)D=D(P)=k(bh2−d2h)F=F(P)=−kqk=(2a+2bh2−e)−1A = 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} A=A(P)=k(h2a+2hc)B=B(P)=k(h2a−2hc)C=C(P)=k(h2b+2hd)D=D(P)=k(h2b−2hd)F=F(P)=−kqk=(h22a+2b−e)−1
如果我们要求a,b≥0,e≤0a,b \geq 0,e\leq0a,b≥0,e≤0,且A,B,C,D∈[0,1],∀P∈D,A+B+C+D≤1A,B,C,D\in [0,1],\forall P \in D,A+B+C+D\leq1A,B,C,D∈[0,1],∀P∈D,A+B+C+D≤1。我们就可以构造以下的随机游动模型:
P(u→u11)=AP(u→u12)=BP(u→u13)=CP(u→u14)=DP(u→u)=1−(A+B+C+D)P(u \rightarrow u_{11}) = A\\ P(u \rightarrow u_{12}) = B\\ P(u \rightarrow u_{13}) = C\\ P(u \rightarrow u_{14}) = D\\ P(u \rightarrow u) = 1-(A+B+C+D) P(u→u11)=AP(u→u12)=BP(u→u13)=CP(u→u14)=DP(u→u)=1−(A+B+C+D)
产生如下的游走路线:
ρP:P=P0→P1→P2...→Pk−1→Q\rho_P :P = P_0\rightarrow P_1 \rightarrow P_2...\rightarrow P_{k - 1}\rightarrow Q ρP:P=P0→P1→P2...→Pk−1→Q
取的随机变量如下:
ζ=∑i=0k−1F(Pi)+f(Q)\zeta = \sum_{i = 0}^{k-1}F(P_i)+f(Q) ζ=i=0∑k−1F(Pi)+f(Q)
由大数定律:
u(P)∼1N∑j=1Nζju(P) \sim \frac{1}{N}\sum_{j = 1}^N\zeta_j u(P)∼N1j=1∑Nζj

四.问题的求解

假设我们要求得方程为以下形式：
14∂2u∂x2+18∂2u∂y2+12∂u∂x+12∂u∂y=x+2y+1,P∈D:x2+2y2≤2\frac{1}{4}\frac{\partial^2u}{\partial x^2} +\frac{1}{8}\frac{\partial^2u}{\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 \leq2 41∂x2∂2u+81∂y2∂2u+21∂x∂u+21∂y∂u=x+2y+1,P∈D:x2+2y2≤2
且边界条件为:
u(Γ)=2,Γ=∂Du(\Gamma) = 2,\Gamma = \partial D u(Γ)=2,Γ=∂D
这个方程的解析解是:
u=x2+2y2u = x^2+2y^2 u=x2+2y2
现求该方程在u(0.1,0.1)u(0.1,0.1)u(0.1,0.1)处的数值计算的代码及值和第一次的随机游走过程。

现在我们求步长的要求值：h≤0.5h\leq 0.5h≤0.5,取h=0.05h = 0.05h=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:NsumV = 0;xValue = xPoint;yValue = yPoint;P = [xValue ,yValue];Ws = 1;%以下是游走过程while(1)if (xValue)^2+2*(yValue)^2>=2P = [P;xValue,yValue];break;endA = (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 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-1Ws = Ws*W(P(j,1),P(j,2));sumV = sumV + Ws*Z(P(j,1),P(j,2));endzeta(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.710.710.71,理论值为0.750.750.75。