傅科摆模拟
Foucault Pendulum simulation
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
我正在尝试编写巴黎原始傅科摆的代码。我的代码似乎可以正常工作，但到目前为止，我只能得到右下方的图形，"the flower" 进化。因此，我不断地改变我的参数来得到左边的图形，但我做不到。
我采用了安装在巴黎的傅科摆的参数，L=67，angular地球速度=2*pi/86400，纬度为 49/180*pi。
我的初始条件如代码中所写。我尝试了一系列参数来改变我的所有初始条件、我的纬度和 angular 速度，但我无法获得所需的结果。
我使用了 Foucault 微分方程如下：我用有限差分法（比 Runge-Kutta 更简单）对它们进行编码，用其中心有限差分代替二阶导数。而一阶是向后有限差分。到那时，我通过在两个方程中隔离 x(i+1) 和 y(i+1) 来构建我的循环。
我的代码对 h（=推导步骤）、地球 angular 速度和纬度（这是正常的）等参数非常敏感。我试图改变大范围的参数，从大的 h 步到小的 h 步，再到最小和高纬度、初始条件……等等，但我永远无法得到我更需要的左图。
左边的可以做什么？
通过将地球自转加速 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
...
我在一张图表中绘制了第 2 和第 3 列，在第二张图表中绘制了第 4 和第 5 列。
有一件事可能是错误的，这取决于您如何管理不同的步长，以及对现实世界示例的物理观察。通过数组的初始化，您意味着在远离中心的 x 方向上的初始速度约为 0.9/0.01=90 [m/s]。要获得不同步长的兼容结果，您需要调整 x(2) 的计算。然而，在图中，绘图从零速度的点开始。您可以通过设置 x(2)=x(1)=1 实现第一顺序。由于使用的积分方法也是一阶的，这样就够了。
对于第二点，注意可以使用复数坐标z=x+iy将系统写成
z'' = -w^2*z - 2*i*E*z', E = Omega*sin(theta)
这是一个常系数线性ODE，它的解是
z(t) = exp(-i*E*t) * (A*cos(w1*t)+B*sin(w1*t)), w1 = sqrt(w^2+E^2)
这描述了频率为 w1 的钟摆运动，其平面以频率 E 顺时针旋转。大自转有T=2*pi/E个周期，期间会出现w1*T/(2*pi)=w1/E个钟摆摆动。
现在插入您的数字，w=sqrt(g/L)=0.383 和 E=2*pi*sin(49°)/86400=5.49e-05，因此基本上 w1=w。每次完整旋转的钟摆周期数为 w/E=6972，因此您可以在图中看到一个密集填充的圆圈。或者，如果只绘制了几个循环，则为非常窄的双楔形。由于每个循环需要 2*pi/w=16.4 [s]，并且积分达到步长 0.01 的 1000 步，在原图中，您可以预期向前摆动和向后摆动的一部分。
为了更真实，将初始速度设置为零，也就是说，钟摆被带到它的起始位置然后松手。同时将时间增加到 30 [s] 以在图中有一个以上的钟摆周期。
从这里我们可以看出解是收敛的，想象一下，它们是线性收敛的。
要获得像引用图片中那样的情节，需要 w/E 的一小部分，算上摆动，它必须在 15 左右。请注意，您无法在地球上的任何地方使用真实比例的钟摆获得此比率。因此设置 w=pi、E=pi/16 并使用一阶方法对 15 时间单位进行积分。
即使对于钟摆周期中 40 点的最小步长，这也会非常快地恶化。
为了获得更好的结果，在一阶导数近似中使用中心差分将局部截断阶数提高到下一个更高的阶数。
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)
复数除法也可以很容易地在轨迹的实分量中进行，
! 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
现在为了尊重增加的顺序，初始点也需要有更多的精度，再次使用零初始速度，这给出了二阶泰勒展开
z(2) = z(1) - 0.5*w^2*z(1)*dt^2
在一阶方法中给出偏离和结构恶化结果的所有步长现在在二阶方法中给出视觉上相同、结构稳定的结果。
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
我正在尝试编写巴黎原始傅科摆的代码。我的代码似乎可以正常工作，但到目前为止，我只能得到右下方的图形，"the flower" 进化。因此，我不断地改变我的参数来得到左边的图形，但我做不到。 我采用了安装在巴黎的傅科摆的参数，L=67，angular地球速度=2*pi/86400，纬度为 49/180*pi。 我的初始条件如代码中所写。我尝试了一系列参数来改变我的所有初始条件、我的纬度和 angular 速度，但我无法获得所需的结果。
我使用了 Foucault 微分方程如下：我用有限差分法（比 Runge-Kutta 更简单）对它们进行编码，用其中心有限差分代替二阶导数。而一阶是向后有限差分。到那时，我通过在两个方程中隔离 x(i+1) 和 y(i+1) 来构建我的循环。
我的代码对 h（=推导步骤）、地球 angular 速度和纬度（这是正常的）等参数非常敏感。我试图改变大范围的参数，从大的 h 步到小的 h 步，再到最小和高纬度、初始条件……等等，但我永远无法得到我更需要的左图。
左边的可以做什么？
通过将地球自转加速 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
...
我在一张图表中绘制了第 2 和第 3 列，在第二张图表中绘制了第 4 和第 5 列。
有一件事可能是错误的，这取决于您如何管理不同的步长，以及对现实世界示例的物理观察。通过数组的初始化，您意味着在远离中心的 x 方向上的初始速度约为 0.9/0.01=90 [m/s]。要获得不同步长的兼容结果，您需要调整 x(2) 的计算。然而，在图中，绘图从零速度的点开始。您可以通过设置 x(2)=x(1)=1 实现第一顺序。由于使用的积分方法也是一阶的，这样就够了。
对于第二点，注意可以使用复数坐标z=x+iy将系统写成
z'' = -w^2*z - 2*i*E*z', E = Omega*sin(theta)
这是一个常系数线性ODE，它的解是
z(t) = exp(-i*E*t) * (A*cos(w1*t)+B*sin(w1*t)), w1 = sqrt(w^2+E^2)
这描述了频率为 w1 的钟摆运动，其平面以频率 E 顺时针旋转。大自转有T=2*pi/E个周期，期间会出现w1*T/(2*pi)=w1/E个钟摆摆动。
现在插入您的数字，w=sqrt(g/L)=0.383 和 E=2*pi*sin(49°)/86400=5.49e-05，因此基本上 w1=w。每次完整旋转的钟摆周期数为 w/E=6972，因此您可以在图中看到一个密集填充的圆圈。或者，如果只绘制了几个循环，则为非常窄的双楔形。由于每个循环需要 2*pi/w=16.4 [s]，并且积分达到步长 0.01 的 1000 步，在原图中，您可以预期向前摆动和向后摆动的一部分。
为了更真实，将初始速度设置为零，也就是说，钟摆被带到它的起始位置然后松手。同时将时间增加到 30 [s] 以在图中有一个以上的钟摆周期。
从这里我们可以看出解是收敛的，想象一下，它们是线性收敛的。
要获得像引用图片中那样的情节，需要 w/E 的一小部分，算上摆动，它必须在 15 左右。请注意，您无法在地球上的任何地方使用真实比例的钟摆获得此比率。因此设置 w=pi、E=pi/16 并使用一阶方法对 15 时间单位进行积分。
即使对于钟摆周期中 40 点的最小步长，这也会非常快地恶化。
为了获得更好的结果，在一阶导数近似中使用中心差分将局部截断阶数提高到下一个更高的阶数。
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)
复数除法也可以很容易地在轨迹的实分量中进行，
! 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
现在为了尊重增加的顺序，初始点也需要有更多的精度，再次使用零初始速度，这给出了二阶泰勒展开
z(2) = z(1) - 0.5*w^2*z(1)*dt^2
在一阶方法中给出偏离和结构恶化结果的所有步长现在在二阶方法中给出视觉上相同、结构稳定的结果。