Foucault摆模拟

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

我正在尝试在巴黎编写原始的福柯摆锤。我的代码似乎可以正常工作，但是到目前为止，我只能得到下面正确的图形“花”的演变。因此，我不断更改参数以获得左图，但无法这样做。我采用了安装在巴黎的福柯摆锤的参数，L = 67，地球角速度= 2 * pi / 86400，纬度为49/180 * pi。我的初始条件如代码中所写。我尝试了一系列参数来改变我的所有初始条件，纬度和角速度，但无法获得所需的结果。

我使用了如下的Foucault微分方程：我用有限差分方法（比Runge-Kutta更简单）对它们进行编码，方法是用其中心有限差分替换二阶导数。一阶是后向有限差分。到那时，我通过隔离两个方程中的x（i + 1）和y（i + 1）来建立循环。

我的代码对诸如h（=推导步长），地球角速度和纬度（正常）之类的参数非常敏感。我试图将较大范围的参数从较大的h步更改为较小的步，更改为最小和高纬度的初始条件...等等，但是我始终无法获得我宁愿需要的左侧图形。

怎样才能得到剩下的一个？

Answer 1

通过将地球的自转速度提高了120倍，并使模拟运行了32次摆，我就获得了两个图表。另外，我注意到Euler集成会增加系统性能，从而导致不良结果，因此我恢复为标准RK4实现。

这是我用来解决此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

生成的文件results.csv对我来说看起来像这样（用于检查）
              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
   ...