3.6.2. Trapezoidal rule and the implicit midpoint rule
So far, we have used different time stepping method, but there are obviously more possibilities. Let us focus on systems of first order ODEs only, which can be written in general as
with the vector of unknowns \(\vec{U}\) and the rhs \(\vec{F}(\vec{U},t)\). Of course, as usual in numerical simulations, \(\vec{U}\) is discretized in time, i.e. we have the time steps \(\vec{U}^{(0)}\), \(\vec{U}^{(1)}\), \(\ldots\), \(\vec{U}^{(n-1)}\) already computed and we want to compute the step \(\vec{U}^{(n)}\). With this notation, an approximation of the time derivative can be written as
Here, \(k\) is the order of the time stepper, e.g. \(k=1\) for "BDF1" and \(k=2\) for "BDF2", and \(w_i\) are the time stepper weights. For example, in "BDF1", the weights are \(w_0=-w_1=1/\Delta t^{(n)}\), which gives
where \(\Delta t^{(n)}\) is the time step we want to take to obtain \(\vec{U}^{(n)}\).
Until now, we have not addressed the right hand side, i.e. \(\vec{F}(\vec{U})\). In fact, if we use any variable bound by var() in the code, it always means the value corresponding to the time step \(n\), i.e. the right hand side up to now was always fully implicit. This means e.g. for "BDF1" time stepping that the discretized equations read
where \(t^{(n)}=t^{(n-1)}+\Delta t^{(n)}\) is the time we are currently solving for. Indeed, this is the implicit Euler integration.
However, there are time integration methods, where also previous values of \(\vec{U}\), i.e. \(\vec{U}^{(n-1)}\) etc., are used in the rhs. One particular example is the trapezoidal rule (sometimes also called Crank-Nicolson method). This one reads
In pyoomph, the lhs can easily be obtained by writing partial_t(...,scheme="BDF1"), but how to evaluate the right hand side appropriately? As mentioned before, var() statements are always evaluated at the current time step \(n\) we are currently solving for, so how to access the values of \(\vec{U}\) at previous time steps?
The answer is the function evaluate_in_past(), which takes an expression as argument and evaluate it previous time steps. Without an argument, evaluate_in_past(expr) evaluates expr at the time step \(n-1\). A second optional argument gives the offset, e.g. evaluate_in_past(expr,2) evaluates the expression expr at the time step \(n-2\) and, beyond that, e.g. evaluate_in_past(expr,0.5) is equal to 0.5*expr+0.5*evaluate_in_past(expr,1), i.e. with fractional values as second argument, one can linearly interpolate between two time steps.
This means that the trapezoidal rule (3.7) in residual form in pyoomph can be written as partial_t(...,scheme="BDF1")-evaluate_in_past(...,0.5), where we have cast the equation already in residual form, i.e. all terms are on the lhs and the rhs is zero.
As an example, let us - once again - consider the harmonic oscillator, but this time integrated with the trapezoidal rule ("TPZ"). The code for the equation class now reads
class HarmonicOscillator(ODEEquations):
def __init__(self,*,omega=1):
super(HarmonicOscillator,self).__init__()
self.omega=omega
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,scheme="BDF1")+evaluate_in_past(self.omega**2*y,0.5))*testfunction(dot_y)
residual += (partial_t(y,scheme="BDF1")-evaluate_in_past(dot_y,0.5)) * testfunction(y)
self.add_residual(residual)
It can be seen that indeed (3.7) is recovered by using scheme="BDF1" for the lhs and evaluate_in_past(...,0.5) for the rhs.
When using the trapezoidal rule for advancing in time, one also has to be consistent with the definition of the total energy, i.e. the argument which is passed to the ODEObservables class (as in the previous section) now has to take the time derivative of \(y\) in the kinetic energy via the "BDF1" scheme, whereas the value of \(y\) in the potential energy has to be evaluated at the average between the time steps \(n\) and \(n-1\), i.e. via evaluate_in_past(...,0.5):
#Calculate the total energy. Important to also stick to the convention: BDF1 derivative and evaluate_in_past(...,0.5)
y=var("y")
total_energy=1/2*partial_t(y,scheme="BDF1")**2+1/2*evaluate_in_past(self.omega*y,0.5)**2
eqs+=ODEObservables(Etot=total_energy) # Add the total energy as observable
Of course, this time stepping can be generalized, which is the so-called \(\theta\)-method:
Obviously, \(\theta=0\) coincides with "BDF1" and \(\theta=1/2\) with the trapezoidal rule. To blend between these methods or any \(0\leq \theta <1\), one can just adjust the second argument of evaluate_in_past(...,theta). The case \(\theta=1\) corresponds to the explicit (forward) Euler method, which is highly unstable and inaccurate.
Another method is the midpoint rule ("MPT"), which reads
If \(\vec{F}\) is linear, the midpoint rule and the trapezoidal rule are identical, but for nonlinear \(\vec{F}\) it is not. The midpoint rule can be used well for symplectic integration, i.e. it is expected that the total energy of the harmonic oscillator is quite well conserved with this method. As explained above, the function evaluate_in_past(...,0.5) gives the trapezoidal rule, not the midpoint rule. To obtain the midpoint rule, one has to use the function evaluate_at_midpoint() instead.
Warning
If you use any history values, e.g. by evaluate_in_past() or evaluate_at_midpoint(), the compilation of the generated code is slower, i.e. you might experience some additional waiting time at the beginning. The reason is that one additional routine has to be generated and compiled to solve for steady solutions in that case, where all history values are identical to the current value. More on steady solutions can be found in Section 3.11.2.