1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

#! /usr/bin/env python3
# vim:fenc=utf-8
import numpy as np
import matplotlib.pyplot as plt

# 定义蒙特卡洛积分函数
def monte_carlo_integration(func, a, b, n=1000):
    x_samples = np.random.uniform(a, b, n)
    y_samples = func(x_samples)
    # 计算平均高度
    average_height = np.mean(y_samples)
    # 计算积分值
    integral = average_height * (b - a)
    return integral

# 定义被积函数
def f(x):
    return x**2

# 解析解计算积分函数
def analytic_integral_f(x,a):
    return (x**3)/3-(a**3)/3

# 使用蒙特卡洛方法求解积分
a, b = -5, 5  # 积分区间
x_points = np.linspace(a, b, 100)  # 用于绘图的x点
monte_carlo_values = []

for i in range(len(x_points)):
    # 对每个x点计算从a到该点的积分
    integral_value = monte_carlo_integration(f, a, x_points[i], n=20000)
    monte_carlo_values.append(integral_value)

# 计算解析解
analytic_values = analytic_integral_f(x_points,a)

# 绘制积分后的函数图像
plt.figure(figsize=(10, 8))
plt.plot(x_points, monte_carlo_values, label='Monte Carlo Integral of $x^2$', color='blue', linestyle='--')
plt.plot(x_points, analytic_values, label='Analytical Integral of $x^2$', color='red')

# 添加标题和标签
plt.title('Comparison between Monte Carlo and Analytical Integration of $x^2$')
plt.xlabel('x')
plt.ylabel('Integral Value')

# 显示图例
plt.legend()

# 显示网格
plt.grid(True)

# 显示图像
plt.show()