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.
