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
--
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.
