一.问题的提出
a ( P ) ∂ 2 u ∂ x 2 + b ( P ) ∂ 2 u ∂ y 2 + c ( P ) ∂ u ∂ x + d ( P ) ∂ u ∂ y + e ( P ) u ( P ) = q ( P ) , p ∈ D u ( Q ) = f ( Q ) , Q ∈ Γ = ∂ D a(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 ) ≤ 0 a(P),b(P)>0,e(P)\leq 0 a(P),b(P)>0,e(P)≤0求 u ( P ) , p ∈ D u(P),p\in D u(P),p∈D的数值解。
二.问题的分析
将上面的式子差分化之后得到的结果如下:
u
=
−
h
2
2
W
Z
+
W
4
(
a
+
b
)
[
(
2
a
+
h
c
)
u
11
+
(
2
a
−
h
c
)
u
12
+
(
2
b
+
h
d
)
u
13
+
(
2
b
−
h
d
)
u
14
]
,
P
∈
D
u
(
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
)
,
u
1
j
=
u
(
P
1
j
)
,
Z
=
Z
(
P
)
=
q
a
+
b
,
W
=
W
(
P
)
=
[
1
−
h
2
e
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),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
≤
2
m
i
n
P
∈
D
a
(
P
)
m
a
x
P
∈
D
∣
c
(
P
)
∣
,
h
≤
2
m
i
n
P
∈
D
b
(
P
)
m
a
x
P
∈
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
→
u
11
)
=
2
a
1
j
+
h
c
1
j
4
(
a
1
j
+
b
1
j
)
P
(
u
→
u
12
)
=
2
a
1
j
−
h
c
1
j
4
(
a
1
j
+
b
1
j
)
P
(
u
→
u
13
)
=
2
b
1
j
+
h
d
1
j
4
(
a
1
j
+
b
1
j
)
P
(
u
→
u
14
)
=
2
b
1
j
−
h
d
1
j
4
(
a
1
j
+
b
1
j
)
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
→
P
1
→
P
2
.
.
.
→
P
k
−
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
ζ为:
ζ
=
−
h
2
2
∑
m
=
0
k
−
1
(
∏
i
=
0
m
W
i
)
Z
m
+
(
∏
i
=
0
k
−
1
W
i
)
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
)
∼
1
N
∑
j
=
1
N
ζ
j
u(P) \sim \frac{1}{N}\sum_{j = 1}^N\zeta_j
u(P)∼N1j=1∑Nζj
三.模型的推广
若上面的模型能化成如下的推广形式:
u
=
A
u
11
+
B
u
12
+
C
u
13
+
D
u
14
+
F
,
P
∈
D
u
(
Q
)
=
f
(
Q
)
,
Q
∈
∂
D
u = 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
(
a
h
2
+
c
2
h
)
B
=
B
(
P
)
=
k
(
a
h
2
−
c
2
h
)
C
=
C
(
P
)
=
k
(
b
h
2
+
d
2
h
)
D
=
D
(
P
)
=
k
(
b
h
2
−
d
2
h
)
F
=
F
(
P
)
=
−
k
q
k
=
(
2
a
+
2
b
h
2
−
e
)
−
1
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}
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
≤
0
a,b \geq 0,e\leq0
a,b≥0,e≤0,且
A
,
B
,
C
,
D
∈
[
0
,
1
]
,
∀
P
∈
D
,
A
+
B
+
C
+
D
≤
1
A,B,C,D\in [0,1],\forall P \in D,A+B+C+D\leq1
A,B,C,D∈[0,1],∀P∈D,A+B+C+D≤1。我们就可以构造以下的随机游动模型:
P
(
u
→
u
11
)
=
A
P
(
u
→
u
12
)
=
B
P
(
u
→
u
13
)
=
C
P
(
u
→
u
14
)
=
D
P
(
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
=
P
0
→
P
1
→
P
2
.
.
.
→
P
k
−
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
=
0
k
−
1
F
(
P
i
)
+
f
(
Q
)
\zeta = \sum_{i = 0}^{k-1}F(P_i)+f(Q)
ζ=i=0∑k−1F(Pi)+f(Q)
由大数定律:
u
(
P
)
∼
1
N
∑
j
=
1
N
ζ
j
u(P) \sim \frac{1}{N}\sum_{j = 1}^N\zeta_j
u(P)∼N1j=1∑Nζj
四.问题的求解
假设我们要求得方程为以下形式:
1
4
∂
2
u
∂
x
2
+
1
8
∂
2
u
∂
y
2
+
1
2
∂
u
∂
x
+
1
2
∂
u
∂
y
=
x
+
2
y
+
1
,
P
∈
D
:
x
2
+
2
y
2
≤
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
,
Γ
=
∂
D
u(\Gamma) = 2,\Gamma = \partial D
u(Γ)=2,Γ=∂D
这个方程的解析解是:
u
=
x
2
+
2
y
2
u = x^2+2y^2
u=x2+2y2
现求该方程在
u
(
0.1
,
0.1
)
u(0.1,0.1)
u(0.1,0.1)处的数值计算的代码及值和第一次的随机游走过程。
现在我们求步长的要求值: h ≤ 0.5 h\leq 0.5 h≤0.5,取 h = 0.05 h = 0.05 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.71 0.71,理论值为 0.75 0.75 0.75。