python如何利用蒙特卡洛平均值法求定积分

一、提出问题：

python中如何利用蒙特卡洛平均值法求定积分？

二、解决方法

（1）基本理论与操作说明

1、蒙特卡洛 (Monte Carlo) 求定积分概述
蒙特卡洛方法也称统计模拟方法、随机抽样技术,是基于“随机数”、概率统计理论为基础的数值计算方法。蒙特卡洛定积分主要思想就是均匀分布生成的随机数,将积分符号转化为求和，从而实现快速求解目的。主要有三种方法：随机投点法、平均值法、重要抽样法。
2、平均值法求定积分
计算过程如图1

其数学公式为：

3、定积分值误差检验
方根误差（Root Square Error）是观测值与真值（理论值）偏差平方的平方根，是用来衡量观测值同真值之间的偏差。如图2，随采样数量的增加，误差逐渐降低。在达到一定数量（200）次采样之后，误差变化不大。说明采样次数提高到一定值后，对运算精度变化影响较少。进而采取改变抽样法来减小方差，从而提高积分计算的精度。（此不在本次运演范围）

三、举例说明蒙特卡洛算法

1、要求：根据上述基本原理，用python程序验证蒙特卡洛求定积分计算过程
2、具体代码

import numpy as np
import pylab as pl

def fy(x):
    return x**3 + 4*x*np.cos(x)


def dfy(x):
    return 3*x**2+4*np.cos(x)-4*x*np.sin(x)

def monto(x,a,b):
    return (b-a)/len(x)*sum(dfy(x))

def integral(a,b):
    return fy(b) - fy(a)

a,b,n=10,15,1000
vm=np.ones(n)
X=np.linspace(a,b,n)
y=np.ones(n)
vi=integral(a,b)
for i in range(n):
    x = np.random.uniform(a,b,i+1)
    vm[i] = monto(x,a,b)
    y[i] = dfy(X[i])
    
pl.subplot(212)
pl.title(" RSE ")
pl.plot(range(n), np.sqrt((vi-vm)**2) )
pl.xlabel("Monte Carlo sampling point")

pl.subplot(211)
pl.title("function curve")
pl.plot(X,y,label='$y=3*x^2+4*cos(x)-4*x*sin(x)$')
pl.legend(loc='upper left')
pl.show()

3、运行结果

以上就是本文所有内容，希望能帮到大家！！！