Here is one way to do it:

https://gist.github.com/moorepants/858503aa180df60a7829

But solve is returning an empty list. I'm not sure why that isn't working.


Jason
moorepants.info
+01 530-601-9791

On Mon, Feb 16, 2015 at 10:28 AM, Ondřej Čertík <[email protected]>
wrote:

> Hi Petr,
>
> On Mon, Feb 16, 2015 at 2:56 AM, Petr Baudis <[email protected]> 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.
>
> >
> > 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))
> >         # -> []
>
> 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
>
> >
> >     (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
> >
> > --
> > 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 [email protected].
> > To post to this group, send email to [email protected].
> > 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.
>
> --
> 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 [email protected].
> To post to this group, send email to [email protected].
> 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/CADDwiVDaF6Q%2BO1RTJRRiEfDiG23bJicPZFCXwq3CQXJfys2XZw%40mail.gmail.com
> .
> 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 [email protected].
To post to this group, send email to [email protected].
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/CAP7f1AiTadyQyDdO%2B3tRGmKG3KaRY%3Dc7v4aPWxnRrqT24RAQdg%40mail.gmail.com.
For more options, visit https://groups.google.com/d/optout.

Reply via email to