Multivariate Taylor approximation in sympy

Question

I aim to write a multidimensional Taylor approximation using sympy, which

• uses as many builtin code as possible,
• computes the truncated Taylor approximation of a given function of two variables
• returns the result without the Big-O-remainder term, as e.g. in sin(x)=x - x**3/6 + O(x**4).

Here is what I tryed so far:

Approach 1

Naively, one could just combine the series command twice for each variable, which unfortunately does not work, as this example shows for the function 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)

Approach 2

Based on this post I first wrote a 1D taylor approximation:

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)])

Checking it with 1D examples works fine

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)

However, if I know do a chained call for doing both expansions in x and y, sympy hangs up

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

Does somebody know how to fix this or achieve it in an alternative way?

Answer 1

You can use expr.removeO() to remove the big O from an expression.

---

Oneliner: expr.series(x, 0, 3).removeO().series(y, 0, 3).removeO()

Answer 2

multivariate taylor expansion

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

test for some function

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)