蒙特卡法洛求π

蒙特·卡罗方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。与它对应的是确定性算法。蒙特·卡罗方法在金融工程学，宏观经济学，计算物理学（如粒子输运计算、量子热力学计算、空气动力学计算）等领域应用广泛

用蒙特卡法洛求π

先贴代码

import math
import random

# 正方形的面积: 2R * 2R = 4 * R * R 
# 内接圆的面积: π * R * R
# π = 内接圆的面积 / 正方形的面积 * 4

circle_inside = 0
circle_outside = 0
count = 0
while True:
    count += 1
    x = random.random() * 2
    y = random.random() * 2
    if x ** 2 + y ** 2 > 4:
        circle_outside += 1
    else:
        circle_inside += 1
    if circle_outside != 0:
        pi = circle_inside / (circle_inside + circle_outside) * 4
        print("π约为:{}".format(pi))
        diff = math.pi - pi
        if abs(diff) < 0.0000000001:
            print("一共进行了{}次模拟，π的值约为:{},误差值为:{}".format(count, pi, diff))
            break

代码运行结果：

....................省略......................
π约为:3.141592516653196
π约为:3.1415923024323424
π约为:3.1415923609659737
π约为:3.1415924194995966
π约为:3.1415924780332114
π约为:3.1415925365668187
π约为:3.1415925951004176
π约为:3.141592653634009
一共进行了14665210次模拟，π的值约为:3.141592653634009,误差值为:-4.4215742178721484e-11

为了方便更深刻的理解和记忆，下面用可视化程序演示

import math
import random

# 正方形的面积: 2R * 2R = 4 * R * R
# 内接圆的面积: π * R * R
# π = 内接圆的面积 / 正方形的面积 * 4
import pygame

pygame.init()
width = 800
height = 800
width_half = width // 2
height_half = height // 2
screen = pygame.display.set_mode((width, height))

circle_inside = 0
circle_outside = 0
count = 0
dots = []
running = True
while running:
    events = pygame.event.get()
    for event in events:
        if event.type == pygame.QUIT:
            running = False
            break
    count += 1
    x = random.randint(0, width)
    y = random.randint(0, height)

    b = (x - width_half) ** 2 + (y - height_half) ** 2 > height_half ** 2
    if b:
        circle_outside += 1
    else:
        circle_inside += 1
    dots.append(((x, y), b))
    if circle_outside != 0:
        pi = circle_inside / (circle_inside + circle_outside) * 4
        diff = math.pi - pi
        if abs(diff) < 0.000001:
            print("一共进行了{}次模拟，π的值约为:{},误差值为:{}".format(count, pi, diff))
            break

    screen.fill((255, 255, 255))
    pygame.draw.circle(screen, (255, 255, 0), (width_half, height_half), width_half, 1)
    for dot, bt in dots:
        pygame.draw.circle(screen, (255, 0, 0) if not bt else (0, 0, 0), dot, 1)
    pygame.display.update()

pygame.quit()

图片