欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  IT编程

python求积分面积的几个方法

程序员文章站 2022-09-21 08:58:43
示例已知积分公式如下求[0.5,5]上积分,即求下图阴影部分面积根据积分公式求源函数等于:则确切解等与F(5)-F(0.5)=3.9002072872864524当不知道源函数时使用以下方法可以求得积分面积首先定义函数def func(x): return np.cos(np.pi) * np.exp(-x) + 11、通过迭代计算每个步长矩形面积求和 x = np.linspace(0.5, 5, 1000) dx = (5-...

示例

已知积分公式如下

python求积分面积的几个方法

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

python求积分面积的几个方法

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

python求积分面积的几个方法

则确切解等与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、使用蒙特卡洛模拟估算积分面积

python求积分面积的几个方法

原理是已知区间【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

相关标签: python