蒙特卡洛方法求定积分

#依旧是dsd的概率论大作业

一、题目描述

1、蒙特卡洛方法：也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。布丰投针实验是蒙特卡洛方法的起源。

蒙特卡洛方法主要是构造或描述一个概率过程，对于本身具有随机行的问题，正确描述过程，对于不是随机过程的问题，比如本次要计算的定积分，就必须认为构造一个概率过程，它的某些参量正好是所要求问题的解。

随机数是我们运用蒙特卡洛方法的基础，虽然一般利用计算机产生的是伪随机数，但经过统计检验，其与真正随机数具有相近的性质，可以用作真实的随机数。

2、问题描述：

题干问题要求用蒙特卡洛方法求定积分，也就是求解：

用微积分思想可知，求定积分就是求函数所围成的图形下的面积。利用这个方式构建一个随机过程。

假设向坐标轴上图示长方形区域（即积分区域）扔小点，小点落在随机的位置，大量随机试验之后落在长方形内函数曲线之下区域的点与总点数之比，可以近似地看做是积分区域与长方形区域的面积之比，而长方形面积易得，故可以得到定积分的值。

即定积分面积=阴影部分的随机数与正方形面积中总的随机数之比。

由以上思想我们可以编写代码，进行大量随机试验。

二、代码

利用python实现上述过程。

假设原函数为gx=x3+2x

则导函数为fx=3x2+2

实际应用过程中可以根据需求更改函数表达式，代码如下： 

import random  # 引入包
import pyecharts.options as opts
from pyecharts.charts import Scatter


def g(x):  # 定义原函数

    return x * x * x + 2 * x


def func(x):  # 定义导函数

    return 3 * x * x + 2


data1 = []
data2 = []


def monte_carlo_method(a, b):  # 定义蒙特卡洛方法，ab为积分区间

    Sum: float = 0
    y1 = func(a)
    y2 = func(b)
    count = 1000000
    global data1
    global data2
    # y1,y2是函数的上下界
    for i in range(count):
        x = random.uniform(a, b)
        y = random.uniform(0, max(y1, y2))

        # (x,y)是在矩形内随机生成的点
        if y < func(x):  # 判断点是否在积分区间内
            Sum += 1
            data1.append([x, y])
        else:
            data2.append([x, y])  # 将新生成的不属于积分面积的点加入列表2中

    return (Sum / count) * (abs(b - a)) * (abs(max(y1, y2)))


def primitive_function_method(a, b):  # 定义利用原函数计算积分

    return g(b) - g(a)


print("Result calculated by monte-carlo method:")
print(monte_carlo_method(0, 1))
print("Result calculated by primitive function method:")
print(primitive_function_method(0, 1))

x1_data = [d[0] for d in data1]  # 设置xy坐标值
x2_data = [d[0] for d in data2]
y1_data = [d[1] for d in data1]
y2_data = [d[1] for d in data2]

c = (
    Scatter()
    .add_xaxis(xaxis_data=x1_data)  # 设置x坐标
    .add_yaxis(  # 设置y坐标
        series_name="积分区间内",
        y_axis=y1_data,
        symbol_size=2,
        label_opts=opts.LabelOpts(is_show=False),
    )
    .add_xaxis(xaxis_data=x2_data)
    .add_yaxis(
        series_name="积分区间外",
        y_axis=y2_data,
        symbol_size=2,
        label_opts=opts.LabelOpts(is_show=False),
    )

    .set_series_opts()
    .set_global_opts(
        xaxis_opts=opts.AxisOpts(
            type_="value", splitline_opts=opts.SplitLineOpts(is_show=True)
        ),
        yaxis_opts=opts.AxisOpts(
            type_="value",
            axistick_opts=opts.AxisTickOpts(is_show=True),
            splitline_opts=opts.SplitLineOpts(is_show=True),
        ),
        tooltip_opts=opts.TooltipOpts(is_show=False),

    )
    .render("basic_scatter_chart.html")
)

代码中利用了百度的pyecharts库进行图形化绘制。

 三、图形化结果