使用Sympy计算多元函数的泰勒级数

Question

我正在尝试使用SymPy计算一个依赖于三角函数sinc（在这里）的函数的泰勒级数。为了简化问题，我们可以假设我需要泰勒级数的函数是：

f(x1, x2) = sinc(x1) * sinc(x2)

我的问题是，在导入sympy.mpmath时，经常出现错误：

无法从...创建mpf

我尝试使用泰勒级数逼近或另一种解决方案（1号）,但它们似乎都失败了。例如，对于后一种替代方法，该行：

(sinc(x)*sinc(y)).series(x,0,3).removeO().series(y,0,3).removeO()

返回：

无法从x创建mpf

我也尝试将函数定义为表达式和lambda函数。但似乎没有什么作用。

任何帮助都将不胜感激。

- lasofivec

Answer 1

这里是一个用于使用Sympy的多元泰勒级数展开式:

def Taylor_polynomial_sympy(function_expression, variable_list, evaluation_point, degree):
    """
    Mathematical formulation reference:
    https://math.libretexts.org/Bookshelves/Calculus/Supplemental_Modules_(Calculus)/Multivariable_Calculus/3%3A_Topics_in_Partial_Derivatives/Taylor__Polynomials_of_Functions_of_Two_Variables
    :param function_expression: Sympy expression of the function
    :param variable_list: list. All variables to be approximated (to be "Taylorized")
    :param evaluation_point: list. Coordinates, where the function will be expressed
    :param degree: int. Total degree of the Taylor polynomial
    :return: Returns a Sympy expression of the Taylor series up to a given degree, of a given multivariate expression, approximated as a multivariate polynomial evaluated at the evaluation_point
    """
    from sympy import factorial, Matrix, prod
    import itertools

    n_var = len(variable_list)
    point_coordinates = [(i, j) for i, j in (zip(variable_list, evaluation_point))]  # list of tuples with variables and their evaluation_point coordinates, to later perform substitution

    deriv_orders = list(itertools.product(range(degree + 1), repeat=n_var))  # list with exponentials of the partial derivatives
    deriv_orders = [deriv_orders[i] for i in range(len(deriv_orders)) if sum(deriv_orders[i]) <= degree]  # Discarding some higher-order terms
    n_terms = len(deriv_orders)
    deriv_orders_as_input = [list(sum(list(zip(variable_list, deriv_orders[i])), ())) for i in range(n_terms)]  # Individual degree of each partial derivative, of each term

    polynomial = 0
    for i in range(n_terms):
        partial_derivatives_at_point = function_expression.diff(*deriv_orders_as_input[i]).subs(point_coordinates)  # e.g. df/(dx*dy**2)
        denominator = prod([factorial(j) for j in deriv_orders[i]])  # e.g. (1! * 2!)
        distances_powered = prod([(Matrix(variable_list) - Matrix(evaluation_point))[j] ** deriv_orders[i][j] for j in range(n_var)])  # e.g. (x-x0)*(y-y0)**2
        polynomial += partial_derivatives_at_point / denominator * distances_powered
    return polynomial

以下是一个双变量问题的验证，遵循练习和答案中的步骤：https://math.libretexts.org/Bookshelves/Calculus/Supplemental_Modules_(Calculus)/Multivariable_Calculus/3%3A_Topics_in_Partial_Derivatives/Taylor__Polynomials_of_Functions_of_Two_Variables

# Solving the exercises in section 13.7 of https://math.libretexts.org/Bookshelves/Calculus/Supplemental_Modules_(Calculus)/Multivariable_Calculus/3%3A_Topics_in_Partial_Derivatives/Taylor__Polynomials_of_Functions_of_Two_Variables
from sympy import symbols, sqrt, atan, ln

# Exercise 1
x = symbols('x')
y = symbols('y')
function_expression = x*sqrt(y)
variable_list = [x,y]
evaluation_point = [1,4]
degree=1
print(Taylor_polynomial_sympy(function_expression, variable_list, evaluation_point, degree))
degree=2
print(Taylor_polynomial_sympy(function_expression, variable_list, evaluation_point, degree))

# Exercise 3
x = symbols('x')
y = symbols('y')
function_expression = atan(x+2*y)
variable_list = [x,y]
evaluation_point = [1,0]
degree=1
print(Taylor_polynomial_sympy(function_expression, variable_list, evaluation_point, degree))
degree=2
print(Taylor_polynomial_sympy(function_expression, variable_list, evaluation_point, degree))

# Exercise 5
x = symbols('x')
y = symbols('y')
function_expression = x**2*y + y**2
variable_list = [x,y]
evaluation_point = [1,3]
degree=1
print(Taylor_polynomial_sympy(function_expression, variable_list, evaluation_point, degree))
degree=2
print(Taylor_polynomial_sympy(function_expression, variable_list, evaluation_point, degree))

# Exercise 7
x = symbols('x')
y = symbols('y')
function_expression = ln(x**2+y**2+1)
variable_list = [x,y]
evaluation_point = [0,0]
degree=1
print(Taylor_polynomial_sympy(function_expression, variable_list, evaluation_point, degree))
degree=2
print(Taylor_polynomial_sympy(function_expression, variable_list, evaluation_point, degree))

在结果上执行simplify()可能会很有用。

- Bernardo Costa

Answer 2

你不能将mpmath函数与符号化的SymPy对象一起使用。你需要以符号方式定义sinc。

一种方法是将其定义为Python函数，该函数返回等效值（这里sin为sympy.sin）:

def sinc(x):
    return sin(x)/x

或者，您可以为此编写自己的SymPy类。

class sinc(Function):
    @classmethod
    def eval(cls, x):
        if x == 0:
            return 1
        # Should also include oo and -oo here
        sx = sin(x)
        # Return only if sin simplifies
        if not isinstance(sx, sin):
            return sx/x

    # You'd also need to define fdiff or _eval_derivative here so that it knows what the derivative is

后者允许你写sinc(x)并以那样的方式打印，你也可以定义它的自定义行为，比如导数的外观、级数的工作方式等等。这也能让你在0、oo等点上定义正确的值。

但是，为了进行级数计算，仅使用sin(x)/x（即第一种解决方案）就足够了，因为你最终的答案不会有sinc，而且在评估0等点时，级数方法会正确地取极限。

- asmeurer

---

还有一个开放的SymPy问题需要实现sinc https://github.com/sympy/sympy/issues/5051 - asmeurer

谢谢。我担心第一个解决方案不起作用，因为限制为0。但它确实有效。谢谢。 - lasofivec