C++で単振り子の数値シミュレーション

はじめに

記述言語はC++となっていますが、C++らしい書き方はしていないため、ほぼC言語となっています。よってC++自体の勉強をしたい方などにはあまり良い内容とは言えませんので、あらかじめご了承ください。C++らしい書き方はこちらに書き直してくださった方がいますので、比較してみるとよろしいかと思います。

物理シミュレーション自体について投稿したものがありますので、こちらから読んでいただいたほうが良いかもしれません。
https://qiita.com/ysuzuki19/items/ec98c5667ecd5fafa8b1

こちらではオイラー法という手法でプログラムを作成してみました。

コメントで頂いたのですが、細かいところで厳密ではないとのことですので、「取り敢えず何か動かしたい」という初心者の方は試してもらうと良いかと思いますが、「ガッツリ数値シミュレーションを勉強したい」という方には向いていないと思います。

実行環境

ubuntu18.04LTS

シミュレーションする対象について

単振り子（ふつうの振り子）を対象として、摩擦を考慮するものとします。

単振り子の式については「単振り子　運動方程式」などと調べれば出てくるのと少し煩雑になってしまいますので割愛させていただきます。

オイラー法

一言で言うと、「微分方程式を用いた数値解析の手法」です。

イメージとしては、「現在の状況と次の状況を比較する」ことが大切だと感じます。

プログラムの作成と実行

single_pendulum.cpp

# include <iostream>
# include <fstream>
# include <cmath>
using namespace std;
constexpr double THETA_START = 12.0;
constexpr double START_TIME = 0.00;
constexpr double DT = 0.01;
constexpr double END_TIME = 100.00;
constexpr double G = 9.8;
constexpr double L = 10.0;
constexpr double M = 2.0;
constexpr double MU = 0.01;

main(){
  // 時間の刻み幅の逆数を格納する変数
  const int reciprocal_dt = round(1/DT);
  //データ数を格納する変数
  const int data_num = (END_TIME-START_TIME)*reciprocal_dt+1;
  //プロット数をカウントする変数,number of plot 
  int plot_num = 0;
  //時間を格納する配列.array of time
  double *time;
  time = new double[data_num+1];
  //角度θを格納する配列,array of degree
  double *theta;
  theta = new double[data_num+1];
  //角度θの微分を格納する配列,array of degree/secont
  double *theta_dot;
  theta_dot = new double[data_num+1];
  double theta_dot_dot = 0.0;

  theta[0] = THETA_START;
  theta_dot[0] = 0.0;
  for(double t=START_TIME; t<END_TIME+DT; t+=DT){
    time[plot_num] = t;
    theta_dot_dot = - G/L*sin(theta[plot_num]/2.0/M_PI) - MU*theta_dot[plot_num];
    theta_dot[plot_num+1] = theta_dot[plot_num] + theta_dot_dot*DT;
    theta[plot_num+1] = theta[plot_num] + theta_dot[plot_num+1]*DT;
    plot_num ++ ; 
  }//単振り子の方程式

  ofstream fout("single_pendulum.txt");
  for(int i=0; i<plot_num; i++){
    fout << time[i] << ", " << theta[i] << endl;  
  }//出力
  fout.close();
}

実行結果はプロット数が多かったのでグラフとして掲載します。

$ gnuplot

gnuplot> plot "single_pendulum.txt"

このプロットから、振動がだんだん減衰していることがわかります。