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

分手后的思念是犯贱 2023-02-26 02:18 107阅读 0赞

## 一、提出问题： ##

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

## 二、解决方法 ##

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

**1、蒙特卡洛 (Monte Carlo) 求定积分概述**  
蒙特卡洛方法也称统计模拟方法、随机抽样技术,是基于“随机数”、概率统计理论为基础的数值计算方法。蒙特卡洛定积分主要思想就是均匀分布生成的随机数,将积分符号转化为求和，从而实现快速求解目的。主要有三种方法：随机投点法、平均值法、重要抽样法。  
**2、平均值法求定积分**  
计算过程如图1  
![在这里插入图片描述][20200713195436223.png]  
其数学公式为：  
![在这里插入图片描述][20200713195511216.png]  
**3、定积分值误差检验**  
方根误差（Root Square Error）是观测值与真值（理论值）偏差平方的平方根，是用来衡量观测值同真值之间的偏差。如图2，随采样数量的增加，误差逐渐降低。在达到一定数量（200）次采样之后，误差变化不大。说明采样次数提高到一定值后，对运算精度变化影响较少。进而采取改变抽样法来减小方差，从而提高积分计算的精度。（此不在本次运演范围）  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ1MTYxNzY2_size_16_color_FFFFFF_t_70]

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

**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、运行结果**  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ1MTYxNzY2_size_16_color_FFFFFF_t_70 1]

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

[20200713195436223.png]: /images/20230209/1103952f753448c89fc2752fa09407bd.png
[20200713195511216.png]: https://img-blog.csdnimg.cn/20200713195511216.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ1MTYxNzY2_size_16_color_FFFFFF_t_70]: https://img-blog.csdnimg.cn/20200713195651828.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ1MTYxNzY2,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ1MTYxNzY2_size_16_color_FFFFFF_t_70 1]: https://img-blog.csdnimg.cn/20200713195846703.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ1MTYxNzY2,size_16,color_FFFFFF,t_70