使用 numpy [duplicate] 应用特殊 sympy 函数时出错-解网

问：

2个月前关闭。

我正在尝试使用磁场的数值解创建一个图形。对象是螺线管。

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

python numpy sympy 符号数学

这个函数最大的问题是它的值是无限的，当 DiracDelta（0）返回 'DiracDelta（0）' 时，它的值是无限的。考虑到高斯近似的存在，我通过观察“被积”在数学上的行为方式来解决这个问题。每当“integrad”中的“t”为 0 时，此函数就会被覆盖，所以我手动删除了它。

1赞 Davide_sd 9/23/2023 #2

我将使用 SymPy 绘图后端。

我们能做的第一件事就是用 .这将在视觉上创建更长的表达式，但使用时计算速度要快得多。HeavisidePiecewisempmath

我将创建一个 3D 流线图。使用 NumPy 进行评估将失败，因为它不知道是什么。但是使用 mpmath 进行评估会成功。但是，mpmath 不支持矢量化操作，因此计算速度会很慢。DiracDelta

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, 
)

它看起来很混乱，但使用 K3D，您可以快速添加剪裁平面（单击 K3D-面板 -> 控件 -> 剪裁平面 -> 添加新的）。

如果你仔细观察，所有的流线似乎都位于某个子午线平面上。因此，我们可以选择一个 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)$"
)

上一个：如何使用带有 numpy 或 scipy 的正态分布，除了平均值和标准差外，还使用最小值和最大值？

下一个：如何演示给定二维函数的 numpy 梯度使用情况？

使用 numpy [duplicate] 应用特殊 sympy 函数时出错

Error applying special sympy functions with numpy [duplicate]

评论