基于蒙特卡洛法的定积分求解

0.摘要

本文首先介绍了蒙特卡洛法的起源和应用，接着详细推导了基于蒙特卡洛法求解定积分的数学理论公式。然后进行相关的实验研究。本文的实验共有三个，分别为蒙特卡洛法可行性研究，nnn值对结果的影响与积分区间对结果的影响。每一类实验都进行大量实验，并对试验结果进行详细的分析最终得出结论。最后附上相关的实验代码

1.蒙特卡洛法简介

蒙特卡罗法也称统计模拟法、统计试验法。是把概率现象作为研究对象的数值模拟方法。是按抽样调查法求取统计值来推定未知特性量的计算方法。它非常强大和灵活，又相当简单易懂，很容易实现。对于许多问题来说，它往往是最简单的计算方法，有时甚至是唯一可行的方法。它诞生于上个世纪40年代美国的"曼哈顿计划"，蒙特卡罗是摩纳哥的著名赌城，该法为表明其随机抽样的本质而命名，故适用于对离散系统进行计算仿真试验。在计算仿真中，通过构造一个和系统性能相近似的概率模型，并在数字计算机上进行随机试验，可以模拟系统的随机特性。
蒙特卡罗方法适用于两类问题，第一类是本身就具有随机性的问题，第二类是能够转化为概率模型进行求解的确定性问题。蒙特卡罗模拟在金融工程学，宏观经济学，生物医学，计算物理学(如粒子输运计算、量子热力学计算、空气动力学计算)等领域也应用广泛。
本文使用蒙特卡罗法求解积分。当无法精确计算和或积分时，通常可以使用蒙特卡罗采样来近似它。这种想法把和或者积分视作某分布下的期望，然后通过估计对应的平均值来近似这个期望。即令

s=∑xp(x)f(x)=Ep[f(x)]s=\sum_{x}p(x)f(x)=E_{p}[f(x)]s=x∑p(x)f(x)=Ep[f(x)]

或者

s=∫p(x)f(x)=Ep[f(x)]s=\int p(x)f(x)=E_{p}[f(x)]s=∫p(x)f(x)=Ep[f(x)]

ppp是一个关于随机变量xxx的概率分布或者概率密度函数
我们可以通过从ppp中抽取n个样本 x(1),⋅⋅⋅,x(n)x^{(1)},···,x^{(n)}x(1),⋅⋅⋅,x(n)来近似s并得到一个经验平均值

sn^=1n∑i=1nf(x(i))\hat{s_n}=\frac{1}{n}\sum^{n}_{i=1}f(x^{(i)})sn^=n1i=1∑nf(x(i))

并且有

E[sn^]=1n∑i=1nE[f(x(i))]=1n∑i=1ns=sE[\hat{s_n}]=\frac{1}{n}\sum^{n}_{i=1}E[f(x^{(i)})]=\frac{1}{n}\sum^n_{i=1}s=sE[sn^]=n1i=1∑nE[f(x(i))]=n1i=1∑ns=s

容易观察到sn^\hat{s_n}sn^这个估计是无偏的。此外，根据大数定理，如果样本x(i)x^{(i)}x(i)是独立同分布的，那么该均值几乎必然会收敛到真实的期望值，即：

lim⁡n→+∞sn^=s{\lim_{n \to +\infty}}\hat{s_n}=sn→+∞limsn^=s

2.公式推导

这一节我们进行详细的公式推导。
设我们的目标函数为g(x)g(x)g(x),g(x)g(x)g(x)为一连续函数。现在我们要计算g(x)g(x)g(x)在区间 [a,b][a,b][a,b]上的定积分∫bag(x)dx\int^a_bg(x)dx∫bag(x)dx。于是我们采用蒙特卡洛法，向区间[a,b][a,b][a,b]上均匀随机连投nnn个随机点，得到坐标X1,X2,...,X3X_1,X_2,...,X_3X1,X2,...,X3为nnn个独立均匀分布的随机变量。
显然有：

E[g(x)]=1b−a∫bag(x)dxE[g(x)]=\frac{1}{b-a}\int^a_bg(x)dxE[g(x)]=b−a1∫bag(x)dx

由大数定理：

lim⁡n→+∞P∣1n∑i=1ng(Xi)−1b−a∫bag(x)dx∣=1{\lim_{n \to +\infty}}P{|\frac{1}{n}\sum^{n}_{i=1}g(X_i)-\frac{1}{b-a}\int^a_bg(x)dx|}=1n→+∞limP∣n1i=1∑ng(Xi)−b−a1∫bag(x)dx∣=1

即当nnn很大时，算数平均值近似等于统计平均值，即有：

1n∑i=1ng(Xi)≈1b−a∫bag(x)dx\frac{1}{n}\sum^{n}_{i=1}g(X_i)\approx\frac{1}{b-a}\int^a_bg(x)dxn1i=1∑ng(Xi)≈b−a1∫bag(x)dx

在实验中，我们用计算机生成大量在区间[0,1][0,1][0,1]上服从均匀分布的随机变量xxx，接着引入变换：

y=(b−a)x+ay=(b-a)x+ay=(b−a)x+a

使变量在区间[a,b][a,b][a,b]上非常均匀分布。于是我们得到：

1b−a∫bag(x)dx≈1n∑i=1ng(Yi)\frac{1}{b-a}\int^a_bg(x)dx\approx\frac{1}{n}\sum^{n}_{i=1}g(Y_i)b−a1∫bag(x)dx≈n1i=1∑ng(Yi)

3.实验过程及结果分析

本部分共有三个实验，第一个实验用蒙特卡洛法求解多种定积分，并于实际值进行比对，目的是研究蒙特卡洛法的可行性和准确性。第二个实验设置多个nnn值，探究nnn蒙特卡洛随机数数目对求解结果的影响。第三个实验探究积分区间对于所的积分值的精确度的影响。每个实验都对试验结果进行分析并得出最终结论。

实验一：蒙特卡洛法求解定积分

3.1.1 实验一的过程与结果

目标函数g(x)g(x)g(x)	区间[a,b][a,b][a,b]	手动求解所得值	蒙特卡洛法所得值	偏差程度	nnn值
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1093.58	0.4025%	10000
4x3+4cos(x)+4xsin(x)4x^3+4cos(x)+4xsin(x)4x3+4cos(x)+4xsin(x)	[1,10][1,10][1,10]	9963.2759	9822.21	1.4158%	10000
2xex2+1x2xe^{x^2}+\frac{1}{x}2xex2+x1	[1,10][1,10][1,10]	2.6881+e43	2.6196e+43	2.5504%	10000
ex2+ln(x)e^{x^2}+ln(x)ex2+ln(x)	[1,10][1,10][1,10]	-	1.24e+42	-	10000
11+x100\frac{1}{\sqrt{1+x^{100}}}1+x1001	[1,10][1,10][1,10]	-	0.02397	-	10000
15+sin(100x)3\frac{1}{\sqrt{5+sin(100x)^3}}5+sin(100x)31	[1,10][1,10][1,10]	-	4.04484	-	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1093.58	0.4025%	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1091.25	0.6146%	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1093.61	0.4000%	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1090.23	0.7073%	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1086.75	1.2046%	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1085.15	1.1697%	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1102 .12	0.3756%	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1088.55	0.8605%	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1098.06	0.0005%	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1099.88	0.1715%	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1093.86	0.3771%	10000

3.1.2 实验一的分析与结论

由上述实验结果我们可以看出，蒙特卡洛确实可以在一定允许的偏差范围内求出近似的定积分值。然而与实际值进行比较我们可以发现，蒙特卡洛法有一定的偏差，且由于随机变量是变化的，所以每次实验所得的偏差都不同。上一部分给出的图片形象的可视化了对于同一函数，每次蒙特卡洛法所得结果的不同和变化。然而，虽然每一次试验都有一定偏差且偏差不同，偏差的期望值却为0，所以我们可以多次进行蒙特卡罗方法，并将结果取平均值，这样会减小最终的误差，更好地近似所求的定积分值。
下图所示为对于函数g(x)=3x2+2g(x)=3x^2+2g(x)=3x2+2多次实验求得区间[1,10][1,10][1,10]积分值的折线图

下图所示为对于函数g(x)=3x2+2g(x)=3x^2+2g(x)=3x2+2多次实验求得区间[1,10][1,10][1,10]积分值偏差程度的折线图

实验二：n值对试验结果的影响

3.1.1 实验二的过程与结果

目标函数g(x)g(x)g(x)	区间[a,b][a,b][a,b]	手动求解所得值	蒙特卡洛法所得值	偏差程度	nnn值
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1087.09	0.9964%	100
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1108.00	0.9109%	1000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1093.58	0.4025%	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1097.37	0.0568%	100000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1098.61	0.05575%	1000000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1098.23	0.0216%	10000000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	1098	1098.23	0.0216%	100000000

3.1.2 实验二的分析与结论

由上述实验我们可以看出，当随机变量的数目n值不断增加的时候，所得积分值的偏差程度不断减小，即所得的积分值不断精确，最终与实际值几乎相等。然而，当n值非常大的时候，偏差程度减小的幅度也在不断减小，逐渐趋于稳定不变。

下图所示为对于函数g(x)=3x2+2g(x)=3x^2+2g(x)=3x2+2多次改变n值所得区间[1,10][1,10][1,10]积分值的折线图

下图所示为对于函数g(x)=3x2+2g(x)=3x^2+2g(x)=3x2+2多次改变n值所得区间[1,10][1,10][1,10]积分值偏差程度的折线图

3.2 实验三：积分区间的影响

3.2.1 实验三的过程与结果

目标函数g(x)g(x)g(x)	区间[a,b][a,b][a,b]	蒙特卡洛法求解偏差程度	nnn值
3x2+2x3x^2+2x3x2+2x	[1,80][1,80][1,80]	0.5643%	10000
3x2+2x3x^2+2x3x2+2x	[1,40][1,40][1,40]	0.2025%	10000
3x2+2x3x^2+2x3x2+2x	[1,20][1,20][1,20]	1.1966%	10000
3x2+2x3x^2+2x3x2+2x	[1,10][1,10][1,10]	0.4025%	10000
3x2+2x3x^2+2x3x2+2x	[1,5][1,5][1,5]	0.0297%	10000
3x2+2x3x^2+2x3x2+2x	[1,2.5][1,2.5][1,2.5]	0.3904%	10000
3x2+2x3x^2+2x3x2+2x	[1,1.125][1,1.125][1,1.125]	0.0225%	10000
3x2+2x3x^2+2x3x2+2x	[1,1+1/16][1,1+1/16][1,1+1/16]	0.4057%	10000
3x2+2x3x^2+2x3x2+2x	[1,1+1/32][1,1+1/32][1,1+1/32]	0.9623%	10000
3x2+2x3x^2+2x3x2+2x	[1,1+1/64][1,1+1/64][1,1+1/64]	0.5342%	10000

3.2.2 实验三的分析与结论

由上述实验可知，在随机变量数目不变的情况下，区间长度变化时与所得积分值的精度无关。这与我原来的直觉不同，我认为随着区间长度减小，所的积分值会上升且上升到一定程度时增加幅度减小并趋于稳定。而当区间长度增加时，精度不断降低，于是就能得出结论，当区间长度较长时，我们需要增加随机变量数目，即增大n的值。然后实际的试验结果显示，区间长度变化时与所得积分值的精度无明显联系,究其原因可能是因为随机变量数目相较而言已经足够近似出积分值，如果减少n值，可能会有和我直觉一致的结果
下两张图从10次试验所的图中随机挑选出，为对于函数g(x)=3x2+2g(x)=3x^2+2g(x)=3x2+2多次改变积分区间所得积分值偏差程度的折线图

4.代码(JuPyter NoteBook)

In(1)
import math
import random
def f(i,x):if i==0:return 3*x*x+2*xif i==1:return 4*(x**3)+4*math.cos(x)-4*x*math.sin(x)if i==2:return 2*x*math.exp(x**2)+1/xif i==3:return 1/math.sqrt(1+x**100)if i==4:return 1/math.sqrt(5+math.sin(x*100)**3)
def F(i,x):if i==0:return x**3+x**2if i==1:return x**4+4*x*math.cos(x)if i==2:return math.exp(x**2)+math.log(x)else:return x
a=1
b=10for i in range(5):sum=0count=1while count<=10000:sum=sum+f(i,random.uniform(a,b))count=count+1A=(b-a)*(sum/10000)B=F(i,b)-F(i,a)print(i)print("用蒙特卡洛方法计算的定积分：")print(A)print("直接用原函数计算的定积分：")print(B)d=abs(A-B)/Bprint("偏差程度为：")print('percent: {:.4%}'.format(d))In(2)
i=0
sum=0
count=1
t=10000
B=F(i,b)-F(i,a)
A=[]
D=[]
v=0
for m in range(10):sum=0count=1while count<=t:sum=sum+f(i,random.uniform(a,b))count=count+1v=(b-a)*(sum/t)A.append(v)D.append(abs(v-B)/B)
print(A)
print(D)In(3)
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import FuncFormatter
x= np.arange(0,10,0.1)
plt.figure('Output')
plt.plot(D,label='Output',color='b')
#plt.plot(x,1098+x*0,color='red')
def to_percent(temp, position):return '%.2f'%(100*temp) + '%'
plt.gca().yaxis.set_major_formatter(FuncFormatter(to_percent))
plt.show()In（4）
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import FuncFormatter
x=[pow(10,i) for i in range(0,7)]
plt.figure('Output')
plt.plot(x,D,label='Output',color='b')
#plt.plot(x,1098+x*0,color='red')
def to_percent(temp, position):return '%.2f'%(100*temp) + '%'
plt.gca().yaxis.set_major_formatter(FuncFormatter(to_percent))
plt.xscale('log')
plt.show()
plt.figure('Output')
y= np.arange(0,100000000)
x=[pow(10,i) for i in range(0,7)]
plt.plot(x,A,label='Output',color='b')
plt.plot(y,1098+y*0,color='red')
plt.xscale('log')#设置纵坐标的缩放
plt.show()