使用numpy [重复]应用特殊渐近函数时出错

bqujaahr 于 2023-10-19 发布在其他

关注(0)|答案(2)|浏览(89)

这个问题已经有答案了：

Solving an integral with dirac delta using sympy（1个答案）
24天前关闭
我正试图用磁场的数值解来创建一个图形。对象是螺线管。

import numpy as np
import sympy as smp
from scipy.integrate import quad_vec

t, x, y, z = smp.symbols("t x y z")

l = smp.Matrix([0.5 * smp.cos(t), 
                0.5 * smp.sin(t), 
                (t / (600 * smp.pi)) * smp.Heaviside(t - 1) * smp.Heaviside(t)])

r = smp.Matrix([x, y, z])
sep = r-l

integrand = smp.diff(l, t).cross(sep) / sep.norm()**3

dBxdt = smp.lambdify([t, x, y, z], integrand[0])
dBydt = smp.lambdify([t, x, y, z], integrand[1])
dBzdt = smp.lambdify([t, x, y, z], integrand[2])

def B(x, y, z):
    return np.array([quad_vec(dBxdt, 0, 2*np.pi, args=(x, y, z))[0],
                     quad_vec(dBydt, 0, 2*np.pi, args=(x, y, z))[0],
                     quad_vec(dBzdt, 0, 2*np.pi, args=(x, y, z))[0]])

x = np.linspace(-2, 2, 20)
xv, yv, zv = np.meshgrid(x, x, x)

B_field = B(xv, yv, zv)
Bx, By, Bz = B_field

File <lambdifygenerated-8>:2, in _lambdifygenerated(t, x, y, z)
      1 def _lambdifygenerated(t, x, y, z):
----> 2     return (-(y - 0.5*sin(t))*((1/600)*t*DiracDelta(t)*select([less(t, 1),equal(t, 1),True], [0,1/2,1], default=nan)/pi + (1/600)*t*DiracDelta(t - 1)*select([less(t, 0),equal(t, 0),True], [0,1/2,1], default=nan)/pi + (1/600)*select([less(t, 0),equal(t, 0),True], [0,1/2,1], default=nan)*select([less(t, 1),equal(t, 1),True], [0,1/2,1], default=nan)/pi) + 0.5*(-1/600*t*select([less(t, 0),equal(t, 0),True], [0,1/2,1], default=nan)*select([less(t, 1),equal(t, 1),True], [0,1/2,1], default=nan)/pi + z)*cos(t))/(abs(x - 0.5*cos(t))**2 + abs(y - 0.5*sin(t))**2 + abs((1/600)*t*select([less(t, 0),equal(t, 0),True], [0,1/2,1], default=nan)*select([less(t, 1),equal(t, 1),True], [0,1/2,1], default=nan)/pi - z)**2)**(3/2)

NameError: name 'DiracDelta' is not defined

numpy

来源：https://stackoverflow.com/questions/77163188/error-applying-special-sympy-functions-with-numpy

2条答案

按热度按时间

cbjzeqam1#

我将使用SymPy Plotting Backends。
我们可以做的第一件事是用Piecewise重写Heaviside。这将在视觉上创建更长的表达式，但使用mpmath时，它们的计算速度会更快。
我将创建一个3D流线图。NumPy的评估将失败，因为它不知道DiracDelta是什么。但是使用mmpmath的评估会成功。但是，mpmath不支持矢量化操作，因此计算速度会很慢。

from spb import *

# rewrite integrand's Heaviside with Piecewise
a = smp.Integral(integrand[0].rewrite(smp.Piecewise), (t, 0, 2*smp.pi))
b = smp.Integral(integrand[1].rewrite(smp.Piecewise), (t, 0, 2*smp.pi))
c = smp.Integral(integrand[2].rewrite(smp.Piecewise), (t, 0, 2*smp.pi))

graphics(
    vector_field_3d(
        a, b, c, (x, -2, 2), (y, -2, 2), (z, -2, 2),

        # `a, b, c` are complicated expressions, whose
        # latex representation would make the backend fail.
        # let's set an easier label
        label="a", 

        # number of discretization points on each direction.
        # Remember you are discretizing a volume, so you will have
        # n^3 evaluation points. The higher this number the slower
        # the evaluation.
        n=10,

        # evaluates the expressions with mpmath.
        # mpmath appears to be able to deal with `DiracDelta`, and
        # it is usually faster than SymPy at evaluating an expression.
        # however, it is much much slower than NumPy, hence n must be
        # kept low
        modules="mpmath",

        # visualize streamlines (True) or quivers (False)
        streamlines=True,
        stream_kw={
            # request random starting locations for the streamlines
            "starts": True
        }
    ),

    # use the plotting library K3D
    backend=KB, 
)

它看起来很乱，但使用K3 D可以快速添加剪切平面（单击K3 D-面板->控件->剪切平面->添加新平面）。

如果你仔细观察，所有的流线似乎都位于某个子午面上。因此，我们可以选择一个2D平面，例如平面XZ（y=0）：

a = smp.Integral(integrand[0].rewrite(smp.Piecewise).subs(y, 0), (t, 0, 2*smp.pi))
b = smp.Integral(integrand[1].rewrite(smp.Piecewise).subs(y, 0), (t, 0, 2*smp.pi))
c = smp.Integral(integrand[2].rewrite(smp.Piecewise).subs(y, 0), (t, 0, 2*smp.pi))

from spb import *
graphics(
    vector_field_2d(
        a, c, (x, -2, 2), (z, -2, 2),
        label="Magnitude",
        modules="mpmath",
        
        # do not visualize a contour plot with magnitude of the
        # vector field, as it evaluate over 100^2 points, which
        # is too slow for this complicated expression
        scalar=False,
        
        # visualize streamlines (True) or quivers (False)
        streamlines=True,
        
        # streamlines benefit from higher number of discretization points
        n=30
    ),
    backend=MB, aspect="equal", grid=False,
    xlabel="x", ylabel="z", title=r"$\vec{B}(x, 0, z)$"
)