2

我正在尝试绘制函数 sin(x)/x 和它的泰勒近似值。我使用 python 3 和 pyzo - 第一个情节有效，但我在将来自 sympy 模块的系列转换为可以工作的 numpy 表达式时遇到问题。

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from sympy.abc import x

x = np.linspace(-10, 10, 100)
y = np.sin(x)/x #first function
plt.plot(x, y, 'k') #this is working fine

### this is a code that removes the "0(x**something)" part of 
the series at the end

我在这里找到了http://pastebin.com/ZNQakWP7

def series(expr, x, x0, n, removeO=False):
"""
sympy bugs avoided
"""
# expr_series = expr.series(x, x0, n)
# return expr_series.removeO() if removeO else expr_series
expansion = list()
for t in expr.lseries(x, x0):
    p = t.as_coeff_exponent(x)[1]
    if p < n:
        expansion.append(t)
    else:
        break
if not removeO:
    expansion.append(sp.O(x**n))
return sp.Add(*expansion)
### my code continued ####

y_t=series(sp.sin(x)/x,x,0,6,removeO=True)

如果我现在看 y_t，我会得到这个近似值

out: x**4/120 - x**2/6 + 1

现在我尝试将其转换为 numpy 以便像使用第一个函数一样绘制它。

f_t = lambdify(x, y_t,modules=['numpy'])
x = np.linspace(-10, 10, 100) #i do this because x has
#been a symbolic variable before

plt.plot(x, y_t, 'b') #this is where the problem occurs

我得到第一个情节还有第二个错误消息：

 File "<console>", line 1, in <module>
  File "F:\pyzo2013_0_2_2\lib\site-packages\matplotlib\pyplot.py", line 2832, in plot
    ret = ax.plot(*args, **kwargs)
  File "F:\pyzo2013_0_2_2\lib\site-packages\matplotlib\axes.py", line 3998, in plot
    for line in self._get_lines(*args, **kwargs):

我怎样才能实现我的想法来绘制来自 sympy 的东西？我的另一个想法是将 sympy 从系列转换为字符串，然后以某种方式将其解析为 numpy 表达式。我会感谢这里的任何帮助！

4

1 回答 1

4

我认为你的问题是：

plt.plot(x, y_t, 'b') #this is where the problem occurs

应该是这样的：

plt.plot(x, f_t(x), 'b')

f_t是lambdified系列，所以它是一个评估其参数的可调用函数。

我lambdify在以下示例中使用了它，它对我有用：

from sympy.abc import x
from sympy import sin, series
from sympy.utilities.lambdify import lambdify

import numpy as np
import matplotlib.pyplot as plt


func = sin(x)/x
taylor = series(func, n=6).removeO()

evalfunc = lambdify(x, func, modules=['numpy'])
evaltaylor = lambdify(x, taylor, modules=['numpy'])

t = np.linspace(-5.5, 5.5, 100)
plt.plot(t, evalfunc(t), 'b', label='sin(x)/x')
plt.plot(t, evaltaylor(t), 'r', label='Taylor')
plt.legend(loc='best')
plt.show()

这是脚本生成的情节。

在此处输入图像描述

于 2013-11-08T17:24:31.027 回答