3.6.3. Easily selecting time stepping methods for first order time derivatives
The easiest method to switch between all implemented time stepping schemes is the usage of the function time_scheme(). This function can be applied to any mathematical expression and will expand the expression to the time scheme selected via the string argument scheme. This is done in the following with the anharmonic oscillator
So let us implement an equation class AnharmonicOscillator for this as follows
# Anharmonic oscillator by first order system with different time stepping schemes
class AnharmonicOscillator(ODEEquations):
def __init__(self, scheme):
super(AnharmonicOscillator, self).__init__()
self.scheme = scheme
def define_fields(self):
self.define_ode_variable("y")
self.define_ode_variable("dot_y")
def define_residuals(self):
y = var("y")
dot_y = var("dot_y")
residual = (partial_t(dot_y) + y ** 3) * testfunction(dot_y)
residual += (partial_t(y) - dot_y) * testfunction(y)
# Here, we evaluate the chosen scheme just by applying time_scheme(scheme,...)
self.add_residual(time_scheme(self.scheme, residual))
The constructor takes and argument scheme, which is - as usual - stored in a member variable. At the end of define_residuals(), just before adding the residuals to the problem, time_scheme() with the passed scheme is applied on the residual. The problem class and the corresponding code to run all simulations in different output directories reads
class AnharmonicOscillatorProblem(Problem):
def __init__(self, scheme): # Passing scheme here
super(AnharmonicOscillatorProblem, self).__init__()
self.scheme = scheme
def define_problem(self):
eqs = AnharmonicOscillator(scheme=self.scheme)
t = var("time") # Time variable
eqs += InitialCondition(y=1,dot_y=0)
# Calculate the total energy. We also use time_scheme here, e.g. the energy is evaluated by the same scheme as the time stepping
y = var("y")
total_energy = time_scheme(self.scheme, 1/2 * partial_t(y) ** 2 + 1/4 * y ** 4)
eqs += ODEObservables(Etot=total_energy) # Add the total energy as observable
eqs += ODEFileOutput()
self.add_equations(eqs @ "anharmonic_oscillator")
if __name__ == "__main__":
for scheme in {"BDF1", "BDF2", "Newmark2", "MPT", "TPZ", "Simpson", "Boole"}:
with AnharmonicOscillatorProblem(scheme) as problem:
problem.set_output_directory("osci_timestepping_scheme_" + scheme)
problem.run(endtime=100, numouts=200)
Before comparing all time stepping schemes here, let us briefly discuss what the function time_scheme() actually does. Let us consider any generic expression \(G(\partial_t \vec{U},\vec{U},t)\), or in Python some expression G which may contain partial_t(var("U")), var("U") and var("time"). To understand what the application of time_scheme() on G, i.e. time_scheme(scheme,G), actually does, let us first introduce the shorthand notation
i.e. the evaluation of the expression G at the \(k^\text{th}\) history time step, with \(k=0\) meaning the step we are currently solving for and \(k=1\) meaning the values at the last successfully taken step. For a fractional \(k\), the arguments are interpolated linearly, e.g. for \(k=\frac{1}{4}\) we get
The function time_scheme() does two things: depending on the selected scheme, it replaces all occurrences of \(\partial_t\vec{U}\), i.e. all partial_t() calls, by an approximation (cf. (3.6)) suitable for the particular scheme. Then, the expression G is expanded by a linear combination of the current and previous values of \(\vec{U}\) and \(t\), i.e.
where \(g_i\) are the weights of the contributions and \(k_i\) are the corresponding history offsets, which might be also fractional. In the latter case, (3.8) is used. The possible time stepping methods with their approximation of \(\partial_t\vec{U}\) and the used pairs \((g_i,k_i)\) are listed in Table 3.1.
scheme |
\(\partial_t \vec{U}\) replacement |
\((g_i,k_i)\) |
|---|---|---|
“BDF1” |
\(\frac{1}{\Delta t^{(n)}}\left(\vec{U}^{(n)}-\vec{U}^{(n-1)}\right)\) |
\((1,0)\) |
“BDF2” |
\(\frac{1}{\Delta t^{(n)}}\left(\frac{3}{2}\vec{U}^{(n)}-2\vec{U}^{(n-1)}+\frac{1}{2}\vec{U}^{(n-2)}\right)\) (*) |
\((1,0)\) |
“Newmark2” |
(**) |
\((1,0)\) |
“TPZ” |
cf. “BDF1” |
\((\frac{1}{2},0)\), \((\frac{1}{2},1)\) |
“MPT” |
cf. “BDF1” |
\((1,\frac{1}{2})\) |
“Simpson” |
cf. “BDF1” |
\((\frac{1}{6},0)\), \((\frac{2}{3},\frac{1}{2})\), \((\frac{1}{6},1)\) |
“Boole” |
cf. “BDF1” |
\((\frac{7}{90},0)\), \((\frac{16}{45},\frac{1}{4})\), \((\frac{2}{15},\frac{1}{2})\), \((\frac{16}{45},\frac{3}{4})\), \((\frac{7}{90},1)\) |
Note that the scheme for partial_t(...,2)-terms within G will not be adjusted by the application of time_scheme(). These are always calculated via the Newmark-beta method.
Fig. 3.5 Total energy conservation of an anharmonic oscillator with different time stepping schemes. To visualize the impact a rather larger time step of \(\Delta t=0.5\) was taken.
In Fig. 3.5, the total energy of the simulations of the anharmonic oscillator with the different time steps is depicted. The time stepping is quite coarse, so that the differences are easily visible: First of all "BDF1" fails to conserve the energy dramatically and also "BDF2" is not suitable for this coarse time step. "Newmark2" conserves the energy quite acceptable over long time, however, it has considerable problems the first time steps. The reason is that "Newmark2" (and also "BDF2") require two history values to be set. However, the initial condition specified with InitialCondition in the code is independent of the time, i.e. there is no variable \(t\) (\(=\)var("time")) occurring in the expressions we set for the initial condition. In this case, the first time step of "Newmark2" and "BDF2" is evaluated by "BDF1", which does not require any history values except of the initial values at \(t=0\). Alternatively, one can also can supplement the InitialCondition object with the keyword argument degraded_start=False to fill all history values with the passed values, i.e. with y=1 and dot_y=0 here. In that case, or if the initial condition explicitly depends on the time, the first time step is not degraded to "BDF1". One can best circumvent this problem if the analytical solution is known: In that case, one can simply set the initial condition based on the analytical solution, as it was done in Section 3.6.1.
All other methods do not have these problems: they (as also "BDF1") have no requirements of further history values. Furthermore, by explicitly considering the evaluation of \(\vec{U}\) at the last successful time step, these methods are quite accurate and energy-conserving, given the large time step. However, in particular the method "Boole" is quite expensive since it involves a lot of evaluations at different sub-steps. If one reduces the time step, all methods increase in accuracy, but for e.g. "BDF1" it has to be reduced drastically to yield acceptable results. Moreover, if there is substantial dissipation in the system, i.e. conservation is not required, "BDF2" can already give quite good results. However, it is always best to check the time stepping scheme for your particular problem with an analytical solution if possible. If an analytical solution is not at hand, one should at least test whether a halved and a doubled time step influences the results. If so, one should take a smaller time step and repeat this procedure.
As a final note, other well-established methods, as e.g. the Runge-Kutta method of fourth order, is currently not possible in pyoomph. The reason is that it requires the storage of the results of the sub-stages and multiple solves for each time step. This is not implemented yet in pyoomph.