Setup¶
Import basic python and additional godot packages.
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
Import basic godot packages regarding automatic differentiation.
from godot.core import num, tempo
from godot.core import autodif as ad
from godot.core.autodif import bridge as br
Import model packages for propagation.
from godot.model import common, prop
# optionally avoid verbose logging messages
import godot.core.util as util
util.suppressLogger()
Definition of differential equation¶
Here we define a class that implement the right-hand-side for integrating the motion of the pendulum. We parametrize the equation with a constant, omega, which is the gravitational force over the length of the pendulum:
x'' = -omega * sin(x)where x denotes the angular displacement of the pendulum mass from the resting position, and x'' denotes the second derivative of x. This is second order differential equation which we need to write as two first order equations as following:
x' = y
y' = -omega * sin(x)The class has to implement the common.VectorTimeEvaluable interface to be able to interact with the propagator.
class pendulumRHS(common.VectorTimeEvaluable):
def __init__(self, vec, omega):
# We have to call C++ trampoline class constructor explicitly
common.VectorTimeEvaluable.__init__(self)
# Assign member variables
self.__vec = vec
self.__omega = omega
def eval(self, epoch):
# Computes derivative of d/dx sin(w * t) = w * cos(w * t)
x = self.__vec.eval(epoch)
omega = self.__omega.eval(epoch)
return br.concatenate([omega * x[1], -omega * br.sin(x[0])])
Script to setup propagator¶
We need to carefully setup the propagator inputs, initial states and differential equations for each thing(tm) we want to propagate. In this example there is only one vector quantity.
vec = prop.PropagatorVectorTimeEvaluable()
First, we setup the initial conditions, epoch and state. The state is created as a time evaluable object as required by the propagator input.
epoch0 = tempo.Epoch(0.0, tempo.TimeScale.TDB)
x0 = ad.Vector([0.0, 1.0], 'x0')
x0Func = common.ConstantVectorTimeEvaluable(x0)
We setup another parameter for the (linear) frequency of the sine functions, which is defined in 1/s = Hz. We choose a (linear) frequency equivalent to one day.
omega = ad.Scalar(1.0 / 86400.0, 'omega')
omegaFunc = common.ConstantScalarTimeEvaluable(omega)
Finally, we create the RHS object.
rhs = pendulumRHS(vec, omegaFunc)
With these we can create the propagator input that defines the quantity to propagate. We also need to define the relative and absolute propagation error tolerances.
rtol = [1e-9, 1e-9]
atol = [1e-9, 1e-9]
inp = [prop.VectorInput(vec, x0Func, rhs, rtol, atol)]
At this point the propagator can be created using this input. We also define a parameter for the delta in initial epoch.
dt0 = ad.Scalar(0.0, 'dt0')
dt0Func = common.ConstantScalarTimeEvaluable(dt0)
propagator = prop.Propagator('Pendulum', epoch0, dt0Func, inp)
Propagate the differential equation¶
To compute the propagation, we define the final epoch for the propagation. We choose two whole periods for the sine function, and carry out the propagation. We also specify with False flag that we do not wish to compute partial derivatives.
T = num.TwoPi / omega.value()
n = 2
epoch = epoch0 + n * T
propagator.compute(epoch, False)
To retrieve a certain propagated quantity to be evaluated, we use the output from the propagator.
x = vec.eval(epoch)
print("Angular displacement [rad] = ", x[0])
print("Angular velocity [rad/s] = ", x[1])
The obtained value for the angular displacement is not zero, because the RHS is non-linear in x. If the sine function on the RHS had been approximated by x, the final state would indeed be zero since we had propagated for two (linear) periods.
Plotting the solution¶
To get a full picture of the evolution of the angle, we can generate data and plot it for the range of the propagation. We will define a function to generate the solution as we will reuse this later.
def generatePlotData(pro):
grid = pro.range().createGrid(60.0)
t = [e.mjd() for e in grid]
x = np.asarray([vec.eval(e) for e in grid])
return t, x
We plot it using matplotlib.
t, x = generatePlotData(propagator)
plt.plot(t, x[:, 0], label = 'Angular displacement')
plt.plot(t, x[:, 1], label = 'Angular velocity')
plt.xlabel('Time (days)')
plt.legend()
plt.grid()
plt.show()
The plot shows that indeed the period of the solution is larger than 2 pi/omega because of the nonlinear RHS.
Propagate including partial derivatives¶
To compute partial derivatives we simplfy have to repropagate with partials flag set to True. The same final state is computed.
propagator.compute(epoch, True)
xepoch = tempo.XEpoch(epoch)
x = vec.eval(xepoch)
print(x)
To get the partial derivatives of the final state, we have to create and XEpoch and use that for the evaluation.
y = common.resolve(x)
print(y)
The resolve() method computes the partial derivatives w.r.t. the leaf, x0. There are mainly two reasons why common.resolve is needed:
- x0 might depend on a large number of other parameters (leafs). To minimise the dimension of the Parameter-transition matrix right-hand-side, they are not considered during propagation, but need to be resolved to after propagation
- One might also be interested in partials of x w.r.t. x(t1) where t1 is between the initial and final propagation epoch (e.g. for OD). The decision on which epoch the user is interested is delayed until the partials are resolved
Changing parameter values¶
To modify the value of a parameter we can simply use the provided interface of the constant time evaluable object.
x0Func.set([0.5, -1.0])
dt0Func.set(86400)
Then, we can simplfy propagate again the solution which will use these changes.
propagator.compute(epoch, False)
t, x = generatePlotData(propagator)
plt.plot(t, x[:, 0], label = 'Angular displacement')
plt.plot(t, x[:, 1], label = 'Angular velocity')
plt.xlabel('Time (days)')
plt.legend()
plt.grid()
plt.show()
