8

Given some f and the differential equation x'(t) = f(x(t)), how do I compute x(n)(t) in terms of x(t)?

For example, given f(x(t)) = sin(x(t)), I want to obtain x(3)(t) = (cos(x(t))2 − sin(x(t))2) sin(x(t)).

So far I've tried

>>> from sympy import diff, sin
>>> from sympy.abc import x, t
>>> diff(sin(x(t)), t, 2)

which gives me

-sin(x(t))*Derivative(x(t), t)**2 + cos(x(t))*Derivative(x(t), t, t)

but I'm not sure how to tell SymPy what Derivative(x(t), t) is and have it figure out Derivative(x(t), t, t), etc. automatically.


Answer:

Here's my final solution based on the answers I received below:

def diff(x_derivs_known, t, k, simplify=False):
    try: n = len(x_derivs_known)
    except TypeError: n = None
    if n is None:
        result = sympy.diff(x_derivs_known, t, k)
        if simplify: result = result.simplify()
    elif k < n:
        result = x_derivs_known[k]
    else:
        i = n - 1
        result = x_derivs_known[i]
        while i < k:
            result = result.diff(t)
            j = len(x_derivs_known)
            x0 = None
            while j > 1:
                j -= 1
                result = result.subs(sympy.Derivative(x_derivs_known[0], t, j), x_derivs_known[j])
            i += 1
            if simplify: result = result.simplify()
    return result

Example:

>>> diff((x(t), sympy.sin(x(t))), t, 3, True)
sin(x(t))*cos(2*x(t))
Community
  • 1
  • 1
user541686
  • 205,094
  • 128
  • 528
  • 886
  • 1
    You should post your answer separately from your question. – user Nov 01 '16 at 18:00
  • @Fermiparadox: It's not really my answer, it's just a description of what I gathered from the other answers. I didn't want it to take away potential votes from existing answers... – user541686 Nov 01 '16 at 18:21
  • 1
    @Mehrdad you should still post it as an answer. It doesn't take away votes (you are allowed to vote on multiple answers). – asmeurer Nov 01 '16 at 19:08
  • @asmeurer: I meant people might vote up my answer instead of theirs because it's better (people often just vote for the best answer) even though it's based on theirs. I don't see the point either. – user541686 Nov 01 '16 at 19:38

2 Answers2

4

Here is one approach that returns a list of all derivatives up to n-th order

import sympy as sp

x = sp.Function('x')
t = sp.symbols('t')

f = lambda x: x**2 #sp.exp, sp.sin
n = 4 #3, 4, 5

deriv_list = [x(t), f(x(t))]  # list of derivatives [x(t), x'(t), x''(t),...]
for i in range(1,n):
    df_i = deriv_list[-1].diff(t).replace(sp.Derivative,lambda *args: f(x(t)))
    deriv_list.append(df_i)

print(deriv_list)

[x(t), x(t)**2, 2*x(t)**3, 6*x(t)**4, 24*x(t)**5]

With f=sp.sin it returns

 [x(t), sin(x(t)), sin(x(t))*cos(x(t)), -sin(x(t))**3 + sin(x(t))*cos(x(t))**2, -5*sin(x(t))**3*cos(x(t)) + sin(x(t))*cos(x(t))**3]

EDIT: A recursive function for the computation of the n-th derivative:

def der_xt(f, n):
    if n==1:
        return f(x(t))
    else:
        return der_xt(f,n-1).diff(t).replace(sp.Derivative,lambda *args: f(x(t)))

print(der_xt(sp.sin,3))

-sin(x(t))**3 + sin(x(t))*cos(x(t))**2

Stelios
  • 5,271
  • 1
  • 18
  • 32
3

Declare f and use substitution:

>>> f = diff(x(t))
>>> diff(sin(x(t)), t, 2).subs(f, sin(x(t)))
-sin(x(t))**3 + cos(x(t))*Derivative(sin(x(t)), t)
Uriel
  • 15,579
  • 6
  • 25
  • 46
  • Substitution doesn't work generally. Try replacing `2` with `4`. – user541686 Oct 29 '16 at 22:47
  • yea. you mean somethiing like automatic conversation of every integral and deifferential? no such thing on sympy. you need to tell it you want the representation substituted by some predefined function, otherwise it will always want to convert to full representation. also, `f` is a function. not some integer. it replaces the core sympy expression in the evaluated expression. – Uriel Oct 29 '16 at 22:53
  • I think substitution should actually work fine if SymPy only used first derivatives (i.e. didn't convert `Derivative(Derivative(x(t),t),t)` into `Derivative(x(t),t,t)` automatically). Is there no way to tell it that either? – user541686 Oct 29 '16 at 23:00
  • they do not specify such way, and i think it goes against their build of the classes. as i tested it, it automatically transform the `Derivative(Derivative(x(t),t),t)` into `Derivative(x(t),t,t)` (when typing it in shell) – Uriel Oct 29 '16 at 23:04