14

In sage it is fairly easy to do a Taylor expansion of an unknown function f(x),

x = var('x')
h = var('h')
f = function('f',x)
g1 = taylor(f,x,h,2)

How can this be done in sympy?

Update

asmeurer points out that this is a feature which will be available soon in sympy from the pull request http://github.com/sympy/sympy/pull/1888. I installed the branch using pip,

pip install -e git+git@github.com:renatocoutinho/sympy.git@897b#egg=sympy --upgrade

However, when I try to calculate the series of f(x),

x, h = symbols("x,h")
f = Function("f")
series(f,x,x+h)

I get the following error,

TypeError: unbound method series() must be called with f instance as first argument (got Symbol instance instead)

Løiten
  • 3,185
  • 4
  • 24
  • 36
Daniel Farrell
  • 9,316
  • 8
  • 39
  • 62
  • You can not. Just use a loop and `diff`. And the function is called `series` not `taylor`. – Krastanov Jun 01 '13 at 09:14
  • 1
    There is a pull request to make this work at https://github.com/sympy/sympy/pull/1888. – asmeurer Jun 02 '13 at 19:45
  • @asmeurer fantastic! This hasn't been merged into the master branch, is it still possible to install this with pip. Or do I need to clone the repo, apply the patch and then build from source? – Daniel Farrell Jun 03 '13 at 10:47
  • @asmeurer I tried to apply https://github.com/sympy/sympy/pull/1888.patch to the latest master copy of sympy using git apply but that failed. What do you recommend? – Daniel Farrell Jun 03 '13 at 13:13
  • The easiest way is to just work off his branch, until that pull request is merged. – asmeurer Jun 03 '13 at 14:24
  • I just updated with what I tried, but the code returns an error. I am not too familiar with sympy but I am probably using the library incorrectly or there is something wrong with the installed branch. – Daniel Farrell Jun 04 '13 at 23:59
  • I`m currently working on a code that does something similar. I plan on providing it to the sympy community when I`m more sure it`s results are accurate. it'll need some adapting though.. – juggler Jan 10 '14 at 18:24
  • @juggler that sounds good, please post the pull commit here, I would be interested in taking a look. – Daniel Farrell Jan 13 '14 at 06:33
  • sure, when I figure out how :-) I'll prioritize first continuing to make sure the results are accurate. I look forward to seeing what it looks like once (if) it's eventually integrated.. – juggler Jan 13 '14 at 14:29
  • still working on it.. (amoung other things) ..so other people don't go about develloping the same thing needlessly, mine takes multiple taylor series involving non-commuting operators, expands them, cross-multiplies them, takes differences, and returns a result to a given order, symbolically, as a pdf.. it's quite quick. – juggler Feb 24 '14 at 17:54
  • @boyfarrell -see above comment.. – juggler Feb 24 '14 at 18:00

2 Answers2

19

As @asmeurer described, this is now possible with

from sympy import init_printing, symbols, Function
init_printing()

x, h = symbols("x,h")
f = Function("f")

pprint(f(x).series(x, x0=h, n=3))

or

from sympy import series
pprint(series(f(x), x, x0=h, n=3))

both returns

                                              ⎛  2        ⎞│                          
                                            2 ⎜ d         ⎟│                          
                                    (-h + x) ⋅⎜────(f(ξ₁))⎟│                          
                                              ⎜   2       ⎟│                          
                ⎛ d        ⎞│                 ⎝dξ₁        ⎠│ξ₁=h    ⎛        3       ⎞
f(h) + (-h + x)⋅⎜───(f(ξ₁))⎟│     + ──────────────────────────── + O⎝(-h + x) ; x → h⎠
                ⎝dξ₁       ⎠│ξ₁=h                2

If you want a finite difference approximation, you can for example write

FW = f(x+h).series(x+h, x0=x0, n=3)
FW = FW.subs(x-x0,0)
pprint(FW)

to get the forward approximation, which returns

                                  ⎛  2        ⎞│
                                2 ⎜ d         ⎟│
                               h ⋅⎜────(f(ξ₁))⎟│
                                  ⎜   2       ⎟│
          ⎛ d        ⎞│           ⎝dξ₁        ⎠│ξ₁=x₀    ⎛ 3    2        2    3                 ⎞
f(x₀) + h⋅⎜───(f(ξ₁))⎟│      + ────────────────────── + O⎝h  + h ⋅x + h⋅x  + x ; (h, x) → (0, 0)⎠
          ⎝dξ₁       ⎠│ξ₁=x₀             2
Løiten
  • 3,185
  • 4
  • 24
  • 36
9

There is no function for this in sympy, but it's rather easy to do it "by hand":

In [3]: from sympy import *
        x, h = symbols('x, h')
        f = Function('f')
        sum(h**i/factorial(i) * f(x).diff(x, i) for i in range(4))

Out[3]: h**3*Derivative(f(x), x, x, x)/6 + h**2*Derivative(f(x), x, x)/2 + h*Derivative(f(x), x) + f(x)

Note that sympy typically works with expressions (like f(x)) and not with bare functions (like f).

Ronan Lamy
  • 577
  • 3
  • 9