Sympy - 绘制矩阵计算的结果

Question

在我的程序中，我使用 sympy 来获取一组常量（r、a、b 和 tau）的解析表达式，然后将 r 的图形绘制为 l 的函数。首先，我这样定义我的常量：

import sympy
import numpy
import matplotlib.pyplot as plt

[z1, z2, z3, k2] = [1, 2, 4, 2 * numpy.pi]
l = sympy.Symbol('l')
r = sympy.Symbol('r')
tau = sympy.Symbol('tau')
base = numpy.linspace(0, 10, 10000)

然后我定义并计算我的矩阵（其中 M2 是 [r, a, b, tau] 的转置）：

M1 = sympy.Matrix([[-1, 1, 1, 0],
               [1, z1/z2, - z1/z2, 0],
               [0, sympy.exp(- 1j * k2 * l), sympy.exp(1j * k2 * l), -1],
               [0, - sympy.exp(- 1j * k2 * l), sympy.exp(1j * k2 * l), z2/z3],
               ])


M3 = sympy.Matrix([[1], [1], [0], [0]])

M2 = (M1 ** -1) * M3

然后我尝试使用：

sympy.plot(M2[0].subs(l, base), base)

这会导致以下错误：

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "C:\Program Files\JetBrains\PyCharm 2020.2.3\plugins\python\helpers\pydev    \_pydev_bundle\pydev_umd.py", line 197, in runfile
pydev_imports.execfile(filename, global_vars, local_vars)  # execute the script
  File "C:\Program Files\JetBrains\PyCharm 2020.2.3\plugins\python\helpers\pydev\_pydev_imps\_pydev_execfile.py", line 18, in execfile
exec(compile(contents+"\n", file, 'exec'), glob, loc)
  File "matrix solver.py", line 47, in <module>
sympy.plot(M2[0].subs(l, base), base)
  File "AppData\Local\Programs\Python\Python39\lib\site-packages\sympy\plotting\plot.py", line 1641, in plot
plot_expr = check_arguments(args, 1, 1)
  File "AppData\Local\Programs\Python\Python39\lib\site-packages\sympy\plotting\plot.py", line 2302, in check_arguments
assert all(isinstance(e, Expr) for expr in exprs for e in expr)
AssertionError

Answer 1

（该帖子已更新，以包含 @OscarBenjamin 富有洞察力的评论。）

.subs(l, base) 不起作用。一般来说，sympy 和 numpy 很难混合。 Sympy 效果最好，所有表达式都包含 sympy 符号、整数和分数。浮点数，由于其不精确的性质，常常使 sympy 很难获得精确的符号解决方案。此外，由于 sympy 是如何集成到 Python 中的，一个烦人的 gotcha 是 Python 将诸如 1/2 之类的表达式转换为 0.5 而不是 sympy 分数。

因此，最好不要在与 sympy 一起使用的代码部分中使用import numpy。为了避免分数问题，教程建议写z1 = sympy.S(1)。 np.pi 是一个浮点数（大约 16 位），而 sympy.pi 是符号常量 pi。将 1j 替换为 sympy.I 让 sympy 象征性地进行计算。

import sympy

[z1, z2, z3, k2] = [sympy.S(1), sympy.S(2), sympy.S(4), 2 * sympy.pi]
l = sympy.Symbol('l')
r = sympy.Symbol('r')
tau = sympy.Symbol('tau')

M1 = sympy.Matrix([[-1, 1, 1, 0],
                   [1, z1 / z2, - z1 / z2, 0],
                   [0, sympy.exp(- sympy.I * k2 * l), sympy.exp(sympy.I * k2 * l), -1],
                   [0, - sympy.exp(- sympy.I * k2 * l), sympy.exp(sympy.I * k2 * l), z2 / z3],
                   ])

M3 = sympy.Matrix([[1], [1], [0], [0]])
M2 = (M1 ** -1) * M3

由于结果包含复数，因此绘图的同情方式将使用sympy.re(M2[0])。 Sympy 的绘图默认使用自适应策略来决定细节级别，这对于具有大量精细细节的绘图可能会很慢。

sympy.plot(sympy.re(M2[0]), (l, 0, 10), adaptive=False)

同样，虚部可以使用sympy.im(M2[0]) 绘制。

要使用 numpy 和 matplotlib，lambdify 可以将 sympy 表达式转换为 numpy 函数：

M2_0_np = sympy.lambdify(l, M2[0])

从此M2_0_np()可以作为一个numpy函数使用。

import numpy
import matplotlib.pyplot as plt

base = numpy.linspace(0, 10, 10000)
plt.plot(M2_0_np(base), base)

绘制函数会发出警告“ComplexWarning：将复数值转换为实数会丢弃虚部”。因此，只绘制了函数的实部：

可以结合实部和虚部创建更多奇特的情节，例如通过颜色。