4.3 The approaches
The overall goal with this report is to investigate how different formulation of the numerical solution can affect the solution, especially with long-time integration.
Different formulation is meant in a broader meaning, the physical problem is still the same, but the numerical method is changed.
The cases for evaluation are:
1. Plain ODE system.
2. Formulation as DAE problem by inclusion of invariant.
3. Postprocessing of solution to ODE system (4.2.1) by use of invariant (4.2.3).
4. Reformulation of model as Hamiltonian system, applied with symplectic integrator.
The numerical methods used to integrate the solution will be of an explicit and implicit kind with fixed stepsizes, the Improved Euler method and the implicit midpoint rule (IMR).
These methods are of order two, and are used to better recognize the effects of the different approaches. Each of the initiatives is applied to calculations with both the explicit and the implicit method.
One could argue that such low order methods should not be used for long integration of this dynamical system (and especially not with fixed stepsize), but in test with use of the Matlabsolverode45, the numerical method was simply too good to see the effects of the initiatives.
The Butcher tableau of the two methods are presented in figure 4.2, since these methods are computed as such (i.e. as RK methods).
0 0
1 1 0
12 1 2 (a) Improved Euler
12 1 2
1
(b) Implicit Midpoint Rule
Figure 4.2: Used numerical schemes
The implicit midpoint rule is a symplectic method (since it is a one stage Gauss-Legendre method) and hence useful to evaluate the reformulation of the system (4.2.2) into a Hamil-tonian system, see section 4.3.4.
4.3.1 Plain ODE system
In this case, nothing is done to the system (4.2.2). Instead the circumstances under which the tests are performed are elaborated.
Evaluation of the numerical solution somehow requires a periodic solution which doesn’t pass too close to the singular points (x, y) = (1−µ,0) or (−µ,0), since a fixed stepsize is used in this report. A more or less classic periodic solution is used, which has the initial conditions:
x(0) = 1.2 x(0) = 0
y(0) = 0 y(0) =−1.04935750983031990726 along with the ratioµand period T
µ= 1/82.45 T = 6.19216933131963970674 . The solution to this IVP can be seen in figure 4.3.
The error of the numerical solution is chosen to be the drift in position in one dimension i.e. when the periodic solution crosses the liney = 0, the error is as
error(k) =|x(0)−xn(kT)| (4.3.1) where xn(kT) is the solution found (interpolated) at y = 0 for the k’th period. The numerical solution isn’t always near the axis at a multiplum of the period, so its the period number that counts.
This error is used to evaluate the numerical method over ’long-time’ integration (10 times a period).
−1.5 −1 −0.5 0 0.5 1 1.5
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
x y
Figure 4.3: Orbit 4.3.2 DAE formulation
A way of including the invariant-’constraint’ into the ODE system, turning it into a DAE system, is by use of Lagrange multipliers. The multiplierλis introduced in (4.2.2) as
y=f(y)−G(y)λ
0 =g(y) (4.3.2)
where g(y) =I(y)−I(y0) and the function G(y) is chosen to be (∂g/∂y)T, which gives an orthogonal projection of the solution onto the constraint manifoldM.
Performing index reduction on the constraint would give λin terms ofy λ= gyf
gyG
giving an ind1 DAE system, for which an appropriate solver could be applied, for ex-ample theMatlabsolver ode15sorode23t.
Unfortunately, a property of the definition of a first integral to the system (4.2.2), is
that ∂I
∂y f = 0
and since ∂I/∂y = gy, the introduction of the multiplier into the ODE system gives nothing.
4.3.3 Postprocessing
By postprocessing we will mean taking a time-step in the solution of the ODE system with some method, and then ’pulling back’ or projecting this estimate onto the constraint manifold. Different ways of doing this are present in the litterature, an approach in [29]
is: having solution value yn−1, advancing one step with an ODE method φh giving an estimate ˜y and then solving
y˜−yn2 = min
yn
g(yn) = 0
which constitutes a nonlinear constrained least squares problem.
This approach was implemented by use of the Matlab built-in function fmincon which
”finds a constrained minimum of a function of several variables” (from fmincon help) by the setup
fmincon(y˜−yn,y˜,[ ],[ ],[ ],[ ],[ ],[ ], g(yn)) whereg(y) is as in (4.3.2).
Note: The use of this function is quite evaluation expensive, implementation of an efficient algorithm set up for this specific problem could probably make this solution form more at-tractive.
Another way of doing postprocessing is to introduce a stabilizing term in the system such that the ODE system becomes
y =f(y)−γF(y) (I(y)−I0) (4.3.3) where the termF is
F =ET
EET−1
, E(y) = ∂I
∂y such that for the system (4.2.2) the stabilization functionF becomes
F(y) =
Note that this stabilization term vanishes on the constraint manifold.
This stabilization is applied after each advancement by the ODE methodφh, such that the scheme becomes
y˜n+1 = φh(yn)
yn+1 = y˜n+1−γF(˜yn+1)(I(˜yn+1)−I0)
The constant γ = 1 is chosen, which should be an optimal value, according to [28] which also holds more information on this stabilization technique.
Since the stabilization term is a vector function, its inclusion doesn’t increase the compu-tational load significantly.
4.3.4 Hamiltonian reformulation
Reformulating the system (4.2.2) into a Hamiltonian form, and applying a symplectic solver to this, could improve the solution with regard to drift in position, since the symplectic method tends to mimic the flow of the Hamiltonian1 system.
A Hamiltonian system is defined as q = ∂H
∂p , p=−∂H
∂q or equivalently, withy=
q,pT
y=J∇H , J =
0 I
−I 0
whereH is the Hamiltonian (scalar) function.
The state variables in the Hamiltonian form of the ODE system (4.2.2) are q=
so that the system is given by q =
with the Hamiltonian function H(q,p) = 12
p21+p22
+q2p1−q1p2− β R1 − µ
R2 , (4.3.6)
withRi(q) as in (4.2.1).
The Hamiltonian system is very close to the Newtonian considered before, the invariantH is actually the same as the Jacobi IntegralI, and the ODE system (4.3.5) differs only with the factor 2 on the derivatives x, y in p. So it can be expected that these two systems behaves much alike.