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

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

0.摘要

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

1.蒙特卡洛法简介

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

$$s=\sum _{x}p(x)f(x)=E_p[f(x)]$$

或者

$$s=\int p(x)f(x)=E_p[f(x)]$$

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

$$\widehat{s_n}=\frac{1}{n}\sum _{i=1}^{n}f(x^{(i)})$$

并且有

$$E[\widehat{s_n}]=\frac{1}{n}\sum _{i=1}^{n}E[f(x^{(i)})]=\frac{1}{n}\sum _{i=1}^{n}s=s$$

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

$$\underset{n\to +\infty }{\lim }\widehat{s_n}=s$$

2.公式推导

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

$$E[g(x)]=\frac{1}{b-a}{\int }_b^{a}g(x)dx$$

由大数定理：

$$\underset{n\to +\infty }{\lim }P|\frac{1}{n}\sum _{i=1}^{n}g(X_i)-\frac{1}{b-a}{\int }_b^{a}g(x)dx|=1$$

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

$$\frac{1}{n}\sum _{i=1}^{n}g(X_i)\approx \frac{1}{b-a}{\int }_b^{a}g(x)dx$$

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

$$y=(b-a)x+a$$

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

$$\frac{1}{b-a}{\int }_b^{a}g(x)dx\approx \frac{1}{n}\sum _{i=1}^{n}g(Y_i)$$

3.实验过程及结果分析

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

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

3.1.1 实验一的过程与结果

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

3.1.2 实验一的分析与结论

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

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

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

3.1.1 实验二的过程与结果

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

3.1.2 实验二的分析与结论

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

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

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

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

3.2.1 实验三的过程与结果

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

3.2.2 实验三的分析与结论

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

4.代码(JuPyter NoteBook)

In(1)
import math
import random
def f(i,x):
    if i==0:
        return 3*x*x+2*x
    if 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/x
    if 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**2
    if 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=10


 
for i in range(5):
    sum=0
    count=1
    while count<=10000:
        sum=sum+f(i,random.uniform(a,b))
        count=count+1
    A=(b-a)*(sum/10000)
    B=F(i,b)-F(i,a)
    print(i)
    print("用蒙特卡洛方法计算的定积分：")
    print(A)
    print("直接用原函数计算的定积分：")
    print(B)
    d=abs(A-B)/B
    print("偏差程度为：")
    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=0
    count=1
    while count<=t:
        sum=sum+f(i,random.uniform(a,b))
        count=count+1
        v=(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()