在sympy中进行多元泰勒逼近

Question

我打算使用 sympy 写一个多元泰勒级数近似，同时：

• 尽可能使用内置代码
• 计算给定两个变量函数的截断泰勒级数近似
• 返回结果不带大O余项，例如 sin(x)=x - x**3/6 + O(x**4)。

这是我迄今为止尝试过的方法:

方法1

朴素地，我们可以对每个变量分别使用series命令进行两次组合，但这不起作用，如下面这个例子所示，用于函数 sin(x*cos(y))：

sp.sin(x*sp.cos(y)).series(x,x0=0,n=3).series(y,x0=0,n=3)
>>> NotImplementedError: not sure of order of O(y**3) + O(x**3)

方法二

根据此帖子，我首先编写了一个一维泰勒近似：

def taylor_approximation(expr, x, max_order):
    taylor_series = expr.series(x=x, n=None)
    return sum([next(taylor_series) for i in range(max_order)])

用1D例子进行检查效果很好

mport sympy as sp
x=sp.Symbol('x')
y=sp.Symbol('y')
taylor_approximation(sp.sin(x*sp.cos(y)),x,3)

>>> x**5*cos(y)**5/120 - x**3*cos(y)**3/6 + x*cos(y)

然而，如果我现在尝试在x和y上执行连锁调用以进行两种扩展，sympy会挂起。

# this does not work
taylor_approximation(taylor_approximation(sp.sin(x*sp.cos(y)),x,3),y,3)

有人知道如何修复这个问题或以替代方式实现它吗？

- flonk

Answer 1

您可以使用expr.removeO()来移除表达式中的大O。

---

一行代码：expr.series(x, 0, 3).removeO().series(y, 0, 3).removeO()

- asmeurer

---

谢谢，这个提示可以让你用一行代码解决问题：expr.series(x,x0=0,n=3).removeO().series(y,x0=0,n=3).removeO()。请将其添加到您的答案中，我会接受它 :) - flonk

我喜欢这个解决方案的简单性。然而，我意识到它是不正确的，因为它包括高阶项。假设你有 f = sqrt(x**2 + y**2 +3*x*y+ 1)，你想在 (x,y) = (0,0) 处展开它。结果的总和包含一个 x^2 y^2 的项，显然不是正确的阶数。对于非对称函数，情况甚至可能更严重。 - kluonk

Answer 2

这里是用于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。请注意，保留了HTML标记。

# 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 3

多元泰勒展开

def mtaylor(funexpr,x,mu,order=1):

    nvars = len(x)
    hlist = ['__h' + str(i+1) for i in range(nvars)]
    command=''
    command="symbols('"+'  '.join(hlist) +"')"
    hvar = eval(command)
    #mtaylor is utaylor for specificly defined function
    t = symbols('t')
    #substitution
    loc_funexpr = funexpr
    for i in range(nvars):
        locvar = x[i]
        locsubs = mu[i]+t*hvar[i]
        loc_funexpr = loc_funexpr.subs(locvar,locsubs)
    #calculate taylorseries
    g = 0
    for i in range(order+1):
        g+=loc_funexpr.diff(t,i).subs(t,0)*t**i/math.factorial(i)

    #resubstitute
    for i in range(nvars):
        g = g.subs(hlist[i],x[i]-mu[i])

    g = g.subs(t,1)    
    return g

测试某些功能

x1,x2,x3,x4,x5 = symbols('x1 x2 x3 x4 x5')
funexpr=1+x1+x2+x1*x2+x1**3
funexpr=cos(funexpr)
x=[x1,x2,x3,x4,x5]
mu=[1,1,1,1,1]
mygee = mtaylor(funexpr,x,mu,order=4)
print(mygee)

- user10586551

---

我认为这个公式在高于1阶的情况下是错误的。例如，对于order=2，Hessian矩阵的非对角线元素没有被考虑在内。请参见：https://math.libretexts.org/Bookshelves/Calculus/Supplemental_Modules_(Calculus)/Multivariable_Calculus/3%3A_Topics_in_Partial_Derivatives/Taylor__Polynomials_of_Functions_of_Two_Variables - Bernardo Costa