【问题标题】:How to do a (trapeze) integration in Python with x^2?如何在 Python 中使用 x^2 进行(空中飞人)集成?
【发布时间】:2020-10-11 13:25:11
【问题描述】:

我的任务是首先与 f(x)=x^2

的 Python 进行集成,然后进行梯形集成
import numpy as np
import matplotlib.pyplot as plt

x = np.arange(-10,10)   
y = x**2

l=plt.plot(x,y)
plt.show(l)

现在我想整合这个函数得到这个:F(x)=(1/3)x^3和图片:

这应该是最后的输出:

有人可以解释一下如何用 python 获得 f(x)=x^2 的反导数 F(x) 吗? 我想通过普通积分和空中飞人积分来做到这一点。对于从(-10 到 10)的梯形积分,步长为 0.01(梯形的宽度)。最后我想在这两种情况下都得到函数 F(x)=(1/3)x^3 。我怎样才能做到这一点?

谢谢你帮助我。

【问题讨论】:

标签: python numpy matplotlib math integration


【解决方案1】:

有两个关键观察:

  • 梯形法则是指数值积分,其输出不是积分函数而是一个数
  • 积分取决于您对F(x) 的定义中未包含的任意常数

考虑到这一点,您可以使用scipy.integrate.trapz() 来定义积分函数:

import numpy as np
from scipy.integrate import trapz


def numeric_integral(x, f, c=0):
    return np.array([sp.integrate.trapz(f(x[:i]), x[:i]) for i in range(len(x))]) + c

或者,更有效的是,使用scipy.integrate.cumtrapz()(从上面进行计算):

import numpy as np
from scipy.integrate import cumtrapz


def numeric_integral(x, f, c=0):
    return cumtrapz(f(x), x, initial=c) 

如下图:

import matplotlib.pyplot as plt


def func(x):
    return x ** 2


x = np.arange(-10, 10, 0.01)
y = func(x)
Y = numeric_integral(x, func)

plt.plot(x, y, label='f(x) = x²')
plt.plot(x, Y, label='F(x) = x³/3 + c')
plt.plot(x, x ** 3 / 3, label='F(x) = x³/3')
plt.legend()

它为您提供所需的结果,但您应该自己指定任意常量。

对于良好的测量,虽然在这种情况下不相关,但请注意 np.arange() 如果与小数步一起使用,则不会提供稳定的结果。通常,人们会改用np.linspace()

【讨论】:

  • 为什么使用 c=0?我可以不说吗?
  • @mathflower 当然你可以忽略它。这只是为了说明正确的数学运算,即积分函数被定义为任意常数,c 在代码中。
  • 我之所以问这个只是因为您手动插入了 x**2 的反导数而不是计算它。所以你已经计算了反导数,但是 c 是常数。我只想得到没有 c 的反导数。
  • @mathflower 我不确定你指的是什么。 numeric_integral() 中的代码一开始就对 f() 一无所知。我只是为了完整性绘制分析结果,仅此而已。 c 的存在具有数学意义。它在这种数值方法中的值仅取决于您使用的实际x。没有任何数学上有效的方法来获得 c 的值,该值始终与您在 分析反导数 中所拥有的值相匹配,该值使用您在微积分中研究的任意函数和 x 范围中的规则找到考虑。
  • @mathflower 换句话说,你想要c = 0,但你得到的是c = F(x0)x0 的下限x。由于您(必须假装)不知道F(x),因此您无法明确计算c
【解决方案2】:

来自 scipy 的 cumtrapz 函数将使用梯形积分提供反导数:

from scipy.integrate import cumtrapz
yy = cumtrapz(y, x, initial=0)

# make yy==0 around x==0 (optional)
i_x0 = np.where(x >= 0)[0][0]
yy -= yy[i_x0]

【讨论】:

    【解决方案3】:

    梯形积分

    import numpy as np
    import matplotlib.pyplot as plt
    
    x = np.arange(-10, 10, 0.1)
    f = x**2
    
    F = [-333.35]
    for i in range(1, len(x) - 1):
        F.append((f[i] + f[i - 1])*(x[i] - x[i - 1])/2 + F[i - 1])
    F = np.array(F)
    
    fig, ax = plt.subplots()
    
    ax.plot(x, f)
    ax.plot(x[1:], F)
    
    plt.show()
    

    这里我应用了理论公式(f[i] + f[i - 1])*(x[i] - x[i - 1])/2 + F[i - 1],而积分是在块中完成的:

    F = [-333.35]
    for i in range(1, len(x) - 1):
        F.append((f[i] + f[i - 1])*(x[i] - x[i - 1])/2 + F[i - 1])
    F = np.array(F)
    

    注意,为了绘制xF,它们必须具有相同数量的元素;所以我忽略了x 的第一个元素,所以它们都有199 元素。这是梯形方法d 的结果:如果你整合一个数组fn 元素,你得到一个数组Fn-1 元素。此外,我将F的初始值设置为-333.35x = -10,这是积分过程中的任意常数,我决定该值是为了在原点附近传递函数。


    分析整合

    import sympy as sy
    import numpy as np
    import matplotlib.pyplot as plt
    
    x = sy.symbols('x')
    f = x**2
    F = sy.integrate(f, x)
    
    xv = np.arange(-10, 10, 0.1)
    fv = sy.lambdify(x, f)(xv)
    Fv = sy.lambdify(x, F)(xv)
    
    fig, ax = plt.subplots()
    
    ax.plot(xv, fv)
    ax.plot(xv, Fv)
    
    plt.show()
    

    这里我通过sympy 模块使用符号数学。集成在块中完成:

    F = sy.integrate(f, x)
    

    注意,在这种情况下,Fx 已经具有相同数量的元素。而且代码更简单。

    【讨论】:

    • 请注意,您对梯形积分的使用不正确。您可以通过绘制x ** 3 / 3 来查看它。
    • 你能解释一下那行吗,拜托:(f[i] + f[i - 1])*(x[i] - x[i - 1])/2 + F[i - 1]
    • @mathflower 这就是x轴和函数之间的梯形面积加上所有previos梯形的总和的公式。
    • 正如@norok2 已经说过的那样,在该公式中,您将空中飞人的主要和次要基础(f[i]f[i - 1])相加,乘以空中飞人的高度(x[i] - x[i - 1])并除以 2,这是秋千面积的公式。最后,将此结果与上一个结果相加F[i - 1]
    • 你是怎么想到 F[-333.35] 的?我可以取另一个值吗?
    猜你喜欢
    • 2021-04-16
    • 1970-01-01
    • 2015-12-08
    • 2020-04-21
    • 2012-07-01
    • 2020-04-26
    • 1970-01-01
    • 2012-08-23
    • 1970-01-01
    相关资源
    最近更新 更多