微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

Scipy 在多个列表上进行数值积分

如何解决Scipy 在多个列表上进行数值积分

我有一个带有跨度步长和相应旋转值的图表。我需要对每个步骤进行数值积分以获得斜率值。我想知道因为 scipy 集成中已经有内置函数,如梯形规则或辛普森规则。如何在没有任何附加功能的情况下在两个数组或数据列表上实现?

enter image description here

version: "3.5"
services:
  rabbitmq:
    image: rabbitmq:3-alpine
    expose:
      - 5672
      - 15672
    volumes:
      - ./rabbit/enabled_plugins:/etc/rabbitmq/enabled_plugins
    labels:
      - traefik.enable=true
      - traefik.http.routers.rabbitmq.rule=Host(`HOST.com`) && PathPrefix(`/rmq`)
      # needed,when you do not have a route "/rmq" inside your container (according to https://stackoverflow.com/questions/59054551/how-to-map-specific-port-inside-docker-container-when-using-traefik)
      - traefik.http.routers.rabbitmq.middlewares=strip-docs
      - traefik.http.middlewares.strip-docs.stripprefix.prefixes=/rmq
      - traefik.http.services.rabbitmq.loadbalancer.server.port=15672
    networks:
      - proxynet

  traefik:
    image: traefik:2.1
    command: --api=true     # Enables the web UI
    volumes:
      - /var/run/docker.sock:/var/run/docker.sock:ro
      - ./traefik/traefik.toml:/etc/traefik/traefik.toml:ro
    ports:
      - 80:80
      - 443:443
    labels:
      traefik.enable: true
      traefik.http.routers.traefik.rule: "Host(`HOST.com`)"
      traefik.http.routers.traefik.service: "api@internal"
    networks:
      - proxynet

预期结果:

import scipy
fraction_of_span = [0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]
rotation = [0.33,1.34,2.62,3.41,3.87,4.02,0] 
result = scipy.trapz(fraction_of_span,rotation,10)

解决方法

如上提议

import scipy
fraction_of_span = [0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]
rotation = [0.33,1.34,2.62,3.41,3.87,4.02,0]
result = scipy.trapz(fraction_of_span,rotation,10)
print(result)
-2.6665

以上已构建解决方案的示例,使用 simp

from scipy.integrate import simps
y = rotation
x = fraction_of_span


result_simps = simps(y,x)
print(result_simps)
2.6790000000000003

请注意,结果非常相似,因为方法不同而略有不同。注意符号应该是正的,因为积分只在正值之间(旋转元素都是正的)

有很好的材料符合你的要求,也许你想看看这里:docs.scipy.org/doc/scipy/reference/tutorial/integrate.html

让我们尝试获得类似的矢量结果。

为此,您可以转到上述函数并修改它们以获得结果。所以我去 https://github.com/numpy/numpy/blob/master/numpy/lib/function_base.py#L4081-L4169 并修改/创建一个新的函数,如下所示:

def trapz_modified(y,x=None,dx=1.0,axis=-1):
    """
    Integrate along the given axis using the composite trapezoidal rule.
    Integrate `y` (`x`) along given axis.
    Parameters
    ----------
    y : array_like
        Input array to integrate.
    x : array_like,optional
        The sample points corresponding to the `y` values. If `x` is None,the sample points are assumed to be evenly spaced `dx` apart. The
        default is None.
    dx : scalar,optional
        The spacing between sample points when `x` is None. The default is 1.
    axis : int,optional
        The axis along which to integrate.
    Returns
    -------
    trapz : float
        Definite integral as approximated by trapezoidal rule.
    See Also
    --------
    sum,cumsum
    Notes
    -----
    Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
    will be taken from `y` array,by default x-axis distances between
    points will be 1.0,alternatively they can be provided with `x` array
    or with `dx` scalar.  Return value will be equal to combined area under
    the red lines.
    References
    ----------
    .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
    .. [2] Illustration image:
           https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
    Examples
    --------
    >>> np.trapz([1,2,3])
    4.0
    >>> np.trapz([1,3],x=[4,6,8])
    8.0
    >>> np.trapz([1,dx=2)
    8.0
    >>> a = np.arange(6).reshape(2,3)
    >>> a
    array([[0,1,2],[3,4,5]])
    >>> np.trapz(a,axis=0)
    array([1.5,2.5,3.5])
    >>> np.trapz(a,axis=1)
    array([2.,8.])
    """
    y = asanyarray(y)
    if x is None:
        d = dx
    else:
        x = asanyarray(x)
        if x.ndim == 1:
            d = diff(x)
            # reshape to correct shape
            shape = [1]*y.ndim
            shape[axis] = d.shape[0]
            d = d.reshape(shape)
        else:
            d = diff(x,axis=axis)
    nd = y.ndim
    slice1 = [slice(None)]*nd
    slice2 = [slice(None)]*nd
    slice1[axis] = slice(1,None)
    slice2[axis] = slice(None,-1)
    try:
        # MODIFIED HERE
        #ret = (d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0).sum(axis)
        ret = d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0
    except ValueError:
        # Operations didn't work,cast to ndarray
        d = np.asarray(d)
        y = np.asarray(y)
        # MODIFIED HERE
        #ret = add.reduce(d * (y[tuple(slice1)]+y[tuple(slice2)])/2.0,axis)
        ret = d * (y[tuple(slice1)]+y[tuple(slice2)])/2.0

    return ret

我们还需要以下库,位于文件/脚本的顶部:

from numpy import diff
from numpy import asanyarray

让我们看看输出:

>>>trapz_modified(y,x=x)
array([0.0835,0.198,0.3015,0.364,0.3945,0.067 ])

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。