Foucault Pendulum simulation

Question

Program Foucault
  IMPLICIT NONE
  REAL,DIMENSION(:),ALLOCATABLE :: t, x,y
  REAL,PARAMETER :: pi=3.14159265358979323846, g=9.81
  REAL           :: L, vitessea, lat, h, omega, beta
  INTEGER :: i , zeta

  zeta=1000

  Allocate(x(zeta),y(zeta),t(zeta))

  L=67.
  lat=49/180*pi
  omega=sqrt(g/L)
  h=0.01

  Do i= 1,zeta

  IF(i==1 .OR. i==2) THEN
    t(1)=0.0
    t(2)=0.0
    x(1)=0.1
    x(2)=1
    y(1)=0.0
    y(2)=0.0
  ELSE

    t(i+1)=real(i)*h

    x(i+1)=(-omega**2*x(i)+2.0*((y(i)-y(i-1))/h)*latang(lat))*h**2+2.0*x(i)-x(i-1)

    y(i+1)=(-omega**2*y(i)-2.0*((x(i)-x(i-1))/h)*latang(lat))*h**2+2.0*y(i)-y(i-1)

  END IF

    WRITE(40,*) t(i), x(i)
    WRITE(60,*) t(i), y(i)
    WRITE(50,*) x(i), y(i)


END DO

Contains

REAL Function latang(alpha)
REAL, INTENT(IN) :: alpha
REAL :: sol

latang=2*pi*sin(alpha)/86400
END FUNCTION

End Program Foucault

I'm trying to code the original Foucault Pendulum in Paris. My code seems to be working but so far, I could only get the below right graphic, "the flower" evolution. Therefore, I changed my parameters constantly to get the left graphic but I couldn't. I took parameters of Foucault Pendulum installed in Paris with L=67, angular velocity of earth =2*pi/86400 and latitude of 49/180*pi. My initial conditions are as written in the code. I tried a way range of parameters varying all of my initial conditions, my latitude and angular velocity but i couldn't get the left desired results.

I used Foucault differential equations as below : i coded them with Finite difference method (more simple than Runge-Kutta) by replacing the 2nd order derivation by its central finite difference. And the first order one by it's backward finite difference. By then, i build my loop by isolating x(i+1) and y(i+1) in both equations.

My code is very sensitive to parameters such as h (=derivation step), earth angular velocity and latitude (which is normal). I tried to change a way big range of parameters from a big h step to a small one, to a minimal and high latitude, initial conditions...etc but i couldn't ever get the left graphic which i rather need.

What could be made to get the left one ?

Answer 1

There is one thing that may be wrong depending on how you manage different step sizes, and an observation on the physics of the real-world example. With the initialization of the arrays, you imply an initial velocity of about 0.9/0.01=90 [m/s] in x direction away from the center. To get compatible results for different step sizes, you would need to adapt the calculation of x(2). However, in the graphs the plot starts from a point with zero velocity. This you can implement to first order by setting x(2)=x(1)=1. As the used integration method is also first order, this is sufficient.

For the second point, note that one can write the system using complex coordinates z=x+iy as

z'' = -w^2*z - 2*i*E*z',   E = Omega*sin(theta)

This is a linear ODE with constant coefficients, the solution of it is

z(t) = exp(-i*E*t) * (A*cos(w1*t)+B*sin(w1*t)),  w1 = sqrt(w^2+E^2)

This describes a pendulum motion of frequency w1 whose plane rotates with frequency E clockwise. The grand rotation has period T=2*pi/E, during which w1*T/(2*pi)=w1/E pendulum swings occur.

Now insert your numbers, w=sqrt(g/L)=0.383 and E=2*pi*sin(49°)/86400=5.49e-05, so that essentially w1=w. The number of pendulum cycles per full rotation is w/E=6972, so that you can expect a densely filled circle in the plot. Or a very narrow double wedge if only a few cycles are plotted. As each cycle takes 2*pi/w=16.4 [s], and the integration goes 1000 steps of step size 0.01, in the plot as it is you can expect a swing forth and part of the swing back.

To be more realistic, set the initial velocity to zero, that is, the pendulum is taken to its start position and then let go. Also increase the time to 30 [s] to have more than one pendulum cycle in the plot.

It from this we can see that the solutions converge, and with some imagination, that they converge linearly.

To get a plot like in the cited images, one needs a much smaller fraction of w/E, counting the swings, it has to be around 15. Note that you can not get this ratio anywhere on earth with a realistically scaled pendulum. So set w=pi, E=pi/16 and integrate over 15 time units using the first order method.

This detoriorates really fast, even for the smallest step size with 40 points in a pendulum cycle.

For a better result, increase the local truncation order to the next higher by using the central difference in the first derivative approximation.

z(i+1) - 2*z(i) + z(i-1) = -w^2*z(i)*dt^2 - i*E*(z(i+1)-z(i-1))*dt

z(i+1) = ( 2*z(i) - z(i-1) - w^2*z(i)*dt^2 + i*E*z(i-1)*dt ) / (1+i*E*dt)

The division by the complex number can also be easily carried out in the real components of the trajectory,

! x(i+1)-2*x(i)+x(i-1) = h^2*(-omega**2*x(i)) + h*earth*(y(i+1)-y(i-1))
! y(i+1)-2*y(i)+y(i-1) = h^2*(-omega**2*y(i)) - h*earth*(x(i+1)-x(i-1))
    t(i) = t(i-1) + h
    cx = (2-(h*omega)**2)*x(i) - x(i-1) - h*earth*y(i-1)
    cy = (2-(h*omega)**2)*y(i) - y(i-1) + h*earth*x(i-1)
    den = 1+(h*earth)**2
    x(i+1) = (cx + h*earth*cy)/den
    y(i+1) = (cy - h*earth*cx)/den

Now to respect the increased order, also the initial points need to have an order of accuracy more, using again zero initial speed, this gives in the second order Taylor expansion

z(2) = z(1) - 0.5*w^2*z(1)*dt^2

All the step sizes that gave deviating and structurally deteriorating results in the first order method now give a visually identical, structurally stable results in this second order method.

Answer 2

I was able to get the two charts, by speeding up the earth's rotation 120× fold, and allowing the simulation to run for 32 swings of the pendulum. Also, I noticed that Euler integration added energy to the system making for bad results, so I reverted to a standard RK4 implementation.

and here is the code I used to solve this ODE:

program FoucaultOde
implicit none    
integer, parameter :: sp = kind(1.0), dp = kind(1d0)

! Constants
real, parameter :: g=9.80665, pi =3.1415926536

! Variables
real, allocatable :: y(:,:), yp(:), k0(:),k1(:),k2(:),k3(:)
real :: lat, omega, h, L, earth, period
real :: t0,x0,y0,vx0,vy0
integer :: i, zeta, f1, swings

! Code starts here
swings = 32
zeta = 400*swings
L = 67
lat = 49*pi/180      
period = 24*60*60    ! period = 86400
earth = (2*pi*sin(lat)/period)*120      !120 multiplier for roation
omega = sqrt(g/L)

allocate(y(5,zeta))    
allocate(yp(5), k0(5),k1(5),k2(5),k3(5))

! make pendulum complete 'swings' cycles in 'zeta' steps
h = swings*2*pi/(omega*zeta)  

t0 = 0
x0 = 0.5    ! Initial displacement
y0 = 0
vx0 = 0
vy0 = 0

! Initial conditions in the state vector Y
Y(:,1) = [t0,x0,y0,vx0,vy0]

do i=2, zeta

    ! Euler method (single step)
    ! Yp = ode(Y(:,i-1))

    ! Runge-Kutta method (four steps)
    k0 = ode(Y(:,i-1))
    k1 = ode(Y(:,i-1) + h/2*k0)
    k2 = ode(Y(:,i-1) + h/2*k1)
    k3 = ode(Y(:,i-1) + h*k2)            
    Yp = (k0+2*k1+2*k2+k3)/6            

    ! Take a step
    Y(:,i) = Y(:,i-1) + h*Yp        
end do

open( newunit=f1, file='results.csv', status = 'replace', pad='no')

! write header
write (f1, '(a15,a,a15,a,a15,a,a15,a,a15)') 't',',', 'x',',','y',',', 'vx',',','vy'

! write rows of data, comma-separated
do i=1, zeta
    write (f1, '(g,a,g,a,g,a,g,a,g)') y(1,i),',',y(2,i),',',y(3,i),',',y(4,i),',',y(5,i)
end do

close(f1)
contains

function ode(Y) result(Yp)
real, intent(in) :: Y(5)
real :: Yp(5), t,px,py,vx,vy,ax,ay

    ! Read state vector Y to component values
    t  = Y(1)
    px = Y(2)
    py = Y(3)
    vx = Y(4)
    vy = Y(5)

    ! Reference paper: 
    ! http://www.legi.grenoble-inp.fr/people/Achim.Wirth/final_version.pdf

    ax = -(omega**2)*px + 2*vy*earth    ! (equation 53)
    ay = -(omega**2)*py - 2*vx*earth    ! (equation 54)

    ! State vector rate. Note, rate of time is aways 1.0
    Yp = [1.0, vx, vy, ax, ay]

end function        

end program FoucaultOde

The resulting file results.csv looks like this for me (for checking)

              t,              x,              y,             vx,             vy
    .000000    ,   5.000000    ,    .000000    ,    .000000    ,    .000000    
   .4105792E-01,   4.999383    ,   .1112020E-06,  -.3004657E-01,   .8124921E-05
   .8211584E-01,   4.997533    ,   .8895339E-06,  -.6008571E-01,   .3249567E-04
   .1231738    ,   4.994450    ,   .3001796E-05,  -.9011002E-01,   .7310022E-04
   .1642317    ,   4.990134    ,   .7114130E-05,  -.1201121    ,   .1299185E-03
   .2052896    ,   4.984587    ,   .1389169E-04,  -.1500844    ,   .2029225E-03
   .2463475    ,   4.977810    ,   .2399832E-04,  -.1800197    ,   .2920761E-03
   .2874054    ,   4.969805    ,   .3809619E-04,  -.2099106    ,   .3973353E-03
   ...

from which I plotted the 2nd and 3rd columns in one chart, and the 4th and 5th for the second chart.