使用蒙塔卡罗方法计算定积分（Python）

目录

蒙特卡罗法

方法一：矩形法

方法二：面积法

---

蒙特卡罗法

        蒙特卡罗法也成为统计模拟法、统计试验法，是把概率现象作为研究对象的数值模拟方法。本文讲述使用蒙特卡罗法计算定积分的两种方法。以y=2x+1为例，求其在区间（0,1）上的定积分。

方法一：矩形法

        在高数课上，老师肯定讲过这一种方法，将函数区域分成若干个矩形（见下图），求这些矩形的面积之和即可得出函数在这片区域的定积分。

假如矩形的宽很小很小，那么矩形的高就可近似认为时f(x1)，其面积就为f(x1)*(x2-x1)，那么所有矩形的面积和为，其中n是矩形的数量，n越大，精度越高。

用Python实现：

import numpy as np


def f(x):
    # 被积函数
    return x * 2 + 1


def Monte_Carlo_Methond(a, b):
    # (a,b)是积分区间
    x = np.random.uniform(a, b, 1000000)
    y = f(x)
    sum = 0
    for i in y:
        sum += i

    return sum / 1000000


print(Monte_Carlo_Methond(0, 1))

方法二：面积法

        这个方法就是我们向一个包围函数的矩形（下图中的黑色矩形）上面投点，然后求出落在函数区域（下图中的蓝线围起来的部分）上的点所占的比例，这个比例再乘以矩形的面积即可得到函数在这个区间上的定积分。

用Python实现：

import random


def f(x):
    # 被积函数
    return x * 2 + 1


def Monte_Carlo_Methond(a, b):
    # (a,b)是积分区间
    sum = 0
    y1 = f(a)
    y2 = f(b)
    # y1,y2是函数的上下界
    for i in range(1000000):
        x = random.uniform(a, b)
        y = random.uniform(0, y2)
        # (x,y)是在矩形内随机生成的点
        if y <= f(x):
            sum += 1
    return (sum / 1000000) * (abs(b - a)) * (abs(max(y1, y2)))


print(Monte_Carlo_Methond(0, 1))