Root Finding algorithm demo¶

This example shows how the root finding algorithm works.

In [1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from godot.core import num, tempo, events

In [2]:

freq = num.TwoPi
f = lambda e: np.cos(e.mjd() * freq)
g = lambda e: np.sin(e.mjd() * freq)

In [3]:

ran = tempo.EpochRange(tempo.Epoch("0.0 TDB"), tempo.Epoch("4.0 TDB"))
grid = ran.createGrid(1e-2 * tempo.SecondsInDay)
t = [e.mjd() for e in grid]
x = [f(e) for e in grid]
y = [g(e) for e in grid]

In [4]:

plt.figure(figsize=(10, 5))
plt.plot(t, x, label="f")
plt.plot(t, y, label="g")
plt.grid()
plt.legend()
plt.show()

In [5]:

tol = 1e-9
event_grid = ran.createGrid(0.3 * tempo.SecondsInDay)
ev = events.generateEventSet(f, g, event_grid, tol)
event_t = [e.mjd() for e in event_grid]
event_x = [f(e) for e in event_grid]
roots = [entry.value().mjd() for entry in ev]

In [6]:

plt.figure(figsize=(10, 5))
plt.plot(t, x, label="f")
plt.plot(t, y, label="g")
for l, u in zip(event_t, event_x):
    plt.plot([l, l], [0, u], '--k.', linewidth=1.0)
plt.plot(roots, [0] * len(roots), 'o', label="roots")
plt.grid()
plt.legend()
plt.show()

In [7]:

f2 = lambda e: f(e) - 0.95
g2 = lambda e: 1.0
x2 = [f2(e) for e in grid]
y2 = [g2(e) for e in grid]
ev = events.generateEventSet(f2, g2, event_grid, tol)
event_x2 = [f2(e) for e in event_grid]
roots = [entry.value().mjd() for entry in ev]

In [8]:

plt.figure(figsize=(10, 5))
plt.plot(t, x2, label="f")
plt.plot(t, y2, label="g")
for l, u in zip(event_t, event_x2):
    plt.plot([l, l], [0, u], '--k.', linewidth=1.0)
plt.plot(roots, [0] * len(roots), 'o', label="roots")
plt.grid()
plt.legend()
plt.show()

In [9]:

ev = events.generateEventSet(f2, g, event_grid, tol)
roots = [entry.value().mjd() for entry in ev]

In [10]:

plt.figure(figsize=(10, 5))
plt.plot(t, x2, label="f")
plt.plot(t, y, label="g")
for l, u in zip(event_t, event_x2):
    plt.plot([l, l], [0, u], '--k.', linewidth=1.0)
plt.plot(roots, [0] * len(roots), 'o', label="roots")
plt.grid()
plt.legend()
plt.show()