Simple Newton Pendulum¶

[1]:

from cosapp.base import System
import numpy as np


class Pendulum(System):
    """Point mass planar pendulum model"""
    def setup(self):
        self.add_inward('L', 1.00, desc='Length of the rod')
        self.add_inward('g', 9.81, desc='Gravitational acceleration')

        self.add_outward('acc', 0., desc="Angular acceleration of the mass")
        self.add_transient('omega', der='acc', desc="Angular velocity")
        self.add_transient('theta', der='omega', desc="Positional angle")

    def compute(self):
        self.acc = -(self.g / self.L) * np.sin(self.theta)


class NewtonPendulum(System):
    def setup(self):
        self.add_child(Pendulum('p1'), pulling=['L', 'g'])
        self.add_child(Pendulum('p2'), pulling=['L', 'g'])

        contact = self.add_event('contact', trigger="p1.theta == p2.theta")
        self.add_event('collision', trigger=contact.filter("p1.omega > p2.omega"))

    def transition(self):
        if self.collision.present:
            p1, p2 = self.p1, self.p2
            # swap angular velocities
            p1.omega, p2.omega = p2.omega, p1.omega

[2]:

from cosapp.drivers import RungeKutta
from cosapp.recorders import DataFrameRecorder

pend = NewtonPendulum("pend")

driver = pend.add_driver(
    RungeKutta(order=3, dt=0.01, time_interval=[0, 10])
)
driver.set_scenario(
    init = {
        "p1.theta" : np.radians(-60),
        "p1.omega" : 2,
        "p2.theta" : np.radians(22),
        "p2.omega" : 0.5,
    },
    values = {
        "L" : 0.5,
        "g" : 9.81,
    },
    # stop = pend.collision.filter('p1.theta < 0'),  # try it!
)

# Recorder
driver.add_recorder(
    DataFrameRecorder(includes=[
        'p?.theta', 'p?.omega',
    ]),
    period = 0.05,
)

# Run simulation & retrieve recorded data
pend.run_drivers()
data = driver.recorder.export_data()
data = data.drop(['Section', 'Status', 'Error code'], axis=1)

[3]:

# Plot results
import plotly.graph_objs as go

traces = [
    go.Scatter(
        x = data['time'],
        y = data['p1.theta'],
        mode='lines', name='Pendulum 1', line=dict(color='red'),
    ),
    go.Scatter(
        x = data['time'],
        y = data['p2.theta'],
        mode='lines', name='Pendulum 2', line=dict(color='blue'),
    ),
]
layout = go.Layout(
    title = "Angular positions",
    xaxis = dict(title="Time"),
    height = 450,
    hovermode = "x",
)
go.Figure(data=traces, layout=layout)

Time driver property recorded_events contains a list of all recorded events. Each element of the list is a named tuple EventRecord(time, events), containing the occurrence time, and the list of cascading events at that time, starting by the primitive event.

[4]:

driver.recorded_events

[4]:

[EventRecord(time=0.3129308816366094, events=[<Event pend.contact>, <Event pend.collision>]),
 EventRecord(time=1.0724999600884106, events=[<Event pend.contact>, <Event pend.collision>]),
 EventRecord(time=1.8314905848560308, events=[<Event pend.contact>, <Event pend.collision>]),
 EventRecord(time=2.5902682511923247, events=[<Event pend.contact>, <Event pend.collision>]),
 EventRecord(time=3.349180640271147, events=[<Event pend.contact>, <Event pend.collision>]),
 EventRecord(time=4.108586734528638, events=[<Event pend.contact>, <Event pend.collision>]),
 EventRecord(time=4.868887485626274, events=[<Event pend.contact>, <Event pend.collision>]),
 EventRecord(time=5.6305622357312926, events=[<Event pend.contact>, <Event pend.collision>]),
 EventRecord(time=6.394215419574377, events=[<Event pend.contact>, <Event pend.collision>]),
 EventRecord(time=7.1606375003784, events=[<Event pend.contact>, <Event pend.collision>]),
 EventRecord(time=7.930878893467985, events=[<Event pend.contact>, <Event pend.collision>]),
 EventRecord(time=8.706314427786635, events=[<Event pend.contact>, <Event pend.collision>]),
 EventRecord(time=9.488609014935344, events=[<Event pend.contact>, <Event pend.collision>])]

Time driver property event_data contains the same fields as the recorder, but only at event times:

[5]:

driver.event_data.drop(['Section', 'Status', 'Error code'], axis=1)

[5]:

	Reference	p1.omega	p1.theta	p2.omega	p2.theta	time
0	t=0.31293088163661	4.787848	0.188669	-1.552723	0.188669	0.312931
1	pend.contact	-1.552723	0.188669	4.787848	0.188669	0.312931
2	t=1.0724999600884	1.682617	-0.118677	-4.831507	-0.118677	1.072500
3	pend.contact	-4.831507	-0.118677	1.682617	-0.118677	1.072500
4	t=1.831490584856	4.855791	0.045391	-1.751205	0.045391	1.831491
5	pend.contact	-1.751205	0.045391	4.855791	0.045391	1.831491
6	t=2.5902682511923	1.757959	0.029126	-4.858203	0.029126	2.590268
7	pend.contact	-4.858203	0.029126	1.757959	0.029126	2.590268
8	t=3.3491806402711	4.838499	-0.102850	-1.702831	-0.102850	3.349181
9	pend.contact	-1.702831	-0.102850	4.838499	-0.102850	3.349181
10	t=4.1085867345286	1.586258	0.173733	-4.798684	0.173733	4.108587
11	pend.contact	-4.798684	0.173733	1.586258	0.173733	4.108587
12	t=4.8688874856263	4.743040	-0.239635	-1.409096	-0.239635	4.868887
13	pend.contact	-1.409096	-0.239635	4.743040	-0.239635	4.868887
14	t=5.6305622357313	1.172499	0.298221	-4.678179	0.298221	5.630562
15	pend.contact	-4.678179	0.298221	1.172499	0.298221	5.630562
16	t=6.3942154195744	4.613127	-0.346796	-0.877849	-0.346796	6.394215
17	pend.contact	-0.877849	-0.346796	4.613127	-0.346796	6.394215
18	t=7.1606375003784	0.526950	0.382039	-4.559355	0.382039	7.160638
19	pend.contact	-4.559355	0.382039	0.526950	0.382039	7.160638
20	t=7.930878893468	4.530440	-0.399599	-0.122997	-0.399599	7.930879
21	pend.contact	-0.122997	-0.399599	4.530440	-0.399599	7.930879
22	t=8.7063144277866	-0.326529	0.393557	-4.540494	0.393557	8.706314
23	pend.contact	-4.540494	0.393557	-0.326529	0.393557	8.706314
24	t=9.4886090149353	4.599304	-0.356116	0.802758	-0.356116	9.488609
25	pend.contact	0.802758	-0.356116	4.599304	-0.356116	9.488609