LoginSignup
0
0
@nakatomoyan(友哉 中村)

【Python】【julia】フーコーの振り子の測定データから地球の自転の角速度を求めてみた

Last updated at Posted at 2023-12-04

概要

　2015 年 4 月 4 日のフーコーの振り子の観測結果を用いて、地球の自転の角速度を有効数字二桁で求めた。
　測定開始からある程度の時間が経過したときの振動面の回転角度と、フーコーの振り子の運動方程式を解いて得られた振動面の回転角度を比較して、最もよく測定値を再現する角速度の値を求めることで、地球の自転の角速度を求めた。
　粘性抵抗の場合の空気抵抗を考慮したとき、地球の自転の角速度の値は四桁目まで理科年表2022 の値と一致した。慣性抵抗の場合の空気抵抗を考慮したとき、地球の自転の角速度の値は二桁目まで理科年表 2022 の値と一致した。
　再現性を示すため、求まる地球の自転の角速度の値の有効数字を増やすため、優位に働く空気抵抗の種類を決定するために、求まる値の有効数字がより多い実験手法による、異なる初期条件での、追加の実験が必要である。

動作環境

  • Windows10(64bit)
  • Python 3.7.7
  • julia 1.8.4

実験の目的

フーコーの振り子を用いて、地球の自転の角速度を求める。

実験原理

　例えば、物体の運動を十分な精度で表す微分方程式（ここではフーコーの振り子の運動方程式）があり、初期条件が課されているとする。このとき、微分方程式を解析的な計算、もしくは数値計算で解けば、各時刻における物体の運動の様子を表すパラメータがわかる。
　ここで、その微分方程式の中に、未知数 (ここでは空気抵抗の比例係数と、地球の自転の角速度) が存在する場合を考える。未知数を変化させたとき、最もよく測定値（ここでは、振動面の回転角度）を再現するときの値 (測定値との差が極小値をとるときの値) を求めることで、未知数を決定できる。
　今回は未知数が二つあるため、まず、理科年表 2022 の地球の自転の角速度の値を代入した状態で、空気抵抗の比例係数を求め、その後、求めた空気抵抗の比例係数を代入して、地球の自転の角速度を求め、最後に、その求めた角速度のときに、最初に求めた空気抵抗の値で極小値をとることを確認した。(測定値を最もよく再現する空気抵抗の比例係数と地球の自転の角速度の組み合わせを決定することが目的である。そのため、最もよく測定値を再現する空気抵抗の比例係数と、最もよく測定値を再現する地球の自転の角速度を別々に求めた後に、それらの値を組み合わせても最もよく測定値を再現するか確認する必要がある。)
　空気抵抗について、粘性抵抗と慣性抵抗の二種類の場合を想定し、それぞれの場合で角速度を求めた。北を y 軸の正の方向、東を x 軸の正の方向とした。振り子の振動中心を原点とした。位置の単位はｍとした。振動面の回転角度は振り子の初期位置と原点を結んだ線分を基準にして、反時計回りを正とした。回転角度の単位は rad とした。

実験方法

　北を y 軸の正の方向、東を x 軸の正の方向とした。振り子の振動中心を原点とした。
　位置の単位はｍとした。振動面の回転角度は振り子の初期位置と原点を結んだ線分を基準にして、反時計回りを正とした。
　回転角度の単位は rad とした。

  1. 空気抵抗を考慮したフーコーの振り子の運動方程式をたてた。
    空気抵抗は、速度の 1 乗に比例する場合速度の 2 乗に比例する場合の、2 パターンを考えた。
  2. 初期条件を課し、そして、未知数（空気抵抗の比例係数）を定め、運動方程式を数値計算で解いた。
  3. 得られた解（振動面の回転角度）と 2015 年 4 月 4 日の実験データの差が、最も小さくなる時の未知のパラメータの値を求めて、未知のパラメータを決定した。

プログラムのアルゴリズム

プログラムのアルゴリズムは、以下の通りである。

  1. 測定値を配列として読み込んだ。
  2. 初期条件を課し、理科年表 2022 に記載されている地球の自転の角速度を代入し、運動方程式を解いた。
  3. 測定値と、計算で求めた値を比較した。
  4. 測定値と最も近い値のときの未知数の値を、空気抵抗の比例係数の値として決定した。
  5. 空気抵抗の比例係数を代入した状態で、初期条件を課し、運動方程式を解いた。
  6. 測定値と最も近い値のときの未知数の値を、地球の自転の角速度の値として決定した。
    1.求めた地球の自転の角速度のときに、空気抵抗の比例係数が極小値をとることを確認した。
    以上がプログラムのアルゴリズムである。

解析に用いたフーコーの振り子の運動方程式は、以下の通りである。

  • 粘性抵抗の場合の空気抵抗を考慮したフーコーの振り子の運動方程式
    $m\frac{{d^2}x}{dt^2} = −\frac{mg}{l}x − 2mωsinθ \frac{dy}{dt} − k(\frac{dx}{dt} )$
    $m\frac{{d^2}y}{dt^2} = −\frac{mg}{l}y − 2mωsinθ \frac{dx}{dt} − k(\frac{dy}{dt} )$

  • 慣性抵抗の場合の空気抵抗を考慮したフーコーの振り子の運動方程式
    $m\frac{{d^2}x}{dt^2} = −\frac{mg}{l}x − 2mωsinθ \frac{dy}{dt} − k(\frac{dx}{dt} )^2$
    $m\frac{{d^2}y}{dt^2} = −\frac{mg}{l}y − 2mωsinθ \frac{dx}{dt} − k(\frac{dy}{dt} )^2$
    x:おもりの x 座標 (m) y:おもりの y 座標 (m) t:時刻 (s) m:おもりの質量 (kg)
    g:弘前における重力加速度 ($m/s^2$) l:振り子のワイヤーの長さ (m)
    ω:地球の自転の角速度 (rad/s) θ:弧度法で表した弘前の緯度 (rad)
    k:おもりに働く空気抵抗の比例係数 (無次元)
    解析に使用した測定値は以下の通りである。
    t_measurement_1=0(s)
    x_measurement_1=-0.011(m)
    y_measurement_1=1.464(m)
    x_d_measurement_1=0(m/s)
    y_d_measurement_1=0(m/s)
    t_measurement_2=6300(s)
    x_measurement_2=-0.1382(m)
    y_measurement_2=0.443(m)
    文献を参考にした他の各物理量の値は以下のとおりである。
    おもりの重さ:50(kg)
    弘前における重力加速度:9.8026121($m/s^2$)
    振り子のワイヤーの長さ:45(m)
    弘前の緯度:北緯 40.5°
    　粘性抵抗の場合の空気抵抗を考慮したフーコーの運動方程式は、解析的に解けたので、解析に解いた解を用いて数値計算を行った。(特性方程式の解を求める部分で数値計算を行っている。)
    　慣性抵抗の場合の空気抵抗を考慮したフーコーの運動方程式は、解析的に解けなかったので、ルンゲクッタ法を用いて、数値計算を行った。刻み幅は 1 である。

プログラム

  • 空気抵抗が粘性抵抗の場合の、地球の自転の角速度を求めるプログラム（使用言語 python）
import numpy as np
from mpmath import pi
from mpmath import exp
from mpmath import sin
from mpmath import sqrt
import matplotlib.pyplot as plt
from multiprocessing import Pool
import mpmath
mpmath.mp.dps=100

# 測定値を代入する
t_measurement_1=0
x_measurement_1=-0.011
y_measurement_1=1.464
x_d_measurement_1=0
y_d_measurement_1=0
t_measurement_2=6300
x_measurement_2=-0.1382
y_measurement_2=0.443

# 捜索範囲を設定する
k_first=1
k_last=6000
k_scale=10**(-5)

# 運動方程式の解を、関数として定義する
def f(k):
    t=t_measurement_2-t_measurement_1
    x_0=x_measurement_1
    x_0_d=x_d_measurement_1
    y_0=y_measurement_1
    y_0_d=y_d_measurement_1
    g=9.8026121
    th=40.5*(pi/180)
    o=7292285*(10**(-11))
    m=50/g
    l=45
    a=-g/l
    b=2*o*sin(th)
    c=0
    d=-k/m
    e=0
    n=1/b
    p=-((2*d)/b)
    q=b+((d**2)/b)-((2*a)/b)
    r=((a*d)/b)+((d*a)/b)
    s=(a**2)/b
    ans=np.roots([n,p,q,r,s])
    lamda_1=ans[0]
    lamda_2=ans[1]
    lamda_3=ans[2]
    lamda_4=ans[3]
    v_1=-(((lamda_1)**3)/(a*b))+(((2*d)/(a*b))*((lamda_1)**2))+(((1/b)-(b/a)-((d**2)/(a*b)))*(lamda_1))-(d/b)
    v_2=-(((lamda_2)**3)/(a*b))+(((2*d)/(a*b))*((lamda_2)**2))+(((1/b)-(b/a)-((d**2)/(a*b)))*(lamda_2))-(d/b)
    v_3=-(((lamda_3)**3)/(a*b))+(((2*d)/(a*b))*((lamda_3)**2))+(((1/b)-(b/a)-((d**2)/(a*b)))*(lamda_3))-(d/b)
    v_4=-(((lamda_4)**3)/(a*b))+(((2*d)/(a*b))*((lamda_4)**2))+(((1/b)-(b/a)-((d**2)/(a*b)))*(lamda_4))-(d/b)
    v_5=((d*a*c)/((b**2)*s))-(e/a)-((d*c)/(a*b))
    d_1=((lamda_3)-(lamda_2))/((lamda_2)-(lamda_1))
    d_2=((lamda_4)-(lamda_2))/((lamda_2)-(lamda_1))
    d_3=(1/((lamda_2)-(lamda_1)))*((lamda_2)*(((a*c)/(b*s))+(x_0))-(x_0_d))
    d_4=((lamda_3)-(lamda_1))/((lamda_1)-(lamda_2))
    d_5=((lamda_4)-(lamda_1))/((lamda_1)-(lamda_2))
    d_6=(1/((lamda_1)-(lamda_2)))*((lamda_1)*(((a*c)/(b*s))+(x_0))-(x_0_d))
    d_7=(((lamda_3)-(lamda_2))*(v_3))/(((lamda_2)-(lamda_1))*(v_1))
    d_8=(((lamda_4)-(lamda_2))*(v_4))/(((lamda_2)-(lamda_1))*(v_1))
    d_9=(1/(((lamda_2)-(lamda_1))*(v_1)))*((lamda_2)*(-(v_5)+(y_0))-(y_0_d))
    d_10=(((lamda_3)-(lamda_1))*(v_3))/(((lamda_1)-(lamda_2))*(v_2))
    d_11=(((lamda_4)-(lamda_1))*(v_4))/(((lamda_1)-(lamda_2))*(v_2))
    d_12=(1/(((lamda_1)-(lamda_2))*(v_2)))*((lamda_1)*(-(v_5)+(y_0))-(y_0_d))
    e_1=(d_1)-(d_7)
    e_2=(d_2)-(d_8)
    e_3=(d_9)-(d_3)
    e_4=(d_4)-(d_10)
    e_5=(d_5)-(d_11)
    e_6=(d_12)-(d_6)
    c_3=((e_5)*(e_3)-(e_2)*(e_6))/((e_5)*(e_1)-(e_4)*(e_2))
    c_4=((e_4)*(e_3)-(e_1)*(e_6))/((e_4)*(e_2)-(e_1)*(e_5))
    c_1=(d_1)*(c_3)+(d_2)*(c_4)+(d_3)
    c_2=(d_4)*(c_3)+(d_5)*(c_4)+(d_6)
    x=(c_1)*exp((lamda_1)*t)+(c_2)*exp((lamda_2)*t)+(c_3)*exp((lamda_3)*t)+(c_4)*exp((lamda_4)*t)-((a*c)/(b*s))
    x_calc=x.real
    return x_calc

def h(k):
    t=t_measurement_2-t_measurement_1
    x_0=x_measurement_1
    x_0_d=x_d_measurement_1
    y_0=y_measurement_1
    y_0_d=y_d_measurement_1
    g=9.8026121
    th=40.5*(pi/180)
    o=7292285*(10**(-11))
    m=50/g
    l=45
    a=-g/l
    b=2*o*sin(th)
    c=0
    d=-k/m
    e=0
    n=1/b
    p=-((2*d)/b)
    q=b+((d**2)/b)-((2*a)/b)
    r=((a*d)/b)+((d*a)/b)
    s=(a**2)/b
    ans=np.roots([n,p,q,r,s])
    lamda_1=ans[0]
    lamda_2=ans[1]
    lamda_3=ans[2]
    lamda_4=ans[3]
    v_1=-(((lamda_1)**3)/(a*b))+(((2*d)/(a*b))*((lamda_1)**2))+(((1/b)-(b/a)-((d**2)/(a*b)))*(lamda_1))-(d/b)
    v_2=-(((lamda_2)**3)/(a*b))+(((2*d)/(a*b))*((lamda_2)**2))+(((1/b)-(b/a)-((d**2)/(a*b)))*(lamda_2))-(d/b)
    v_3=-(((lamda_3)**3)/(a*b))+(((2*d)/(a*b))*((lamda_3)**2))+(((1/b)-(b/a)-((d**2)/(a*b)))*(lamda_3))-(d/b)
    v_4=-(((lamda_4)**3)/(a*b))+(((2*d)/(a*b))*((lamda_4)**2))+(((1/b)-(b/a)-((d**2)/(a*b)))*(lamda_4))-(d/b)
    v_5=((d*a*c)/((b**2)*s))-(e/a)-((d*c)/(a*b))
    d_1=((lamda_3)-(lamda_2))/((lamda_2)-(lamda_1))
    d_2=((lamda_4)-(lamda_2))/((lamda_2)-(lamda_1))
    d_3=(1/((lamda_2)-(lamda_1)))*((lamda_2)*(((a*c)/(b*s))+(x_0))-(x_0_d))
    d_4=((lamda_3)-(lamda_1))/((lamda_1)-(lamda_2))
    d_5=((lamda_4)-(lamda_1))/((lamda_1)-(lamda_2))
    d_6=(1/((lamda_1)-(lamda_2)))*((lamda_1)*(((a*c)/(b*s))+(x_0))-(x_0_d))
    d_7=(((lamda_3)-(lamda_2))*(v_3))/(((lamda_2)-(lamda_1))*(v_1))
    d_8=(((lamda_4)-(lamda_2))*(v_4))/(((lamda_2)-(lamda_1))*(v_1))
    d_9=(1/(((lamda_2)-(lamda_1))*(v_1)))*((lamda_2)*(-(v_5)+(y_0))-(y_0_d))
    d_10=(((lamda_3)-(lamda_1))*(v_3))/(((lamda_1)-(lamda_2))*(v_2))
    d_11=(((lamda_4)-(lamda_1))*(v_4))/(((lamda_1)-(lamda_2))*(v_2))
    d_12=(1/(((lamda_1)-(lamda_2))*(v_2)))*((lamda_1)*(-(v_5)+(y_0))-(y_0_d))
    e_1=(d_1)-(d_7)
    e_2=(d_2)-(d_8)
    e_3=(d_9)-(d_3)
    e_4=(d_4)-(d_10)
    e_5=(d_5)-(d_11)
    e_6=(d_12)-(d_6)
    c_3=((e_5)*(e_3)-(e_2)*(e_6))/((e_5)*(e_1)-(e_4)*(e_2))
    c_4=((e_4)*(e_3)-(e_1)*(e_6))/((e_4)*(e_2)-(e_1)*(e_5))
    c_1=(d_1)*(c_3)+(d_2)*(c_4)+(d_3)
    c_2=(d_4)*(c_3)+(d_5)*(c_4)+(d_6)
    y=(v_1)*(c_1)*exp((lamda_1)*t)+(v_2)*(c_2)*exp((lamda_2)*t)+(v_3)*(c_3)*exp((lamda_3)*t)+(v_4)*(c_4)*exp((lamda_4)*t)+(v_5)
    y_calc=y.real
    return y_calc

def tangent_angle(x_1,y_1,x_2,y_2):
    i=x_1*x_2+y_1*y_2
    n=sqrt((x_1**2)+(x_2**2))*sqrt((y_1**2)+(y_2**2))
    theta=i/n
    return theta

def task():
    plt.figure()
    plt.plot(array_k,array_result.T,'blue',label='Calc')
    plt.title('Relation of k and Error')
    plt.xlabel('k')
    plt.ylabel('Error(rad/s)')
    plt.legend(loc=0)
    plt.savefig("Relation of k and Error1.png")

# 各kにおける回転角度を計算して、結果を配列にする

array_k=np.arange(k_first,k_last,1)
array_k=array_k*k_scale
array_calc=np.empty((array_k.size),dtype=float)

if __name__=="__main__":
    pool=Pool(3)
    array_x_calc=np.array(pool.map(f,array_k))
    array_y_calc=np.array(pool.map(h,array_k))
    for i_k in range(0,array_k.size,1):
        x_calc=array_x_calc[i_k]
        y_calc=array_y_calc[i_k]
        theta_calc=tangent_angle(x_measurement_1,y_measurement_1,x_calc,y_calc)
        array_calc[i_k]=theta_calc

# kを決定する
theta_measurement=tangent_angle(x_measurement_1,y_measurement_1,x_measurement_2,y_measurement_2)
array_measurement=np.full((1,array_k.size),(theta_measurement))
array_result=abs(array_measurement-array_calc)
index_min=np.argmin(array_result)
result=array_k[index_min]
print(result)
task()
  • 空気抵抗が慣性抵抗の場合の、空気抵抗の比例係数を求めるプログラム（使用言語 julia）
using Plots
setprecision(100)

# 測定値を代入する
t_initial = 0
x_measurement_0 = -0.011
x_measurement_1 = 1.464
y_measurement_0 = 0
y_measurement_1 = 0
t_final = 6300
x0_measurement_2 = -0.1382
y0_measurement_2 = 0.443

# 捜索範囲を設定する
N = 6300.0
k_initial = 1
k_final = 1000
k_scale = 10^(-3)

# 関数を定義する
array_k = Array(k_initial:1:k_final)
array_k = k_scale*array_k
array_calc = []

function f1(t, x0, x1, y0, y1)
    return x1
end

function f2(t, x0, x1, y0, y1, i_k)
    g = 9.8
    th = 40.5
    l = 45.0
    o = 7292115.0*(10^(-11))
    k = array_k[i_k]
    m = 50.0
    a = -g/l
    b = 2.0*o*sind(th)
    c = -k/m
    return a*x0 - b*y1 + c*(x1^2)
end

function f3(t, x0, x1, y0, y1)
    return y1
end

function f4(t, x0, x1, y0, y1, i_k)
    g = 9.8
    th = 40.5
    l = 45.0
    o = 7292115.0*(10^(-11))
    k = array_k[i_k]
    m = 50.0
    a = -g/l
    b = 2.0*o*sind(th)
    c = -k/m
    return a*y0 + b*x1 + c*(y1^2)
end

function g(t_initial,t_final,N,x0,x1,y0,y1)
    h=(t_final-t_initial)/N
    t_points=Array(t_initial:h:t_final)
    for i_k in range(start=1,stop=Int(length(array_k)))
        x0=x_measurement_0
        x1=x_measurement_1
        y0=y_measurement_0
        y1=y_measurement_0
        for i in range(start=1,stop=Int(N))
            k_11=h*f1(t_points[i],x0,x1,y0,y1)
            k_12=h*f2(t_points[i],x0,x1,y0,y1,i_k)
            k_13=h*f3(t_points[i],x0,x1,y0,y1)
            k_14=h*f4(t_points[i],x0,x1,y0,y1,i_k)
            k_21=h*f1(t_points[i]+(h/2),x0+((k_11)/2),x1+((k_12)/2),y0+((k_13)/2),y1+((k_14)/2))
            k_22=h*f2(t_points[i]+(h/2),x0+((k_11)/2),x1+((k_12)/2),y0+((k_13)/2),y1+((k_14)/2),i_k)
            k_23=h*f3(t_points[i]+(h/2),x0+((k_11)/2),x1+((k_12)/2),y0+((k_13)/2),y1+((k_14)/2))
            k_24=h*f4(t_points[i]+(h/2),x0+((k_11)/2),x1+((k_12)/2),y0+((k_13)/2),y1+((k_14)/2),i_k)
            k_31=h*f1(t_points[i]+(h/2),x0+((k_21)/2),x1+((k_22)/2),y0+((k_23)/2),y1+((k_24)/2))
            k_32=h*f2(t_points[i]+(h/2),x0+((k_21)/2),x1+((k_22)/2),y0+((k_23)/2),y1+((k_24)/2),i_k)
            k_33=h*f3(t_points[i]+(h/2),x0+((k_21)/2),x1+((k_22)/2),y0+((k_23)/2),y1+((k_24)/2))
            k_34=h*f4(t_points[i]+(h/2),x0+((k_21)/2),x1+((k_22)/2),y0+((k_23)/2),y1+((k_24)/2),i_k)
            k_41=h*f1(t_points[i]+h,x0+(k_31),x1+(k_32),y0+(k_33),y1+(k_34))
            k_42=h*f2(t_points[i]+h,x0+(k_31),x1+(k_32),y0+(k_33),y1+(k_34),i_k)
            k_43=h*f3(t_points[i]+h,x0+(k_31),x1+(k_32),y0+(k_33),y1+(k_34))
            k_44=h*f4(t_points[i]+h,x0+(k_31),x1+(k_32),y0+(k_33),y1+(k_34),i_k)
            x0=x0+(1/6)*((k_11)+2*(k_21)+2*(k_31)+(k_41))
            x1=x1+(1/6)*((k_12)+2*(k_22)+2*(k_32)+(k_42))
            y0=y0+(1/6)*((k_13)+2*(k_23)+2*(k_33)+(k_43))
            y1=y1+(1/6)*((k_14)+2*(k_24)+2*(k_34)+(k_44))
        end
        x_calc=x0
        y_calc=y0
        theta_calc=tangent_angle(x_measurement_0,y_measurement_0,x_calc,y_calc)
        push!(array_calc,theta_calc)
        println(i_k)
    end
    theta_measurement=tangent_angle(x_measurement_0,y_measurement_0,x0_measurement_2,y0_measurement_2)
    array_theta_measurement=fill(theta_measurement,(1,Int(length(array_k))))
    array_result=abs.(array_theta_measurement.-array_calc)
    array_result=array_result[:,1]
    index_min=findmin(array_result)[2]
    result=array_k[index_min]
    plot(array_k,array_result)
    savefig("Foucault_pendulum_5.png")
    return result
end

function tangent_angle(x_1,y_1,x_2,y_2)
    i=BigFloat(x_1*x_2+y_1*y_2)
    n=BigFloat(sqrt((x_1^2)+(x_2^2)))*BigFloat(sqrt((y_1^2)+(y_2^2)))
    theta=BigFloat(i/n)
    return theta
end

#角速度を決定する
g(t_initial,t_final,N,x_measurement_0,x_measurement_1,y_measurement_0,y_measurement_1)

#o=7292115.0*(10^(-11))

解析結果

  • 空気抵抗が粘性抵抗の場合
    空気抵抗の比例係数：0.00371（有効数字二桁）
    地球の自転の角速度：$7.292285 × (10 ^{-5}$)rad/s（有効数字二桁）
  • 空気抵抗が慣性抵抗の場合
    空気抵抗の比例係数：0.116（有効数字二桁）
    地球の自転の角速度：$7.236585 × (10 ^{-5}$)rad/s（有効数字二桁）
    (参考)
  • 理科年表 2022 に記載されている地球の自転の角速度の値：$7.292115 × (10^ {-5}$)rad/s（有効数字七桁）

考察

　地球の自転の角速度を有効数字二桁（実験データの記録に単位と、具体的な実験方法が書かれていなかったため、実際の有効数字が何桁かは不明確である。実験ノートを読み込み、時刻と位置の有効数字は二桁であると推測した。）で求めた。おもりの質量の有効数字は二桁、弘前における重力加速度の有効数字は七桁、振り子のワイヤーの長さの有効数字は二桁、測定場所の緯度の有効数字は三桁として考えた。有効数字が二桁になったのは位置、時刻、おもりの質量、振り子のワイヤーの長さの有効数字が二桁であるためである。地球の自転の角速度の有効数字を増やすために、有効数字の考慮が必要な物理量の個数が少なく、かつ、測定値の有効数字が多い実験手法を検討すべきであると考える。
　今回は空気抵抗として、粘性抵抗の場合の空気抵抗と慣性抵抗の場合の空気抵抗両方で解析を行った。粘性抵抗の場合の空気抵抗を考慮したとき、求めた地球の自転の角速度の値は、四桁目まで理科年表 2022 の値と一致した。慣性抵抗の場合の空気抵抗を考慮したとき、求めた地球の自転の角速度の値は、二桁目まで理科年表 2022 の値と一致した。慣性抵抗の場合の空気抵抗を考慮したフーコーの振り子の運動方程式は、非線形方程式である。非線形方程式は、初期条件がわずかに異なるだけで、時間発展の結果が大きく異なる性質を持つことが知られている。複数の異なる初期条件で、実験、解析を行い、フーコーの振り子に優位に働く空気抵抗を決定する必要があると考える。
　今回は、2015 年 4 月 4 日の観測結果一回分だけを用いて、解析を行った。そのため、今回の実験だけでは実験結果に再現性があることを示すことができない。解析結果の再現性を示すために、追加の実験を行う必要があると考える。

結論

　フーコーの振り子を用いて、地球の自転の角速度を有効数字二桁で求めた。
　粘性抵抗の場合の空気抵抗を考慮した場合、求めた地球の自転の角速度は理科年表 2022 の値と四桁目まで一致した。
　慣性抵抗の場合の空気抵抗を考慮した場合、求めた地球の自転の角速度は理科年表 2022 の値と二桁目まで一致した。
　再現性を示すため、求まる地球の自転の角速度の値の有効数字を増やすため、優位に働く空気抵抗の種類を決定するために、求まる値の有効数字がより多い実験手法による、異なる初期条件での、追加の実験が必要である。

参考文献

・力学　御領潤　日本評論社　第一版第一刷
・理科年表 2022 　国立天文台
・弘前大学　フーコー振り子
http://www.st.hirosaki-u.ac.jp/Foucault/
閲覧日 2022/3/5

0
0
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0