使用积分时计算错误
Wrong calculation when using integrals
对于上下文：我有一些用 Mathcad 做的旧方程，我想用 sympy 重写，所以我确实得到了我应该得到的结果。
但是，考虑到公式：
uu, vv, rr = sp.symbols("u v r")
h2 = sp.Integral(2*sp.exp( (sp.I*uu*rr**2) /2 )*sp.besselj(0,vv*rr)*rr,(rr,0,1))
如果我将其中一个变量（uu 或 vv）代入 0，则计算正确。
f=sp.lambdify((uu),(sp.Abs(h2.subs(vv,0).doit()))**2)
g=sp.lambdify((vv),(sp.Abs(h2.subs(uu,0).doit()))**2)
但是，如果我想使用这两个变量制作 3D 图形，结果是不正确的。
f2 = sp.lambdify((uu,vv),(sp.Abs(h2.doit()))**2)
f2(5,0) 应该得到 0.57，对于 f(5)，但它没有
所以我用一些循环解决了这个问题，而不是同时对两个变量进行 lambdifying :
z_list = []
for j in yy:
f3=sp.lambdify([uu],sp.re((sp.Abs(h2.subs(vv,j).doit())**2)))
for i in xx:
z_list.append(f3(i))
这需要更长的时间来计算，但结果很好。
关于积分或 lambdify 或两者的组合，我是否遗漏了什么？
此外，我需要计算的下一个公式是前面 h2 的积分：
I6 = sp.Integral(vv*sp.Abs(h2)**4,(vv,0,10))
从 11-01 编辑
感谢 Oscar Benjamin 的回答，计算中的问题来自 lambdify 的返回值，returns 只有实部。
我测试了三种情况：
f3=sp.lambdify((uu),h2.subs(vv,0).doit())
f4=sp.lambdify((uu),h2.subs(vv,0).evalf())
f5=sp.lambdify((uu),h2.subs(vv,0))
只有 f3 returns 具有实部和虚部的值（这解释了我之前遇到的其他计算问题）。
(0.23938885764158258+0.7204574462187735j)
0.23938885764158258
0.23938885764158258
可以用其他模块调用 lambdify 来防止这种情况发生吗？
doit 积分函数预期用于以符号方式计算积分：
In [47]: f = Integral(exp(-x**2), (x, 0, oo))
In [48]: f
Out[48]:
∞
⌠
⎮ 2
⎮ -x
⎮ ℯ dx
⌡
0
In [49]: f.doit()
Out[49]:
√π
──
2
In [50]: _.evalf()
Out[50]: 0.886226925452758
当它起作用时它很好，因为一个不涉及整数符号的显式表达式可以更有效地进行数值计算。您可以使用 evalf 进行多精度计算的评估，或者您可以使用 lambdify 并让 numpy/scipy 等在机器精度浮点计算中更快地进行计算。
然而，并不是所有的积分都可以用符号计算，所以 doit 有时会 return 相同的积分：
In [53]: h2 = Integral(2*r*exp(I*r**2*u/2)*besselj(0, r*v), (r, 0, 1))
In [54]: h2
Out[54]:
1
⌠
⎮ 2
⎮ ⅈ⋅r ⋅u
⎮ ──────
⎮ 2
⎮ 2⋅r⋅ℯ ⋅besselj(0, r⋅v) dr
⌡
0
In [55]: h2.doit() # slow
Out[55]:
1
⌠
⎮ 2
⎮ ⅈ⋅r ⋅u
⎮ ──────
⎮ 2
2⋅⎮ r⋅ℯ ⋅besselj(0, r⋅v) dr
⌡
0
这里我们看到积分没有求值。有时，如果积分具有积分变量以外的符号，我们可以通过代入值来获得符号结果。如果将 u 或 v 设置为零，就会发生这种情况：
In [59]: h2.subs({u:0, v:5}).doit()
Out[59]:
2⋅besselj(1, 5)
───────────────
5
In [60]: h2.subs({u:5, v:0}).doit()
Out[60]:
5⋅ⅈ
───
2
2⋅ⅈ⋅ℯ 2⋅ⅈ
- ──────── + ───
5 5
因为被积函数变得更简单，所以可以用符号来计算积分。但是，在这种情况下，如果 u 或 v 都不为零，则无法对积分进行符号计算：
In [61]: h2.subs({uu:2,vv:3}).doit()
Out[61]:
1
⌠
⎮ 2
⎮ ⅈ⋅r
2⋅⎮ r⋅ℯ ⋅besselj(0, 3⋅r) dr
⌡
0
我希望这个特定的积分没有直接的解析公式。我检查了 Wolfram Alpha，它也没有找到符号结果（只有数字结果）：
不可能从积分中找到解析结果是符号学的一个基本事实：不存在可识别的数学函数来表达大多数可能的积分结果。即使在许多情况下，函数确实存在，但并不存在可以找到它们的算法。
所以我们接受符号集成是不可能的，这意味着调用 doit 没有意义，这是迄今为止您使用的代码中最慢的操作。 doit 的目的是获得一个符号结果，然后我们可以更有效地使用它来处理 u 和 v 的许多不同值，但是您却为每个调用 doit u 和 v 的可能值，但每次都没有 return 有用的东西。
如果无法进行符号积分，则需要进行数值积分。对于数值积分，我们首先需要替换 u 和 v 的值，然后我们可以使用数值积分算法。由于我们没有积分的符号结果，因此对于每个可能的参数值，数值积分算法需要 运行 一次。在 SymPy 中有两种直接的方法可以做到这一点：.evalf 或 lambdify（请注意 N(obj) 与 obj.evalf() 和 obj.n() 相同）。
使用 .evalf 速度较慢，因为它使用多精度算法，但这也使它更准确。 returned 结果中的每个数字都可以被合理地认为是正确的（除非有错误...），如果需要，您还可以计算更多数字：
In [64]: %time h2.subs({uu:2,vv:3}).evalf()
CPU times: user 1.35 s, sys: 0 ns, total: 1.35 s
Wall time: 1.35 s
Out[64]: 0.236305280388988 + 0.0146270775320418⋅ⅈ
In [65]: %time h2.subs({uu:2,vv:3}).evalf(50)
CPU times: user 3.64 s, sys: 8 ms, total: 3.65 s
Wall time: 3.65 s
Out[65]:
0.23630528038898808262356997039156552180149071451802 + 0.014627077532041813385103750863443465733130
09416169⋅ⅈ
这比调用 doit 快，但与顶级标准机器精度数字代码相比仍然慢。虽然它非常准确。
为了加快评估速度，您可以使用 lambdify:
In [66]: f2 = lambdify((uu, vv), h2)
In [67]: %time f2(2, 3)
/home/oscar/current/sympy/sympy.git/38venv/lib/python3.8/site-packages/scipy/integrate/quadpack.py:463: ComplexWarning: Casting complex values to real discards the imaginary part
return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 438 µs
Out[67]: 0.2363052803889881
这现在 return 非常快并且给出了准确的结果，除了它只是 return 的实部。我猜 scipy 例程只处理实数积分。我们可以分开做实部和虚部，例如虚部是：
In [91]: %time f2(2, 3)
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 387 µs
Out[91]: 0.01462707753204179
在 lambdify 表达式时使用模块“mpmath”return 完整的复数，而不仅仅是它的实部。
f3 = sp.lambdify([uu,vv],h2,modules='mpmath')
f3(5,0)
mpc(real='0.2393888576415826', imag='0.7204574462187735')
对于上下文：我有一些用 Mathcad 做的旧方程，我想用 sympy 重写，所以我确实得到了我应该得到的结果。
但是，考虑到公式：
uu, vv, rr = sp.symbols("u v r")
h2 = sp.Integral(2*sp.exp( (sp.I*uu*rr**2) /2 )*sp.besselj(0,vv*rr)*rr,(rr,0,1))
如果我将其中一个变量（uu 或 vv）代入 0，则计算正确。
f=sp.lambdify((uu),(sp.Abs(h2.subs(vv,0).doit()))**2)
g=sp.lambdify((vv),(sp.Abs(h2.subs(uu,0).doit()))**2)
但是，如果我想使用这两个变量制作 3D 图形，结果是不正确的。
f2 = sp.lambdify((uu,vv),(sp.Abs(h2.doit()))**2)
f2(5,0) 应该得到 0.57，对于 f(5)，但它没有
所以我用一些循环解决了这个问题，而不是同时对两个变量进行 lambdifying :
z_list = []
for j in yy:
f3=sp.lambdify([uu],sp.re((sp.Abs(h2.subs(vv,j).doit())**2)))
for i in xx:
z_list.append(f3(i))
这需要更长的时间来计算，但结果很好。 关于积分或 lambdify 或两者的组合，我是否遗漏了什么？
此外，我需要计算的下一个公式是前面 h2 的积分：
I6 = sp.Integral(vv*sp.Abs(h2)**4,(vv,0,10))
从 11-01 编辑
感谢 Oscar Benjamin 的回答，计算中的问题来自 lambdify 的返回值，returns 只有实部。
我测试了三种情况：
f3=sp.lambdify((uu),h2.subs(vv,0).doit())
f4=sp.lambdify((uu),h2.subs(vv,0).evalf())
f5=sp.lambdify((uu),h2.subs(vv,0))
只有 f3 returns 具有实部和虚部的值（这解释了我之前遇到的其他计算问题）。
(0.23938885764158258+0.7204574462187735j)
0.23938885764158258
0.23938885764158258
可以用其他模块调用 lambdify 来防止这种情况发生吗？
doit 积分函数预期用于以符号方式计算积分：
In [47]: f = Integral(exp(-x**2), (x, 0, oo))
In [48]: f
Out[48]:
∞
⌠
⎮ 2
⎮ -x
⎮ ℯ dx
⌡
0
In [49]: f.doit()
Out[49]:
√π
──
2
In [50]: _.evalf()
Out[50]: 0.886226925452758
当它起作用时它很好，因为一个不涉及整数符号的显式表达式可以更有效地进行数值计算。您可以使用 evalf 进行多精度计算的评估，或者您可以使用 lambdify 并让 numpy/scipy 等在机器精度浮点计算中更快地进行计算。
然而，并不是所有的积分都可以用符号计算，所以 doit 有时会 return 相同的积分：
In [53]: h2 = Integral(2*r*exp(I*r**2*u/2)*besselj(0, r*v), (r, 0, 1))
In [54]: h2
Out[54]:
1
⌠
⎮ 2
⎮ ⅈ⋅r ⋅u
⎮ ──────
⎮ 2
⎮ 2⋅r⋅ℯ ⋅besselj(0, r⋅v) dr
⌡
0
In [55]: h2.doit() # slow
Out[55]:
1
⌠
⎮ 2
⎮ ⅈ⋅r ⋅u
⎮ ──────
⎮ 2
2⋅⎮ r⋅ℯ ⋅besselj(0, r⋅v) dr
⌡
0
这里我们看到积分没有求值。有时，如果积分具有积分变量以外的符号，我们可以通过代入值来获得符号结果。如果将 u 或 v 设置为零，就会发生这种情况：
In [59]: h2.subs({u:0, v:5}).doit()
Out[59]:
2⋅besselj(1, 5)
───────────────
5
In [60]: h2.subs({u:5, v:0}).doit()
Out[60]:
5⋅ⅈ
───
2
2⋅ⅈ⋅ℯ 2⋅ⅈ
- ──────── + ───
5 5
因为被积函数变得更简单，所以可以用符号来计算积分。但是，在这种情况下，如果 u 或 v 都不为零，则无法对积分进行符号计算：
In [61]: h2.subs({uu:2,vv:3}).doit()
Out[61]:
1
⌠
⎮ 2
⎮ ⅈ⋅r
2⋅⎮ r⋅ℯ ⋅besselj(0, 3⋅r) dr
⌡
0
我希望这个特定的积分没有直接的解析公式。我检查了 Wolfram Alpha，它也没有找到符号结果（只有数字结果）：
不可能从积分中找到解析结果是符号学的一个基本事实：不存在可识别的数学函数来表达大多数可能的积分结果。即使在许多情况下，函数确实存在，但并不存在可以找到它们的算法。
所以我们接受符号集成是不可能的，这意味着调用 doit 没有意义，这是迄今为止您使用的代码中最慢的操作。 doit 的目的是获得一个符号结果，然后我们可以更有效地使用它来处理 u 和 v 的许多不同值，但是您却为每个调用 doit u 和 v 的可能值，但每次都没有 return 有用的东西。
如果无法进行符号积分，则需要进行数值积分。对于数值积分，我们首先需要替换 u 和 v 的值，然后我们可以使用数值积分算法。由于我们没有积分的符号结果，因此对于每个可能的参数值，数值积分算法需要 运行 一次。在 SymPy 中有两种直接的方法可以做到这一点：.evalf 或 lambdify（请注意 N(obj) 与 obj.evalf() 和 obj.n() 相同）。
使用 .evalf 速度较慢，因为它使用多精度算法，但这也使它更准确。 returned 结果中的每个数字都可以被合理地认为是正确的（除非有错误...），如果需要，您还可以计算更多数字：
In [64]: %time h2.subs({uu:2,vv:3}).evalf()
CPU times: user 1.35 s, sys: 0 ns, total: 1.35 s
Wall time: 1.35 s
Out[64]: 0.236305280388988 + 0.0146270775320418⋅ⅈ
In [65]: %time h2.subs({uu:2,vv:3}).evalf(50)
CPU times: user 3.64 s, sys: 8 ms, total: 3.65 s
Wall time: 3.65 s
Out[65]:
0.23630528038898808262356997039156552180149071451802 + 0.014627077532041813385103750863443465733130
09416169⋅ⅈ
这比调用 doit 快，但与顶级标准机器精度数字代码相比仍然慢。虽然它非常准确。
为了加快评估速度，您可以使用 lambdify:
In [66]: f2 = lambdify((uu, vv), h2)
In [67]: %time f2(2, 3)
/home/oscar/current/sympy/sympy.git/38venv/lib/python3.8/site-packages/scipy/integrate/quadpack.py:463: ComplexWarning: Casting complex values to real discards the imaginary part
return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 438 µs
Out[67]: 0.2363052803889881
这现在 return 非常快并且给出了准确的结果，除了它只是 return 的实部。我猜 scipy 例程只处理实数积分。我们可以分开做实部和虚部，例如虚部是：
In [91]: %time f2(2, 3)
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 387 µs
Out[91]: 0.01462707753204179
在 lambdify 表达式时使用模块“mpmath”return 完整的复数，而不仅仅是它的实部。
f3 = sp.lambdify([uu,vv],h2,modules='mpmath')
f3(5,0)
mpc(real='0.2393888576415826', imag='0.7204574462187735')