【C言語】ペンデュラムウェーブの糸の長さを数値計算で求めてみた

概要

　ペンデュラムウェーブの糸の長さを求めた。そのために楕円方程式を数値計算(シンプソン法)で解いた。

動作環境

• Windows10(64bit)

やりたいこと

　ペンデュラムウェーブの糸の長さを計算したい。ペンデュラムウェーブの周期の表式には楕円方程式が含まれる。そこで数値計算で楕円方程式を解いてペンデュラムウェーブの糸の長さを求める。

ペンデュラムウェーブとは

　ペンデュラムウェーブとは周期を規則的にそろえた振り子を複数並べたもののことである。具体的には振り子を時刻0秒に一斉に話した後、振り子たちはバラバラに運動したり、規則的なタイミングで振れ方が揃ったりを繰り返す。

• 実際に制作したペンデュラムウェーブ
  https://twitter.com/i/status/1731532866537030101

プログラム

数値計算ではシンプソン法を用いた。
theta_0で振り出しの角度を指定している。
空気抵抗は無視した。

#include <stdio.h>
#include <math.h>

#被積分関数を定義する
double f(double phi,double theta_0){
return 1.0/(sqrt(1.0-pow(sin(theta_0/2.0),2.0)*pow(sin(phi),2.0)));
}

#シンプソン法で振り子の糸の長さを計算する関数を定義する
double calc_l(double i)
{
const double g = 9.8;　#重力加速度
const double T_c = 60.0; #振り子を振れさせる時間(s)
const double theta_0 = (M_PI/180.0)*50.0; #振り子の振出しの角度
const double phi_start = 0.0;
const double phi_finish = M_PI/2.0;
const double m = 10000.0;　
const double h = (phi_finish - phi_start) / m;　#刻み幅
int j = 0;
double sum = 0.0;

for (j = 0; j < m; j++) {
	double phi = phi_start + j * h;
	sum += (h / 6.0) * (f(phi,theta_0) + 4.0 * f(phi+h/2.0,theta_0) + f(phi + h,theta_0));
}
return (g*pow(T_c,2.0))/(16.0*pow(sum,2.0)*pow(i,2.0));
}

int main(){
int i,l;　#i:60sの間に振り子が往復する回数　l:iに対応した糸の長さ
int i_max = 100;
double array_i[i_max];
double array_l[i_max];
for (i=1;i<=i_max;i++){
    array_i[i-1]=(double)i;
    array_l[i-1]=calc_l((double)i);
}
FILE *fp;　#計算結果をtxtファイルに書き出す
fp = fopen("data_c25_50.txt","w");
for (int l = 0; l <= i_max-1; l++)
{
fprintf(fp,"%f %f\n",array_i[l],array_l[l]);
}
fclose(fp);
return 0;
}