Ask Your Question
1

Taylor Series With Order (Big O)

asked 2017-02-19 23:56:02 +0100

douggard gravatar image

I'm trying to replicate some Mathematica code, and I know very little about Mathematica and Sage.

The Mathematica code creates a Taylor series for f(x) about (x0), with degree 5. I was able to replicate that.

x, x0, dt = var('x, x0, dt')
f = function('f')
t = f(x).taylor(x, x0, 5)

Next it does a substitution.

k1 = f(x0)*dt
s = x0 + k1/2
k2 = t.subs(x=s)

The problem is it that it then removes the higher order terms using the big O notation.

k2 = Expand[g[x] \[CapitalDelta]t /. x -> s] + O[\[CapitalDelta]t]^5

It looks like this is available for series and polynomial rings (no clue what those are) in SageMath, but I haven't figured out how to apply it to the Taylor series expansion.

Is there a way to truncate certain order terms off of an expression?

Is there a way to replicate the Taylor series expression I want with .series() or polynomial rings?

I tried importing sage.rings.big_oh and using Order() but neither seemed applicable to the expression I had. I also made half an attempt to recreate a Taylor series with .series() but didn't quite get what I was hoping for.

tq = 1/factorial(n)*f(x)*x^n
tq.series(n, 5)

Thanks for the help, both my math and SageMath skills are lacking.

edit retag flag offensive close merge delete

2 Answers

Sort by ยป oldest newest most voted
1

answered 2017-02-22 22:18:47 +0100

nbruin gravatar image

While manipulating the expression tree explicitly certainly solves your problem, it might be worthwhile seeing a solution that's a little more general (also: you might want to truncate throughout rather than only at the end, because a lot of time and memory might get devoted to computing the higher order terms that get discarded anyway).

You can get a result by turning the thing into a power series object over SR:

sage: e = x*dt + x*x0*9*dt^2 + x*100*dt^3
sage: R.<DT>=PowerSeriesRing(SR)

Ideally at this point you'd use e(dt=DT) to turn the object into what you want. Unfortunately that doesn't work because of the ambiguity of whether R maps into SR or SR into R. It would be nice if coercion could figure out SR mapping into R here, but currently it doesn't. But we can do it manually without having to go all the way to expression tree manipulation:

sage: E=sum(c*DT^i for c,i in e.coefficients(dt)); E
x*DT + 9*x*x0*DT^2 + 100*x*DT^3

Now we have an object that understands big-Oh:

sage: E+O(DT^3)
x*DT + 9*x*x0*DT^2 + O(DT^3)
edit flag offensive delete link more
0

answered 2017-02-22 01:10:38 +0100

douggard gravatar image

Decided to just write this myself. Found the expr2tree() function here.

x, x0, dt = var('x, x0, dt')

def expr2tree(expr): 
    if expr.operator() is None: 
        return expr 
    else: 
        return [expr.operator()]+map(expr2tree, expr.operands())

def tree2expr(tree):
    args = []
    for ea in tree[1:]:
        if type(ea) is list:
            args.append(tree2expr(ea))
        else:
            args.append(ea)
    return tree[0](*args)

def remove(tree, var, exp):
    if type(tree) != list:
        return False

    if not hasattr(tree[0], '__name__'):
        return False

    if tree[0].__name__ == 'pow':
        if tree[1] == var and tree[2] >= exp:
            return True

    for ea in tree[1:]:
        if remove(ea, var, exp):
            return True

    return False

def truncate_terms(expr, term):
    orig = expr
    expr = expr.expand()
    tree = expr2tree(expr)

    # check is addition
    if tree[0].__name__ != 'add_vararg':
        return expr

    # expand term (eg. dt^5)
    term_tree = expr2tree(term)
    if term_tree[0].__name__ != 'pow':
        return expr

    var, exp = term_tree[1:3]

    args = []
    for ea in tree[1:]:
        if not remove(ea, var, exp):
            args.append(ea)

    return tree2expr([tree[0]] + args)

Test:

e = x*dt + x*x0*9*dt^2 + x*100*dt^3
truncate_terms(e, dt^3)

Output:

9*dt^2*x*x0 + dt*x
edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2017-02-19 23:56:02 +0100

Seen: 1,265 times

Last updated: Feb 22 '17