As can be seen from the results, the angular momentum, which should be constant, is not constant, and even increases during the transition manoeuvre. I've checked my code up and down, but could not find any mistakes. I'm guessing it has something to do with the ode45 function, but I'm not sure. Does anyone here have a clue about what's going on?

$\begingroup$2‰ drift over the simulation, with initial conditions given to four significant figures, probably down to rounding errors.$\endgroup$
– JCRMApr 16 at 9:58

2

$\begingroup$I would say regardless of the initial conditions, the angular momentum should remain constant. So then there must be something in the ode45 algorithm that causes this drift.$\endgroup$
– woeterbApr 16 at 10:06

$\begingroup$"ode45 is a versatile ODE solver and is the first solver you should try for most problems. However, if the problem is stiff or requires high accuracy, then there are other ODE solvers that might be better suited to the problem."$\endgroup$
– JCRMApr 16 at 13:44

2 Answers
2

I rewrote your matlab into Python (which is really easy to do and that's not a coincidence!) and used the default SciPy ODE integrator. Checking info['mused'] I see that it always used the Adams (non-stiff) algorithm, so stiffness isn't the problem.

I would guess that you should start reading your Matlab documentation for ode45 and set your error tolerance lower. I am guessing that ode45 is what I would call RK45 and is a variable step size 4/5 order Runge-Kutta. It should be fine but if you don't specify a tolerance, it will choose one for you.

It's possible that between the time your reference was published and when you ran your example, Matlab changed some defaults or the author may have used some and forgotten to tell you, or you didn't notice them mentioned earlier in what you've been reading.

Long story short, it works for me just fine.

The initial angular momentum is not 3000 using the initialization values, but is actually 3000.0045which is where the simulaion starts. It drifts to 3000.001 after 1000 seconds. That's a change of 0.003. Deviation from an even 3000 (or 3e3) is shown in the plot.

I can improve that a little bit, to a total drift of only about 0.001 by inserting rtol=1E-10 or rtol=1E-11 in the call to ODEint.

$\begingroup$Thanks @uhoh, this confirms what I was suspecting. Good to know there's nothing wrong with the code. All in for Python, when I'm not working with Simulink :)$\endgroup$
– woeterbApr 16 at 16:44

$\begingroup$@woeterb yep! If you find a way to coax ode45 to give you what you need, then post another answer here. If it really solves your problem you can also accept it. The goal is to provide the best information possible to future users. Have fun!$\endgroup$
– uhohApr 16 at 16:49

1

$\begingroup$I posted a new answer, which provides the actual solution for the MATLAB example. This indeed really solves the question. Thanks for your efforts.$\endgroup$
– woeterbApr 16 at 17:30

When adjusting the relative and absolute error tolerances, better results are attained. For example, the options variable in the MATLAB code can be defined as follows:

options = odeset('RelTol',1e-10,'AbsTol',1e-10);

Also, to let the initial angular momentum be 3000, instead of 3000.0045 (as was pointed out in this answer), the equation for calculating the angular momentum can be solved for the initial minor axis spin rate, as follows:

$\begingroup$Very nice answer and great news! I would stick to specifying rtol only. Absolute tolerance or atol is awkward to use when you have numbers with physical units. rtol is all you need. Also, see if you can find a way to let the vertical axis of angular momentum autoscale so that you can always see how big the deviation is. By fixing it to [2999 to 3001] you can't see what's going on.$\endgroup$
– uhohApr 16 at 17:39

$\begingroup$@uhoh Turns out when I specify RelTol only, accuracy is only down to 10^(-2). Specifying AbsTol as well gives 10^(-6) accuracy. You're right about the scaling to show the deviation, but the goal was to reproduce the figure in the article. Therefore, fixing the axis suffices.$\endgroup$
– woeterbApr 16 at 17:51

1

$\begingroup$Try an examination of the specific values you use. The effect of the two will be different because one is normalized and the other isn't. It's not enough to say you used or didn't use the parameter, you need to optimize the value of the parameter. They will both have an impact, and if you just set them equal one may have an impact at a different value than the other, but there's no reason to set them to the same value. Good luck, over and out!$\endgroup$
– uhohApr 16 at 22:27

1

$\begingroup$Oh, except feel free to accept your own answer! It best describes the solution to your question.$\endgroup$
– uhohApr 16 at 22:44