数值概率算法和蒙特卡罗算法

设计伪随机数发生器并数值近似求解

运用python编程,实现了一个伪随机数发生器[1],并计算的近似值以及所给积分的近似值,具体结果如图1。横坐标表示随机投点数,上图纵坐标表示对应的真实值,下图纵坐标表示所需的时间。从图1(左)中分析得到,当随机投点数为pow(10,3)时,使用该伪随机数发生器计算的近似值已经比较接近圆周率pi的真实值,并且从时间消耗图来看,对应的时间消耗并不大;从图1(右)中分析得到,当随机投点数为pow(10,4)时,其所计算的近似值已经接近给定积分(详见代码部分)的真实值,所对应的时间消耗也不大。从图1左右两图均可看出,只有当随机投点数大于等于pow(10,6)时,其时间消耗开始大幅增加。
在这里插入图片描述在这里插入图片描述
图1 圆周率pi与给定积分的近似值以及对应的时间消耗

python代码如下:

import time
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

seed = 65536
start, end = 0, int(1e8)

def Rand(start, end, seed=9999):
    a = 32310901
    b = 1729
    r_old = seed
    m = end - start
    while True:
        r_new = int((a * r_old + b) % m)
        yield r_new
        r_old = r_new

r = Rand(start, end, seed)

def pi(n):
    sum = 0
    for i in range(n):
        x = next(r) / end
        y = next(r) / end
        if (x * x + y * y <= 1):
            sum = sum + 1
    return 4 * sum / n

def fun(x):
    return 3 * x ** 2 - 4 * x + 3

def definite_integration(n):
    a, b = 1, 4
    data = fun(np.linspace(a, b, num=1000, endpoint=True))
    M, L = int(max(data)), int(min(data))
    # print(M,L)
    c, d = (M - L) * (b - a), L * (b - a)
    sum = 0
    for i in range(n):
        x = next(r) / end
        y = next(r) / end
        fun_new = 1 / (M - L) * (fun(a + (b - a) * x) - L)
        if (y <= fun_new):
            sum = sum + 1
    return c * (sum / n) + d

def figure_1():
    t0, result0 = [], []
    numbers = [pow(10, i) for i in range(1, 8)]
    for num in numbers:
        start = time.clock()
        result0.append(pi(num))
        t0.append(time.clock() - start)

    plt.figure()
    plt.subplot(211)
    plt.hlines(np.pi, 1, len(numbers), colors='b')
    plt.plot(list(range(1, len(t0) + 1)), result0, linestyle='--', marker='*', color='r')
    # plt.xlabel('随机投点数(10^x)')
    plt.ylabel('数值计算值')

    plt.subplot(212)
    plt.plot(list(range(1, len(t0) + 1)), t0)
    plt.xlabel('随机投点数(10^x)')
    plt.ylabel('所需要的时间(秒)')
    plt.savefig("result_1.png")
    # plt.show()

def figure_2():
    t1, result1 = [], []
    numbers = [pow(10, i) for i in range(1, 8)]
    for num in numbers:
        start = time.clock()
        result1.append(definite_integration(num))
        t1.append(time.clock() - start)

    plt.figure()
    plt.subplot(211)
    plt.hlines(42, 1, len(numbers), colors='b')
    plt.plot(list(range(1, len(t1) + 1)), result1, linestyle='--', marker='*', color='r')
    # plt.xlabel('随机投点数(10^x)')
    plt.ylabel('数值计算值')

    plt.subplot(212)
    plt.plot(list(range(1, len(t1) + 1)), t1)
    plt.xlabel('随机投点数(10^x)')
    plt.ylabel('所需要的时间(秒)')
    plt.savefig('result_2.png')
    # plt.show()

if __name__ == '__main__':
    figure_1()
    figure_2()

参考文献:

[1] Retrieved May 10, 2019, from https://blog.csdn.net/qq_39022311/article/details/83477469.

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值