皆さんは「運動方程式」「ニュートン方程式」という言葉を聞いたことがありますか?高校で物理学を勉強している人なら一度は聞いたことがあるでしょう。高校物理を学んでいない人でも耳にしたことがあるかもしれません。
この「運動方程式」「ニュートン方程式」は、読んで字の如く、物体の運動を記述(説明)することのできる方程式で、1680年代にかの有名なイギリスの物理学者、アイザック・ニュートンによって作られた方程式です。この方程式のすごいところは、私たちの身の回りの物理(例えば投げた野球ボールの軌道)から天体の運動までを統一的に理解できるということです。
そんな運動方程式は次のような形をしています。
$$F=ma$$この単純な1つの式が、さまざまな物理を説明することができるのです。
ではこの「運動方程式」を作っている$m,a,F$について具体的に見てみましょう。
まず、右辺に現れている$m$は、考えている物体の質量です。この量は秤(はかり)を用いれば簡単に知ることができますね。左辺にある$F$は、その物体に働いている力の大きさです。中学の理科の授業で力を表す単位として「ニュートン(N)」を習ったことがあると思います。その力を$F$と書きます。
では右辺にある$a$は何でしょうか。これは「加速度」と呼ばれる量で、「ある時間の間に、物体の速度がどれだけ変化しますか」という量です。例えば皆さんが歩いているときに、他の人に背中を押されたとしましょう。すると皆さんの歩く速度は少し早くなりますよね。つまり、力を加えられたことによって加速度が生じているのです。
もっと正確にいうと、加速度$a$は、ある時間$\Delta t$の間に速さが$\Delta v$だけ変化するとき
$$a=\frac{\Delta v}{\Delta t}$$という式で与えられます。最初の速さが$2\ \mathrm{m/s}$で、$0.1$秒間背中を押されたことによって速さが$3\,\mathrm{m/s}$に変化したとしましょう。ここでの速さの変化は$\Delta v=+1\ \mathrm{m/s}$となります。時間の変化は$\Delta t=0.1\ \mathrm{s}$ですから、加速度は$a=\Delta v/\Delta t = 10\ \mathrm{m/s^2}$となります。
しかし、速度や加速度という量は目に見えるものではありません。座標(物体の位置)は目に見えてわかりやすいのに…そんな速度や加速度ですが、実は座標(と時間)を用いて書くことができます。これで加速度を使わずに運動方程式を書くことができますね!
まず、速度は「ある時間の間に、物体の座標がどれだけ変化するか」という量なので
$$v=\frac{\Delta x}{\Delta t}$$とかけます。$x$というのは物体の座標で、一般に時間に依存して異なる値を持ちます。その意味で「$x$は時間$t$の関数」という見方もできるでしょう。
これらのことを用いると、運動方程式は次のように書くことができそうです。
$$m\frac{\Delta (\Delta x)}{(\Delta t)^2}=F$$加速度$a$が少し気持ちの悪い形になってしまいました。これについて少し詳しくみてみましょう。$\Delta x$というのは$x$の変化量でした。つまり$\Delta(\Delta x)$は「$x$の変化量が、時間変化に伴ってどれだけ変化したか」という量に対応します。
これを具体的に考えてみます。ある時刻$t_0$において、物体が$x_0$という座標に存在しています。$-\Delta t$時間前の座標は$x_-$でした。この間の速度は
$$v_-=\frac{x_0-x_-}{\Delta t}$$となります。また、$\Delta t$時間後の座標が$x_+$になるとしましょう。すると、ここでの速度は
$$v_+=\frac{x_+-x_0}{\Delta t}$$です。この$v_-,v_+$という速度は、それぞれ時刻$t_0-\Delta t/2\,,t_0+\Delta t/2$の速度に対応します($t_0\pm \Delta t$ではないことに気を付けてください)。このことから、速度の変化、加速度は次のように書くことができます。
$$a=\frac{v_+-v_-}{\Delta t}=\frac{x_+-2x_0+x_-}{(\Delta t)^2}$$これは時刻$t_0$での加速度です。したがって、運動方程式は
$$m\frac{x_+-2x_0+x_-}{(\Delta t)^2}=F(t_0)$$と表されます($F(t_0)$は時刻$t_0$での力の大きさという意味です)。これによって、加速度を用いて表されていた運動方程式$ma=F$を、時刻と座標を用いた形に書き換えることができました!
先ほどまでは時間$\Delta t$の間に座標や速度がどれぐらい変わるかを考えました。ですが、このような方法で求めた速度・加速度は、もちろん$\Delta t$の選び方に依存してしまいます。同じ物体の運動を見ているのに、人によって速度や加速度の値が変わったら困りますよね。
なので$\Delta t$をとても小さく(無限に小さく)とった時の速度・加速度を考えます。
すでに述べましたが、座標$x$は時刻$t$の関数と見ることができます。そのため、時刻$t$から$t+\Delta t$までの座標の変化を$\Delta t$で割った量、すなわち速度$\tilde v$は
$$\tilde v=\frac{x(t+\Delta t)-x(t)}{\Delta t}$$とかけます。今は$\Delta t$を$0$に近づけたいので、$\Delta t\to 0$という極限を取りましょう。すると、これは「微分」そのものです!
$$\lim_{\Delta t\to 0}\frac{x(t+\Delta t)-x(t)}{\Delta t}=\frac{dx(t)}{dt}=v(t)$$これが「速度」(高校物理では「瞬間の速度」とも言います)の定義です。
同じように考えてみると、(瞬間の)加速度は速度の時間微分になります。
$$a(t)=\frac{dv(t)}{dt}=\frac{d^2x(t)}{dt^2}$$このことから、運動方程式は
$$m\frac{d^2x(t)}{dt^2}=F(t)$$という方程式になります。このように微分が含まれている方程式のことを「微分方程式」と言います。今の場合では、力の関数$F(t)$がわかっているときに、物体の軌道を表す関数$x(t)$がこの方程式を満たしている、関数$x(t)$は何ですか?という問題になるのです。
しかし一般的に、この微分方程式を正確に解く($x(t)$を求める)ことはできません(できるケースはすこしだけあります)。つまり、力の関数がわかっても、運動方程式を眺めているだけでは物体がどのような運動をするのか知ることはできないのです!
そのために「数値計算」という技術が非常に重要になってくるのです。
では、運動方程式を解いて物体の運動を知るためにはどうすればよいでしょうか。そのために、すでに求めた次のような運動方程式の形を用います。
$$m\frac{x_+-2x_0+x_-}{(\Delta t)^2}=F(t_0)$$ここで$x_0$はある時刻$t_0$における物体の座標で、$x_+\,,x_-$はそれぞれ時刻$t_0+\Delta t\,,t_0 -\Delta t$での座標でした。さてここで運動方程式を次のように変形してみましょう。
$$x_+=(\Delta t)^2\frac{F(t_0)}{m}+2x_0-x_-$$この式が意味していることは何でしょうか?それは、「時刻$t_0-\Delta t$と$t_0$での座標、時刻$t_0$において物体にはたらいている力の大きさ$F(t_0)$がわかっていれば、右辺を計算することによって、時刻$t_0+\Delta t$における物体の座標$x_+$が決まる」ということです。ここで得られた$x_+$をさらに用いることで、さらに先の時刻$t_0+2\Delta t$での座標も求めることができます。
このようにして、「最初の座標」「その次の時刻での座標」の2つが分かれば、運動方程式から物体の座標を次々と求めることができるのです。
ここでは具体例として、物体の自由落下を考えてみましょう。地球上で自由落下する物体にはたらく力は
$$F=-mg$$と書くことができます(符号が$-$なのは、力が下向きにはたらいているからです)。ここで$g$は重力定数$g=9.8$です。運動方程式は
$$x_+=-(\Delta t)^2g+2x_0-x_-$$となりますので、これを用いて数値計算で物体の運動を計算します。
東京スカイツリー(高さ634m)から、初速度0で物体を落とす際は、
$$x_0=634, x_1=634$$となります。$x_0$は最初の座標、$x_1$は$\Delta t$だけ時間が進んだ時の座標です。初速度を与えないという条件から、$x_0=x_1=634$としています。
次に載せているプログラムは、時間幅$\Delta t=0.1\ \mathrm{s}$として、自由落下を計算するプログラムです。▶︎ボタンを押すと計算が始まります。計算が終わると、物体の座標(高さ)と速度が時間によってどのように変化するかを示したグラフが現れます(計算を実行しなくても、すでにプログラムの後に計算結果が表示されています)。
自由落下の運動の公式
$$x(t)=h-\frac{1}{2}gt^2$$のグラフも載せていますので、数値計算で求められた結果と比較してみてください。
import sys
from matplotlib import pyplot as plt
print(sys.version_info)
g = 9.8
def free_fall_scheme(h,v,dt,t_max):
x = [h,h+v*dt]
for t in range(2,t_max):
x_nxt = -pow(dt,2)*g + 2*x[-1] - x[-2]
if x_nxt < 0: x_nxt = 0
x.append(x_nxt)
return x
def free_fall_main():
dt = 0.1
t_max = int(12/dt)
h = 634
v = 0
t_ax = [dt*i for i in range(t_max)]
x_num = free_fall_scheme(h,v,dt,t_max)
x_exact = [h-0.5*g*pow(t,2) if h-0.5*g*pow(t,2) > 0 else 0 for t in t_ax]
v_num = [(x_num[i+1]-x_num[i])/dt for i in range(len(x_num)-1)]
v_exact = [(x_exact[i+1]-x_exact[i])/dt for i in range(len(x_exact)-1)]
fig = plt.figure(figsize=(15,5))
ax1 = plt.subplot(1,2,1)
ax1.set_xlabel("time [s]")
ax1.set_ylabel("heiht [m]")
ax1.plot(t_ax[:-1], x_num[1:], c="r", lw=0.5, label="Numerical")
ax1.plot(t_ax[:-1], x_exact[1:], c="b", lw=0.5, label="Exact")
ax1.legend(loc = 'lower left')
ax1.grid()
ax2 = plt.subplot(1,2,2)
ax2.set_xlabel("time [s]")
ax2.set_ylabel("velocity [m/s]")
ax2.plot(t_ax[:-1], v_num, c="r", lw=0.5, label="Numerical")
ax2.plot(t_ax[:-1], v_exact, c="b", lw=0.5, label="Exact")
ax2.legend(loc = 'lower left')
ax2.grid()
plt.show()
free_fall_main()
sys.version_info(major=3, minor=7, micro=13, releaselevel='final', serial=0)
グラフ上の"Numerical"(赤色)は数値計算を用いた結果、"Exact"(青色)は自由落下の公式から得られる結果を示しています。赤色の線は青色の線とほとんど一致しており、数値計算を用いて運動の様子を比較的精度良く計算することができたことがわかります。
なお、今回はスカイツリーから物体を落とすケースを考えましたが、物体が地面に着いた後の跳ね返りなどは考慮していません。そのため、左のグラフでは時刻$11$秒後以降のデータが$0\mathrm{[m]}$になっており、右のグラフでは速度が$-110\mathrm{[m/s]}$から急激に$0\mathrm{[m/s]}$に変わっています。
1章では、私たちの身の回りで起こっている物理は、「運動方程式」という1つの式で説明できることを見ました。しかし、比較的単純な形をしている運動方程式で記述されるはずの私たちの身の回りの世界は、そう簡単に理解することができないほど複雑な現象で満ち溢れています。これはどういうことなのでしょうか。
その大きな理由の一つは、私たちの世界が「複数の物体が集まってできている」からです。私たちの身の回りに溢れているさまざまな「物質」は、すべて100数種類の「原子」からできています。そしてその「原子」は、「原子核」とその周りを回る「電子」からできています。「原子核」もまた、いくつかの「陽子」「中性子」からできていますし、「陽子」「中性子」は「クォーク」と呼ばれる「素粒子」からできています。つまり、私たちが普段目にしている「物質」は、膨大な量の小さな粒子が集まってできているのです。(「素粒子」などについての詳しい説明は、こちらの記事をご覧ください!)
一つ一つの粒子を記述する理論は単純でも、多くの粒子が集まると非常に複雑で多様な現象が現れます。これが、物理法則が単純な形をしているのに、私たちの身の回りにある現象が複雑である理由です。
1つの粒子を説明する理論を知っていても、粒子が複数集まった時に生まれる性質を直ちに理解することはできません。このような問題を「多体問題」と言います。私たち原子核理論研究室は、比較的単純な理論から生まれる、多体系の複雑な現象を理解することを目標に、日々研究をおこなっています。
ここでは多体問題をより身近に感じてもらうために、2つのトピックについてお話しします。
どちらのトピックも当研究室で扱っている問題ではありませんが、多体問題とはどのようなものか、実感していただけると幸いです。
ここからは振り子が2つ繋がった「二重振り子」(下図)の運動方程式を計算して、初期状態が僅かに違うだけで最終的な結果が大きく変わってしまう現象(カオス)を見てみましょう。
※このパートでは、最後に時間がかかる計算を行います(3分程度)。最初に「3.2.1 二重振り子を計算するためのプログラム」に載せているプログラムを実行させてから、説明を読むことをお勧めします。
二重振り子を考える前に、まずは単振り子の運動について考えてみましょう。振り子の振れ角を$\theta$とすると、運動方程式は
$$ml\frac{d^2\theta}{dt^2}=-mg\sin\theta$$となります。高校物理では、$\theta$が十分に小さく$\sin\theta\fallingdotseq\theta$と近似できる場合の運動方程式
$$ml\frac{d^2\theta}{dt^2}=-mg\theta$$を扱うことがほとんどでしょう。この場合には、運動方程式を解くと周期$ T=2\pi\sqrt{\frac{l}{g}}$で振動する解(単振動)が得られます。この周期は振り子の重さ、振れ幅には依存しません。このことを「振り子の等時性」と言います。
しかし、これはあくまでも「$\theta$が十分に小さい」という仮定のもとで成り立つ性質です。では、一般の$\theta$の場合はどうなるでしょうか。この計算は大学生でも難しいのですが、結果は次のようになります。
$$T=4\sqrt{\frac{l}{g}}\int^{\pi/2}_0d\phi\frac{1}{\sqrt{1-\sin^2\frac{\theta_0}{2}\sin^2\phi}}$$なお、$\theta_0$は振り子の最大角です。余談ですが、ここで現れている積分は「第一種完全楕円積分」と呼ばれる積分です。
では、高校で習う振り子の周期$T=2\pi\sqrt{l/g}$と、厳密に計算した周期がどれほど一致しているのか、計算してみましょう。
次に載せているプログラムを用いると、さまざまな$\theta_0$での周期$T$を求めることができます($l=1$に固定してあります)。▶︎ボタンを押して計算してみてください(先ほどと同様に、すでに結果は表示されています)。
from matplotlib import pyplot as plt
from scipy.special import ellipk
import math
pi = 3.1415926536
g = 9.8
l = 1
theta = [t for t in range(1,45,5)]
T_approx = [2*pi*pow(l/g,1/2) for _ in theta]
T_exact = [4*pow(l/g,1/2)*ellipk(math.sin(t*pi/180/2)**2) for t in theta]
fig = plt.figure(figsize=(10,5))
plt.ylim([0,2.5])
plt.plot(theta,T_approx, lw=1, c="r", label="approximation")
plt.plot(theta,T_exact, lw=1, c="b", label="exact")
plt.legend()
plt.xlabel(r"Maximum angle $\theta_0$")
plt.ylabel(r"Period $T$")
plt.grid()
plt.show()
$\theta_0$が大きくなるにつれて、"approximation" ($\theta_0$が小さいと近似した時の結果)よりも"exact"(近似を用いてない場合の結果)の方が大きくなることがわかります。ですが差があると言っても、$\theta_0=40^\circ$の時の差は$0.0733$程度。近似で求めた周期はおおよそ$T=2$ですので、大きく異なっているわけではなさそうです。
ここからは本題であった二重振り子を、「1. はじめに 〜物体の運動を数値計算で解く?〜」で説明した方法で計算していきます(高精度の計算が必要なので、全く同じ方法ではありません)。
運動を計算するためには運動方程式が必要です。単振り子の運動方程式は簡単でしたが、二重振り子の運動方程式はどうでしょうか。なお、二重振り子には2つの変数$\theta_1\,,\theta_2$がありますので、運動方程式は2つ出てきます。
その形はというと…
$$(m_{1}+m_{2})l_{1}\frac{d^2\theta_{1}}{dt^2}+m_{2}l_{2}\frac{d^2\theta_{2}}{dt^2}\cos(\theta _{1}-\theta _{2})+m_{2}l_{2}\left(\frac{d\theta_{2}}{dt}\right)^{2}\sin(\theta _{1}-\theta _{2})+(m_{1}+m_{2})g\sin \theta _{1}=0$$$$ l_{1}l_{2}\frac{d^2\theta_{1}}{dt^2}\cos(\theta _{1}-\theta _{2})+l_{2}^{2}\frac{d^2\theta_{2}}{dt^2}-l_{1}l_{2}\left(\frac{d\theta_{1}}{dt}\right)^{2}\sin(\theta _{1}-\theta _{2})+gl_{2}\sin \theta _{2}=0$$難しすぎる!この方程式を直接解くことは不可能です(数学的にも)。振り子の数が2つに増えただけで、運動を理解するのは難しくなってしまうのです。
ですが、数値計算を用いて運動を計算することはできます。次に載せているプログラムは、この複雑な運動方程式を「4次精度Runge-Kutte法」という方法で計算するためのプログラムです。
分量が多いですが、1つずつ▶︎ボタンを押して進んでいってください。
ここに書いてるプログラムは、計算用のプログラム(言語:C++)と表示用のプログラム(言語:python)です。
左に出ている▼ボタンを押してセルを非表示にした後、「3個のセルが非表示」と書かれている左側にある▶︎ボタンを押 してください。
!g++ --version
g++ (Ubuntu 7.5.0-3ubuntu1~18.04) 7.5.0 Copyright (C) 2017 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%%writefile double_pendulum.cpp
#include <cmath>
#include <ctime>
#include <iostream>
#include <iomanip>
#include <string>
#include <unistd.h>
#include <valarray>
#define TIME 50000
#define N 2
using std::string;
using std::to_string;
using std::valarray;
const double PI = 3.141592653589793;
// angle function
// d2(ang[0])/dt2 = f1
// d2(ang[1])/dt2 = f2
// ↓
// dang0 := d(ang[0])/dt = fd1
// dang1 := d(ang[1])/dt = fd2
// d(dang0)/dt = f1
// d(dang1)/dt = f2
void credit(void){
std::cout << "Double Pendulum Calculation" << std::endl;
std::cout << "Runge-Kutte with 4-th Order Accuarcy" << std::endl;
std::cout << "Created : 2022.06.09" << std::endl;
std::cout << "Credit : Nuclear Theory Group, Department of Physics, Graduate School of Science, Tohoku University" << std::endl;
std::cout << std::endl;
}
// time evolution functions
double f1(valarray<double> ang, valarray<double> l, valarray<double> m, double g){
double A = m[1]*l[1]*sin(ang[0]-ang[1])*pow(ang[3],2)+(m[0]+m[1])*g*sin(ang[0]);
double B = -l[0]*sin(ang[0]-ang[1])*pow(ang[2],2)+g*sin(ang[1]);
return (
A-B*m[1]*cos(ang[0]-ang[1])
)/(
(m[1]*pow(cos(ang[0]-ang[1]),2)-(m[0]+m[1]))*l[0]
);
}
double f2(valarray<double> ang, valarray<double> l, valarray<double> m, double g){
double A = m[1]*l[1]*sin(ang[0]-ang[1])*pow(ang[3],2)+(m[0]+m[1])*g*sin(ang[0]);
double B = -l[0]*sin(ang[0]-ang[1])*pow(ang[2],2)+g*sin(ang[1]);
return (
A*cos(ang[0]-ang[1])-B*(m[0]+m[1])
)/(
((m[0]+m[1])-m[1]*pow(cos(ang[0]-ang[1]),2))*l[1]
);
}
double fd1(valarray<double> ang, valarray<double> l, valarray<double> m, double g){
return ang[2];
}
double fd2(valarray<double> ang, valarray<double> l, valarray<double> m, double g){
return ang[3];
}
// to enable to access the functions like array
double (* const angleFuncArr[])(valarray<double> ang, valarray<double> l, valarray<double> m, double g){
fd1,fd2,f1,f2,
};
// convert the angle to euclid postion
// the return is string-formated-data
string euclid_position(valarray<double> ang, valarray<double> l){
valarray<double> x,y;
x.resize(N+1), y.resize(N+1);
x[0] = 0.0, y[0] = 0.0;
string text = "";
for(int i = 0; i < N; i++){
text += to_string(x[i]) + "," + to_string(y[i]) + ",";
x[i+1] = x[i] + l[i]*sin(ang[i]);
y[i+1] = y[i] - l[i]*cos(ang[i]);
}
text += to_string(x[N]) + "," + to_string(y[N]);
return text;
}
// generate random string
string gen_random(const int len) {
std::srand( time(NULL) );
static const char alphanum[] =
"0123456789"
"ABCDEFGHIJKLMNOPQRSTUVWXYZ"
"abcdefghijklmnopqrstuvwxyz";
string tmp_s;
tmp_s.reserve(len);
for (int i = 0; i < len; ++i) {
tmp_s += alphanum[rand() % (sizeof(alphanum) - 1)];
}
return tmp_s;
}
int main(void){
credit();
// generate file id
string ID = gen_random(16);
// open file
FILE *fp;
if( N == 2 ){
string filename = "double_pendulum_" + ID + ".txt";
fp = fopen(filename.c_str(), "w");
std::cout << "> The data file is " << filename << std::endl;
}else{
std::cout << "> Ops! 2 != N case is currently not supported!" << std::endl;
return 0;
}
// global parameters
const double g = 9.8;
const double dt = 0.001;
int i, j, t;
// set angles
// angle[TIME][2*number of angles]
valarray<valarray<double>> angle;
angle.resize(TIME);
for(t = 0; t < TIME; t++){ angle[t].resize(2*N); }
// initial angle
double theta0, phi0;
std::cout << "> Enter the init angle theta_1[deg] : ";
std::cin >> theta0;
angle[0][0] = PI*theta0/180.0;
std::cout << "> Enter the init angle theta_1[deg] : ";
std::cin >> phi0;
angle[0][1] = PI*phi0/180.0;
angle[0][2] = 0.0, angle[0][3] = 0.0; // velocity
// length of arms and mass
// taking all identical
valarray <double> l(N);
l = 1.0;
valarray <double> m(N);
m = 1.0;
fprintf(
fp, "{'l1':%.5f,'l2':%.5f,'m1':%.5f,'m2':%.5f,'dt':%.5f,'TIME':%d,'ang0':%.15f,'dang0':%.15f,'ang1':%.15f,'dang1':%.15f}\n"
,l[0],l[1],m[0],m[1],dt,TIME,theta0,phi0,0.0,0.0
);
// RK4 scheme
// k1~k4 + k0
// k0: dummy k ( to calculate ks with single function )
valarray<valarray<double>> k;
k.resize(5);
for(i = 0; i < 5; i++){ k[i].resize(2*N); }
k[0] = 0.0;
// weight
valarray <double> w = {1.0,0.5,0.5,1.0};
for(t = 0; t < TIME-1; t++){ //for each time
for(i = 1; i <= 4; i++){ //calculate k1~k4
for(j = 0; j < 2*N; j++){ //calculate for each angle freedom
k[i][j] = dt*angleFuncArr[j](k[i-1]*w[i-1]+angle[t],l,m,g);
}
}
// time evoluation
angle[t+1] = angle[t] + (k[1]+2.0*k[2]+2.0*k[3]+k[4])/6.0;
// convert angle -> euclidian position
string dat = euclid_position(angle[t],l);
// write to the file
fprintf(fp,"%s\n",dat.c_str());
}
return 0;
}
Overwriting double_pendulum.cpp
%matplotlib inline
import os,sys
import json
from IPython.display import Image
from matplotlib import (
pyplot as plt,
animation
)
FLAME = 1000
def loading(path,delimiter=","):
try:
print(f"> Loading {path} ...")
with open(path, "r", encoding="utf-8") as fp:
param = json.loads(fp.readline().replace("\'","\""))
data = [list(map(float,l.split(delimiter))) for l in fp.readlines()]
except Exception as e:
print(e)
print("ERROR: could not load the file. Check if the file name is correct or file exists.")
sys.exit(1)
return param,data
def drawing_double(xydata):
L = len(xydata)
N = len(xydata[0])
dat = xydata[::L//FLAME]
fig = plt.figure(figsize=(7,7))
ax1 = fig.add_subplot(1,1,1)
txt_title = ax1.set_title('')
ax1.set_xlabel("X")
ax1.set_ylabel("Y")
ax1.set_xlim(-(N//2+0.5),N//2+0.5)
ax1.set_ylim(-(N//2+0.5),N//2+0.5)
ax1.grid()
imgs = []
colors = ["r","g","b","black"]
for d in dat:
img = []
for j in range(N//2):
txt_title.set_text(f"time = {j*0.001}")
img.append(plt.scatter(d[2*j],d[2*j+1],marker=".",c=colors[j],s=40))
imgs.append(img)
ani = animation.ArtistAnimation(fig, imgs, interval=50)
plt.close()
return ani
def graphics_double(filename):
param,data = loading(filename)
print("> Drawing the Graph ...")
ani = drawing_double(data)
print("> Visualizing the Graph ... ( It takes about 3 min. )")
gifname = filename.replace(".txt",".gif")
ani.save(gifname,writer="pillow")
print("> Opening the gif file ...")
return gifname
ここからは一個づつ実行していきましょう。
計算を実行する準備です。
!g++ -o double_pendulum double_pendulum.cpp
次に計算を行います。初期条件として2つの角度の値を入力してください。
$\theta_1$の値は
> Enter the init angle theta_1[deg] :
の後に、$\theta_2$の値は> Enter the init angle theta[deg]
の後にそれぞれ入力します。
> The data file is ~~.txt
に記載されているファイル名は、データを表示するために必要ですのでコピーしておいてください。
なお、今回の計算は$l_1=l_2=1.0 \mathrm{[m]}\,, m_1=m_2=1 \mathrm{[kg]}$として計算を行なっています。
!./double_pendulum
Double Pendulum Calculation Runge-Kutte with 4-th Order Accuarcy For the Open Cumpus 2022 Created 2022.06.09 Credit : Nuclear Theory Group, Department of Physics, Graduate School of Science, Tohoku University > The data file is double_pendulum_T0XzXPMVYysJOJsp.txt > Enter the init angle theta[deg] : 35 > Enter the init angle phi[deg] : 135
計算が終了したので、データをグラフにします。先ほど表示されていたデータファイル名("doublependulum***.txt")を
gif = graphics_double("double_pendulum_***.txt")
の括弧の中に入れてください。
例)ファイル名が"double_pendulum_TohokuOpenCumps2022.txt"なら、
gif = graphics_double("double_pendulum_TohokuOpenCumps2022.txt")
※クォーテーションマークを忘れずに入力してください!
▶︎ボタンを押して実行すると、二重振り子のアニメーションが作成され、表示されます。
※この処理は3~4分程度かかります。
gif = graphics_double("double_pendulum_T0XzXPMVYysJOJsp.txt")
Image(gif, format='png')
> Loading double_pendulum_T0XzXPMVYysJOJsp.txt ... > Drawing the Graph ... > Visualizing the Graph ... ( It takes about 3 min. ) > Opening the gif file ...