Forming and solving simple physics equation systems

132 views
Skip to first unread message

Petr Baudis

unread,
Feb 16, 2015, 12:22:57 PM2/16/15
to sy...@googlegroups.com
Hello!

This is going to be a rather long post, so I'll first explain what I'm
trying to do - then list all the essentially unsuccessful ways I tried
so far.

Basically, my long-term goal is to build a database of high-school level
physics equations (dynamics, energy, electromagnetism, waves, ...) and
then be able to answer queries like "what is A given X, Y, Z?" with
formulas (or even numerical results). I need to make this completely
automatic, without any human input on which equations are relevant or
how to massage them to apply them to a problem.

However, as an initial toy problem, I'm simply considering:

Let us release a mass object at height $x_0$. Given gravity
uniform acceleration $g$, what will its landing velocity be?

Basically, I'm looking for a way to supply SymPy some basic kinematic
equations and perform a query such that I get something like

sym.Eq(v, sym.sqrt(2 * g * x_0))

on output.


I have not been very successful so far, though! (Using SymPy-0.7.6.)
Either (i) I'm doing various trival novice mistakes, (ii) there are some bugs
or easy-to-add missing features in SymPy preventing me to do this, (iii) what
I want to do is fundamentally hard problem, or (iv) SymPy is in too early
development stage and there is a more appropriate symbolic computation system
to use for this right now. I'd appreciate any advice!


Here's what I tried so far:

(A) Most naive and least flexible approach - define functions s(t) and v(t),
build some simple equations around them and then ask for definition
of v(t) given that s(t) = x_0:

import sympy as sym
import sympy.physics.mechanics as me
from sympy import init_printing

x_0, t, g = sym.symbols('x_0 t g')
ua = sym.symbols('ua', positive=True) # acceleration is function of time, uniform acceleration is not
s, v = sym.symbols('s v', cls=sym.Function) # functions of t
initial = {v(0): 0, s(t): x_0}

velocity_ds_eq = sym.Eq(v(t), s(t) / t)
unifaccel_dv_eq = sym.Eq(ua, (v(t) - v(0))/2/t)
eqs = [eq.subs(initial) for eq in [velocity_ds_eq, unifaccel_dv_eq]]
print(sym.solve(eqs, v(t), implicit=True))
# -> []

Here's the kick - if I apply three rather arbitrary changes to the
last code segment, I *do* get a solution I'm after (thanks to
@vramana for helping me out with this!):

velocity_ds_eq = sym.Eq(v(t) * t, s(t))
unifaccel_dv_eq = sym.Eq(ua * t, (v(t) - v(0))/2)
eqs = [eq.subs(initial) for eq in [velocity_ds_eq, unifaccel_dv_eq]]
print(sym.solve(eqs, [v(t), t], implicit=True))
# -> [(-sqrt(2)*ua*sqrt(x_0/ua), -sqrt(2)*sqrt(x_0/ua)/2), (sqrt(2)*ua*sqrt(x_0/ua), sqrt(2)*sqrt(x_0/ua)/2)]

However, I don't understand why is it required to do either of these
changes - that is, solving for [v(t),t] (hmm, maybe I do have an idea
but not sure) and especially converting /t to *t in both of the
equations (doing that only in the first one will yield
{v(t): x_0/t, t: x_0/(2*t*ua)} - I guess I could deal with that).
The latter two changes give me an impression of a brittle approach
that might break randomly in more complex cases than the toy problem.

(Of course I would prefer to represent velocity and acceleration as
vectors etc., but I got stuck much sooner than that!)


(B) Inspired by "Taming math and physics using SymPy", representing
velocity as integrated acceleration and position as integrated
velocity. I originally aimed to make both functions of time:

import sympy as sym
import sympy.physics.mechanics as me
from sympy import init_printing

t, a, v_0, x_0 = sym.symbols('t a v_0 x_0')
v, x = sym.symbols('v x', cls=sym.Function) # functions of (t)
# initial = {x_0: 100, v_0: 0, a: -9.8}
initial = {v_0: 0}

eqv = sym.Eq(v(t), v_0 + sym.integrate(a, (t, 0,t)))
eqx = sym.Eq(x(t), x_0 + sym.integrate(v(t), (t, 0,t)))
eql = sym.Eq(x(t), 0)
eqs = [eq.subs(initial) for eq in eqx,eqv,eql]
print(sym.solve(eqs, [v(t), t], implicit=True))
# -> []

Is my expectation that the sym.solve() call should yield the
solution I'm after correct?


(C) More literal to "Taming math and physics using SymPy", representing
vi and xi by simple expressions instead of functions:

import sympy as sym
import sympy.physics.mechanics as me
from sympy import init_printing

t, a, v, v_0, x_0 = sym.symbols('t a v v_0 x_0')
initial = {v_0: 0}

vi = v_0 + sym.integrate(a, (t, 0,t))
xi = x_0 + sym.integrate(vi, (t, 0,t))
xeq = sym.Eq(xi, 0)
veq = sym.Eq(vi, v)
print(sym.solve([xeq, veq], [v,t], implicit=True))
# -> {v: (-a*(a*t**2 + 2*x_0)/2 + v_0**2)/v_0, t: -(a*t**2/2 + x_0)/v_0}

Basically, I'm not sure how to phrase my query - especially how to
obtain a solution for v which is independent of t. I thought that is
what the exclude= parameter of sym.solve() is for, but using it has
no effect.

This works, but xeq is completely arbitrary invention, which is why
such a solution is of little use to me (it's not clear to me how to
come up with such equations in a generic way):

xeq = sym.Eq((v*v), ((v*v).expand() - 2*a*x).simplify())
print(sym.solve(xeq, v)[1].subs(initial))
# -> sqrt(2)*sqrt(-a*x_0)


(D) Starting with position function x(t) and expressing velocity and
acceleration as time derivations of that. This was actually my
first approach (when I still hoped things would be easy ;).

I tried to leverage sympy.physics.mechanics for this. I spent
quite a lot of time going through various examples and docs but
I think I still didn't fully get it. I can define the point and
references, describe velocity and acceleration, and solve the
differential equation for x:

import sympy as sym
import sympy.physics.mechanics as me
from sympy import init_printing

N = me.ReferenceFrame('N') # world reference frame
o = me.Point('o') # reference point (origin) for any world positions
o.set_vel(N, 0)
x_0, g = sym.symbols('x_0 g')

r = me.Point('r')
x = me.dynamicsymbols('x')
r.set_pos(o, x * N.y)

t = sym.symbols('t')
r_v = sym.Derivative(x, t)
r.set_vel(N, r_v * N.y)

r_aeq = me.dot(g * N.y - r.acc(N), N.y)
r_xeq = sym.dsolve(r_aeq, x, ics={x.subs(t, 0): x_0, x.diff(t).subs(t, 0): 0})
# -> x(t) == C1 + C2*t + g*t**2/2
# Huh, ics parameter did not work:
r_xeq = r_xeq.replace('C1', x_0).replace('C2', 0)
# -> x(t) == x_0 + g*t**2/2

# x(t) == 0 on landing
# (I would've hoped that [r_xeq, x] would work instead of this:)
r_landing_eq = r_xeq.rhs

# This step is a little arbitrary for my taste already:
r_teq = sym.solve(r_landing_eq, t)
# -> [-sqrt(2)*sqrt(-x_0/g), sqrt(2)*sqrt(-x_0/g)]
r_teq = r_teq[1]

But now I'm stuck - I need to re-express r_v using r_aeq and r_teq.
My few naive attempts mostly using (d)solve failed. How to proceed?


I'd be glad for any assistance or advice. I hope someone finished
reading up to here... ;-)

--
Petr Baudis
If you do not work on an important problem, it's unlikely
you'll do important work. -- R. Hamming
http://www.cs.virginia.edu/~robins/YouAndYourResearch.html

Ondřej Čertík

unread,
Feb 16, 2015, 1:29:00 PM2/16/15
to sympy
Hi Petr,

On Mon, Feb 16, 2015 at 2:56 AM, Petr Baudis <pa...@ucw.cz> wrote:
> Hello!
>
> This is going to be a rather long post, so I'll first explain what I'm
> trying to do - then list all the essentially unsuccessful ways I tried
> so far.
>
> Basically, my long-term goal is to build a database of high-school level
> physics equations (dynamics, energy, electromagnetism, waves, ...) and
> then be able to answer queries like "what is A given X, Y, Z?" with
> formulas (or even numerical results). I need to make this completely
> automatic, without any human input on which equations are relevant or
> how to massage them to apply them to a problem.

I think this might work to some extent, but I suspect some manual
tweaking might be necessary.
I think you have to specify two unknowns to solve for, if you have a
system of two equations, e.g. change v(t) -> [v(t), t], as you did in
your second try:

In [15]: print solve(eqs, [v(t), t])
[(-sqrt(2)*sqrt(ua)*sqrt(x_0), -sqrt(2)*sqrt(x_0)/(2*sqrt(ua))),
(sqrt(2)*sqrt(ua)*sqrt(x_0), sqrt(2)*sqrt(x_0)/(2*sqrt(ua)))]

Then it works.

Why do you use "implicit=True"? It seems it doesn't work with it.

>
> Here's the kick - if I apply three rather arbitrary changes to the
> last code segment, I *do* get a solution I'm after (thanks to
> @vramana for helping me out with this!):
>
> velocity_ds_eq = sym.Eq(v(t) * t, s(t))
> unifaccel_dv_eq = sym.Eq(ua * t, (v(t) - v(0))/2)
> eqs = [eq.subs(initial) for eq in [velocity_ds_eq, unifaccel_dv_eq]]
> print(sym.solve(eqs, [v(t), t], implicit=True))
> # -> [(-sqrt(2)*ua*sqrt(x_0/ua), -sqrt(2)*sqrt(x_0/ua)/2), (sqrt(2)*ua*sqrt(x_0/ua), sqrt(2)*sqrt(x_0/ua)/2)]
>
> However, I don't understand why is it required to do either of these
> changes - that is, solving for [v(t),t] (hmm, maybe I do have an idea
> but not sure) and especially converting /t to *t in both of the
> equations (doing that only in the first one will yield
> {v(t): x_0/t, t: x_0/(2*t*ua)} - I guess I could deal with that).
> The latter two changes give me an impression of a brittle approach
> that might break randomly in more complex cases than the toy problem.

I don't think multiplying by "t" has anything to do with the result,
see my comment above.

After we find a solution that you like, I can look at the other
problems you mentioned.

Ondrej
> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/20150216095625.GQ6082%40machine.or.cz.
> For more options, visit https://groups.google.com/d/optout.

Jason Moore

unread,
Feb 16, 2015, 1:54:32 PM2/16/15
to sy...@googlegroups.com
But solve is returning an empty list. I'm not sure why that isn't working.

Jason Moore

unread,
Feb 16, 2015, 2:10:33 PM2/16/15
to sy...@googlegroups.com
I had a sign error. The code works now.

Btw, sympy.physics.mechanics is likely overkill for high school physics. You can use it to solve these kinds of problem, but the overhead is high.

Petr Baudis

unread,
Feb 16, 2015, 11:52:28 PM2/16/15
to sy...@googlegroups.com


On Tuesday, February 17, 2015 at 3:29:00 AM UTC+9, Ondřej Čertík wrote:
> (A) Most naive and least flexible approach - define functions s(t) and v(t),
>     build some simple equations around them and then ask for definition
>     of v(t) given that s(t) = x_0:
>
>         import sympy as sym
>         import sympy.physics.mechanics as me
>         from sympy import init_printing
>
>         x_0, t, g = sym.symbols('x_0 t g')
>         ua = sym.symbols('ua', positive=True)  # acceleration is function of time, uniform acceleration is not
>         s, v = sym.symbols('s v', cls=sym.Function)  # functions of t
>         initial = {v(0): 0, s(t): x_0}
>
>         velocity_ds_eq = sym.Eq(v(t), s(t) / t)
>         unifaccel_dv_eq = sym.Eq(ua, (v(t) - v(0))/2/t)
>         eqs = [eq.subs(initial) for eq in [velocity_ds_eq, unifaccel_dv_eq]]
>         print(sym.solve(eqs, v(t), implicit=True))
>         # -> []

I think you have to specify two unknowns to solve for, if you have a
system of two equations,

Oh, I see. I thought that it's rather because t is parameter of v(t). So there is a hard requirement
that len(f) == len(symbols)? (Maybe sym.solve() should just throw an exception if that's not the
case then?)

Basically, this requirement is a bit troublesome for me. Basically, I know that:
  - I'm interested in v(t)
  - I don't want to solve for ua and x_0 (these are parameters)
So I don't know that I want to also solve for t, and I don't care about this solution.  Here it might
be simple to just walk the expression and take note that t is the only other possible variable,
but I may also have other completely irrelevant equations for v(t), with other variables. 

What I would have *thought* would work is sym.solve(eqs, exclude=[ua, x_0]), but that returns

> NotImplementedError: could not solve -2*t*ua + v(t)

and it turns out the issue here is that the implicit list of free symbols does not return v(t).
Is this a bug? Or does it make sense not to include undefined functions in free symbols list
in some situations?

In my pretty limited experiments, the issue could have been fixed by unconditionally appending
the auto-generated list of free symbols to the passed list of symbols:
    print(sym.solve(eqs, [v(t), ua, x_0, t]))
    # -> [{v(t): -sqrt(2)*sqrt(ua)*sqrt(x_0), t: -sqrt(2)*sqrt(x_0)/(2*sqrt(ua))}, {v(t): sqrt(2)*sqrt(ua)*sqrt(x_0), t: sqrt(2)*sqrt(x_0)/(2*sqrt(ua))}]

Does this sound sensible?  (It would probably have to be done only in _solve_system(),
specifically after elimination of unrelated equations.)

Then I could in addition pass exclude parameter listing all the parameter symbols for which
I really don't want a solution.

e.g. change v(t) -> [v(t), t], as you did in
your second try:

In [15]: print solve(eqs, [v(t), t])
[(-sqrt(2)*sqrt(ua)*sqrt(x_0), -sqrt(2)*sqrt(x_0)/(2*sqrt(ua))),
(sqrt(2)*sqrt(ua)*sqrt(x_0), sqrt(2)*sqrt(x_0)/(2*sqrt(ua)))]

Then it works.

Why do you use "implicit=True"? It seems it doesn't work with it.

Oh. I'm a bit dumbfounded here, I'm pretty sure I tried that too, but I must have missed this. Thanks!
Obviously I didn't like to use implicit=True (if I understand what it does from the documentation)
but I must have made some error in my experiments without it.
Note that based on the insights above, I tried to replace this with

         print(sym.solve(eqs, [v(t), t, x(t)]))

but that didn't help.

Thanks!
  Petr

Petr Baudis

unread,
Feb 17, 2015, 5:58:36 AM2/17/15
to sy...@googlegroups.com
Oh, thanks a lot! I didn't think about approaching it this way, and it makes a lot of sense.
I'll play with this more, but it's a big step forward for me!

Petr Baudis

unread,
Feb 17, 2015, 11:41:50 PM2/17/15
to sy...@googlegroups.com
On Tue, Feb 17, 2015 at 02:58:35AM -0800, Petr Baudis wrote:
>
> On Tuesday, February 17, 2015 at 3:54:32 AM UTC+9, Jason Moore wrote:
> >
> > Here is one way to do it:
> >
> > https://gist.github.com/moorepants/858503aa180df60a7829
>
> Oh, thanks a lot! I didn't think about approaching it this way, and it
> makes a lot of sense.
> I'll play with this more, but it's a big step forward for me!

I have rewrote your solution a bit to lend better to autogeneration from
predicate form that I use as an intermediate form:

import sympy as sym
from sympy import init_printing
t = sym.symbols('t', real=True)
eqs = []

g = sym.symbols('g', real=True, positive=True)

x_r, v_r = sym.symbols('x_r, v_r', cls=sym.Function) # functions of time
x0_r, v0_r = sym.symbols('x0_r, v0_r', real=True)
accel_r = sym.Eq(x_r(t).diff(t, 2), -g)
position_r = sym.dsolve(accel_r, x_r(t)).subs({'C1': x0_r, 'C2': v0_r})
eqs.append(position_r)
velocity_r = sym.dsolve(accel_r.subs({x_r(t).diff(t): v_r(t)}), v_r(t)).subs({'C1': v0_r})
eqs.append(velocity_r)

h = sym.symbols('h', real=True)
eqs.append(position_r.subs({x_r(t): h, t: 0}))
eqs.append(velocity_r.subs({v_r(t): 0, t: 0}))

t_e2 = sym.symbols('t_e2', real=True)
eqs.append(position_r.subs({x_r(t): 0, t: t_e2}))

print(eqs)
print(sym.solve(eqs, v_r(t)))

So one curious observation I have is that even if sym.solve() is given
multiple equations and only a single variable, it does yield a solution
in this case. I'm still confused about what the exact API contract of
sym.solve() is. Working "sometimes" leaves me unsure about what I can
or cannot rely on...

Of course, the much more important problem for me is that this will
again generate just a formula of v_r(t) depending on t, instead of the
other parameters, and I see no easy way to ask SymPy for a solution
without t in it. Jason's solution is to first compute the time:

# Remove position_r, velocity_r from eqs
sol = sym.solve(eqs, t_e2, x0_r, v0_r, dict=True)
final_velocity = velocity.subs({t: tf}).subs(sol[0])

but that's the problem - I need to know I need to compute the time
first, and it's not clear what kind of rules to figure out what to
compute I will need with different combinations of equations. (In other
words, I would have hoped a CAS would have this figured out already. :-)

(The other thing of also including x0_r, v0_r in solve set is not so
bad, that might be fixed by a simple semi-automated substitution step.)

Right now I'm inclined to conclude that what I want to do *is* rather
hard and I'd have to come up with some new algorithms to deal with this.
(Or likely take a look at some other CASes first if any handle this case
for me already.)

Aaron Meurer

unread,
Feb 18, 2015, 12:02:16 PM2/18/15
to sy...@googlegroups.com
Sadly, there is none. Part of this is due to algorithmic limitations
(solve() is basically a bunch of heuristics, there are few guarantees
that it will find a solution if it exists). But a lot of it is just
poor design. We are trying to make the design better with the new
solveset module.

Aaron Meurer

>
> Of course, the much more important problem for me is that this will
> again generate just a formula of v_r(t) depending on t, instead of the
> other parameters, and I see no easy way to ask SymPy for a solution
> without t in it. Jason's solution is to first compute the time:
>
> # Remove position_r, velocity_r from eqs
> sol = sym.solve(eqs, t_e2, x0_r, v0_r, dict=True)
> final_velocity = velocity.subs({t: tf}).subs(sol[0])
>
> but that's the problem - I need to know I need to compute the time
> first, and it's not clear what kind of rules to figure out what to
> compute I will need with different combinations of equations. (In other
> words, I would have hoped a CAS would have this figured out already. :-)
>
> (The other thing of also including x0_r, v0_r in solve set is not so
> bad, that might be fixed by a simple semi-automated substitution step.)
>
> Right now I'm inclined to conclude that what I want to do *is* rather
> hard and I'd have to come up with some new algorithms to deal with this.
> (Or likely take a look at some other CASes first if any handle this case
> for me already.)
>
> --
> Petr Baudis
> If you do not work on an important problem, it's unlikely
> you'll do important work. -- R. Hamming
> http://www.cs.virginia.edu/~robins/YouAndYourResearch.html
>
> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/20150218044144.GI6082%40machine.or.cz.

Ondřej Čertík

unread,
Feb 18, 2015, 1:57:45 PM2/18/15
to sympy
Hi Petr,
If you find some other CAS that can do this better, please let us know
and post here your solution. As Aaron said, we would like to improve
our solver module.

Thanks,
Ondrej

>>
>> --
>> Petr Baudis
>> If you do not work on an important problem, it's unlikely
>> you'll do important work. -- R. Hamming
>> http://www.cs.virginia.edu/~robins/YouAndYourResearch.html
>>
>> --
>> You received this message because you are subscribed to the Google Groups "sympy" group.
>> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
>> To post to this group, send email to sy...@googlegroups.com.
>> Visit this group at http://groups.google.com/group/sympy.
>> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/20150218044144.GI6082%40machine.or.cz.
>> For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6LcGaCDCJWcvHwsrTq39nEuYY8xopbYLesa839bHk6MGA%40mail.gmail.com.

Petr Baudis

unread,
Feb 19, 2015, 7:35:27 AM2/19/15
to sy...@googlegroups.com
On Wed, Feb 18, 2015 at 11:01:54AM -0600, Aaron Meurer wrote:
> On Tue, Feb 17, 2015 at 10:41 PM, Petr Baudis <pa...@ucw.cz> wrote:
> > So one curious observation I have is that even if sym.solve() is given
> > multiple equations and only a single variable, it does yield a solution
> > in this case. I'm still confused about what the exact API contract of
> > sym.solve() is. Working "sometimes" leaves me unsure about what I can
> > or cannot rely on...
>
> Sadly, there is none. Part of this is due to algorithmic limitations
> (solve() is basically a bunch of heuristics, there are few guarantees
> that it will find a solution if it exists). But a lot of it is just
> poor design. We are trying to make the design better with the new
> solveset module.

Thanks! Knowing this is quite valuable for me.

> >> Right now I'm inclined to conclude that what I want to do *is* rather
> >> hard and I'd have to come up with some new algorithms to deal with this.
> >> (Or likely take a look at some other CASes first if any handle this case
> >> for me already.)
>
> If you find some other CAS that can do this better, please let us know
> and post here your solution. As Aaron said, we would like to improve
> our solver module.

I will for sure - my plan is to take a hard look at Maple and a brief
look at Mathematica. If I won't be successful, I might come back to
SymPy and try to help improve the solveset module (I can't promise
I will have the time though, this was supposed to be just a small piece
of puzzle in a larger project - http://21robot.org/).

Thanks,

Petr Baudis

Ondřej Čertík

unread,
Feb 19, 2015, 4:43:54 PM2/19/15
to sympy
Are you trying to do all of this:

http://21robot.org/research_activities/math/

? If so, that's a *major* project. There is Wolfram Alpha
(http://www.wolframalpha.com/), that has some understanding
of a human language, e.g.:

http://www.wolframalpha.com/input/?i=integrate+cos^2x+from+0+to+2pi&lk=4&num=1

But in my experience once things get a little more complicated, it
can't understand it either, e.g.:

http://www.wolframalpha.com/input/?i=Solve%5B{v%5Bt%5D+%3D%3D+x0%2Ft%2C+a+%3D%3D+v%5Bt%5D%2F%282*t%29}%2C+{v%5Bt%5D%2C+t}%5D

The same expression in Mathematica returns the correct answer (as in SymPy).

Ondrej

Petr Baudis

unread,
Feb 20, 2015, 11:57:18 AM2/20/15
to sy...@googlegroups.com
Actually, I'm working on physics questions (not sure why these aren't
listed on the homepage). The nature of the questions is not always
the same, but SAT questions are reasonable examples of what we are
trying to solve, e.g.:

http://papers.xtremepapers.com/SAT/SAT%20II%20Success%20Physics.pdf.pdf

The pipeline should convert text to logical form (set of predicates
that describe the situation and the question) and currently the logical
form is transformed to Modelica statements and the situation is
simulated. I'm trying to explore an alternate avenue of converting the
predicates to a system of equations and solving them symbolically.

Ondřej Čertík

unread,
Feb 20, 2015, 12:16:50 PM2/20/15
to sympy
I see, thanks for the clarification. To me that sounds like a major
project. Definitely keep us posted how it goes.

Ondrej
Reply all
Reply to author
Forward
0 new messages