2023-01-04

英語を極める！

こんにちは、おととしです！

今回は2023年の最初の記事ということで、今年の目標を書いていきたいと思います！

今年の目標は、「英語を極める」です！

このブログは数値シミュレーションをやっていきたいから書いてきたのですが、自分の学びを発信していく場としてやっていくことにしました。
数値シミュレーションを続けるにしても、今の時代、多くの情報が英語でもたらされます。日本語に翻訳されたものが悪いというわけではなく、一次情報を知っておきたいという好奇心があるのです。
最初に英語を極めてしまって、今後の発信のために役立てていきたいと考えています。

私の英語の実力を示すものは、3年前に受けたTOEIC L&R750点くらいです。当時は普段から英語を使っていましたが、今となってはほとんど使っていません。
きっと点数も下がっていることと思います。直近の目標は決まっていませんが、最終的な目標は英検1級かなとぼんやり思っています。
そのための本も早速購入してみました。

【音声DL付】英検1級合格マップ

作者:Jun
アルク

Amazon

今のところ、英語を使えなくて困っていることは正直ありません。しかし、昨年、会社の指示で国際学会を聴講することがありました。
そのときに私は発表の内容さえまともに聴き取れないのに、英語を流暢に話し、活発に議論を行う日本人の研究者が何人もいました。
これから会社を支える新製品を開発しなけらばならない部署として、最新の情報を手に入れることはとても重要です。
これまでなあなあにしてきた部分でしたが、自分は英語を使えないのだから、きちんと勉強しなければいけないのだと考えるようになりました。
どのように勉強していくかは未定ですが、普段から英語に触れ、慣れていくようにしたいと思います。

2022-12-19

針金内の温度分布の時間変化

シミュレーション

こんにちは。おととしです。
針金内の温度分布について、与えられた境界条件から定常状態となった時の温度分布については、以前やってみました。

【よくあるやつ】温度分布シミュレーション - あめつちひらく

今度は、針金の両端の温度のみ固定し、徐々に定常状態へ遷移していく様子をシミュレーションしてみたいと思います。
考えるのは、両端をそれぞれ0℃と100℃に固定した針金です。時間変化を考えるには、まず針金をいくつかの細かい部分に分割します。

次にそれぞれの針金の部分について、隣り合ったところからどのように熱が伝わってくるのかを考えます。
点Cへ、点Aおよび点Bから伝わる熱は次の式のようになると考えます。

これは距離あたりの温度の差に対して、温度の伝わる度合いを示す比例定数kをかけたものとなっています。左からは熱が流出し、右から熱が流入すると考えています。
今度はこれを温度変化と対応させます。熱が伝わってきたときの温度の上がり方は、cΔxに比例すると考えられます。
cは比熱に相当する係数です。ここでΔxを掛けるのは、針金の質量は長さに比例するからです。以上のことを踏まえると、Δtの時間が経った後の点Cにおける温度について、次の式が成り立ちます。

変形すると、

と書くことができます。この式をすべての分割に対して行うことで、針金内の温度分布変化を計算することができるのです。
実際にプログラムを実行すると、次のような図が得られます。

# Const Definition
d=21			# split
a=0			# start point
b=1			# end point
x=np.linspace(a,b,d)

itr=500			# iterations
dt=0.01			# Time interval
dx=x[1]			# X interval
k=0.1			# Temperature Rise Rate
c=1			# Heat Capacity
r=k*dt/(c*dx**2)

u=np.zeros((d,itr,))

# Boundary Condition
u[0,0]=0
u[-1,0]=100
t=0

# Matrix Definition
C=np.zeros((d,d,))
C[0,0]=1
C[-1,-1]=1
for i_row in range(1,d-1):
    C[i_row,i_row-1]=r
    C[i_row,i_row]=1-2*r
    C[i_row,i_row+1]=r

# Calculation
for i in range(itr-1):
    u[:,i+1]=np.dot(C,u[:,i])

# Final Plot
fig=plt.figure(figsize=(4,4))
def update(i, x, u, dt):
    if i !=0:
        plt.cla()

    t=i*dt
    plt.plot(x, u[:,i], color = 'orange')
    plt.title('t= %.1f' %t)

# Animation
def animation_play(x, u, dt, itr):
	params={
		'fig':fig,
		'func':update,
		'fargs':(x, u, dt),
		'interval':30,
		'frames':itr,
		'repeat':False,
	}
	ani=anm.FuncAnimation(**params)
	plt.show()

animation_play(x, u, dt, itr)

kやcには本来、物理的な意味があるはずですが、ここではよくわからないまま使っています。
今回行ったシミュレーションは陽解法と呼ばれる種類のものです。t+Δtの様子を、tのときの値のみを使って表そうとしています。これに対し、t+Δtの様子をtのときの値だけでなく、t+Δtのときの値も用いて計算する方法もあるようです。今回の陽解法では、計算条件によってはうまくシミュレーションが行えません。具体的にいうと、途中で1-2rという数がでてきます。この値が負になるとうまくいかなくなります。ですので、rが0.5を超えないようにdt、dx、k、cを設定しなければならないという縛りがあるのです。

時間変化を計算するので、針金内の温度を何度も計算することになります。その回数をitrで指定しています。uという「箱」を作っておいて、計算するたびに結果を書き込んでいきます。Calculationの部分が終わった段階で、itr回の計算はすべて終了しています。

アニメーションについて、図はfigへと引き渡されます。この時にupdateと呼ばれる関数が呼ばれます。関数は30msごとに呼び出され、引数としてx、u、dtを渡します。呼び出しはitr回行われ、itr回の呼出し後の繰り返しは行いません。

プロットについて、animation_playによって関数updateが呼び出されると、引数としてi、x、u、dtが渡されます。iはアニメーションの関数によって指定されているもので、呼び出した回数（0, 1, 2, ..., itr-1)を表しています。updateが呼び出されるたびに一度プロット画面(fig)をリセットし、プロットするということを行います。

2022-12-10

二次元の温度分布

シミュレーション

以前の記事では、一次元の温度分布をシミュレーションしました。
【よくあるやつ】温度分布シミュレーション - あめつちひらく

次は二次元の温度分布をやりたくなるというものです。
というわけで、自分のできる範囲でやってみました。
参考にしたのは前回と同じく、次の書籍です。

数値シミュレーション入門 (Computer Science Library)

作者:河村哲也
サイエンス社

Amazon

シミュレーションするのは、4辺を一定の温度に保った平板内の温度分布です。

やることはほとんど変わらないと思っています。二次元のものをどのように表し、計算できるようにするかを考えなければなりません。
まずは、平板を次のように分割することにします。dは分割の数、mはdの2乗です(m=d^2)。pythonで計算することを想定し、スタートを0にして描いています。

Pにおける温度は、周囲のA、B、C、Dの温度の影響を受けて、それら4点の平均の温度となると考えます。式に表すと次のようになります。

これは注目しているそれぞれの点について成り立ちます。ですので、平板内のすべての点において当てはめ、連立方程式をたてることができます。平板の4辺はそれぞれ一定の温度に保たれているので、次のような式ができます。

上の式では上手く表せず、申し訳ないです。今回は平板の温度を考えるのですが、計算の都合上、1列に並びかえています。uの部分が平板を分割したそれぞれの部分の温度を表しています。4辺の温度は一定なので、周囲には影響されず、80℃や40℃、0℃のままとなります。平板の内部は周囲の温度の影響を受けるので、行列に-4が表れています。つまり、上の式の1行1行は分割された小平板1つ1つに対応していて、その小平板が4辺にある時は一定の値となり、4辺でないときは周囲の温度の影響を受けるということを表しているのです。この計算を行うと、次のような図が得られます。

図の左下（0のところ）で0℃となっており、10℃刻みで線が引かれています。一番上が80℃になっていることが確認できます。この図を出力するために書いたのが下のコードとなります。

import numpy as np
import matplotlib.pyplot as plt

# Const. Definition
d=25			# split
m=d**2
u0=np.zeros((m,1,))

a=0			# start point
b=1			# end point
x=np.linspace(a,b,d)
y=np.linspace(a,b,d)
X,Y=np.meshgrid(x,y)

# Boundary Condition
for i in range(m):
    if i<d-1:
        u0[i]=80
    elif i>m-d-1:
        u0[i]=0
    elif i%d==0:
        u0[i]=0
    elif i%d==d-1:
        u0[i]=40

# Matrix Definition
A=np.zeros((m,m,))
for i_row in range(m):
    if i_row<d-1:
        A[i_row,i_row]=1
    elif i_row>m-d-1:
        A[i_row,i_row]=1
    elif i_row%d==0:
        A[i_row,i_row]=1
    elif i_row%d==d-1:
        A[i_row,i_row]=1
    else:
        A[i_row,i_row]=-4
        A[i_row,i_row-1]=1
        A[i_row,i_row+1]=1
        A[i_row,i_row-d]=1
        A[i_row,i_row+d]=1

B=np.linalg.inv(A)

# Solve
u=np.dot(B,u0)

# Final Plot
fig=plt.figure(figsize=(3,3))
u_Reshape=np.reshape(u,[d,d])
u_UpsideDown=np.flipud(u_Reshape)
plt.contour(X, Y, u_UpsideDown)
plt.show()

まず、numpyとmatplotlibをインポートします。

次に定数の定義を行います。dが分割数で、1辺をいくつに分割するかを表します。私のPCでは、100を超えると目に見えて計算が遅くなりました。50では余裕でした。今回は平板なので、その2乗の数の計算が行われます。分割数は小さくて良いならなるべく小さい数にしたいところです。平板のx方向とy方向をそれぞれ分割し、後にプロットするためのメッシュをここで定義しておきます。

次に、境界条件を付与します。u0には温度の初期状態を入力します。境界すなわち平板の4辺をそれぞれの条件で表しています。i＜d-1は一番上の辺(80℃)、i＞m-d-1は一番下の辺(0℃)を表します。i%d==0では、dで割り切れる場所、ここでは左の辺(0℃)を表しています。i%d==d-1では、dの倍数より1小さい場所、ここでは右の辺(40℃)を表しています。このとき、頂点の部分を何℃の境界にするかによって微妙にシミュレーション結果が異なるので注意が必要です。

次に行列の定義を行います。初期条件のときと同様、対応する小平板が4辺にあるかどうかで形が異なります。4辺にあるかどうかの条件はu0のときと同様で、そのときA(i_row,i_row)=1となります。4辺にないとき、A(i_row,i_row)=-4で、その上下左右が1となります。上下は列の値に±1をした要素、左右は列の値に±dをした要素が対応します。逆行列をとって左から掛けることで温度分布uが求まります。

最後に、プロットを行います。このままでは直線になっているので、reshape関数で正方形に戻します。また、上下が逆転してしまったので、上下を戻す関数を使っています。この辺はpythonに詳しくないのでもっと良い方法があるかと思います。メッシュの場合は、contour関数を使ってプロットを行います。

2次元になると急に計算量も手間も増えたので大変でした。次は非定常状態から平衡状態となる過程を描写していきたいと思います。

2022-11-27

結果報告！甲種危険物取扱者

資格

こんにちは！

先日11月3日に甲種危険物取扱者の試験を受けてきました。結果が11月15日に発表され、手元に結果通知書が到着したのでその結果を書こうかと思います。

早速ですが、結果は合格でした！！！！

法令、物理・化学、性質・消火それぞれ60%以上で合格なので、思ったより良い点数だなという感想です。

試験問題はやはり法令が一番難しいと感じました。これまで全く聞いたことがないことばかりだったので、ひたすら参考書を覚えるだけでした。勉強する中で、ガソリンスタンドの色分けとか表示とか、身の周りに法令のもと行われているものがあるのだと現実と結びつけることができたのが楽しかったです。実際、この資格を取ったからには建物の構造の決まりとか行政手続きとかそういったことをする際に役立てたいなと思います。

物理・化学については、大学の専攻が化学だったので比較的簡単でした。ただ、危険物の資格ということで燃焼に関する内容が求められ、そこはなじみが無かったので重点的に勉強しました。静電気除去の方法や、分解燃焼など燃焼の仕組み、引火点・燃焼点・発火点の違い、消火剤の種類と効果などこの試験独自の視点が必要だなと感じました。実際に受けた試験問題では、めんどくさい系の計算問題があったり、燃焼とは直接関係のないベンゼンの性質であったり、大学受験をしている気持ちでした。

性質・消火については、とても大変でした。第1-6類まで全ての危険物が対象であり、名称と消火方法とか特有の性質とか参考書をひたすら読みました。勉強時間の半分くらいはこのセクションだった気がします。それでも、何類に分類されるのかとそれぞれの消火方法の特徴をしっかり把握しておけば、なぜこっちの危険物には水が使えないのかとか、この危険物だけハロゲン消火剤が使えないのかとか、だんだんわかるようになってきました。ひたすら読んでいく中で法則性と例外を覚えてなんとか本番に間に合ったといった感じでした。

試験にむけた勉強は8月半ばくらいから始めたので、３か月間くらいでしょうか。もともと化学の知識があったということもあって、法令と個々の危険物の性質を覚えるくらいでした。危険物の性質についても、もともとの知識が助けになってくれました。１日に電車の往復１時間くらいの勉強と参考書についていた問題を解くくらいで、比較的とりやす印象でした。以前からとりたいと思っていた資格なので、今回合格できてとても嬉しいです！これからもどんどん資格取っていきたいです！！

amzn.to

2022-10-09

【よくあるやつ】温度分布シミュレーション

シミュレーション

こんにちは。おととしです。数値シミュレーションといえばこれかな、というのがあるのでやってみました。それは、熱が伝わっていくのをシミュレーションするものです。今回、数値シミュレーションというものがどういうものかよくわかっていないので、本を1冊購入しました。この本をもとにやっていきたいと思います。コードもそれっぽいものを考えて書いてきます。ただ、この本には考え方が書いてあるのみで、実際のコードがどのようなものかは書いてありません。ですので、この記事で紹介するコードは正解であるとは限りません。

数値シミュレーション入門 (Computer Science Library)

作者:河村哲也
サイエンス社

Amazon

最初にやるのは、針金（１次元）の温度分布です。一端を0℃、もう一方の端を100℃に保って放置すると、最終的にどのような温度分布になるのかというものです。

数値シミュレーションでは、針金をいくつかに区切り、それぞれの点においてどのような温度となるかを計算していきます。上の図でいうと、点Cの温度は隣り合う点Aおよび点Bの温度の影響を受けて決定されます。これをプログラムしていきます。

用いる式は次の式です。点Cの温度は、点Aと点Bの温度の平均であるという式です。

これは注目している1つの点に対して成り立ちます。これを針金内のすべての点について当てはめると、連立方程式ができます。また、針金の両端は0℃と100℃に固定されています。これらのことをまとめると、次のような行列を用いた式ができます。

あとは、逆行列を求めて、両辺に左から掛けると終了です。

import numpy as np
import matplotlib.pyplot as plt
# Const Difinition
num=11      # split
a=0         # start point
b=10        # end point
x=np.linspace(a,b,num)

u0=np.zeros((num,1,))

# Initial state
u0[-1]=100

# Matrix Definition
A=np.zeros((num,num,))
for i_row in range(1,num-1):
    A[i_row,i_row-1]=1
    A[i_row,i_row]=-2
    A[i_row,i_row+1]=1
A[0,0]=1
A[num-1,num-1]=1

B=np.linalg.inv(A)

# Solve
u=np.dot(B,u0)

fig=plt.figure(figsize=(3,3))
plt.plot(x, u)
plt.show()

この結果は次のように出てきます。

温度は針金内に線形に分布していることがわかります。

2022-09-27

甲種危険物の演習問題解いてみた

資格

こんにちは。おととしです。
2022年11月3日に迫っている甲種危険物取扱者の試験に向けて、問題集に載っていた演習問題を解いてみました。

らくらく突破　甲種危険物取扱者　合格テキスト＋問題集　第2版

作者:飯島晃良
技術評論社

Amazon

法令：15問、物理および化学：10問、危険物の性質と消火方法：20問となっています。
試験時間が2時間30分とあったので時間も計ってやってみました。

結果は･･･
法令：6/15、物理および化学：8/10、危険物の性質と消火方法：12/20でした。
狙ったように合格ラインギリギリで驚きです。もっと勉強しなければいけません。

やってみた感想としては、まず物理と化学は余裕でした。ただ、燃焼に特化した問題が多く、慣れが必要だと思いました。一般的な物理と化学を知っているだけでは、有利ではありますが簡単には解けない印象です。
法令については、これは覚えるしかないなと思いました。指定数量や製造所等の種類とそれぞれに関わる数字がたくさんあって、系統立ててうまく覚えていこうと思います。
危険物の性質と消火方法は、意外と個別の危険物についてきかれることが多く、これまた覚えなければいけないなと思いました。

テキストを一通り読んでから演習問題を解いたのですが、答えに自信があるものは少なかったです。どこに注力して勉強すればよいか方向性がわかったので、ここからは弱点を重点的に覚えていこうと思います。
めざせ、一発合格！！！

2022-09-27

二重振り子のシミュレーション

シミュレーション

こんにちは。おととしです。
前回、二重振り子の運動方程式の計算が終わったので、次にプログラムをしていきます。

ラグランジアン使ってみた - あめつちひらく

プログラミングは少しだけやったことがあって、学生の頃にFortranとMTALABを用いていました。ただ、どちらもたくさんの人に使われている印象はないので、今回はPythonでやっていこうと思います。とはいっても、そもそもプログラミングは趣味程度しかやったことがないので、参考程度にご覧いただくと幸いです。

プログラムの流れとしては、
１．定数と変数の宣言
２．初期状態
３．θの更新式の定義
４．数値シミュレーションの更新式の定義（ルンゲ・クッタ法）
５．計算
６．プロット
というのを考えています。前回やったのは、３のθの更新式の定義だけです。

まず、定数と変数の宣言を行います。

import matplotlib.pyplot as plt
import matplotlib.animation as anm
import numpy as np

# Const Definition
g=9.8		# m/s2
m1=0.1		# kg
m2=0.1		# kg
l1=0.55		# m
l2=0.45		# m
mu=(m1+m2)/m2

dt=0.01
itr=1000 	# Iteration Times

theta1=[0]*itr
theta2=[0]*itr
V1=[0]*itr
V2=[0]*itr

X1=[0]*itr
Y1=[0]*itr
X2=[0]*itr
Y2=[0]*itr

gは重力加速度、m1とm2は2つのおもりのそれぞれの質量、l1とl2は2つの振り子のそれぞれの紐の長さです。
muは次の計算のためにここで定義しておきます。
dtは時間の区切りで、細かいほどシミュレーションした時の動きが滑らかになります。
itrは計算の終点で、多いほど長い時間のシミュレーションを行えますが、計算にも時間がかかります。
theta1からY2までは計算結果を保持するための「箱」を用意しています。

次に、初期状態を定義しておきます。

# Initial State
theta1[0]=np.pi/180*120
theta2[0]=np.pi/180*120
V1[0]=0
V2[0]=0
X1[0]=l1*np.sin(theta1[0])
Y1[0]=l1*np.cos(theta1[0])
X2[0]=X1[0]+l2*np.sin(theta2[0])
Y2[0]=Y1[0]+l2*np.cos(theta2[0])
t=0

二重振り子は最初の重りの位置でどのような挙動をするか変化します。今回は角度を120°にしています。
V1とV2はθの時間微分を表しています。
X1からY2までは、計算して得られるθの値を直交座標に変換しています。

次に、θの更新式の定義を行います。

# Function Definition
def F_V1(g, mu, l1, l2, theta1, theta2, V1, V2):
	D_theta=theta1-theta2
	F_V=(g*(np.sin(theta2)*np.cos(D_theta)-mu*np.sin(theta1))\
	-(l1*V1*V1*np.cos(D_theta)+l2*V2*V2)*np.sin(D_theta))\
	/(l1*(mu-(np.cos(D_theta))**2))
	return F_V

def F_V2(g, mu, l1, l2, theta1, theta2, V1, V2):
	D_theta=theta1-theta2
	F_V=(g*mu*(np.sin(theta1)*np.cos(D_theta)-np.sin(theta2))\
	+(mu*l1*V1*V1+l2*V2*V2*np.cos(D_theta))*np.sin(D_theta))\
	/(l2*(mu-(np.cos(D_theta))**2))
	return F_V

ここは、前回の記事で求めた式を間違いのないように書いたものとなります。

次に、今回のプログラムの肝となる、ルンゲ・クッタ法の記述となります。

# RK Definition
def RK4(g, mu, l1, l2, theta1, theta2, V1, V2, dt, func1, func2):
	k1_x1=V1
	k1_x2=V2
	k1_v1=func1(g, mu, l1, l2, theta1, theta2, V1, V2)
	k1_v2=func2(g, mu, l1, l2, theta1, theta2, V1, V2)

	k2_x1=V1+dt*k1_v1/2
	k2_x2=V2+dt*k1_v2/2
	k2_v1=func1(g, mu, l1, l2, theta1+dt*k1_x1/2, theta2+dt*k1_x2/2, V1+dt*k1_v1/2, V2+dt*k1_v2/2)
	k2_v2=func2(g, mu, l1, l2, theta1+dt*k1_x1/2, theta2+dt*k1_x2/2, V1+dt*k1_v1/2, V2+dt*k1_v2/2)

	k3_x1=V1+dt*k2_v1/2
	k3_x2=V2+dt*k2_v2/2
	k3_v1=func1(g, mu, l1, l2, theta1+dt*k2_x1/2, theta2+dt*k2_x2/2, V1+dt*k2_v1/2, V2+dt*k2_v2/2)
	k3_v2=func2(g, mu, l1, l2, theta1+dt*k2_x1/2, theta2+dt*k2_x2/2, V1+dt*k2_v1/2, V2+dt*k2_v2/2)

	k4_x1=V1+dt*k3_v1
	k4_x2=V2+dt*k3_v2
	k4_v1=func1(g, mu, l1, l2, theta1+dt*k3_x1, theta2+dt*k3_x2, V1+dt*k3_v1, V2+dt*k3_v2)
	k4_v2=func2(g, mu, l1, l2, theta1+dt*k3_x1, theta2+dt*k3_x2, V1+dt*k3_v1, V2+dt*k3_v2)

	theta1+=(k1_x1+2*k2_x1+2*k3_x1+k4_x1)*dt/6
	theta2+=(k1_x2+2*k2_x2+2*k3_x2+k4_x2)*dt/6
	V1+=(k1_v1+2*k2_v1+2*k3_v1+k4_v1)*dt/6
	V2+=(k1_v2+2*k2_v2+2*k3_v2+k4_v2)*dt/6
	return theta1, theta2, V1, V2

数値シミュレーションでは、ルンゲ・クッタ法と呼ばれる方法がよく用いられているようです。運動方程式のように、微分方程式の形で表されているものの、解析解を求めるのが難しい場合に用いられます。漸化式のようになっていて、前のステップの値を用いて運動の傾きを近似し、次の位置を近似的に求めています。私も方法として用いているだけで、原理や精度については完璧には理解していません。

次に、実際に計算を行います。

# Calculation
for i in range(itr-1):
	theta1[i+1], theta2[i+1], V1[i+1], V2[i+1]=RK4(g, mu, l1, l2, theta1[i], theta2[i], V1[i], V2[i], dt, F_V1, F_V2)

	# Convert to XY
	X1[i+1]=l1*np.sin(theta1[i+1])
	Y1[i+1]=l1*np.cos(theta1[i+1])

	X2[i+1]=X1[i+1]+l2*np.sin(theta2[i+1])
	Y2[i+1]=Y1[i+1]+l2*np.cos(theta2[i+1])

ここまでに関数を定義していたので、計算の記述はシンプルになります。

最後に、プロットしていきます。

# Final Plot
fig=plt.figure(figsize=(6,6))
def update(i, dt, X1, Y1, X2, Y2):
	if i !=0:
		plt.cla()

	if i<150:
		len=i
	else:
		len=150
	T_x1=[0]*len
	T_y1=[0]*len
	T_x2=[0]*len
	T_y2=[0]*len

	for k in range(len):
		T_x1[k]=X1[i-len+k]
		T_y1[k]=Y1[i-len+k]
		T_x2[k]=X2[i-len+k]
		T_y2[k]=Y2[i-len+k]
	
	t=i*dt
	plt.plot([0, X1[i], X2[i]],[0, Y1[i], Y2[i]], color='k')	# String
	plt.plot(0,0,marker='.', markersize=10, color='k')		# Origin
	plt.plot(T_x1,T_y1, color='r')					# m1 Trace
	plt.plot(T_x2,T_y2, color='b')					# m2 Trace
	plt.plot(X1[i],Y1[i], marker='.', markersize=20, color='r')	# m1
	plt.plot(X2[i],Y2[i], marker='.', markersize=20, color='b')	# m2
	plt.xlim(-1.1,1.1)
	plt.ylim(1.2,-1.0)
	plt.title('t= %.1f s' %t)

# Animation
def animation_play(itr, dt, X1, Y1, X2, Y2):
	params={
		'fig':fig,
		'func':update,
		'fargs':(dt, X1, Y1, X2, Y2),
		'interval':100,
		'frames':itr,
		'repeat':False,
	}
	ani=anm.FuncAnimation(**params)
	plt.show()

animation_play(itr, dt, X1, Y1, X2, Y2)

アニメーションまでできるというので、調べてやってみました。また、質点の軌跡を描いていた例があって、キレイだなと思ったのでマネしてみました。

ここまでやってみて、この程度のシミュレーションであれば、思ったより簡単にできたという印象です。ただ、世の中に存在する物理現象は二重振り子だけではありません。それらを網羅していてはどれだけリソースがあっても足りないと感じました。いちいち運動方程式を立てて計算するのではなく、もっと違うアプローチがありそうだと考えています。今度はそっちの方法でいろいろシミュレーションできたらいいな、と思います。

コード全文：

import matplotlib.pyplot as plt
import matplotlib.animation as anm
import numpy as np

# Const Definition
g=9.8		# m/s2
m1=0.1		# kg
m2=0.1		# kg
l1=0.55		# m
l2=0.45		# m
mu=(m1+m2)/m2

dt=0.01
itr=1000 	# Iteration Times

theta1=[0]*itr
theta2=[0]*itr
V1=[0]*itr
V2=[0]*itr

X1=[0]*itr
Y1=[0]*itr
X2=[0]*itr
Y2=[0]*itr

# Initial State
theta1[0]=np.pi/180*120
theta2[0]=np.pi/180*120
V1[0]=0
V2[0]=0
X1[0]=l1*np.sin(theta1[0])
Y1[0]=l1*np.cos(theta1[0])
X2[0]=X1[0]+l2*np.sin(theta2[0])
Y2[0]=Y1[0]+l2*np.cos(theta2[0])
t=0

# Function Definition
def F_V1(g, mu, l1, l2, theta1, theta2, V1, V2):
	D_theta=theta1-theta2
	F_V=(g*(np.sin(theta2)*np.cos(D_theta)-mu*np.sin(theta1))\
	-(l1*V1*V1*np.cos(D_theta)+l2*V2*V2)*np.sin(D_theta))\
	/(l1*(mu-(np.cos(D_theta))**2))
	return F_V

def F_V2(g, mu, l1, l2, theta1, theta2, V1, V2):
	D_theta=theta1-theta2
	F_V=(g*mu*(np.sin(theta1)*np.cos(D_theta)-np.sin(theta2))\
	+(mu*l1*V1*V1+l2*V2*V2*np.cos(D_theta))*np.sin(D_theta))\
	/(l2*(mu-(np.cos(D_theta))**2))
	return F_V

# RK Definition
def RK4(g, mu, l1, l2, theta1, theta2, V1, V2, dt, func1, func2):
	k1_x1=V1
	k1_x2=V2
	k1_v1=func1(g, mu, l1, l2, theta1, theta2, V1, V2)
	k1_v2=func2(g, mu, l1, l2, theta1, theta2, V1, V2)

	k2_x1=V1+dt*k1_v1/2
	k2_x2=V2+dt*k1_v2/2
	k2_v1=func1(g, mu, l1, l2, theta1+dt*k1_x1/2, theta2+dt*k1_x2/2, V1+dt*k1_v1/2, V2+dt*k1_v2/2)
	k2_v2=func2(g, mu, l1, l2, theta1+dt*k1_x1/2, theta2+dt*k1_x2/2, V1+dt*k1_v1/2, V2+dt*k1_v2/2)

	k3_x1=V1+dt*k2_v1/2
	k3_x2=V2+dt*k2_v2/2
	k3_v1=func1(g, mu, l1, l2, theta1+dt*k2_x1/2, theta2+dt*k2_x2/2, V1+dt*k2_v1/2, V2+dt*k2_v2/2)
	k3_v2=func2(g, mu, l1, l2, theta1+dt*k2_x1/2, theta2+dt*k2_x2/2, V1+dt*k2_v1/2, V2+dt*k2_v2/2)

	k4_x1=V1+dt*k3_v1
	k4_x2=V2+dt*k3_v2
	k4_v1=func1(g, mu, l1, l2, theta1+dt*k3_x1, theta2+dt*k3_x2, V1+dt*k3_v1, V2+dt*k3_v2)
	k4_v2=func2(g, mu, l1, l2, theta1+dt*k3_x1, theta2+dt*k3_x2, V1+dt*k3_v1, V2+dt*k3_v2)

	theta1+=(k1_x1+2*k2_x1+2*k3_x1+k4_x1)*dt/6
	theta2+=(k1_x2+2*k2_x2+2*k3_x2+k4_x2)*dt/6
	V1+=(k1_v1+2*k2_v1+2*k3_v1+k4_v1)*dt/6
	V2+=(k1_v2+2*k2_v2+2*k3_v2+k4_v2)*dt/6
	return theta1, theta2, V1, V2

# Calculation
for i in range(itr-1):
	theta1[i+1], theta2[i+1], V1[i+1], V2[i+1]=RK4(g, mu, l1, l2, theta1[i], theta2[i], V1[i], V2[i], dt, F_V1, F_V2)

	# Convert to XY
	X1[i+1]=l1*np.sin(theta1[i+1])
	Y1[i+1]=l1*np.cos(theta1[i+1])

	X2[i+1]=X1[i+1]+l2*np.sin(theta2[i+1])
	Y2[i+1]=Y1[i+1]+l2*np.cos(theta2[i+1])

# Final Plot
fig=plt.figure(figsize=(6,6))
def update(i, dt, X1, Y1, X2, Y2):
	if i !=0:
		plt.cla()

	if i<150:
		len=i
	else:
		len=150
	T_x1=[0]*len
	T_y1=[0]*len
	T_x2=[0]*len
	T_y2=[0]*len

	for k in range(len):
		T_x1[k]=X1[i-len+k]
		T_y1[k]=Y1[i-len+k]
		T_x2[k]=X2[i-len+k]
		T_y2[k]=Y2[i-len+k]
	
	t=i*dt
	plt.plot([0, X1[i], X2[i]],[0, Y1[i], Y2[i]], color='k')	# String
	plt.plot(0,0,marker='.', markersize=10, color='k')		# Origin
	plt.plot(T_x1,T_y1, color='r')					# m1 Trace
	plt.plot(T_x2,T_y2, color='b')					# m2 Trace
	plt.plot(X1[i],Y1[i], marker='.', markersize=20, color='r')	# m1
	plt.plot(X2[i],Y2[i], marker='.', markersize=20, color='b')	# m2
	plt.xlim(-1.1,1.1)
	plt.ylim(1.2,-1.0)
	plt.title('t= %.1f s' %t)

# Animation
def animation_play(itr, dt, X1, Y1, X2, Y2):
	params={
		'fig':fig,
		'func':update,
		'fargs':(dt, X1, Y1, X2, Y2),
		'interval':100,
		'frames':itr,
		'repeat':False,
	}
	ani=anm.FuncAnimation(**params)
	plt.show()

animation_play(itr, dt, X1, Y1, X2, Y2)

参考：
ルンゲ＝クッタ法 - Wikipedia
2重8足振り子 - YouTube