The strategy is to start out by implementing the simple plane pendulum and find out how to handle the singularity. Then the rest of the joints can be appended one at a time.
For each joint a new singularity has to be handled by the solver. Hopefully it will be the first two joints that gives the most dificulties. Once the singularities of the double pendulum are handled it should not be a problem writing the code for handling the rest of the singularities. Numerical problems might arise though, as the number of singularities is raised.
The simple pendulum, in reduced form, is given by:
Note that it is not formulated as a system with a mass matrix. The standard ODE-solvers of MatLab can all be applied to systems formulated with a mass matrix, but none of them are applicable to this problem because the mass matrix is singular. Some of them can handle an ODE-system with a singular mass matrix (i.e. a DAE) but it has to be index one and apperently the index of this system is greater than one at the singularity. It is not possible to useMatLabs event detection in handling the singularity since it cannot detect events where the mass matrix is singular. Therefore it is nescesarry to write a method that can handle the singularity.
A first approach is to use a simple method with constant stepsize. This will make the risk of hitting the singularity smaller than when the stepsize is variable, in which case the stepsize will be controlled by the size of the error. When the solution approaches the singularity the error grows and the stepsize is reduced and thus it almost certain that the singularity is hit. With constant stepsize it should be possible to step over the singularity.
We try with the trapezoidal rule:
yn+1=yn+h
2(fn+fn+1) (3.5.1)
where the ODE is given byy=f(y) andfn=f(yn). The method is implicit so we have to use iterations to take a step along the solution curve. It turns out that fixed point iterations are very fast, when the mass is not too close to the singularity. When|x|>0.15 and the stepsize is h = 10−3 it takes only iterations to obtain an accuracy of size 10−6. When
|x|<0.15 the stepsize has to be smaller other wise the fixed point iteration diverge quite fast. In fact to avoid divergence of the fixed point iteration when close to the singularity one has to keep decreasing the stepsize, and this will make it impossible to step over the singularity.
Newton iterations might then be interesting. When solving the nonlinear equation y=φ(y)
using Newtons method the iteration scheme is given by [?]:
yv+1=yv− our caseφis the righthand side of equation (3.5.1) and we are solving foryn+1. Therefore the scheme is:
Thus it is not a problem that the Jacobi matrix is singular because it is subtracted from the unit matrix before inversion.
After stepping over the singularity we use the continuous extension [14] of trapezoidal rule to determine the exact time passing the singularity. In this manner we can avoid drifting.
The solution obtained when stepping over the singularity is probaly not very accurate, and by going back a step and calculating the time of crossing using a continuous version of the trapezoidal rule it should be possible to get a more accurate value of the dynamical variables right after the crossing.
The continuous extension is given by:
y(tn+θh) =yn+h viewing the trapezoidal rule as a two-stage Runge-Kutta method. They are thus given by:
k1 =fn
k2 =fn+1
Note that we have already made the step so therefore k2 is known. We know that the singularity occurs when x = 0 therefore we can obtain the parameter value θ by solving the scalar equation:
where the number 1 in parentheses refers to the first element of the vector. The solution, θ∗ to this equation reveals the time for the occurence of the singularity. By settingyn+1 = y(tn+ (θ∗ +ε)h) we get a starting point after the singularity which is very close to the singularity and this should give us a better solution curve.
A last comment in this section is concerneing the stepsize. Allthough stated in the begin-ning it is not nescesarry to keep stepsize completely constant. When crossing the singularity it is important to have quite small stepsize, but away from the singularity there is no reason for this. As long as there is a lower bound for the stepsize there is no problem in using differentstepsizes in different domains. In this implementation the stepsize is not deter-mined by the error but by the size of the x-coordinate. We operate with three different stepsizes.
Figure 3.3: The trajectories for the pendulum with different initial conditions without at-temps to step over singularity. Thus the stepsize is not changed as the singularity is ap-proached and the fixed point iterartions fail at some point. Before this happens the mass stay on a circular trajectory. The small circle marks the starting point.
The numerical results are not very extensive. The singularity in the simple pendulum has been passed but there has been made no continuous extension. The double pendulum has not been implemented.
There are good news though. When the solution curve is not close to the singularity it