A little question

[me7pako]

Hello there, I have a little question :
could you share with us the theoretical part, I mean how you derived the equations and considered the Euler method.

I am interesting of studying its theory to try and implement it myself.

Thanks again

[mapellidario]

Hello!

TL;DR

As i mention in the README, I wrote this code years ago, before attending classes about computational physics, approximations of ODEs and PDEs. I had a quick look at the code and I would do things a bit differently today. I would use the runge-kutta method, I would be more explicit when dropping constants, I would clean up the code, it looks awful. I would also not use https://root.cern to plot the result, I would use pyplot and interactive figures in a jupyter notebook.

Some more info.

From the code, I could dig out what my choices at the time were. The main loop over time calls a single function that updates the status of the pendulum, here

inverted-pendulum/pend.cpp

Line 99 in 0e73456

asta.GravStabVert(dt,i,A,w,theta,perno) ;

I made multiple such functions, describing different models.

• inverted-pendulum/bar.cc

  Line 162 in 0e73456

  void bar::Gravity (double deltaStep)

  describes a simple pendulum, integrating over time with finite difference model [3]

• inverted-pendulum/bar.cc

  Line 194 in 0e73456

  void bar::GravFriction (double deltaStep, double k)

  describes a simple pendulum with friction, using this model [1], integrating over time with euler method [2]

• inverted-pendulum/bar.cc

  Line 206 in 0e73456

  void bar::GravStabVert (double deltaStep, int i, double A, double w, double phi, vect2d perno )

  modelling the stabilization of an inverted pendulum, with simple euler method [2]

• inverted-pendulum/bar.cc

  Line 223 in 0e73456

  void bar::GravStabVert1 (double deltaStep, int i, double A, double w, double phi, vect2d perno )

  modelling the stabilization of an inverted pendulum, with finite difference method [3]

As you see, I used Euler method or finite difference method pretty much interchangeably without even mentioning it in a comment. today I would be a bit more explicit, and I would use a runge-kutta method.

In the stabilized inverted pendulum, there are a couple little assumptions that are not obvious from the code, nor described in any comment:

• I assume that the pivot of the pendulum moves on a vertical line with an oscillating sinusoidal motion, here

  inverted-pendulum/bar.cc

  Line 209 in 0e73456

  double arg = w*i*deltaStep + phi ;

  I compute the angle of this sinusoid, here

  inverted-pendulum/bar.cc

  Line 211 in 0e73456

  vect2d traslazione (0., A*sin(arg));

  I move the pivot to the correct place. At every step, i place the pivot at (0,0), then I translate it to where it should be.
• Do not ask me why in the main, I call this function passing the theta variable as phi argument. phi supposed to be a phase in the sinusoidal motion of the pivot, theta is the starting angle of the pendulum. I consider the angle 0 to be pointing downwards, to a theta angle of 3.14 radians means that the pendulum is pointing upwards. phi and theta should not mix, bad dario of the past. However, changing the initial phi should not have any qualitative effect on the simulation.
• When computing the effect of the moving pivot on to the pendulum, I sum the pivot acceleration to gravity

  inverted-pendulum/bar.cc

  Line 214 in 0e73456

  double accTheta = A*w*w*sin(arg)*sin(theta_p) - 9.806*sin(theta_p) ;

  . Note that if the pivot at time i*deltaStep is at A*sin(w * i*deltaStep), its second derivative is -A*w*w*sin(w * i*deltaStep) [4]. Do not ask me why I did not put the minus sign there, bad dario of the past. In any case, it would just mean that the effect of the pendulum is not in phase with the position of the pivot, this should have not have any qualitative efffect on the simulation. Moreover, here i dropped the mass of the pendulum, the length of the pendulum and where along the length of the pendulum its barycenter resides. I needed to manually tune A (amplitude of the sinusoidal motion of the pivot) and w (angular velocity of the sinusoidal motion of the pivot).

I hope this helps a bit!

[1] https://en.wikipedia.org/wiki/Drag_(physics)#Very_low_Reynolds_numbers:_Stokes'_drag
[2] https://en.wikipedia.org/wiki/Euler_method
[3] https://www.physics.udel.edu/~bnikolic/teaching/phys660/numerical_ode/ode.html
[4] https://www.wolframalpha.com/input?i=second+derivative+of+sin%28w*x%29