あらきけいすけのメモ帳

あらきけいすけの雑記帳2

Runge-Kutta-Gill法のPythonコード

4段4次のルンゲ=クッタ法で記憶領域の量が少ないRunge-Kutta-Gill法*1のコードをPythonで書いてみたのでメモ*2。関数rk4()のインターフェースは古典的ルンゲ=クッタ法のPythonプログラムをWikiPediaの記述と同じ変数名で書いてみた - あらきけいすけのメモ帳に合わせてある。解いている微分方程式\ddot{y}=-y, すなわち\frac{d}{dt}\left(y(t),\dot{y}(t)\right)=\left(\dot{y}(t),-y(t)\right), 初期条件は\left(y(0),\dot{y}(0)\right)=\left(\cos0,-\sin0\right)=(1,0), 解は\left(y(t),\dot{y}(t)\right)=\left(\cos t,-\sin t\right).

import numpy

def rk4 ( y, t, h, f ):
  # Gill's low strage method is implemented.
  p= numpy.zeros_like(y); q= numpy.zeros_like(y)
  p=f(t,y); p=p*h*(1./2.); rt2=2**-.5; y+=p; q=p*2
  p=f(t,y); p=p*h*(1-rt2) - q*(1-rt2); y+=p; q=p*2 + q*rt2
  p=f(t,y); p=p*h*(1+rt2) - q*(1+rt2); y+=p; q=p*2 - q*rt2
  p=f(t,y); p=p*h*(1./6.) - q*(2./6.); y+=p;
  return y, t+h

def RHS(t,y):
  return numpy.array([y[1],-y[0]])

y=numpy.array([1.0,0.0]); t=0; h=0.001
for i in range(100):
  y, t =rk4( y, t, h, RHS )
  print(y[0]-numpy.cos(t))

*1:ストレージは3個(このコードでは y, p, q)で、典型的なRunge-Kutta法より係数がちょっと複雑。[参考]ルンゲ・クッタ・ギル法
4次の公式の最小ストレージは2個らしい: David I.Ketcheson, "Runge–Kutta methods with minimum storage implementations", J. Comput. Phys. Vol. 229, Issue 5, 1 March 2010, Pages 1763-1773 https://doi.org/10.1016/j.jcp.2009.11.006 http://www.davidketcheson.info/assets/papers/2010_LSRK_postprint.pdf.

*2:意外にPythonでのチートシート的なページが無い。fortranなどの他言語に移植することを少しだけ意識した書き方にしたつもり。

Raspberry Pi の初回起動でやっておくべきこと

授業のための覚書
ラズパイ(3B, 3B+, zero w, 4B を想定)のGPIOを使って外付けのA/DコンバータやRCサーボモータを動かすことを主目的としている。初回起動ではどうしても本体にキーボード、マウス、ディスプレイの接続をしてOSの設定をしなくてはいけないので「次からはLAN経由でPCから操作するための最小限の設定」のつもり*1

  1. [Raspberry Piの設定]→[インターフェース]で[カメラ]から[リモートGPIO]まで一応、全部、[有効]にしておく。とくに[SSH]と[VNC]はリモートログインに必要。*2
  2. [Raspberry Piの設定]→[システム]の[解像度](RPi4では[設定]→[Screen Configuration]→[ウィザードのウィンドウ部分で右クリック])を1024x768以上にしておく。さもないとリモートからVNCでデスクトップ環境を表示したときに、[Raspberry Pi の設定]のウィザードが画面からはみ出して、非常に操作しづらくなる。(学生にGUIで設定を指導するときに面倒になる)
    • (Raspberry Pi Zero W Rev 1.1) sudo raspi-config → 7 Advanced Options → A5 Resolution
  3. Raspberry Pi Zero W/4 B の場合にはUSBガジェット機能を有効にする*3
    1. /boot/config.txt の最後尾に "dtoverlay=dwc2" と書かれた行を追加する。
    2. /boot/cmdline.txt の "rootwait" と "quiet" の間に " modules-load=dwc2,g_ether" を書き加える。
    3. 母艦のWindows機側の設定はRaspberry Pi Zero W(RPi4BもOK)をUSBガジェットにしてUSBケーブル1本でWindows 10マシンから操作する - あらきけいすけのメモ帳
  4. 次の起動でWiFi経由でリモートログインしたい場合には、WiFiへの接続(ネットワークキーの入力)を済ませる。*4

*1:アカウント管理まわりのことをまだ書いていない。

*2:OpenCVでUSBカメラの画像を表示させたいならVNCでのリモートログインは必須。xrdpでうまくいかなかった経験がある。

*3:この作業はRPi3Bでディスクを作成して、RPi0wに積み替えることも考えて書いている。

*4:次回の接続では(同じサブネットに1台だけつながっているという前提で)ホスト名 "raspberrypi" で接続し、デフォルトのユーザ名とパスワードでログインする。

古典的ルンゲ=クッタ法のPythonプログラムをWikiPediaの記述と同じ変数名で書いてみた

教育用の覚書
Python だと汎用の Runge-Kutta 法のルーチンが6行で書けてしまう*1。以下のコードではルンゲクッタ法ルーチン rk4() を書き換えることなく、1階, 2階, 3階の常微分方程式(ordinary differential equation(s), ODE)を解いている*2*3。関数 rk4() で使っている変数名はRunge–Kutta methods - Wikipediaに揃えた。C言語*4数値計算をコーディング/教育する気力が一気に失われてしまったwwwwww
次のコンプリートコードは次の微分方程式積分して、あらかじめ求めておいた解の関数の値と比較している*5

  1. class getExpt():
    [1階]\displaystyle\frac{dy}{dt}=y (比較する解:y=e^t);
  2. class getSint():
    [2階]\displaystyle\frac{d^2y}{dt^2}=-y (比較する解:\displaystyle\left(y,\frac{dy}{dt}\right)=(\sin t,\cos t));
    [1階連立ODE] \displaystyle\frac{d}{dt}\Big(\underbrace{y}_{y_0},\underbrace{\frac{dy}{dt}}_{y_1}\Big)=\underbrace{\left(\frac{dy}{dt},-y\right)}_{f(t,(y_0,y_1))}
  3. class getExptSint():
    [3階]\displaystyle\frac{d^3y}{dt^3}+\frac{d^2y}{dt^2}+\frac{dy}{dt}+y=0
    (比較する解:\displaystyle\left(y,\frac{dy}{dt},\frac{d^2y}{dt^2}\right)=(e^{-t}+\sin t,-e^{-t}+\cos t,e^{-t}-\sin t));
    [1階連立ODE] \displaystyle\frac{d}{dt}\bigg(\underbrace{y}_{y_0},\underbrace{\frac{dy}{dt}}_{y_1},\underbrace{\frac{d^2y}{dt^2}}_{y_2}\bigg)=\underbrace{\left(\frac{dy}{dt},\frac{d^2y}{dt^2},-y-\frac{dy}{dt}-\frac{d^2y}{dt^2}\right)}_{f(t,(y_0,y_1,y_2))}
  4. class spring():
    [2階]単振動(調和振動子)の問題 \displaystyle m a=- k x, すなわち
    位置:\displaystyle x(t), 速度:\displaystyle v(t)=\frac{dx}{dt}, 加速度:\displaystyle a(t)=\frac{d^2x}{dt^2} より\displaystyle \frac{d^2x}{dt^2}=-\frac{k}{m}x
    (比較する解:(x(t),v(t))=(\sin\omega t,\omega\cos\omega t), \displaystyle\omega=\sqrt{\frac{k}{m}});
    [1階連立ODE] \displaystyle\frac{d}{dt}\Big(\underbrace{x}_{y_0},\underbrace{v}_{y_1}\Big)=\underbrace{\left(v,-\frac{k}{m}x\right)}_{f(t,(y_0,y_1))}

初期条件は比較する解に t=0 を代入して作っている。

import numpy

# classical 4-th order Runge-Kutta method subroutine
# variables letters are the same as those in
# https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
def rk4 ( y, t, h, f ):
    k1 = h * f( t, y )
    k2 = h * f( t + h/2, y + k1/2 )
    k3 = h * f( t + h/2, y + k2/2 )
    k4 = h * f( t + h, y + k3 )
    return y + ( k1 + 2*k2 + 2*k3 + k4 )/6, t + h

# each of the following four classes has
#   __init__(): constructor,
#   rhsODE(t,y): f(t,y) of the 1st order ODE dy/dx=f(t,y),
#   ref(t): a reference solution to the ODE.
class getExpt():
    def __init__ (self):
        print(self.__class__,' y(t)=exp(t)')
    def rhsODE ( self, t, y ):
        return y
    def ref ( self, t ):
        return numpy.exp(t)

class getSint():
    def __init__ (self):
        print(self.__class__,' y(t)=sin(t)')
    def rhsODE ( self, t, y ):
        return numpy.array([y[1],-y[0]])
    def ref ( self, t ):
        return numpy.array([numpy.sin(t),numpy.cos(t)])

class getExptSint():
    def __init__ (self):
        print(self.__class__,' y(t)=exp(-t)+sin(t)')
    def rhsODE ( self, t, y ):
        return numpy.array([y[1],y[2],-y[0]-y[1]-y[2]])
    def ref ( self, t ):
        et=numpy.exp(-t); st=numpy.sin(t);
        return numpy.array([et+st,-et+numpy.cos(t),et-st])

class spring():
    def __init__ (self,mass=1,springConstant=1):
        self.mass= mass
        self.springConstant= springConstant
        m= mass
        k= springConstant
        self.angularFrequency= (k/m)**0.5
        print(self.__class__,' linear spring motion m=',m,' k=',k)
    def rhsODE ( self, t, y ):
        k= self.springConstant
        m= self.mass
        x= y[0]
        v= y[1]
        a= -(k/m)*x
        return numpy.array([v,a])
    def ref ( self, t ):
        w= self.angularFrequency
        x= numpy.sin(w*t)
        v= w*numpy.cos(w*t)
        return numpy.array([x,v])

if __name__=="__main__":

    # choose class to e[x]ecute
    #x= getExpt();
    #x= getSint();
    #x= getExptSint();
    x= spring(mass=2,springConstant=3);

    # set [ref]erence [sol]ution function
    #   and [f]unction appears in the [r]ight [h]and [s]ide of
    #   [O]rdinary [D]ifferential [E]quation
    sol= x.ref; f= x.rhsODE;

    # set parameters
    h = 0.01;   # step-size
    t = 0;      # initial time
    y = sol(t); # initial condition (y[,dy/dt[,d^2y/dt^2]])_t=0
    tEnd = 0.1; # terminal time
    print(t,y,sol(t),y-sol(t))

    # integration and comparison to reference solution
    while t < tEnd:
        y, t = rk4( y, t, h, f )
        print(t,y,sol(t),y-sol(t))

*1:ルンゲクッタ python - Google 検索で検索すると上位に上がるページに意外と汎用性のあるコードが無い。

*2:ルンゲクッタ法は1階の連立常微分方程式を解く計算スキームなので、n階の微分方程式はn個の関数 \boldsymbol{y}=(y,y',...,y^{(n-1)}) に対する連立方程式に書き換える。

*3:Pythonなので関数 rk4() には識別子さえ渡せばいいので、任意サイズの配列 y で計算ができる。例えば3次元流体をスペクトル法で解くとか(配列のメモリアロケーションどうなっとるか知らんけどwwwwww)。

*4:C言語だと構造体へのポインタ
 str *rk4( str *y, double t, double h, (str*)f(double,*str) ) { str *k1, *k2, *k3, *k4; ...; return k?; }
で宣言すればいいのかな?

*5:上掲のプログラムは[4]の調和振動子を解いている。[1][2][3]の方程式を解くためには x= の行を適宜、書き換えればよい。

GCC Developer Lite (GDL)をインストールした 64bit Windows 10 PC 上で gcc をコマンドプロンプトから使うための設定方法

授業用の書きかけの覚書

GDL(バージョン2.5)インストール記事はこちら*1
GCC Developer Lite (GDL)を64bitのWindows 10にインストールし hello world プログラムを作成する - あらきけいすけのメモ帳
H8などのマイコンのプログラミングツールとして GCC Developer lite を 64bit Win10 PC にフルインストールしてしまった場合に、ふつーの Windows x64 (Windows 64bit)プログラミング(C, C++, fortran95*2)もできちゃうので、使うためのセットアップを記す;↓新記述[2022.9.23]
メモ帳などのテキストエディタを使って、コマンド

set PATHC2=\PROGRA~2\BestTech\GCCDEV~1\GCC\x64\bin
set PATHC3=\PROGRA~3\BestTech\GCCDEV~1\GCC\x64\bin
set PATH=%PATHC2%;%PATHC3%;%PATH%
gcc -std=c99 -Wall -lm %~n1.c

をコピペして build.bat というファイル名で保存し、C言語のプログラムのソースコードファイル [なんとか].c を build.bat にドラッグアンドドロップする。*3

↓旧記述
GCC Developer Lite 内の GCC へのパス
C:\Program Files (x86)\BestTech\GCC Developer Lite\GCC\x64\bin
環境変数 Path に書き加える*4
と、コマンドプロンプトを起動して、命令 gcc hello.c で実行ファイル a.exe を作ってくれる。GDLの[ツール(T)]⇒[コンパイルオプション(O)]⇒[設定リスト]の[x64 (Windows)]をいちいち設定するより簡単。

環境変数 Path の設定方法

  1. [スタート(窓マーク)]⇒(スタートメニューが開く)
  2. ⇒[設定(メニューの左端下の歯車マーク)]⇒([Windowsの設定]ウィンドウが開く)
  3. ⇒[設定の検索]枠に「環境変数」と入力⇒(検索結果が現れる)
  4. ⇒[環境変数を編集]を選択⇒([環境変数]のウィザードが現れる)
  5. ⇒[[Username]のユーザー環境変数(U)]内の[Path]を選び、[編集(E)...]をクリックする⇒([環境変数名の編集]ウィザードが開く)
  6. ⇒[新規(N)]をクリック⇒(環境変数の最後の行の下に入力用の枠が現れる)
  7. ⇒その枠に"C:\Program Files (x86)\BestTech\GCC Developer Lite\GCC\x64\bin"と書き込み*5⇒[環境変数名の編集]ウィザードの[OK]をクリック⇒[環境変数]ウィザードの[OK]をクリック

*1:2020.7.1に追記

*2:GCC, the GNU Compiler Collection - GNU Project

*3:PATH2 はGDL 2.5以前のパス、PATH3 は GDL 2.7 のパス。GDL 2.6 は短命だったので、知らない。

*4:新記述:パス名を8.3形式で書けばよい。GDLの吐いたバッチファイルを読めば一目瞭然であった。(2022.9.23)
旧記述→参考:組込みエンジニアのブログ gcc: error: CreateProcess:No such file or directory の対処法 にも書いてあるように、バッチファイルに Path を追加して書いても 、ファイルをコンパイルするときに gcc: error: CreateProcess: No such file or directory が出てしまう。どうもディレクトリ名の " (x86)"の部分が悪さをしているようである。

*5:"C:\Program Files (x86)\"は GCC Developer Lite を64bit Windows PC にインストールしたときのディレクトリ名。32bit OS の場合は""C:\Program Files\"である。

JavaScript/Google Apps Script/Excel/Spreadsheet 覚え

GASでフォルダ名, ファイル名を指定して, Excel ファイルをスプレッドシートに変換する

Google Drive 上の excel と spreadsheet を Google Application Script で操作する練習

function testConvertExcel2Spreadsheet() {
  nameFolder= '[folder name]'
  nameExcel= '[file name].xlsx'
  idSpreadsheet= convertExcel2Spreadsheet( nameFolder, nameExcel )
  Browser.msgBox(idSpreadsheet)
}

function convertExcel2Spreadsheet( nameFolder, nameExcel ) {
  // {String} nameFolder で最初に見つかったフォルダの
  // {String} nameExcel で最初に見つかった Excel ファイルを
  // Google スプレッドシートに変換して、
  // そのスプレッドシートの id を戻す。
  if ( ! nameExcel.match(/.xls/) ) return 0
  idFolder= DriveApp.getFoldersByName(nameFolder).next().getId()
  idExcel= DriveApp.getFolderById(idFolder).getFilesByName(nameExcel).next().getId()
  fileExcel= DriveApp.getFileById(idExcel)
  idSpreadsheet= Drive.Files.insert(
    options = {
      title: fileExcel.getName()
     ,mimeType: MimeType.GOOGLE_SHEETS
     ,parents: [{id: idFolder}]
    }
    ,fileExcel.getBlob()
  ).id
  return idSpreadsheet
}

参考にしたサイト
GAS入門 - DriveAppクラスリファレンス - Qiita
【GAS】GoogleDrive上のExcelをGoogleスプレッドシートに一括変換 - logicoffee プログラミング勉強日記

2019年7月28日3時31分頃, M6.5, 三重県南東沖, 深さ420km, プレート境界に沿って揺れてる

震源三重県沖だが、震度分布がプレート境界に沿って連なっていて、震源から遠い関東から東北の太平洋側が大きく揺れているのが、とても興味深かったのでメモしておく。Hi-netはM6.3, depth:400kmとなっていた。「三重県南東沖でM6.5の深発地震 津波の心配なし 宮城県で震度4「異常震域」」との記事を見たので調べた。震源域はしょっちゅう深いところで地震が帯状に起きているところなので、地球規模では「平常運転」ではある。今回も「大地震」である*1
weathernews.jp

2019年7月28日3時31分頃 M6.5 深さ420km
気象庁のデータ
2019年7月28日3時31分頃 M6.3 深さ400km
防災科学技術研究所のデータ