Python 实现定积分与二重定积分的操作

1.概述

最近项目需要使用程序实现数学微积分，最初想用java实现，后来发现可用文档太少，实现比较麻烦，后来尝试使用python实现，代码量较少，主要有sympy与scipy两种实现方式，本文主要记录scipy的实现方式。

2.内容

2.1 所求函数

2.2 python代码

# 引入需要的包
import scipy.integrate
from numpy import exp
from math import sqrt
import math

# 创建表达式
f = lambda x,y : exp(x**2-y**2)

# 计算二重积分：（p:积分值，err:误差）
# 这里注意积分区间的顺序
# 第二重积分的区间参数要以函数的形式传入
p,err= scipy.integrate.dblquad(f, 0, 2, lambda g : 0, lambda h : 1)
print(p)

2.3 注意问题

1. exp尽量使用numpy的exp

2. 注意积分区间参数的顺序

3. 第二重积分的区间参数要以函数的形式传入

补充：python实现求解积分

例子 1：

假设有随机变量 x，定义域 X，其概率密度函数为 p(x)，f(x) 为定义在 X 上的函数，目标是求函数 f(x) 关于密度函数 p(x) 的数学期望 。

蒙特卡洛法根据概率分布 p(x) 独立地抽样 n 个样本 x1,x2,…..xn，得到近似的 f(x) 期望为：

其实这个的理解就是要求一个拥有概率密度的函数期望值

期望=积分（每个点的密度函数*每个点的价值函数）

例子 2：

假设我们想要求解 h(x) 在 X 上的积分：

我们将 h(x) 分解成一个函数 f(x) 和一个概率密度函数 p(x) 的乘积，进而又将问题转换为求解函数 f(x) 关于密度函数 p(x) 的数学期望 ：

这里的Ep(x)是相当于把整个分布当时了概率分布，即总发生概率为1.

这里，f(x) 表示为 ，则有：

更一般的，假设我们想要求解 ，熟悉积分的同学肯定已经知道答案为 ，那么如何用采样的方法来得到这个值呢？

令 ，0<x<10，那么 。

下面是代码：

'''import random
num=1000000
sum=0
for i in range(0,num):
    x=random.uniform(0,10)
    sum+=x*x*10
sum/=1000000
print(sum)'''
import random
numSamples=10000
samples=[random.uniform(0,10)for _ in range(numSamples)]
f_samples=[10*sample**2 for sample in samples]
result=1/10000.0*sum(f_samples)
print(result)

result=333.10527012455066

random.uniform(x,y)表示在[x,y)之间生成一个 实数

对于复杂的 h(x)，这种方法计算起来显然就更加方便了（特别是忘记积分怎么算的同学）。

蒙特卡洛方法其实就是利用大数定理通过大量统计来算出最后的值。

到这里为止，我们简单的介绍了蒙特卡洛方法，但是依旧没有提到要怎么利用复杂的概率密度函数进行采样。

接下来我们来看一下接受-拒绝法（accept-reject sampling method），它也是蒙特卡洛法中的一种类型适用于不能直接抽样的情况。

以上为个人经验，希望能给大家一个参考，也希望大家多多支持我们。

赞 (0)

相关推荐

  • python 求定积分和不定积分示例

    求f(x) = sin(x)/x 的不定积分和负无穷到正无穷的定积分 sin(x)/x 的不定积分是信号函数sig ,负无穷到正无穷的定积分为pi import math import numpy as np import matplotlib.pyplot as plt from sympy import * #用于求导积分等科学计算 def draw_plot_set():#设置画图格式 ax = plt.gca() #改变坐标轴位置 ax.spines['right'].set_color

  • python、Matlab求定积分的实现

    python求定积分 计算 from sympy import * x = symbols('x') print(integrate(sin(2*x)/(1+x**2), (x, 0, 3))) sympy库中integrate函数 integrate(f, (x, lower_bound, upper_bound)) # f-函数,x-变量,lower_bound-下限,upper_bound-上限 但是发现求不出来,如果是sin(2*x)就可以,为什么? matlab求定积分 syms x

  • python编程通过蒙特卡洛法计算定积分详解

    想当初,考研的时候要是知道有这么个好东西,计算定积分...开玩笑,那时候计算定积分根本没有这么简单的.但这确实给我打开了一种思路,用编程语言去解决更多更复杂的数学问题.下面进入正题. 如上图所示,计算区间[a b]上f(x)的积分即求曲线与X轴围成红色区域的面积.下面使用蒙特卡洛法计算区间[2 3]上的定积分:∫(x2+4*x*sin(x))dx # -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt

  • Python 实现定积分与二重定积分的操作

    Python 实现定积分与二重定积分的操作

    1.概述 最近项目需要使用程序实现数学微积分,最初想用java实现,后来发现可用文档太少,实现比较麻烦,后来尝试使用python实现,代码量较少,主要有sympy与scipy两种实现方式,本文主要记录scipy的实现方式. 2.内容 2.1 所求函数 2.2 python代码 # 引入需要的包 import scipy.integrate from numpy import exp from math import sqrt import math # 创建表达式 f = lambda x,y

  • Python实现对PPT文件进行截图操作的方法

    本文实例讲述了Python实现对PPT文件进行截图操作的方法.分享给大家供大家参考.具体分析如下: 下面的代码可以为powerpoint文件ppt进行截图,可以指定要截取的幻灯片页面,需要本机安装了powerpoint,可以指定截图的大小分辨率 import os import comtypes.client def export_presentation(path_to_ppt, path_to_folder): if not (os.path.isfile(path_to_ppt) and

  • Python文件夹与文件的相关操作(推荐)

    最近在写的程序频繁地与文件操作打交道,这块比较弱,还好在百度上找到一篇不错的文章,这是原文传送门,我对原文稍做了些改动. 有关文件夹与文件的查找,删除等功能 在 os 模块中实现.使用时需先导入这个模块, 导入的方法是: import os 一.取得当前目录 s = os.getcwd() # s 中保存的是当前目录(即文件夹) 比如运行abc.py,那么输入该命令就会返回abc所在的文件夹位置. 举个简单例子,我们将abc.py放入A文件夹.并且希望不管将A文件夹放在硬盘的哪个位置,都可以在A

  • python条件变量之生产者与消费者操作实例分析

    python条件变量之生产者与消费者操作实例分析

    本文实例讲述了python条件变量之生产者与消费者操作.分享给大家供大家参考,具体如下: 互斥锁是最简单的线程同步机制,面对复杂线程同步问题,Python还提供了Condition对象.Condition被称为条件变量,除了提供与Lock类似的acquire和release方法外,还提供了wait和notify方法.线程首先acquire一个条件变量,然后判断一些条件.如果条件不满足则wait:如果条件满足,进行一些处理改变条件后,通过notify方法通知其他线程,其他处于wait状态的线程接到

  • Python双精度浮点数运算并分行显示操作示例

    Python双精度浮点数运算并分行显示操作示例

    本文实例讲述了Python双精度浮点数运算并分行显示操作.分享给大家供大家参考,具体如下: #coding=utf8 def doubleType(): ''''' Python中的浮点数是双精度浮点数,可以用十进制或科学计数法表示. 实际精度依赖于机器架构和创建Python解释器的编译器. 浮点数值通常都有一个小数点和一个可选的后缀e(大写或小写,表示科学计数法). 在e和指数之间可以用正(+)或负(-)表示指数的正负(正数可以省略符号) ''' (one,two,three,four,fiv

  • Python中列表与元组的乘法操作示例

    本文实例讲述了Python中列表与元组的乘法操作.分享给大家供大家参考,具体如下: 直接上code吧,还可以这么玩儿 列表乘法: li=[1,] li=li*3 print(li) out: [1, 1, 1] 元组乘法: >>> t=(1,2) >>> t*3 (1, 2, 1, 2, 1, 2) 但字典,集合不能这么玩 例如: >>> dict1={'k1':1,'k2':2} >>> dict1*2 #报错 Traceback

  • Python中文件的读取和写入操作

    从文件中读取数据 读取整个文件 这里假设在当前目录下有一个文件名为'pi_digits.txt'的文本文件,里面的数据如下: 3.1415926535 8979323846 2643383279 with open('pi_digits.txt') as f: # 默认模式为'r',只读模式 contents = f.read() # 读取文件全部内容 print contents # 输出时在最后会多出一行(read()函数到达文件末会返回一个空字符,显示出空字符就是一个空行) print '

  • Python实现的列表排序、反转操作示例

    本文实例讲述了Python实现的列表排序.反转操作.分享给大家供大家参考,具体如下: 排序: 使用sorted方法和列表的sort方法: sorted方法适用范围更广,sort方法只有列表有. li = [{'a':'23'}, {'a':'12'}] def sort_fun(mp): s = mp['a'] return int(s) print(sorted(li, key = sort_fun, reverse = True)) #这会返回一个排好序的列表,原列表不变. print(li

  • Python反射和内置方法重写操作详解

    本文实例讲述了Python反射和内置方法重写操作.分享给大家供大家参考,具体如下: isinstance和issubclass isinstance(obj,cls)检查是否obj是否是类 cls 的对象,类似 type() class Foo(object): pass obj = Foo() isinstance(obj, Foo) issubclass(sub, super)检查sub类是否是 super 类的派生类 class Foo(object): pass class Bar(Fo

随机推荐