日々の雑記帳

日々感じたこと、学んだことの備忘録

Pythonによる計算物理（３）

Study Python

本稿から具体的な数値計算の手法についての話に入る。 まずは、古典力学の基本方程式であるNewtonの運動方程式を数値的に解く手法から。

古典力学と常微分方程式

  
  m \frac{d^2 x(t)}{dt^2}=F(x(t),v(t),t)

ここで、 F(t,x(t),v(t))は外部から質点に作用する力である。
t：時間、 x(t)：質点の位置、 v(t)：質点の速度である。

  • 一般に、 Fは時刻 tと位置 x(t)および v(t)に依存する。
  • たとえば、重量の場合は F=mgとなる。
  • 強制振動の場合は、 F=F_0cos(\omega t)。
  • バネの場合は、 F=-kx(t)。
  • 空気抵抗の場合は、 F=-bv(t)。

といったように、作用する力の種類によって、 Fの関数形が決まる。

2階常微分方程式の変形

再び、Newtonの運動方程式に戻る。

 m\frac{d^2x(t)}{dt^2}=F(x(t),v(t),t)

2階の常微分方程式を数値計算するには、新たに速度 v(t)=\frac{dx(t)}{dt}を導入することで、以下のような２つの1階の常微分方程式に分けると、数値解法を適用しやすくなる。


\begin{align}
\frac{dv(t)}{dt} &= \frac{1}{m}F(t,x(t),v(t)) \\
\frac{dx(t)}{dt} &= v(t)
\end{align}

常微分方程式は、一般にN個の成分 y_i(t)(i=1,...,N)それぞれに対する1階の常微分方程式として、以下のような形で表せる。

 \frac{dy_i(t)}{dt} = f_i(t,y_1(t), ..., y_N(t))

ベクトルの形で書くと以下となる。

 \frac{d\boldsymbol{y}(t)}{dt}=\boldsymbol{f}(t,\boldsymbol{y}(t))

ここで、ベクトル表記 \boldsymbol{y}=(y_1,...,y_N)を用いている。
Newtonの運動方程式を一般的なベクトル表記の形式に当てはめると、以下のようになる。
- 1次元の場合　 \boldsymbol{y}=(x, v)と定義できる。
- 2次元の場合　 \boldsymbol{y}=(x, y, v_x, v_y)と定義できる。
- 3次元の場合　 \boldsymbol{y}=(x, y, z, v_x, v_y, v_z)と定義できる。

常微分方程式の数値解法

基本的な考え方として、位置 \boldsymbol{y}、時間 tを離散化する。
時刻 t=t_nでの位置ベクトル \boldsymbol{y}(t_n)=\boldsymbol{y_n}とおく。
\boldsymbol{y}を既知として、次の時刻[tex: t{n+1}]での位置ベクトル[tex: \boldsymbol{y}{n+1}]を以下のように逐次的に求める。

 \boldsymbol{y}_{n+1}=\boldsymbol{y}_n + \Delta\boldsymbol{y}

ここで、 \Delta\boldsymbol{y}を見積もる方法として、次のようなものがある。
- 陽解法 - 陰解法 - 予測子-修正子法

陽解法

陽解法は、もっとも直接的な方法で、代表的なのがオイラー法である。

<基本的な考え方> - t=t_nにおける関数 \boldsymbol{y}(t)の接線で \boldsymbol{y}(t)を近似。 - 接線は、 \boldsymbol{y}(t_n+h)を t_n=tの周りでテイラー展開し、[tex: h2]以上の項を無視することで得られる。

式で書くと以下となる。

 \boldsymbol{y}(t_n+h)=\boldsymbol{y}(t_n)+f\frac{d\boldsymbol{y}(t_n)}{dt}+O(h^2)

Euler法（中点法）

オイラー法の考え方をより精度を高いものにした手法として、中点法がある。
以下に中点法の考え方と模式図を示す。

  • 傾き k_1で中間点 t=t_n+h/2まで進み、そこで新たに傾き k_2を見積もる。
  • 傾き[yex: k_2]を使って、元の出発点から終点 t=t_n+hまで進める。

式で書くと、次のようになる。

 
\begin{align}
k_1 &= \boldsymbol{f}(t,\boldsymbol{y}_n) \\
k_2 &= \boldsymbol{f}(t+\frac{h}{2}, \boldsymbol{y}_n+\frac{h}{2}k_1) \\
\boldsymbol{y}_{n+1} &= \boldsymbol{y}_n + hk_2
\end{align}

ポイント

  • 新たな傾き k_2は t=t_n+\frac{h}{2}で計算するが、 \boldsymbol{y}_{n+1}の計算は、元の出発点 t=t_nに戻って、 k_2を使って計算する。
  • 1回の時間発展において、関数 \boldsymbol{f}(t,\boldsymbol{y}(t))を2回計算することで、精度の次数が1つ向上する。

ルンゲ・クッタ法

精度と実装の簡潔さのバランスが取れていて最もよく使われる手法が4次のルンゲ・クッタ法（Runge-Kutta method）である。
4次のルンゲ・クッタ法の詳細な考え方は省略するが、以下の式で与えられる。


\begin{align}
\boldsymbol{y}_{n+1} &= \boldsymbol{y}_n + \frac{h}{6}(k_1+2k_2+2k_3+k_4) \\
k_1 &= \boldsymbol{f}(t,\boldsymbol{y}_n) \\
k_2 &= \boldsymbol{f}(t+\frac{h}{2}, \boldsymbol{y}_n+\frac{h}{2}k_1) \\
k_3 &= \boldsymbol{f}(t+\frac{h}{2}, \boldsymbol{y}_n+\frac{h}{2}k_2) \\
k_4 &= \boldsymbol{f}(t+\frac{h}{2}, \boldsymbol{y}_n+\frac{h}{2}k_3) \\
\end{align}

Pythonによる計算物理（２）

Python Study

Matplotlib基本事項メモ

ここでは、数値計算結果をグラフなどの可視化の際に便利なMatplotlibライブラリの基本的事項を記しておく。

Matplotlib

Matplotlibはグラフ描画用のライブラリである。数値計算の結果を簡易的に確認するだけではなく、論文やプレゼンにも使える見栄えのよいグラフ作成までできる。
Matplotlibを使用するには、次のようにpyplotをインポートする。

import matplotlib.pyplot as plt

２つのインターフェース

Matplotlibには以下２つのインターフェース（呼び出し方）が存在する。
- オブジェクト指向インターフェース（object-orinted interface） - pyplotインターフェース（state-based interface） 以下に２つのインターフェースでのMatplotlibよりグラフを描画するコード例を示す。







import numpy as np
x = np.linspace(-2, 2, 101)    # -2から2までを100等分した配列を生成
y = np.tanh(x)                 # 各x点におけるtanh(x)を計算





オブジェクト指向インターフェースを使ってグラフを描くコードは以下となる。



fig, ax = plt.subplots()       # まずFigureオブジェクトとAxesオブジェクトを生成
ax.plot(x, y)                  # axesオブジェクトにプロット
ax.set_xlabel('x labael')      # x軸にラベルをつける
ax.set_ylabel('y labael')      # y軸にラベルをつける
fig.show()                     # グラフを表示





/var/folders/yt/sj88vtyd5d1d_gdrrqdq469h0000gn/T/ipykernel_29949/1623430655.py:5: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()                     # グラフを表示



[f:id:AKH1116:20240131011919p:plain]




次に、同じ図をpyplotインターフェースで描くコードを以下に示す。



plt.plot(x, y)
plt.xlabel('x labael')
plt.ylabel('y labael')
plt.show()








上記のオブジェクト指向インターフェースとpyplotインターフェースを比べると、pyplotインターフェースの方がシンプルに書けるが、オブジェクト指向インターフェースの方がより手の込んだグラフを作成できるので、こちらで書くことを推奨する。　　
たとえば、オブジェクト指向インターフェースを利用すると、１つの図の中に複数のグラフを描くことができる。



1つの図に複数のグラフを描く方法



1つのFigureの中に複数のAxesオブジェクトを作り、各Axesオブジェクトに対してグラフを描けば良い。

複数のAxesオブジェクトを生成するには、plt.subplots()の引数に、縦と横に並べるグラフの数を与える。



x = np.linspace(0, 1, 51) * np.pi * 2
y1 = np.sin(x)
y2 = np.cos(x)
y3 = x
y4 = x*x
fig, ax = plt.subplots(2, 2)
ax[0,0].plot(x, y1)
ax[0,1].plot(x, y2)
ax[1,0].plot(x, y3)
ax[1,1].plot(x, y4)
plt.show()











    
    
  



    
  

  




      
      
    
  
    
      
        

  

    

  
    

    
      
    
      

  

  Pythonによる計算物理（１）



  
  

  

    
    Python
    
    Study
    
  



  

  




    


    

  
    
現在、以下の書籍で数値計算（特に物理に関する）におけるPythonの活用を学習している。

Pythonによる計算物理
Pythonによる計算物理
Amazon



まずその１として、Pythonで数値計算コードを書く際に威力を発揮するNumpyライブラリごく基本的な事項のメモである。
※VScode上でJupyternoteを作成しそれをGithub Gistへ登録、はてなブログのGist貼り付け機能を使っているが、Jupytetnote貼り付け部分は別のwindowで表示されるようで、長いJupyternoteだとスクロールする必要があり、ちょっと見にくい・・・






Numpy基本事項メモ


    
    
  



    
  

  




      
      
    
  
    
      
        

  

    

  
    

    
      
    
      

  

  MacでC言語の開発環境を整備する



  
  

  

    
    プログラミング
    
    C言語
    
  



  

  




    


    

  
    
目的



この記事は、MacでC言語の開発環境を整備した際のメモである。



コンパイラのインストール



MacではXcodeをインストールすることで、Mac専用のC言語開発環境であるclangが自動的にインストールされる。  xcodeのインストールについては、ターミナルを開いて以下のコマンドを打つことでインストールされる。



xcode-select --install




エディタの準備



プログラムを書くエディタとしてはVisual Studio Code（VS code）があれば十分であろう。



VScodeにC言語開発環境のプラグインをインストール



拡張機能：マーケットプレイスより以下のC/C++をインストールする。

　　



VScodeでプログラムを書き、コンパイル



たとえば、以下のようなプログラムをVScodeで書き、「HelloWorld.c」で保存する。



#include <stdio.h>
int main(void){
     puts("Hello, World!");
}




VScode内にターミナルを表示して（上のメニューバーのターミナルをクリックするとウィンドウ内に現れる）以下のコマンドを打つことでコンパイルを実行する。



gcc -o HelloWord HelloWorld.c




  



プログラムを実行する



VScode内のターミナルで以下のコマンドを実行することで、プログラムを実行する。実行後、ターミナルに「Hello, World!」と表示されれば成功。



./HelloWorld

    
    
  



    
  

  




      
      
    
      
        

  

    

  
    

    
      
    
      

  

  2024元旦



  
  

  

    
    雑感
    
  



  

  




    


    

  
    
函館旅行から大晦日に帰ってきて、年越しは昨年同様自宅でゆっくり過ごした。

旅行先での年越しも非日常感があり、良いものであるが、やっぱり自宅が一番落ち着く。



天気が良かったので、湘南海岸まで散歩に出かけて、帰り途中に近くの龍口寺で初詣をしてきた。
天気も良かったこともあり、富士山もはっきりと見えた。



そうこうしているうちに、夕方に北陸地方で震度7の地震が発生したとのニュース速報が入り、ほとんどの正月特番番組から地震に関するニュース報道に切り替わった。東日本大震災ほどではないが、津波も発生していた。元旦から大きな地震で被害が出ているのは何とも心が痛む。
2024年最初から何とも不吉なスタートであるが、自然災害など自分がコントロールできない事象なので、せめて外部環境に合わせて対応できるように柔軟な姿勢を持ち続けていきたい。



    
    
  



    
  

  




      
      
    
  
    
      
        

  

    

  
    

    
      
    
      

  

  2023年末函館旅行



  
  

  

    
    旅行
    
  



  

  




    


    

  
    
2023年の仕事納めは28日で、翌日の29日〜31日の2泊3日の旅程で函館へ行ってきた。
29日は妻が仕事納めで、就業後直で東京駅19:20の最終の新幹線で函館へ向かい、函館の宿泊ホテルにチェックインしたのは日付が回っていたこともあり、初日は移動のみ。29日の東京駅は、ちょうど帰省ラッシュのピークのようで、かなりの混雑ぶりだった。



函館では、朝市で海鮮を食べ、八幡坂、元町周辺の教会など洋風建築を見学、五稜郭の見学、夜景鑑賞といったほぼ王道の観光ルートを回った。

函館は、見どころがコンパクトにまとまっていて、1日でも回れるから短い旅程でも満足できた。



Youtubeで青函連絡船がまだ通っていた頃の函館駅の様子の動画を見たことがあるが、その頃はかなり人通りが多く、かなり活気がある様子が窺える。
今回も観光客もいてそれなりの人はいたが、外国人が多い印象を受けた。確かに、函館駅はまだ新しく、洗練されていて駅前も広く余裕があるのだが、青函連絡船がまだ現役だった頃の動画を見ると、昔の昭和感（レトロ感）万歳の方が賑やかさと元町や八幡坂といった洋風な雰囲気が混ざった感じがどことなく魅力的に思えるのである。



青函連絡船が退役したのは、確か1987年で小学校に上がったばかりで乗る機会がなかった。現在のように新幹線で便利に行けるのもいいが、連絡船みたいな旅情がある旅行もしてみたいと思った。



名物のイカが品薄で活イカの刺身食べれなかったのが少し残念だったが、函館はやはり、魚介類が美味しかった。

イカが獲れるのは冬場では2月が多いようで、次はイカが食べれる時期にリベンジしたい。

    
    
  



    
  

  




      
      
    
  
    
      
        

  

    

  
    

    
      
    
      

  

  2023年夏 金沢旅行



  
  


  

  




    


    

  
    
8/11(金)〜8/13(日)の２泊3日の旅程で金沢へ行ってきた。
移動手段は車。圏央道〜関越道〜上信越道〜北陸道の片道約500kmという久しぶりの長距離ドライブであった。



初日は、藤沢の自宅を7時頃に出発したが、途中渋滞にはまったこともあり、金沢に着いたのは夕方の5時頃。この日は、北陸の鮮魚を堪能できる小料理屋で食事を楽しんだ。



2日目は、主に金沢市内を観光した。近江市場で朝昼兼用で海鮮寿司を食した。その後、東茶屋街へ移動し、昔ながらの町屋が並ぶ街並みを楽しんだ。昼間ということもあり、かなり暑く、途中町屋カフェのような場所で休憩した。雰囲気も相まってか町屋カフェで飲んだアイスコーヒーや抹茶カフェが冷たく美味しかった。



同じ古都でも京都や鎌倉とはやはり違う。あまり混雑がひどくないということもあるが、金沢はどこか上品な雰囲気が漂う街という感じがした。
やはり加賀百万石と言われることもあり、文化や芸術が発達して余裕があるのだろうか。



3日目（最終日）は、能登半島をドライブということで、能登島へ行った。
以前、確か2015年（もう８年前か・・・）に金沢に旅行した時も行ったが、その時の海や島へ渡る橋からの景色が記憶に残っていて、今回も行ってみたいと思っていた。幸い、天気も良くて、海の青さが強調されていて、８年前の景色も思い出された。
また、能登島の海沿いの道路にある洋食屋（ハイカラ食堂）で昼食を取ったが、ここも８年前に行った場所だった。海が見えるところで、景色がよかった。



台風が近づいていて、天気が心配だったが、滞在中は良い天気が続き、また食事が美味しく、良い旅行であった。
国内で１カ所を旅行するには、２泊３日程度がちょうど良い期間であると思う。