[WIP] Object model for vibrational concepts

[moorepants]

My current thought is that if we want these materials to standout and be unique in the way vibrations is taught that we do need to develop an object model to go along with the text. Here are some ideas to start with.

We need at least two types of parameters: time varying (SymPy undefined functions of time) and constant (SymPy symbols).

NOTE : Does it make sense to assign a value to a time varying parameter? If so, does the value of time also need to be indicated?

>>> from sympy import Eq
>>> from resonance import meter, second, Parameter
>>> x = Parameter('x', units=meter, description='position')
>>> x
x(t) : position in m
>>> v = Parameter('v', value=2.0, units=meter / second, description='speed')
>>> v
v(t) = 2.0 : speed in m/s

There is the concept of a system's state. The state tells you what the numerical values are at a specificed time and should be mutable. As time changes we can update the state values. It is desirable to preserve the order of the state variables, so that matrix operations work nicely in subsequent code.

>>> s = State([Eq(x, 1.0), Eq(v, 2.0)], time=0.1, units=second)
>>> s
|x(t)|   |1.0| : position in m
|    | = |   |
|v(t)|   |2.0| : speed in m/s
t = 0.1 s
>>> s.x.sym
x(t)
>>> s.x.val
1.0
>>> s.time.val
0.1
>>> s.time.sym
t
>>> s.vals
[1.0, 2.0]
>>> s.vals_as_array
array([1.0, 2.0])
>>> s.syms
[x(t), v(t)]
>>> s.syms_as_matrix
Matrix([[x(t)],
        [v(t)]])

We also need an object that contains a monotonic collection of state values. Maybe a subclass of a Pandas DataFrame is appropriate (NOTE : I've had little luck successfuly subclassing a DataFrame in the past.)

>>> times = linspace(0, 5, num=6)
>>> st = TimeFrame(times, s)
t  x    v
0  1.0  2.0
1  1.0  2.0
2  1.0  2.0
3  1.0  2.0
4  1.0  2.0
5  1.0  2.0
>>> st.values
array([[1.0, 1.0, 1.0, 1.0, 1.0],
       [2.0, 2.0, 2.0, 2.0, 2.0]])
>>> st.plot()  # plots a line chart versus time

A row would then be a State instance.

The students will have to specify the equations of motion manually in symbolic form. This forces them to design the model versus some kind of approach that assembles models from Springs, Masses, etc. This would be the result of Lagrange's method in our case.

>>> m = ConstantParameter('m', value=1.0, units=kilogram, description='mass')
>>> m
m = 1.0 : mass in kg
>>> c = ConstantParameter('c', value=1.0, units=newton * second / meter, description='viscous damping coefficient')
>>> k = ConstantParameter('k', value=1.0, units=newton / meter, description='stiffness')
>>> eom = Eq(m * x.diff(t, 2) + c * x.diff(t) + k * x, 0)
>>> eom
c*x'(t) + k*x(t) + m*x''(t) = 0

There is also the concept of a model that contains the mathematical equations that attempt to describe reality.

>>> m = Model(odes=eom, aux_var=Eq(v, x.diff(t)))
>>> m.state
|x(t)|   |0.0| : position in m
|    | = |   |
|v(t)|   |2.0| : speed in m/s
t = 0.1 s
>>> m.constants
[m, c, k]
>>> m.inputs  # these would be any time varying forcing parameters
[]

>>> m.first_order_eom
[x'(t)] = [v(t)                ]
[v'(t)]   [-c*v(t)/m - k*x(t)/m]
>>> m.eval_first_order(x=1, v=4, t=1.5)
|x(t)|   |1.0| : position in m
|    | = |   |
|v(t)|   |4.0| : speed in m/s
t = 1.5 s

I'd like the students to be able to use custom integrators, even one they might write themselves.

>>> def euler_integrate(rhs, init_cond, times):
...     state_traj = np.zeros((len(times), *init_cond.shape))
...     state_traj[0] = init_cond
...     for i, t in enumerate(times[:-1]):
...         h = t - times[i - 1]
...         state_traj[i + 1] = states_traj[i] + h * rhs(t, states_traj[i])
...     return state_traj
>>> m.integrator = euler_integrator

# simulate the sytem given some monotonic time values
>>> m.evolve_state([0, 1, 2, 3])
array([[?, ?, ?, ?],
       [?, ?, ?, ?]])
# maybe this should return a StateEvolution object

In many cases, we'd like to linearize a non-linear model about an equilibrium point. If the symbolic form of the model i

>>> theta = Parameter('theta', value=0.0)
>>> omega = Parameter('omega', value=0.0)
# passing in a list to State preserves the state order
>>> s = State([theta, omega])
>>> s.theta_val
0
>>> s.theta_sym
theta
>>> s.state_vals
|0|
|0|
>>> s.state_syms
[theta]
[omega]
>>> s.state_vec
|θ| = |0|
|ω|   |0|
>>> eom = Eq(omega.diff() + g / l * sin(theta), 0)
>>> pendulum_model = Model(eom)
# linearize the system (if not already linear) about an equilibrium state,
# returns a LinearModel
>>> lin_m = pendulum_model.linearize(
>>> lin_m.eom
omega' + g / l * theta = 0
>>> lin_m.first_order_form()
[theta'(t)] = [omega(t)      ]
[omega'(t)]   [-g / l * theta]
>>> lin_m.eval_first_order_form(s)
[0]
[0]

The coefficients of the canoncial form can be accessed:

>>> lin_m.M
1
>>> lin_m.C
0
>>> lin_m.K
g / l
>>> lin_m.A
[0,      1]
[-g / l, 0]
>>> lin.m.evolve(StateEvolution(...))
[0, 1, 2, 3]
[0, 4, 5, 6]

Frequency analysis of a linear system.

>>> lin_m.natural_frequency()
>>> lin_m.damping_ratio()

Time series. Allen uses a TimeSeries object to hold values with a time index. I'm not that fond of indexing Series and DataFrames using floating point time values but maybe some hack to use an integer index for microseconds from epoch would work?

[moorepants]

@ixjlyons Here are some starting ideas. I'll keep playing around with this. My only worry is that we hide too many things that we'd want the students to learn. I think trying to create the notebooks first with "raw" code will show us what might be useful to abstract away, but the above stuff shows some of the things I've thought about for a long time, i.e. how to expose the symbolics and numerics in a useful way.

[ixjlyons]

This looks really cool, but I agree that it's hard to tell if it hides too much. One possibility is to implement the high-level stuff (e.g. a linearize method) but withhold showing it to them until after working through it using lower level methods -- essentially "unlocking" features only after they've learned how they work.

Starting out with pure numpy/scipy/sympy sounds like a plan. I think we'll quickly start to see patterns and pain-points that are ok to abstract away.

[moorepants]

Yeah, there are at least two ways:

Student does something the hard way (e.g. linearize an ode on paper) and then show them a tool that does it for them.

OR

Students learn high level concept, e.g. "a linear version of the model only predicts the motion at small deviations from a nominal state", give them the model with the lmitations in mind, let them sovle some problems with it, and then maybe they will ask "how do you make a linear model?".

[moorepants]

Closing this as we have a simpler object model. Just looked again here and there are some cool ideas, but mostly out of the scope of pass 1 for resonance.