Question about diff() and Derivative()

[Gilbert Gede]

Hi, 

I was trying to implement some functionality for PyDy for this year's GSoC, and was looking for some advice.  

In dynamics problems, you usually have time-varying quantities, like generalized coordinates, speeds, and accelerations.  Often, you want to take the partial derivative of an expression with respect to the time derivative of one of these quantities.  This come up when using Lagrange's Method (or Kane's Method).  It's described to some degree here:  

http://en.wikipedia.org/wiki/Lagrangian_mechanics

https://gist.github.com/1005937

In Lagrange's Method, you end up taking the partial derivative of the energy with respect to the time derivative of a generalized coordinate.  I'm trying to figure out a way to make this work in PyDy/SymPy. Derivative won't take in anything but a Symbol.  

The only idea I have come up with is to extend Symbol and write my own .diff() method for it which returns a new symbol representing the time differentiation of the original extended Symbol.  Once my new object is inside a Mul or Add sympy object, then my .diff() method is no longer called.  

Can anyone give some insight into how I could get this desired behavior, taking the derivative of an expression wrt a time-differentiated symbol, to work in a way consistent with existing SymPy behavior?  Thanks.  

-Gilbert

[Aaron S. Meurer]

This has actually been discussed quite a bit before (a lot of people want to use Lagrangians).  You can search the mailing list.  From what I've seen, you will either have to write your own custom diff routine or do clever substitution of functions and derivatives with symbols.  I don't think I've ever seen anyone suggest extending Symbol to hold a time derivative, which is essentially just a more formal way of doing the substation method.  It might work.

Aaron Meurer

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sympy?hl=en.

[Gilbert gede]

Yeah, I had read some of them. I had already thought of writing my
own Diff method or something and do substitution with it, but was
hoping to have the functionality I want work like standard SymPy
operations. That's what I've been trying to do with my PyDy classes;
make them work more like you would expect other SymPy objects to.
I have read through Derivative() and diff(), and couldn't really find
a way to make them do what I want (like I said about my symbol
extension no longer having its methods called once it is inside a
SymPy add or mul). I guess what I was hoping for was input on whether
I could make Derivative do what I want with my extended Symbol, as I
couldn't really see how. But if writing my own Diff method is the
only option, there's not much I can do then.

Thanks,
-Gilbert

[Gilbert gede]

I guess I should have asked this as well; is it considered bad to
write diff() and replace the current Sympy diff() within my code? Or
is that OK?

-Gilbert

[Aaron Meurer]

Can you just get what you want by overriding _eval_derivative()?

Aaron Meurer

[Gilbert gede]

You're talking about Symbol._eval_derivative? I tried that within my
extended Symbol class. It returns 0 or 1 testing self == symbol. I
tried making some changes to it, but I don't think I can use it. I
think it doesn't even get called unless you do something like:
t = timevaryingsymbols('t')
Derivative(2+3*t,t)
I think only when t is both part of (expr, symbols, ...) within
Derivative's __new__ definition does t's _eval_derivative() method get
called.

-Gilbert

[Aaron Meurer]

I'm not sure I follow. If I create a test _eval_derivative method to
Symbol, like

diff --git a/sympy/core/symbol.py b/sympy/core/symbol.py
index 635c2c5..e1d36b1 100644
--- a/sympy/core/symbol.py
+++ b/sympy/core/symbol.py
@@ -24,6 +24,11 @@ class Symbol(AtomicExpr, Boolean):

"""

+ def _eval_derivative(self, s):
+ from sympy.core.function import Derivative
+ print 'test'
+ return Derivative(self, s)
+
is_comparable = False

__slots__ = ['is_commutative', 'name']

Then I get

In [1]: diff(x, x)
test
Out[1]:
d
──(x)
dx

So it's definitely being called. It should be the same for any
subclass of Symbol.

The idea, for what you want to do, would be to make x.diff(t) return
some kind of object instead of 0 (i.e., make x kind of like x(t), but
the object returned would be a Symbol subclass, so that it can be
passed to diff()).

Aaron Meurer

[Ronan Lamy]

Le vendredi 03 juin 2011 à 19:16 -0700, Gilbert gede a écrit :
> You're talking about Symbol._eval_derivative? I tried that within my
> extended Symbol class. It returns 0 or 1 testing self == symbol. I
> tried making some changes to it, but I don't think I can use it. I
> think it doesn't even get called unless you do something like:
> t = timevaryingsymbols('t')
> Derivative(2+3*t,t)
> I think only when t is both part of (expr, symbols, ...) within
> Derivative's __new__ definition does t's _eval_derivative() method get
> called.

That seems to be a recent "optimisation" from commit 2361dd86. You
should revert this to the old behaviour: call
expr._eval_derivative(symbol) in all cases.

[Gilbert gede]

Yeah, if you try diff(x,y), the method is not called.
Ronan, I'm not sure how to do what you're describing? How would I
call that?

-Gilbert

[Aaron Meurer]

I see. Yeah, that's because of the optimization built in to
Derivative.__new__ that does diff(expr, x) => 0 if x is not in
expr.free_symbols. This should probably be factored out to just be a
default _eval_derivative action in Expr, so that subclasses can
override it. On the other hand, if your object is going to be
something akin to x(t) or x'(t), then t should be in the free_symbols.

Aaron Meurer

[Renato Coutinho]

I'm probably missing something basic, but why not use functions x(t),
v(t) etc. in the first place?

Renato

[Aaron Meurer]

The problem is that with the Euler-Langrange formulas, you need to be
able to do something like expr.diff(f(t)), which is not allowed by
SymPy because f(t) is not a Symbol. You also need to be able to do
expr.diff(f(t).diff(t)). So the idea is to have something that acts
like f(t) or f(t).diff(t) but looks and acts like a Symbol to
Derivative.

Aaron Meurer

[Renato Coutinho]

On Sat, Jun 4, 2011 at 12:13 AM, Aaron Meurer <asme...@gmail.com> wrote:
> The problem is that with the Euler-Langrange formulas, you need to be
> able to do something like expr.diff(f(t)), which is not allowed by
> SymPy because f(t) is not a Symbol.  You also need to be able to do
> expr.diff(f(t).diff(t)).  So the idea is to have something that acts
> like f(t) or f(t).diff(t) but looks and acts like a Symbol to
> Derivative.

I see. But then there needs to be some way to indicate when a
derivative is partial or total, right? Because dx/dt = x'(t) but
\partial x/\partial t = 0. A subclass of Symbol can handle the first
overloading _eval_derivative and free_symbols, and keeping "internal
variables" on which it depends (like time), The second is harder
though.

Renato

[Aaron Meurer]

I don't understand. How is ∂x/∂t != dx/dt?

Aaron Meurer

[Tim Lahey]

On Fri, Jun 3, 2011 at 8:05 PM, Gilbert Gede <gilbe...@gmail.com> wrote:
> Can anyone give some insight into how I could get this desired behavior,
> taking the derivative of an expression wrt a time-differentiated symbol, to
> work in a way consistent with existing SymPy behavior?  Thanks.
> -Gilbert
>

I implemented Euler-Lagrange equations in Maple. I've put the code for
doing it below. Basically, I create symbols for each of the functions
and substitute these. Then, I take the derivative, and revert back.
This code is quite fast in Maple (at the time I originally wrote the
code, it was more than 2 orders of magnitude faster than Maple's
equivalent) and I think most of the functionality necessary is in
Python and SymPy. It should be fairly straightforward to convert to
SymPy. This code handles any number of variables (as a list), but it
only does the traditional Euler-Lagrange equation, not some of the
others that can arise due to Calculus of Variations.

I hope this helps.

Cheers,

Tim.

---
# Calculate the Mass (d/dt(dL_dqv)) and Stiffness (dL_dq) contributions to the
# Euler-Lagrange equation. Note that the Stiffness contribution does
not include the negative sign.
EulerLagrange := proc(Lagrangian::anything, variables::list)
local num_list, qv_name, vel_var, qv_subs, qv_unsubs, Lagrange_subs1,
Lagrange_subs2, dL_dqv1, dL_dqv2, dL_dqv, dL_dqvt, dL_dq, dL_dq1,
dL_dq2, dL_dq3, q_name, q_subs, q_unsubs:

# create a list of indices from 1 to the number of variables used in
the formulation
num_list := [seq(i,i=1..nops(variables))]:

# Define a list of generalized velocity and position variables
qv_name := map2(cat,qv,num_list):
q_name := map2(cat,q,num_list):

# Equate the time derivatives of the system variable to the
generalized velocities
# also define the reverse mapping
vel_var := map(diff,variables,t):
qv_subs := zip(equate,vel_var,qv_name):
qv_unsubs := zip(equate,qv_name,vel_var):

# Equate the generalized positions to the system variables and define
the reverse mapping
q_subs := zip(equate,variables,q_name):
q_unsubs := zip(equate,q_name,variables):

# Convert the Lagrangian to the generalized position and velocity variables
Lagrange_subs1 := subs(qv_subs,Lagrangian):
Lagrange_subs2 := subs(q_subs,Lagrange_subs1):

# Differentiate the Lagrangian with respect to the generalized
velocities and positions
dL_dqv1 := map2(diff,Lagrange_subs2,qv_name):
dL_dq1 := map2(diff,Lagrange_subs2,q_name):

# Revert back to the system variables
dL_dq2 := map2(subs,qv_unsubs,dL_dq1):
dL_dqv2 := map2(subs,qv_unsubs,dL_dqv1):
dL_dqv := map2(subs,q_unsubs,dL_dqv2):
dL_dq := map2(subs,q_unsubs,dL_dq2):
dL_dqvt := map(diff,dL_dqv,t):

# Return the two components of the Euler-Lagrange Equation
return (dL_dqvt, dL_dq):
end proc:

--
Tim Lahey
PhD Candidate, Systems Design Engineering
University of Waterloo
http://about.me/tjlahey

[Renato Coutinho]

On Sat, Jun 4, 2011 at 12:49 AM, Aaron Meurer <asme...@gmail.com> wrote:
> I don't understand.  How is ∂x/∂t != dx/dt?

Sorry, I thought about one thing and wrote another. I was thinking
something like L(v, x, t) = v**2/2 - x**2/2, then take dL/dt = v *
v'(t) + x * x'(t), but ∂L/∂t = 0. I know one should explicit in a
partial derivative which other variables are to be held constant, but
this "convention" was used when I took classes on mechanics, and it
really make things shorter to write -- and anyway, how woud you write
(∂L/∂t)_{v, x} in sympy?

Renato

[Gilbert gede]

Yeah, that could work.
I was hoping that I could make it work in a way which was hidden to
the user; that just using diff would work. Writing a method (or
overriding one) to do this pre- and post-substitution is the only
solution that I currently understand how to do; but again, it is
undesirable.

-Gilbert


On Jun 3, 9:09 pm, Tim Lahey <tim.la...@gmail.com> wrote:

[Renato Coutinho]

On Sat, Jun 4, 2011 at 1:18 AM, Gilbert gede <gilbe...@gmail.com> wrote:
> Yeah, that could work.
> I was hoping that I could make it work in a way which was hidden to
> the user; that just using diff would work.  Writing a method (or
> overriding one) to do this pre- and post-substitution is the only
> solution that I currently understand how to do; but again, it is
> undesirable.

Using a class that overrides some methods will only be visible to
users when creating the variables, it could be something like "x =
DynamicSymbol(x, t)" instead of "x = Symbol(x)". All the rest would be
just the same, I believe. Is this too "user-annoying"?

Renato

[Gilbert gede]

Right now, I made a timevaryingsymbol class; it's diff is indeed
rewritten and works. The problem is once you put that inside an add
or mul, it's diff method no longer gets called.


On Jun 3, 9:28 pm, Renato Coutinho <renato.couti...@gmail.com> wrote:

[Renato Coutinho]

On Sat, Jun 4, 2011 at 1:44 AM, Gilbert gede <gilbe...@gmail.com> wrote:
> Right now, I made a timevaryingsymbol class; it's diff is indeed
> rewritten and works.  The problem is once you put that inside an add
> or mul, it's diff method no longer gets called.

I think this can be solved overloading the free_symbols method in your
class, making it always include the time variable, besides itself.

Renato

[Vinzent Steinberg]

On 4 Jun., 06:44, Gilbert gede <gilbertg...@gmail.com> wrote:
> Right now, I made a timevaryingsymbol class; it's diff is indeed
> rewritten and works.  The problem is once you put that inside an add
> or mul, it's diff method no longer gets called.

What does it differently than Function('name')(t)?

The most obvious way I see to implement this using current sympy
functionality is using Function('...')(t) and using subs() when taking
partial derivatives.

Vinzent

[Chris Smith]

It would help me understand if you gave an example input/output pair so I could see what you were wanting to do. I did look at the gist but didn't see the output that you wanted. Would idiff be of any use?

    >>> help(idiff)

    Help on function idiff in module sympy.geometry.util:

    idiff(eq, y, x, dep=None)

        Return dy/dx assuming that y and any other variables given in dep

        depend on x.

        >>> from sympy.abc import x, y, a

        >>> from sympy.geometry.util import idiff

        >>> idiff(x**2 + y**2 - 4, y, x)

        -x/y

        >>> idiff(x + a + y, y, x)

        -1

        >>> idiff(x + a + y, y, x, [a])

        -D(a, x) - 1

[Gilbert gede]

Adding Symbol('t') to the free_symbols property of my class worked!
Here's what it can do now:

In [1]: run kinematics.py

In [2]: x = TVS('x')

In [3]: y = TVS('y')

In [4]: t = Symbol('t')

In [5]: diff(2 * x**2 + 4 + y,t)
Out[5]: 4*x*xdt + ydt

In [7]: diff(diff(2 * x**2 + 4 + y,t),TVS('xdt'))
Out[7]: 4*x


TVS is timevaryingsymbol, athough I think Renato's DynamicSymbol
sounds better.
I think this is exactly what I wanted; Thanks!

-Gilbert

[Vinzent Steinberg]

On 4 Jun., 19:48, Gilbert gede <gilbertg...@gmail.com> wrote:
> Adding Symbol('t') to the free_symbols property of my class worked!
> Here's what it can do now:
>
> In [1]: run kinematics.py
>
> In [2]: x = TVS('x')
>
> In [3]: y = TVS('y')
>
> In [4]: t = Symbol('t')
>
> In [5]: diff(2 * x**2 + 4 + y,t)
> Out[5]: 4*x*xdt + ydt
>
> In [7]: diff(diff(2 * x**2 + 4 + y,t),TVS('xdt'))
> Out[7]: 4*x
>
> TVS is timevaryingsymbol, athough I think Renato's DynamicSymbol
> sounds better.
> I think this is exactly what I wanted; Thanks!

You can also make it pretty print unicode with dots above the symbol:

>>> print u"x\u0307"
ẋ
>>> print u"x\u0308"
ẍ

Vinzent