示例

已知积分公式如下

求[0.5，5]上积分，即求下图阴影部分面积

根据积分公式求源函数等于：

则确切解等与F(5)-F(0.5)=3.9002072872864524

当不知道源函数时使用以下方法可以求得积分面积

首先定义函数

def func(x):

return np.cos(np.pi) * np.exp(-x) + 1

1、通过迭代计算每个步长矩形面积求和

x = np.linspace(0.5, 5, 1000)

dx = (5-0.5)/1000

y = func(x)

area = np.sum(y * dx)

结果等于：3.8994262124710404

2、通过scipy模块quad直接求积分

area, error = quad(func, 0.5, 5)

quad的结果是 3.9002072872864515

3、使用蒙特卡洛模拟估算积分面积

原理是已知区间【0.5,5】上矩形面积，在矩形范围内随机打点，统计在阴影部分点的占比，从而估算出阴影部分面积。

k = 0

for i in range(N):

r_x = np.random.uniform(low=0.5, high=5, size=1)

r_y = np.random.uniform(low=0, high=1, size=1)

c_v = func(r_x)

if r_y <= c_v:

k += 1

area = k / N * 4.5

蒙特卡洛的结果是 3.9393

所有代码如下

# @Time : 2020/8/26

# @Author : 大太阳小白

# @Software: PyCharm

# @blog：https://blog.csdn.net/weixin_41579863

import numpy as np

from scipy.integrate import quad

from matplotlib import pyplot as plt

def func(x):

return np.cos(np.pi) * np.exp(-x) + 1

def drawing_question(fig):

x = np.linspace(0, 6, 1000)

y = func(x)

fig = plt.figure()

ax1 = fig.add_subplot()

ax1.axis([np.min(x), np.max(x), 0, np.max(y)]) # 坐标范围

ax1.plot(x, y, label="$cos(π)e^{-x}+1$") # 画曲线，带图示

ax1.fill_between(x, y1=y, y2=0, where=(x >= 0.5) & (x <= 5), # 填充积分区域

facecolor='blue', alpha=0.2)

ax1.text(0.5 * (0.7 + 4), 0.4, r"$\int_{0.5}^5(cos(π)e^{-x}+1)\mathrm{d}x$",

horizontalalignment='center', fontsize=14) # 增加说明文本

ax1.legend() # 显示图示

plt.show()

def area_sum():

x = np.linspace(0.5, 5, 1000)

dx = (5-0.5)/1000

y = func(x)

area = np.sum(y * dx)

return area

def area_quad():

area, error = quad(func, 0.5, 5)

return area

def montecarlo(N=5000):

k = 0

for i in range(N):

r_x = np.random.uniform(low=0.5, high=5, size=1)

r_y = np.random.uniform(low=0, high=1, size=1)

c_v = func(r_x)

if r_y <= c_v:

k += 1

area = k / N * 4.5

return area

def montecarlo_show(N=5000):

k = 0

x0 = []

x1 = []

for i in range(N):

r_x = np.random.uniform(low=0.5, high=5, size=1)

r_y = np.random.uniform(low=0, high=1, size=1)

c_v = func(r_x)

if r_y <= c_v:

k += 1

x0.append([r_x, r_y])

else:

x1.append([r_x, r_y])

x = np.linspace(0, 6, 1000)

y = func(x)

x0 = np.array(x0)

x1 = np.array(x1)

fig = plt.figure()

ax1 = fig.add_subplot()

ax1.axis([np.min(x), np.max(x), 0, np.max(y)])

ax1.plot(x, y, label="$cos(π)e^{-x}+1$")

ax1.scatter(x0[:, 0], x0[:, 1])

ax1.scatter(x1[:, 0], x1[:, 1])

ax1.fill_between(x, y1=y, y2=0, where=(x >= 0.5) & (x <= 5),

facecolor='blue', alpha=0.2)

ax1.text(0.5 * (0.7 + 4), 0.4, r"$\int_{0.5}^5(cos(π)e^{-x}+1)\mathrm{d}x$",

horizontalalignment='center', fontsize=14) # 增加说明文本

ax1.legend() # 显示图示

plt.show()

if __name__ == '__main__':

# drawing_question()

original = lambda x: np.exp(-x) + x

print('期望值是：', original(5)-original(0.5))

print('基于矩形面积求和的结果是', area_sum())

print('quad的结果是', area_quad())

print('蒙特卡洛的结果是', montecarlo())

montecarlo_show(N=500)

本文地址：https://blog.csdn.net/weixin_41579863/article/details/108239755

希望与广大网友互动？？

点此进行留言吧！