「プログラミングと数値計算」カテゴリーアーカイブ

弾道計算のコード(Excel)

Excelのマクロ(VBA)で弾道計算を行うコードを作成しました。

https://slpr.sakura.ne.jp/qp/supplement_data/BBullet_course_excel/BBbullet_course_ver1.41.xlsm
上記プログラムの利用の際は連絡先と著作権についてをご覧ください。

以下の微分方程式を計算します。

エクセルシート”bullet”
\(
\begin{align}
m\frac{d^2 \vec{r}}{dt^2}&=
-mg\vec{k}-\frac{1}{2}C_d \rho \pi R^2 |\vec{V}|\vec{V}
-C_l \frac{4}{3}\pi R^3 2\rho |\vec{\omega}| |\vec{V}|\frac{\vec{V}\times\vec{\omega}}{|\vec{V}\times\vec{\omega}|} \\
\frac{d\mathbf{\omega}}{dt}&=N/I \\
N&=\frac{\rho C_f}{2}R\frac{4\pi R^2}{2}\left(|v_{up}|v_{up}+|v_{down}|v_{down}\right)
\end{align}
\)
回転の減衰は本来、弾道計算(BB弾)の理論で導いた通り積分を行わなければなりませんが、簡単にするためにBB弾を真四角の箱として近似して導いています。
この近似は本来の積分値と1パーセント程度の違いしか出てこないので、BB弾の計算には非常に良い近似です。
回転の減衰を取り入れたくない場合、入力パラメータの”回転の減衰”項目をNOにすることでそのように設定が出来ます。
ゼロインを固定した時の弾道の探索機能は実装していないので出来ません。
Excelの設定項目では、Fortranで実装したコードから求めたパラメータをあらかじめ用意しておきました。
設定されたものと書いてあっても、高精度に積分計算をして求めた結果と上の近似をもちいた本エクセルの結果は違くなります。数10cm程度のずれは出てきますので、その点は注意してください。

エクセルシート”barrel”
\(
\begin{align}
m\frac{d^2}{dt^2}x(t)&=\left[P(t)-P_0\right]S_{BB}-\frac{1}{2}C_d \rho \pi R^2 |v(t)|v(t) \\
m_s\frac{d^2}{dt^2}x_0(t)&=-k\left[x_0(t)-x_B-l\right]-[P(t)-P_0]S_s -F_f\label{bbin2} \\
\frac{d}{dt}V(t)&=v(t)S_{\text{barrel}}-v_0(t)S_s+c'(S_{barrel}-\pi R^2)\mathrm{sgn}(P(t)-P_0)\sqrt{\frac{2}{\rho}|P(t)-P_0|}
\end{align}
\)

内部は以下の構造をしていると想定しています。詳細はバレル内部でのBB弾の方程式をご覧ください。

計算手法はどちらも陽的4次ルンゲクッタ法を用いています。

使い方


  1. ファイルのダウンロード
  2. 最新
    https://slpr.sakura.ne.jp/qp/supplement_data/BBullet_course_excel/BBbullet_course_ver1.41.xlsm

    をクリックでダウンロード、もしくは
    右クリックして、”名前を付けてリンク先を保存”からダウンロードしてください。

  3. 計算パラメータの入力
  4. 入力する場所は、色付きのセルです。

    • 弾道計算について

    • 出力結果は大きく分けて2つの項目に分かれています。
      1つ目は定数、射出高さに着弾、地面に着弾の項目は重要と思われるパラメータの抜粋、2つ目は弾道の時間ごとのデータです。
      グラフは、2つ目のデータを元に表示しているだけですので、必要に応じて追加、変更を行ってください。

    • バレル内部の計算について

    • 出力結果は大きく分けて2つの項目に分かれています。
      1つ目は定数、内圧が大気圧に等しい時に射出される時と指定したバレル長におけるパラメータの抜粋、2つ目はBB弾、ピストンの位置、速度の時間ごとのデータです。
      グラフは、2つ目のデータを元に表示しているだけですので、必要に応じて追加、変更を行ってください。

    手入力で入力することも可能ですが、参考設定と書かれたボタンを押すことで、予め設定されたパラメータを読み込ませることもできます。

  5. 計算開始
  6.   入力パラメータを入力した後、計算開始ボタンを押すことで計算が開始されます。

注意点

常識的なパラメータの範囲内でしか動作確認をしておりません。
この結果はあくまで一例です。本シミュレーションを信用して設計・組立を行い、本シミュレーションと結果が違う事があります。その場合私は責任を負いませんので、その点をご了承して使用してください。
BB弾で出来るからと言って、銃弾等の計算はそれなりの信用しかできないことを考慮してください。

追記


バレル内部の計算において、抜弾抵抗の考慮はしていません。すなわちホップが全く掛かっていない状況の計算です。
もし、ホップを掛けた状況を再現したいのであれば、以下のような考えで出来るでしょう。

まず、ホップを掛けた状態でBB弾をバレル内部から引き出す際に必要な力をどうにかして実験的に求めます。
このホップを掛けた際に必要な力\(F\)は、内部圧力\(P\)とBB弾の面積\(S_{BB}\)を用いて
\(
F=PS_{BB}
\)
と書けるはずです。すなわち、内圧がこの力に達した時にBB弾が射出されるとして考え、その内部圧力に達した際にBB弾の運動方程式を解き始める、という方法で可能なはずです。

参考文献

[1] 弾道計算の理論 -シキノート https://slpr.sakura.ne.jp/qp/bullet-course/
[2] BB弾の弾道学 -ハイパー道楽 http://www.hyperdouraku.com/colum/ballistics/index.html

極値を求めるニュートン法

ニュートン法に関するお話です。
”ニュートン法”と呼ばれる方法は
2種類(1.初期値近傍の極値を求めるニュートン法、2. ゼロ点を求めるためのニュートン法
があるようですが、ここでは1. 初期値近傍の極値を求める方のニュートン法についてのお話です。

まとめ

ベクトル\(\mathbf{x}=(x_1, x_2, \cdots, x_N)\)で指定される実数値関数\(f(\mathbf{x})\)の極値を求める事を考えます。
ここで、初期値近傍の点\(\mathbf{x}=\mathbf{x}^{(0)}\)が分かっている時、その近傍にある\(f(\mathbf{x})\)の極小値を与える\(\mathbf{x}=\mathbf{a}\)は適当な反復計算を行い、
\begin{equation}
\mathbf{a}=\lim_{k\to\infty}\mathbf{x}^{(k)}
\end{equation}
で与えられます。極値を求めるニュートン法を用いる場合、\(\mathbf{x}^{(k)}\)は漸化式として与えられ、
\begin{align}
\mathbf{x}^{(k+1)}&=\mathbf{x}^{(k)}+\mathbf{d}^{(k)} \\
\mathbf{d}^{(k)}&=-\left[H f\bigl(\mathbf{x}^{(k)}\bigr)\right]^{-1}\cdot \left[\nabla f\bigl(\mathbf{x}^{(k)}\bigr)\right]
\end{align}
として逐次的に与えられます。ここで、\(H f(\mathbf{x}), \nabla f\bigl(\mathbf{x}\bigr)\)は関数\(f(\mathbf{x})\)のヘッセ行列、傾き(Gradient)を表し、
\begin{equation}
\hat{H}f=
\begin{pmatrix}
\frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} &\cdots &\frac{\partial^2 f}{\partial x_1 \partial x_n} \\
\frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} &\cdots & \frac{\partial^2 f}{\partial x_2 \partial x_n} \\
\vdots & \vdots &\ddots & \vdots \\
\frac{\partial^2 f}{\partial x_n \partial x_1} & \frac{\partial^2 f}{\partial x_n \partial x_2} &\cdots &\frac{\partial^2 f}{\partial x_n^2} \\
\end{pmatrix}
,~~\nabla f=
\begin{pmatrix}
\frac{\partial f}{\partial x_1} \\
\frac{\partial f}{\partial x_2} \\
\vdots \\
\frac{\partial f}{\partial x_n}
\end{pmatrix}
\end{equation}
で与えられます。

特に\(f\)が\(x, y\)の2変数関数\(f(x,y)\)であるならば、
\begin{align}
\begin{pmatrix}
x^{(k+1)} \\
y^{(k+1)}
\end{pmatrix}
&=
\begin{pmatrix}
x^{(k)} \\
y^{(k)}
\end{pmatrix}
-\frac{1}{f_{xx}f_{yy}-f_{xy}f_{yx}}
\begin{pmatrix}
f_{yy} & -f_{xy} \\
-f_{yx} & f_{xx}
\end{pmatrix}
\begin{pmatrix}
f_x \\
f_y
\end{pmatrix}
\end{align}
ここで、
\begin{equation}
\hat{H}f
=
\begin{pmatrix}
f_{xx} & f_{xy} \\
f_{yx} & f_{yy}
\end{pmatrix}
,~~
\nabla f=
\begin{pmatrix}
f_x \\
f_y
\end{pmatrix}
\end{equation}
です。ここで、\(f_p\equiv\frac{\partial f}{\partial p}, ~f_{pq}\equiv\frac{\partial^2 f}{\partial p\partial q}, (p, q=x,y)\)と書きました。
微分を差分で近似するとします。\(x, y\)方向の刻み幅を\(h=h_x=h_y\)とすればそれぞれの偏微分は点\((x,y)\)近傍の7点を用いて
\begin{align}
f_x &= \frac{f(x+h, y) – f(x-h, y)}{2h} +O(h^2)\\
f_y &= \frac{f(x, y+h) – f(x, y-h)}{2h} +O(h^2)\\
f_{xx} &= \frac{f(x+h, y) – 2f(x, y)+ f(x-h, y)}{h^2} +O(h^2)\\
f_{xy} &=f_{yx}= \frac{1}{2h^2}\bigl[f(x+h, y) + f(x-h, y) + f(x, y+h) \\
&\hspace{1em}+ f(x, y-h)-2f(x,y)-f(x+h,y+h)-f(x-h,y-h)\bigr]+O(h^2) \\
f_{yy} &= \frac{f(x, y+h) – 2f(x, y)+ f(x, y-h)}{h^2} +O(h^2)
\end{align}
と書けます[2]。

導出


詳しくは述べませんが、1次元の問題における極小値を求めるニュートン法くらいは導いておきましょう。
多次元の場合は拡張すれば良く、考え方は変わりません。多次元では他の文献、例えば[1]を参考にして下さい。

関数\(f(x)\)を\(x=x^{(k)}\)周りでテイラー展開を行えば、
\begin{equation}
f(x)=f\bigl(x^{(k)}\bigr)+f’\bigl(x^{(k)}\bigr)\bigl(x-x^{(k)}\bigr)+\frac{1}{2}f^{\prime\prime}\bigl(x^{(k)}\bigr)\bigl(x-x^{(k)}\bigr)^2+O(\Delta^3)
\end{equation}
となります。ここで、\(\Delta=|x-x^{(k)}|\)とおきました。
今、欲しい解\(a\)は関数\(f(x)\)の極値なので、
\(\frac{df}{dx}=0\)を満たす\(x\)を探せばよいのです。よって、両辺を\(x\)で微分して、
\begin{align}
\frac{df(x)}{dx}=f’\bigl(x^{(k)}\bigr)+f^{\prime\prime}\bigl(x^{(k)}\bigr)\bigl(x-x^{(k)}\bigr)+O(\Delta^2)
\end{align}
を得ます。これが\(0\)になる\(x\)が解\(a\)なので、
\begin{equation}
\left. \frac{df(x)}{dx}\right|_{x=a}=f’\bigl(x^{(k)}\bigr)+f^{\prime\prime}\bigl(x^{(k)}\bigr)\bigl(a-x^{(k)}\bigr)+O(\Delta^2) = 0
\end{equation}
となります。今、\(O(\Delta^2)\)を無視すれば得られる解\(a\)は近似解\(a\approx x^{(k+1)}\)となるので、
\begin{equation}
f’\bigl(x^{(k)}\bigr)+f^{\prime\prime}\bigl(x^{(k)}\bigr)\bigl(x^{(k+1)}-x^{(k)}\bigr) = 0 \\
\end{equation}
を解いて
\begin{equation}
x^{(k+1)} = x^{(k)}-\frac{f’\bigl(x^{(k)}\bigr)}{f^{\prime\prime}\bigl(x^{(k)}\bigr)}
\end{equation}
となります。…ここまで来てしまうとニュートン・ラフソン法と同じですね。関数値がゼロになる点を探すのか、導関数値がゼロになる点を探すのかの違いなだけです。
多次元だと行列になってしまいシンプルではない形になり難しく見えますが、1次元では非常に単純です。

問題は1階(傾き)、2階導関数(ヘッセ行列)をどのように求めるかです。厳密な値が無い場合は近似的にしか得るしかありません。基本的にこれらの計算量は膨大になります。
最も簡単な方法は微分を差分に置き換えてしまうことですが、差分は余り良い方法とは言えません。

厳密なヘッセ行列は、例えば自動微分の考えで得ることが出来ます。例えば、Hyper-dual numbersによる二階偏微分の計算を参考にして得ることが出来ます。

また、引数が複素数で計算結果が実数で返ってくる関数\(f(z)\)の最小値を求めたいのであれば、\(z=x+iy\)として、\(f(x,y)\)の二変数関数としてその極致を探せば良いです。

2次元のプログラム


ここでは2次元の問題に限り、傾き、ヘッセ行列を差分で近似した場合のプログラムを記述します。

下のプログラムは関数
\begin{equation}
f(x,y)=x e^{-(x^2+y^2)/2}
\end{equation}
の、点\((x_0, y_0)=(-1.2, -0.3)\)近傍にある極値を求めるプログラムです。
グラフで示せば、以下の通りになります。
この関数の極値は黒点、初期値は赤点で示しています。

以下に実行例を示します。

$ gfortran main.f90
$ ./a.out
  -1.2000000000000000      -0.29999999999999999      -0.55840071716917605    
  -1.0000000016663813        3.4150669573473641E-013 -0.60653065971263342      4
$

1行目1,2,3列目はそれぞれ\(x,y\)の初期値, 関数\(f(x,y)\)の値を示しています。
また、2行目1,2,3,4列目は求まった\(x,y\)の値, 関数\(f(x,y)\)の値、そして収束に至るまでに行った繰り返しの回数を示しています。

参考文献


[1]塩浦昭義, 数理計画法 第13回
[2]Abramowitz and Stegun, “Handbook of Mathematical Functions”, http://people.math.sfu.ca/~cbm/aands/page_883.htm, http://people.math.sfu.ca/~cbm/aands/page_884.htm

反復計算における収束判定について

収束判定


漸化式
\begin{equation}
x^{(k+1)}=x^{(k)}+\Delta^{(k)}
\end{equation}
に従って、\(\lim_{k\to\infty}x^{(k)}=a\)が適当な定数に収束する問題を考えます。
絶対誤差\(\varepsilon_A\)、相対誤差\(\varepsilon_R\)とすると、数値計算を終了する時は、以下の不等式が満たされる時にすると良いです。
\begin{equation}
|x^{(k+1)}-x^{(k)}|\lt \varepsilon_A + \varepsilon_R(|x^{(k)}|+|x^{(k+1)}|)
\end{equation}

判定式の意味


上記判定式の意味を考えましょう。上の判定式は絶対誤差による評価と相対誤差による評価をまとめて評価する式になっています。

相対誤差が重要な場合

本当の解が\(a=100\)であるとします。その時、\(x^{(k)}=100.3, x^{(k+1)}=100.2\)であったと仮定すると、大雑把に計算して
\begin{align}
|x^{(k+1)}-x^{(k)}| &\lt \varepsilon_A + \varepsilon_R(|x^{(k)}|+|x^{(k+1)}|) \\
0.1 &\lt \varepsilon_A + \varepsilon_R\times 100
\end{align}
となります。右辺は本来2倍ですが、2倍はそんなに重要ではないので無視しています。もし仮に、\(\varepsilon\equiv\varepsilon_A=\varepsilon_R\)としていれば、
\begin{align}
0.1 &\lt 100\varepsilon_R
\end{align}
と近似できます。以上から、本当の解が1より十分大きい場合、相対誤差で収束を評価する式になっていると言えます。
つまり上記の問題設定の場合、\(\varepsilon_R\gt 10^{-3}\)にしてあれば”収束した”として判定します。

絶対誤差が重要な場合

本当の解が\(a=0.01000\)であるとします。その時、\(x^{(k)}=0.01002, x^{(k+1)}=0.01001\)であったと仮定すると、大雑把に計算して
\begin{align}
|x^{(k+1)}-x^{(k)}| &\lt \varepsilon_A + \varepsilon_R(|x^{(k)}|+|x^{(k+1)}|) \\
0.00001 &\lt \varepsilon_A + \varepsilon_R\times 0.01
\end{align}
となります。右辺は本来2倍ですが、2倍はそんなに重要ではないので無視しています。もし仮に、\(\varepsilon\equiv\varepsilon_A=\varepsilon_R\)としていれば、
\begin{align}
0.00001 &\lt 0.01\varepsilon_A
\end{align}
と近似できます。以上から、本当の解が1より十分小さい場合、絶対誤差で収束を評価する式になっていると言えます。
つまり上記の問題設定の場合、\(\varepsilon_A\gt 10^{-3}\)にしてあれば”収束した”として判定します。

参考


[1]杉山耕一朗, OBOROの収束判定条件の設定

動かない壁に対する束縛運動と反射

動かない壁に対する束縛運動と反射を考えます。
例えば、初め跳ねてた運動が、壁に沿って動く運動に変化する、という状況です。
あんまり見たことが無いので、面白そうだと思いました。

束縛されている時と反発するときは、式(1),(2)によってあらわされます。それらは
壁に束縛されている場合

壁と反発する場合

です。壁と反発する場合、反発後の速度は式(3)に沿って動きます。

※式(1), (2)では壁が時間依存しても良いように定式化しています。この定式化は恐らく正しいです。また、本稿に載せているプログラムも壁が時間依存しても良いように作成していますが、動く壁の場合、プログラムではうまく計算が出来ません。

ここで、\(e\)は反発係数、\(\mathbf{n}\)は壁の法線ベクトルであり、

\(\nabla\)はナブラ演算子であり

と与えられます。また、\(\hat{H}f\)は関数\(f\)のヘッセ行列であり,

と与えられます。

定式化や数値計算手法の詳細は以下のページを参照してください。

壁との反発と束縛運動の定式化

質点と壁との反発を表す運動方程式
束縛条件下の運動 – 自由度がうまく落とせない運動

数値計算手法

ルンゲ=クッタ法の説明と刻み幅制御
Hyper-dual numbersによる二階偏微分の計算
ゼロ点を探す(二分法、挟み撃ち法、Anderson-Björk法、Brent法、Newton法、Steffensen法)


解放・束縛判定


ここでいう”解放”とは、束縛されていなく、反射を繰り返している状態を表します。また、”束縛”は壁に沿って運動している状態です。

前提として、壁を通り抜けることは無いと考えます。

すなわち、時刻\(t=t_0\)で位置\(\mathbf{r}=\mathbf{r}_0\)の時、任意の時刻\(t\)について

が満たされるとします。
もし、\(f(\mathbf{r}_0,t_0) =0\)ならば判断がつかないため、プラスかマイナスはこちらから与えます。

解放→解放判定、解法→束縛判定


解放状態から壁によって単に反発する場合

関数\(f(\mathbf{r},t)\)の符号が変化した時、壁との反発を考えます。
壁の法線方向の速度が十分に大きい場合、壁と反発し、そうでない場合、壁に束縛されると考えます。
その前提の元、

を満たす\(t=t_i, \mathbf{r}=\mathbf{r}_i\)を見つけます。
その後、式(6)に従い、速度を変化させます。

数値計算的には、衝突直前の時刻を採用します。すなわち、
ゼロ点を探す際にある範囲\(t_a \leq t_i \leq t_b\)で挟み込んでいくのですが、\(t_b\)は壁を越えてしまうので採用しません。

解放状態から壁沿いに束縛される場合

もし、壁の法線方向の速度が十分に小さい場合(ある閾値を下回った場合)、壁に束縛されると考えます。
この時、壁の法線方向の速度はゼロに変更します。
すなわち、速度は時刻\(t=t_i\)で

を持ちますが、\(v_{\parallel}=0\)にしてから、束縛運動に移行するということです。
これは、\(\mathbf{v}\)と壁の法線方向のベクトル\(\mathbf{e}_n\)の内積を取ることで得られます。
また、束縛された瞬間(\(t=t_i\))の束縛力\(C(t_i)\)を計算し、その符号を記録しておきます。
束縛力\(C(t)\)は、

です。
この符号が変化した瞬間が壁からの束縛が無くなる時なので、そのために記録します。

束縛→解放判定

束縛力\(C(t)\)の符号と\(C(t_i)\)の符号が変わるまで、式(1)に従い、時間発展させます。
すなわち、

を満たす\(t=t’, \mathbf{r}=\mathbf{r}’\)を見つけます。
その後、式(2)に従い運動します。
式(2)の運動では束縛力は働かないので、符号は自然と初期条件の符号と同一になる(※)。

※この条件はあまり良くありません。この判別方法のせいで、壁が時間依存している場合、束縛力が働いていない一瞬で質点が壁を超えてしまいます。プログラム自体は壁は時間依存しても良いことになっていますが、この条件分岐は上手く動きません。以下に示す計算プログラムは、質点が束縛されている場合に壁が時間依存しなければ正しいです。

プログラム


プログラムは以下のリンク先においておきます。

https://slpr.sakura.ne.jp/qp/supplement_data/lag_ver1.0.tar.gz

適当なディレクトリで展開し、lag_ver1.0というディレクトリに移動してから以下のコマンドで実行できます。

 $ sh cm.sh
 $ ./a.out
&INPUT
 MASS=  1.0000000000000000     ,
 G=  1.0000000000000000     ,
 TA=  0.0000000000000000     ,
 TB=  20.000000000000000     ,
 NT=        101,
 ELS= 0.59999999999999998     ,
 RX0= -1.0000000000000000     ,
 RY0=  1.0000000000000000     ,
 VX0=  3.0000000000000000     ,
 VY0=  0.0000000000000000     ,
 RKTOL=  1.0000000000000000E-008,
 ZRTOL=  1.0000000000000000E-008,
 TRTOL=  1.0000000000000000E-008,
 REGION=          1,
 /
           0
 $ gnuplot

    G N U P L O T
    Version 4.6 patchlevel 4    last modified 2013-10-02
    Build System: Linux x86_64

    Copyright (C) 1986-1993, 1998, 2004, 2007-2013
    Thomas Williams, Colin Kelley and many others

    gnuplot home:     http://www.gnuplot.info
    faq, bugs, etc:   type "help FAQ"
    immediate help:   type "help"  (plot window: hit 'h')

Terminal type set to 'wxt'
gnuplot> load "anime.plt"
gnuplot>

動かすと以下のような動画が描画されます。

デフォルトでは
ポテンシャル\(V=mgy\)(サブルーチンfp2d)
壁\(f(x,y)=x^2+y^2-4\)(サブルーチンfw2d)
に書かれています。
摩擦、空気抵抗は入っていません。唯一、反発係数(els)だけがinputファイルの中に書かれています。

初期条件の

rx0  = -2d0,  ! Initial condition
ry0  =  2d0,  !     position and velocity
vx0  =  1d0,  !  \mathbf{r} = (rx, ry)
vy0  =  0d0,  !  \mathbf{v} = (vx, vy)

だけを上の通り変更すると以下のような振る舞いをします。

確かめ


確かめを行います。
重力\(g\)の下で、質量\(m\)の物体が、半径\(r\)の円形の壁の内側に沿ってボールが進み、壁からの抗力が無くなり、壁から離れる状況を考えます(下の図を参照)。

円形の壁に沿っている時、垂直抗力\(N\)は

と書けます。エネルギー保存則より、

が成り立っています。ここで、\(v_0\equiv v(t=0)\)と置きました。
垂直抗力\(N\)がゼロになる点が壁から離れる条件ですので、\(v_0\)を用いて

と書けます。初速度が分かっている時、壁から離れる角度は

です。もしくは、壁から離れる角度が分かっている場合、初速度は

と与えられます。重力加速度, 半径を\(g=1, r=2\)とし、\(y=1\)と壁との交点、すなわち\(\theta=2\pi/3\)の場合、初速度は

です。実際に本稿のプログラムを動かしますと、\(y=1\)でちょうど離れていることが確認できます。

ここで、青線は壁に沿って動いて運動している状態であり、赤線は壁から離れている運動している状態です。

剛体に対する散乱

ポテンシャルを無くし(g=0)、弾性散乱(els=1)を考えると
古典的な散乱問題的なものもできます。

Cubic補間

多項式補間です。

Cubic補間 (1次元)


4点を厳密に通る3次多項式によって関数を補間します。
概要は以下の図の通りです。
補間したい点よりも小さいデータ点を2点、大きいデータ点を2点使って補間します。
いわば小さい区間で区切ったラグランジュ補間です。
補間は,既知の係数\(a_{i}\)を用いて関数
\(\displaystyle
g(x)=\sum_{i=0}^{3}a_{i}x^i
\)

誤差は\(O(h^4)\)です。

Fortranプログラムはこちら。

program main
  implicit none
  integer::i,N,M
  double precision::x,f,df,df2,s
  double precision,allocatable::xdata(:),fdata(:)
  double precision::pi=dacos(-1d0)
  double precision,external::cubic
 
  N=30
  M=100
  allocate(xdata(0:N),fdata(0:N))
  xdata=0d0
  fdata=0d0
 
  do i=0,N
     xdata(i)=dble(i)*0.1d0*pi
     fdata(i)=sin(xdata(i))
     write(10,*)xdata(i),fdata(i)
  enddo
 
  ! Cubic-spline interpolation given position as point
  do i=0,M
     x=dble(i)*0.03d0*pi-1d0
     write(11,'(2e20.7e2)')x,cubic(x,N,xdata,fdata)
  enddo
 
  stop
end program main

double precision function cubic(x,N,x0,f0)
  implicit none
  integer,intent(in)::N
  double precision,intent(in)::x,f0(0:N),x0(0:N)
 
  integer::i,i0,i1,i2,i3
  double precision::tx  
  double precision::a,b,c,d,p,q,r,s,t,u

  tx = x-x0(0)
  i1 = 0
  do i=1,N-2
     tx = x-x0(i)
     if(tx.gt.0d0)then
        i1 = i
     else
        exit
     endif
  enddo
  if(i1.eq.0)i1=1
  i0=i1-1
  i2=i1+1
  i3=i1+2

  a = x-x0(i0)
  b = x-x0(i1)
  c = x-x0(i2)
  d = x-x0(i3)  
 
  p = x0(i1)-x0(i0)
  q = x0(i2)-x0(i1)
  r = x0(i3)-x0(i2)

  s = x0(i2)-x0(i0)
  t = x0(i3)-x0(i1)

  u = x0(i3)-x0(i0)

  cubic = a*c*( d*f0(i1)/(p*q) + b*f0(i3)/(u*r) )/t &
        - b*d*( c*f0(i0)/(p*u) + a*f0(i2)/(q*r) )/s
 
  return
end function cubic

Bicubic補間 (2次元)


16点を厳密に通る多項式によって関数を補間します。
補間したい点を囲むように存在する16点を用いて補間します。
補間は,既知の係数\(a_{i,j}\)を用いて関数
\(\displaystyle
g(x,y)=\sum_{i=0}^{3}\sum_{j=0}^{3}a_{i,j}x^i y^j
\)

によって行います[1]。誤差は\(O(h^4)\)です。

Fortranプログラムはこちら。

program main
  implicit none
  integer::i,j
  integer::Nx,Ny
  double precision::x,y
  double precision::xa,xb,ya,yb
  double precision,allocatable::x0(:),y0(:),z0(:,:)
  double precision,external::f,ip16

  Nx=20
  Ny=20
  allocate(x0(0:Nx), y0(0:Ny), z0(0:Nx,0:Ny))

  xa=-3d0
  xb= 3d0
  ya=-3d0
  yb= 3d0
   
  ! generate references points
  do i=0,Nx
     x0(i)=dble(i)*(xb-xa)/Nx+xa
  enddo
  do j=0,Ny
     y0(j)=dble(j)*(yb-ya)/Ny+ya
  enddo
  do i=0,Nx
     do j=0,Ny
        z0(i,j)=f(x0(i),y0(j))
        write(10,*)x0(i),y0(j),z0(i,j)
     enddo
     write(10,*)
  enddo
 
  ! Return interpolated results
  do i=0,50
     do j=0,50
        x=i*0.03d0-1d0
        y=j*0.03d0-1d0
        write(11,*)x,y,ip16(x,y,Nx,Ny,x0,y0,z0)
     enddo
     write(11,*)
  enddo

  stop
end program main

function f(x,y)
  implicit none
  double precision::f
  double precision,intent(in)::x,y
 
  f=x*y*exp(-x*x-y*y)
 
  return
end function f

double precision function ip16(x,y,Nx,Ny,x0,y0,z0)
  implicit none
  integer,intent(in)::Nx,Ny
  double precision,intent(in)::x,y,x0(0:Nx),y0(0:Ny),z0(0:Nx,0:Ny)
  !
  ! Bicubic interpolation
  !   x0,y0 are sorted by ascending order,
  !   suppose x0 & y0 are equal interval.
  !   z0(x0,y0)
  !
  integer::i,j
  integer::i0,i1,i2,i3
  integer::j0,j1,j2,j3
  double precision::tx,u,u2,u3
  double precision::ty,v,v2,v3
  double precision::a00,a01,a02,a03,a10,a11,a12,a13
  double precision::a20,a21,a22,a23,a30,a31,a32,a33
  double precision::p(0:3,0:3)
 
  if(Nx.le.3.or.Ny.le.3)then
     write(6,*)" *** Error at ip16"
     stop
  endif
 
  tx = x-x0(0)
  i1 = 0
  do i=1,Nx-2
     tx = x-x0(i)
     if(tx.gt.0d0)then
        i1 = i
     else
        exit
     endif
  enddo
  if(i1.eq.0)i1=1
  i0=i1-1
  i2=i1+1
  i3=i1+2

  ty = y-y0(0)
  j1 = 0
  do j=1,Ny-2
     ty = y-y0(j)
     if(ty.gt.0d0)then
        j1 = j
     else
        exit
     endif
  enddo
  if(j1.eq.0)j1=1
  j0=j1-1
  j2=j1+1
  j3=j1+2

 
  p(0,0) = z0(i0,j0)
  p(0,1) = z0(i0,j1)
  p(0,2) = z0(i0,j2)
  p(0,3) = z0(i0,j3)
  p(1,0) = z0(i1,j0)
  p(1,1) = z0(i1,j1)
  p(1,2) = z0(i1,j2)
  p(1,3) = z0(i1,j3)
  p(2,0) = z0(i2,j0)
  p(2,1) = z0(i2,j1)
  p(2,2) = z0(i2,j2)
  p(2,3) = z0(i2,j3)
  p(3,0) = z0(i3,j0)
  p(3,1) = z0(i3,j1)
  p(3,2) = z0(i3,j2)
  p(3,3) = z0(i3,j3)
   
  a00 = p(1,1)
  a01 = -0.5d0*p(1,0) + 0.5d0*p(1,2)
  a02 = p(1,0) - 2.5d0*p(1,1) + 2*p(1,2) - 0.5d0*p(1,3)
  a03 = -0.50d0*p(1,0) + 1.50d0*p(1,1) - 1.50d0*p(1,2) + 0.5d0*p(1,3)
  a10 = -0.50d0*p(0,1) + 0.50d0*p(2,1)
  a11 =  0.25d0*p(0,0) - 0.25d0*p(0,2) - 0.25d0*p(2,0) + 0.25d0*p(2,2)
  a12 = -0.50d0*p(0,0) + 1.25d0*p(0,1) - p(0,2) + 0.25d0*p(0,3) &
       + 0.50d0*p(2,0) - 1.25d0*p(2,1) + p(2,2) - 0.25d0*p(2,3)
  a13 =  0.25d0*p(0,0) - 0.75d0*p(0,1) + 0.75d0*p(0,2) - 0.25d0*p(0,3) &
       - 0.25d0*p(2,0) + 0.75d0*p(2,1) - 0.75d0*p(2,2) + 0.25d0*p(2,3)
  a20 = p(0,1) - 2.5d0*p(1,1) + 2.d0*p(2,1) - 0.5d0*p(3,1)
  a21 = -0.50d0*p(0,0) + 0.5d0*p(0,2) + 1.25d0*p(1,0) - 1.25d0*p(1,2) &
       - p(2,0) + p(2,2) + 0.25d0*p(3,0) - 0.25d0*p(3,2)
  a22 =  p(0,0) - 2.5d0*p(0,1) + 2.0d0*p(0,2) - 0.50d0*p(0,3) &
       - 2.5d0*p(1,0) + 6.25d0*p(1,1) - 5d0*p(1,2) + 1.25d0*p(1,3) &
       + 2.0d0*p(2,0) -  5.0d0*p(2,1) + 4d0*p(2,2) - p(2,3) &
       - 0.50d0*p(3,0) + 1.25d0*p(3,1) - p(3,2) + 0.25d0*p(3,3)
  a23 = -0.50d0*p(0,0) + 1.50d0*p(0,1) - 1.50d0*p(0,2) + 0.5d0*p(0,3) &
       + 1.25d0*p(1,0) - 3.75d0*p(1,1) + 3.75d0*p(1,2) - 1.25d0*p(1,3) &
       - p(2,0) + 3*p(2,1) - 3*p(2,2) + p(2,3) &
       + 0.25d0*p(3,0) - 0.75d0*p(3,1) + 0.75d0*p(3,2) - 0.25d0*p(3,3)
  a30 = -0.50d0*p(0,1) + 1.50d0*p(1,1) - 1.50d0*p(2,1) + 0.5d0*p(3,1)
  a31 =  0.25d0*p(0,0) - 0.25d0*p(0,2) - 0.75d0*p(1,0) + 0.75d0*p(1,2) &
       + 0.75d0*p(2,0) - 0.75d0*p(2,2) - 0.25d0*p(3,0) + 0.25d0*p(3,2)
  a32 = -0.50d0*p(0,0) + 1.25d0*p(0,1) - p(0,2) + 0.25d0*p(0,3) &
       + 1.50d0*p(1,0) - 3.75d0*p(1,1) + 3*p(1,2) - 0.75d0*p(1,3) &
       - 1.50d0*p(2,0) + 3.75d0*p(2,1) - 3*p(2,2) + 0.75d0*p(2,3) &
       + 0.50d0*p(3,0) - 1.25d0*p(3,1) + p(3,2) - 0.25d0*p(3,3)
  a33 =  0.25d0*p(0,0) - 0.75d0*p(0,1) + 0.75d0*p(0,2) - 0.25d0*p(0,3) &
       - 0.75d0*p(1,0) + 2.25d0*p(1,1) - 2.25d0*p(1,2) + 0.75d0*p(1,3) &
       + 0.75d0*p(2,0) - 2.25d0*p(2,1) + 2.25d0*p(2,2) - 0.75d0*p(2,3) &
       - 0.25d0*p(3,0) + 0.75d0*p(3,1) - 0.75d0*p(3,2) + 0.25d0*p(3,3)
 
  ! Parallel translation u,v [0:1]
  u  = (x-x0(i1))/(x0(i2)-x0(i1))
  v  = (y-y0(j1))/(y0(j2)-y0(j1))
  u2 = u*u
  u3 = u2*u
  v2 = v*v
  v3 = v2*v
  ip16=  (a00 + a01 * v + a02 * v2 + a03 * v3)      &
       + (a10 + a11 * v + a12 * v2 + a13 * v3) * u  &
       + (a20 + a21 * v + a22 * v2 + a23 * v3) * u2 &
       + (a30 + a31 * v + a32 * v2 + a33 * v3) * u3

  return
end function ip16

補間 (2次元)


Bicubic補間の低次元バージョンです。
補間したい点を囲むように存在する4点を用いて補間します。
一応Four point Formulaと名づけられています[2]。
補間は,既知の係数\(a_{i,j}\)を用いて関数
\(\displaystyle
g(x,y)=\sum_{i=0}^{1}\sum_{j=0}^{1}a_{i,j}x^i y^j
\)

によって行います[2]。誤差は\(O(h^2)\)です。

function ip4(x,y,Nx,Ny,x0,y0,z0)
  implicit none
  integer,intent(in)::Nx,Ny
  double precision,intent(in)::x,y,x0(1:Nx),y0(1:Ny),z0(1:Nx,1:Ny)
  double precision::ip4
  !
  ! 2D Interpolation, 4-point formula,
  !           xy equal distance grid
  !
  integer::i,ix0,ix1,iy0,iy1
  double precision::tx,ty,dx,dy,p,q

  tx=1d100
  ix0=1
  do i=2,Nx
     if(abs(x-x0(i)).le.tx)then
        tx=abs(x-x0(i))
        ix0=i
     endif    
  enddo
  tx=1d100
  ix1=1
  do i=2,Nx
     if(i.ne.ix0)then
        if(abs(x-x0(i)).le.tx)then
           tx=abs(x-x0(i))
           ix1=i
        endif
     endif
  enddo

  ty=1d100
  iy0=1
  do i=2,Ny
     if(abs(y-y0(i)).le.ty)then
        ty=abs(y-y0(i))
        iy0=i
     endif    
  enddo
  ty=1d100
  iy1=1
  do i=2,Ny
     if(i.ne.iy0)then
        if(abs(y-y0(i)).le.ty)then
           ty=abs(y-y0(i))
           iy1=i
        endif
     endif
  enddo
 
  dx=x0(ix1)-x0(ix0)
  dy=y0(iy1)-y0(iy0)
  p=(x-x0(ix0))/dx
  q=(y-y0(iy0))/dy
 
  ip4=(1d0-p)*(1d0-q)*z0(ix0,iy0) &
     +p*(1d0-q)*z0(ix1,iy0) &
     +(1d0-p)*q*z0(ix0,iy1) &
     +p*q*z0(ix1,iy1)
   
  return
end function ip4

参考文献


[1]Cubic interpolation
[2]Abramowitz and Stegun.
Handbook of Mathematical Functions.

LU分解による連立一次方程式の解法

数値計算ライブラリLapackを使って線形連立一次方程式
\(
\begin{equation}
A\mathbf{x}=\mathbf{b}
\end{equation}
\)

をLU分解を用いて数値的に解くプログラムを載せます。

詳しくは説明しませんが、逆行列を用いて解く方法と比べ、
行列が疎である時(行列要素にゼロが多い時)に比較的計算量が減らせます。
また、逆行列を求める際に生じる桁落ちの問題を回避することが出来ます。

解法


連立一次方程式
\(
\begin{equation}
A\mathbf{x}=\mathbf{b}
\end{equation}
\)

を行列\(A\)のLU分解を利用して解きます。
ここで、未知なのはベクトル\(\mathbf{x}\)、既知なのは行列\(A\)とベクトル\(\mathbf{b}\)です。

\(A\)がLU分解されているとします。
すなわち、行列\(A\)に適当な操作を施して、

\(
\begin{align}
A=LU
\end{align}
\)

と書けているとします。ここで、行列\(L, U\)は下三角行列、上三角行列であり、
\(
\begin{align}
L=
\left(
\begin{array}{*4{>{\displaystyle}c}}
1&0&\cdots &0 \\[1em]
l_{2,1}&1&\cdots &0 \\[0.5em]
\vdots&\vdots&\ddots&\vdots \\[0.3em]
l_{s,1}&l_{s,2}&\cdots &1
\end{array}
\right) ,~~~
U=
\left(
\begin{array}{*4{>{\displaystyle}c}}
u_{1,1}&u_{1,2}&\cdots &u_{1,s} \\[1em]
0&u_{2,2}&\cdots &u_{2,s} \\[0.5em]
\vdots&\vdots&\ddots&\vdots \\[0.3em]
0 & 0 &\cdots &u_{s,s}
\end{array}
\right)
\end{align}
\)

という行列です。
もし、行列\(A\)がLU分解されていれば、

\(
\begin{align}
& A\mathbf{x}=\mathbf{b} \\
& LU\mathbf{x} = \mathbf{b} \\
& L\mathbf{w} = \mathbf{b}~~~…(*)
\end{align}
\)

と書けます。ここで、\(\mathbf{w} \equiv U\mathbf{x}\)と置きました。
まず、式(*)を解いて、\(\mathbf{w}\)を求めます。これは、下三角行列に対する問題なので、前進代入を用いて簡単に解くことが出来ます。

続いて、
\(
U\mathbf{x}=\mathbf{w}
\)

を後退代入を利用して解\(\mathbf{x}\)を求めます。

以上がLU分解を利用して連立一次方程式を解く方法です。
数値計算のルーチンは大きく2つのステップに分かれており、

  1. 行列\(A\)のLU分解を求める(ルーチン dgetrf)
  2. LU分解された行列\(A\)を用いて連立一次方程式を解く(ルーチン dgetrs)

というステップです。
これは、例えば問題
\(
\begin{equation}
A\mathbf{x}=\mathbf{b}
\end{equation}
\)


\(
\begin{equation}
A\mathbf{x}=\mathbf{b’}
\end{equation}
\)

の右辺だけが変わる複数の問題を解きたい場合なんかに便利です。この問題の場合、行列\(A\)は変わらないため、LU分解した結果を両方の問題に流用することが出来るのです。

プログラム


プログラムは以下の通りです。
いつもと同じように、ワーク配列を減らす為だけに導入したルーチンを挟んでいます。

program main
  implicit none
  integer::i,N
  double precision,allocatable::a(:,:),b(:),x(:)
  integer,allocatable::ipiv(:)
 
  N=3
  allocate(a(1:N,1:N))
  allocate(ipiv(1:N),b(1:N),x(1:N))
  a(1,1:N)=(/1d0,3d0,5d0/)
  a(2,1:N)=(/0d0,3d0,1d0/)
  a(3,1:N)=(/6d0,2d0,5d0/)
  b(1:N)=(/33d0,10d0,66d0/)

  do i=1,N
     write(6,'(3f10.5)')a(i,1:N)
  enddo

  call LUfact(N,a,ipiv)
  call axbsolve(N,a,ipiv,b,x)

  write(6,*)"---- solution ----"
  do i=1,N
     write(6,'(2f10.5)')x(i)
  enddo
 
  stop
end program main

subroutine LUfact(N,LU,ipiv)
  implicit none
  integer,intent(in)::N
  double precision,intent(inout)::LU(1:N,1:N)
  integer,intent(out)::ipiv(1:N)  

  ! LU factorization for lapack
  ! Overwrite matrix LU by factorized LU matrix
 
  integer::m,lda,info

  m=N
  lda=N
  call dgetrf(m,N,LU,lda,ipiv,info)
  if(info.ne.0)then
     write(6,*)"**error at LUlapack"
     stop
  endif
 
  return
end subroutine LUfact

subroutine axbsolve(N,LU,ipiv,b,x)
  implicit none
  integer,intent(in)::N
  integer,intent(in)::ipiv(1:N)
  double precision,intent(in)::LU(1:N,1:N),b(1:N)
  double precision,intent(out)::x(1:N)
  !
  ! Solve simultaneous equations by Lapack
  !
  double precision,allocatable::bl(:,:)
  integer::nrhs,lda,ldb,info
  character(1)::trans  
 
  nrhs=1
  trans='N'
  allocate(bl(1:N,1:nrhs))  
  bl(1:N,1)=b(1:N)
  lda=N
  ldb=N
  call dgetrs(trans,N,nrhs,LU,lda,ipiv,bl,ldb,info)
  if(info.ne.0)then
     write(6,*)"**error at LUlapack"
     stop
  endif
 
  x(1:N)=bl(1:N,1)  

  return
end subroutine axbsolve

上のプログラムをlapackと一緒にコンパイルして動かすと、端末上で

> ./a.out
   1.00000   3.00000   5.00000
   0.00000   3.00000   1.00000
   6.00000   2.00000   5.00000
 -------------
   7.00000
   2.00000
   4.00000
>

という結果が得られるかと思います。
これは、連立一次方程式
\(
\begin{align}
\left(
\begin{array}{ccc}
1& 3& 5 \\
0& 3& 1 \\
6& 2& 5
\end{array}
\right)
\left(
\begin{array}{ccc}
x \\
y \\
z
\end{array}
\right)
=
\left(
\begin{array}{ccc}
33 \\
10 \\
66
\end{array}
\right)
\end{align}
\)
を解いた結果です。
解析解(\(x=7, y=2, z=4\))は、例えばwolframを利用して求められます。
{{1,3,5},{0,3,1},{6,2,5}}.{x,y,z}=={33,10,66} –wolfram alpha

質点と壁との反発を表す運動方程式

質点が壁に衝突し、反発することを数式で表現します。

  1. 壁の定義
  2. 運動方程式の導出
  3. プログラム
  4. 実行結果
  5. 参考文献

壁の定義


壁とは、壁の法線方向に平行な質点の速度成分を反転させるデバイスである。
と定義します。

運動方程式の導出


まず議論を簡単にする為に一次元の運動について考え、その後多次元の運動について定式化をしていきます。

1次元の壁との衝突を表す運動方程式


ポテンシャル\(V(x)\)保存力下の質点の運動を考えます。
壁が
\(f(x,t)=0\)
で表現されているとします。
壁との衝突では位置は連続、速度は不連続な振る舞いをすると考えると、壁に衝突した時、力はデルタ関数の振る舞いをしていることが予想できます。
よって、未知の定数\(c\)を用いて、

と書くことが出来ます。
ここで、時刻\(t=t_i\)は、
\(
f(x(t),t)=0~~\to~~t=t_i
\)

を満たす時刻(壁との衝突時刻)です。\(f(x,t)\)の時間微分は、

と書くことが出来ます。
\(c\)を定めるために式(2)の両辺を、時刻\(t=t_j\)周りを微小時間\(\Delta\)で積分します。すると

を得ます。ここで、保存力の時間変化は連続であることを仮定します。すなわち、

がいかなる時刻\(t\)で成立するとします。計算を進めると、

という結果が導かれます。ここで、\(\delta_{i,j}\)はクロネッカーのデルタを表します。式(6)を式(2)の右辺に代入し、\(c\)を消去すると、運動方程式

を得ます。式(7)はこのままでは解くことが出来ません。なぜなら、未来の時刻の速度\(v(t_i+0)\)が含まれているからです。これをどうにかするには、新たな条件式、すなわち反発前後の条件式が必要になります。

壁と衝突する場合、質点の衝突前後の速度\(v(t_i-0), v(t_i+0)\)と、壁の衝突前後の速度\(v_\text{w}(t_i-0), v_\text{w}(t_i+0)\)の間には、反発係数\(e\)を用いて

の関係があると予想します。今、壁の質量が無限大であり、衝突前後で速度変化がない場合を考えると、壁の速度は、質点との衝突によって変化しないと考えます。すなわち

が成立していると考えます。すると、衝突後の質点の速度は、\(e\)を用いて

と書けます。よって、式(7)に代入して、

を得ます。

最後に、位置\(x(t_i)\)における壁の速度\(v_\text{w}(t_i)\)を考えましょう。衝突の前後のごく短い時間では、質点の動きは壁に追従すると考えます。すなわち、

が成り立っているとします。衝突時には質点の位置\(x\)と衝突位置は同じであることを注記します。式(12)を書き換えると、衝突が起こる時刻\(t=t_i\)の周りで

が成立しています。このことから、衝突前後のごく短い時間では

が成立します。よって、

から、壁の速度

を得ます。よって、壁との衝突を記述する運動方程式は

となります。

多次元の壁との衝突を表す運動方程式




と質点が衝突することを考えます。運動方程式は

であり、1次元の場合と同様に時刻\(t=t_j\)周りの微小時間で積分して

を得ます。式(20)に代入すれば運動方程式

を得ます。次の節で述べる結果(多次元の反発)を先取りすると、式(22)は

と変形することが出来ます。ここで、\(e\)は反発係数、\(\mathbf{v}_\text{w}(t)\)は壁の速度ベクトル、\(\mathbf{n}\)は壁の法線ベクトルであり、\(\mathbf{n}\)は

と表されます。
実際に数値計算を行う際には、衝突時刻\(t_i\)と位置\(\mathbf{r}=\mathbf{r}_i\)を求めた後、

に従って速度ベクトルを変更すると良いと思います(根拠は特にありません)。

多次元の反発



衝突時刻\(t=t_i\), 点\(\mathbf{r}(t_i)\)で質点が\(f(\mathbf{r},t)=0\)で表される壁に衝突することを考えます。
衝突前、後の質点の速度ベクトルをそれぞれ\(\mathbf{v}(t_i-0), \mathbf{v}(t_i+0),\)と置きます。
位置\(\mathbf{r}_i\)における壁の単位法線ベクトルを\(\mathbf{n}, (\mathbf{n}^2=1)\)と書くと、

と書くことが出来ます。ここで、\(c\)は未知の定数でこれから定めていきます。
衝突では壁の法線方向の速度成分のみが変化すると仮定しているので、

が成立します。ここで、\(e\)は反発係数で\(0\leq e\leq 1\)である(\(e=0\):完全非弾性衝突、\(e=1\):完全弾性衝突。
また、\(v_\perp, v_{\text{w}\perp}\)は\(\mathbf{v}, \mathbf{v}_\text{w}\)の、壁の法線方向の成分であり、

と書きあらわすことが出来ます。
式(27)の\(c\)を求めましょう。式(27)の両辺に\(\mathbf{n}^\mathsf{T}\)を掛けて

のように得ます。よって衝突後の質点の速度は

と表すことが出来ます。壁が\(f(\mathbf{r},t)=0\)を満たす線と表現されていれば、時刻\(t=t_i\)で\(\mathbf{n}\)は

と求める事が出来ます。
続いて壁面の速度\(\mathbf{v}_\text{w}(t)\)を考えます。1次元の場合と同様に衝突前後のごく短い時間では質点の動きは壁に追従すると予想します。すなわち、

が成立すると考えます。ここで、式(35)の\(\mathbf{r}(t)\)は質点の位置ベクトルです。すなわち、衝突が起こる時刻\(t=t_i\)の周りで

が成立しており、衝突のごく短い時間では

が得られます。よって

を得ますので、壁の速度

を得ることが出来ます。

プログラム


2次元平面の運動に対するFortran90のプログラムはこちらです。
時間発展は、刻み幅制御陽的ルンゲクッタ法
衝突時刻を求める際の根の探索には、Anderson-Björk’s法を用いています。

壁の関数(サブルーチン fw)とその偏微分(サブルーチン pwf)は手で入力しています。

壁との衝突は、関数の”符号が変わった時”で判定しているので、法線の方向はどちらでも構いません。
すなわち、
\(
f(x,y,t)=\pm g(x,y,t)
\)

は同じです。

プログラム

gnuplot用のスクリプト

実行結果


プログラムを実行した結果です。
適当な壁を定義して実行しています。

重力なしの運動、動かない壁

楕円の焦点から放たれた質点の軌跡

\(
\begin{align}
x(0)&=4,~~x'(0)=(\text{const}) \\
y(0)&=0,~~y'(0)=(\text{const}’)
\end{align}
\)


\(
\begin{gather}
f(x,y,t)=\left(\frac{x}{l_x}\right)^2+\left(\frac{y}{l_y}\right)^2-1=0,\\
l_x=5,~l_y=3
\end{gather}
\)

重力下の運動、動かない壁


\(
\begin{align}
x(0)&=0,~~x'(0)=2 \\
y(0)&=0,~~y'(0)=7
\end{align}
\)


\(
\begin{align}
f(x,y,t)=y – \sin(x)=0
\end{align}
\)

重力下の運動、動く壁


\(
\begin{align}
x(0)&=4,~~x'(0)=0 \\
y(0)&=4,~~y'(0)=5
\end{align}
\)


\(
\begin{gather}
f(x,y,t)=y – \sin(kx+\phi_x)\sin(\omega t+\phi_t)=0,\\
k=1,~w=1,~\phi_x=-\pi/2,~\phi_t=0
\end{gather}
\)

参考文献


3. 壁との衝突 -物理学の見つけ方

4. 動く壁との衝突 -物理学の見つけ方

Hyper-dual numbersによる二階偏微分の計算

Hyper-dual numbersと呼ばれる、実数を拡張した考えを導入すると導関数が計算できます。

あらかじめ、引数がHyper-dual numbersである時の数々の関数の定義を実装しておけば、その関数の組み合わせで作られる関数の一階導関数、二階導関数のほぼ厳密な答えを得ることが出来ます。
この考えに従う関数の微分方法は、Forward型の自動微分と呼ばれます。

  1. Dual number
    1. Dual numberと導関数
  2. Hyper-dual numbers
    1. Hyper-dual numbersの演算
    2. Hyper-dual numbersのプログラム
    3. ヘッセ行列
  3. 参考文献

Dual number


Dual numberという考えがあります[2,3,5]。

これは実数を拡張する、という考えで複素数に似た考え方です。

良く知られている実数の拡張方法の一つは、複素数です。
通常の実数に\(i^2=-1\)を満たす\(i\)という数を付加するのが複素空間です。

拡張の方法は何も\(i\)だけではありません。
例えば、\(\epsilon^2=0, (\epsilon\ne 0)\)を満たす\(\epsilon\)という数を付加することもできます。

この\(\epsilon\)は複素数ではありません
複素平面上でこの性質を満たす数は無いことからも、この\(\epsilon\)を追加するということは新しい方向への数の拡張です。
実際、複素平面上で\(\epsilon^2=0, (\epsilon\ne 0)\)を満たす数があるのかを探しても
\(
\begin{align}
(a+ib)^2=a^2-b^2+i2ab=0
\end{align}
\)

であるので、これを満たすのは\(a=b=0\)しか存在せず、複素平面上の数ではないことが分かります。

そんな新しい方向\(\epsilon\)の平面で定義された数をDual number(二重数、または双対数)と呼びます。

二重数\(a\)は実数部と非実数部から構成されており、
\(
a=a_0+a_1\epsilon
\)

と書くことが出来ます。ここで、\(a_0, a_1\)は実数であり、 \(\epsilon\)は、虚数単位\(i\)に倣って二重数単位とでも名付けておきましょう。

Dual numberと導関数


二重数は関数の導関数と大きな関係がある事を示しましょう。
テーラー展開を行います。
関数\(f\)の\(x\)周りの展開は刻み幅\(\Delta\)を用いて
\(\displaystyle
f(x+\Delta)=f(x)+f'(x) \Delta+\frac{1}{2!}f^{\prime\prime}(x)\Delta^2+\cdots
\)

と記述することが出来ます。もし仮に\(\Delta\)が純非実数だとしましょう。すなわち、\(\Delta=h\epsilon\)とします。ここで、\(h\)は実数です。テーラー展開の式に代入すれば、
\(
\begin{eqnarray}
f(x+h\epsilon)&=& f(x)+f'(x) h\epsilon+\frac{1}{2!}f^{\prime\prime}(x){h}^2\epsilon^2+\cdots \\
&=& f(x)+f'(x) h\epsilon
\end{eqnarray}
\)

となるわけです。注記しますが、上の式は\(h\)の1次で打ち切っているのではなく、厳密にイコールが成り立っています。
この式が言っているのは、関数\(f(a), (a\text{は二重数})\)を計算し、その非実数部を\(h\)で割ると関数の導関数になっているということです。

二重数の非実数部を取り出す関数\(\text{Dual}\)を定義します。
すると、導関数は
\(\displaystyle
f'(x)=\frac{1}{h}\text{Dual}(f(x+h\epsilon))
\)

として得られます。
言葉で書けば、二重数空間に拡張した\(f(x+h\epsilon)\)を計算すると、その非実数部に導関数が現れる、ということです。

二重数のプログラムはJeffrey Fikeさんによる[4]にありますので、そちらをご参考ください。

Hyper-dual numbers


さて、ここまでで二重数の概念を簡単に説明し、新しい実数の拡張を行いました。
Dual Numberのままでは高階導関数は得られません。なぜなら、導関数の二次以降は二重数の性質\(\epsilon^2=0\)によってゼロになるからです。

そこで、若干工夫します。
Dual Numberを二種類用意することで二階微分を得ることが出来ます[2,3]。
すなわち、
\(\begin{gather}
a=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2\\
\epsilon_1^2=\epsilon_2^2=(\epsilon_1\epsilon_2)^2=0,~\epsilon_1\ne 0,~~\epsilon_2\ne 0,~~\epsilon_1\epsilon_2 \ne 0
\end{gather}
\)

と実数を拡張します。この様に拡張した数\(a\)をHyper-dual numbersと呼びます[2,3]。

Hyper-dual numbersの演算


Hyper-dual numbersである
\(
\begin{align}
a&=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2\\
b&=b_0+b_1\epsilon_1+b_2\epsilon_2+b_3\epsilon_1\epsilon_2
\end{align}
\)

を用意します。和、積、商はそれぞれ
\(
\begin{align}
a+b&=(a_0+b_0)+(a_1+b_1)\epsilon_1+(a_2+b_2)\epsilon_2+(a_3+b_3)\epsilon_1\epsilon_2\\
ab&=a_0b_0+(a_0b_1+a_1b_0)\epsilon_1+(a_0b_2+a_2b_0)\epsilon_2+(a_0b_3+a_1b_2+a_2b_1+a_3b_0)\epsilon_1\epsilon_2\\
\frac{1}{a}&=\frac{1}{a_0}-\frac{a_1}{a_0^2}\epsilon_1-\frac{a_2}{a_0^2}\epsilon_2+\left(-\frac{a_3}{a_0^2}+\frac{2a_1a_2}{a_0}\right)\epsilon_1\epsilon_2
\end{align}
\)
と定義されます[2]。任意の関数は
\(
f(x)=f(x_0)+x_1f'(x_0)\epsilon_1+x_2f'(x_0)\epsilon_2
+\left(x_3f'(x_0)+x_1x_2f^{\prime\prime}(x_0)\right)\epsilon_1\epsilon_2
\)

として計算することが出来ます。なので、この結果から\(\epsilon_1\epsilon_2\)の係数として関数の二階微分が得られます。

Hyper-dual numbersのプログラム


実際にプログラムを組みましょう。
Hyper-dual numbersの型を持つ変数はFortran90では定義できません。
なので、構造体を利用して自分で型と、演算を定義します。

Fortran90のプログラムは以下の通りになるかと思います。
基本的な四則演算、基本的な初等関数のHyper-dual numbersの演算をモジュールとして書いています。

上のモジュールと下のメインプログラムを一緒にコンパイルすることにより、関数の二階微分が得られます。

例として
2変数関数
\(\displaystyle
f(x,y)=\frac{\text{ln}(xy^2)e^x}{\sqrt{\sin^3{x}+\sin^3{y}}}
\)

の偏微分
\(\begin{align}
\frac{\partial f}{\partial x},~~\frac{\partial f}{\partial y},~~\frac{\partial^2 f}{\partial x\partial y}
\end{align}
\)

を得ることを考えます。
\(
f(x+1\epsilon_1,y+1\epsilon_2)
\)

を計算すると
\(
f(x+1\epsilon_1,y+1\epsilon_2)=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2
\)

のように実数係数\(a_0, a_1, a_2, a_3\)が得られます。
すると、
\(
\begin{align}
\frac{\partial f}{\partial x}=a_1 \\
\frac{\partial f}{\partial y}=a_2 \\
\frac{\partial^2 f}{\partial x\partial y} =a_3
\end{align}
\)

として偏微分が得られます。
ちなみに、二階微分が欲しい場合は
\(
f(x+1\epsilon_1+1\epsilon_2, y)
\)

を考えると
\(
\begin{align}
\frac{\partial f}{\partial x}=a_1 \\
\frac{\partial f}{\partial x}=a_2 \\
\frac{\partial^2 f}{\partial x^2} =a_3
\end{align}
\)

として得られます。\(a_1,~a_2\)はどちらを採用しても構いません。

プログラムでは変数の型Hyperdualを持つ入力変数を\(\text{xH,yH}\),出力を\(\text{rH}\)と置いています。

program main
  use Hyperdualmod
  implicit none
  type(Hyperdual)::xH,yH,rH
 
  xH%x0 = 0.3d0 ! real part
  xH%x1 = 1d0   ! unreal part \epsilon_1
  xH%x2 = 0d0   ! unreal part \epsilon_2
  xH%x3 = 0d0   ! unreal part \epsilon_1\epsilon_2

  yH%x0 = 0.4d0 ! real part
  yH%x1 = 0d0   ! unreal part \epsilon_1
  yH%x2 = 1d0   ! unreal part \epsilon_2
  yH%x3 = 0d0   ! unreal part \epsilon_1\epsilon_2

  write(6,'(4f23.16)')xH%x0,xH%x1,xH%x2,xH%x3
  write(6,'(4f23.16)')yH%x0,yH%x1,yH%x2,yH%x3  
 
  rH = log(xH*yH**2)*exp(xH)/sqrt(sin(xH)**3+cos(yH)**3)
  !rH = asin(2d0*xH)*acos(yH)/atan(xH*yH)
  !rH = xH**yH
 
  write(6,'(4f23.16)')rH%x0,rH%x1,rH%x2,rH%x3  

  stop
end program main

ヘッセ行列


二階偏微分の計算が出来たので、ヘッセ行列が簡単に計算できます。
ルーチンを作れば、以下の通りになるかと思います。
下のプログラムは3変数関数
\(\displaystyle
f(x,y,z)=\exp(xy)\tan(z)
\)

の\(x=-2,~y=3,~z=1\)におけるヘッセ行列を計算します。

program main
  use Hyperdualmod
  implicit none
  integer::N
  double precision::x,y,z,f
  double precision,allocatable::nabla(:),Hesse(:,:)
  double precision,allocatable::w(:)
  external::func

  N=3

  allocate(nabla(1:N),Hesse(1:N,1:N))
  nabla=0d0
  Hesse=0d0
  allocate(w(1:N))
  w=0d0

  x = -2d0
  y =  3d0
  z =  1d0
  w(1) = x
  w(2) = y
  w(3) = z
  call Hessian(N,w,func,nabla,Hesse)

  f=exp(x*y)*tan(z)
  write(6,'(1e24.16)')nabla(1)
  write(6,'(1e24.16)')y*f
  write(6,'(1e24.16)')nabla(2)
  write(6,'(1e24.16)')x*f
  write(6,'(1e24.16)')nabla(3)
  write(6,'(1e24.16)')exp(x*y)/(cos(z)**2)
  write(6,'(3e24.16)')Hesse(1,1:3)
  write(6,'(3e24.16)')f*y**2, f*(1+x*y), f*y/(cos(z)*sin(z))
  write(6,'(3e24.16)')Hesse(2,1:3)
  write(6,'(3e24.16)')f*(1+x*y), x**2*f, f*x/(cos(z)*sin(z))
  write(6,'(3e24.16)')Hesse(3,1:3)
  write(6,'(3e24.16)')f*y/(cos(z)*sin(z)), f*x/(cos(z)*sin(z)), f*2d0/(cos(z)**2)

  stop
end program main

subroutine func(N,x,f)
  use Hyperdualmod
  implicit none
  integer::N
  type(Hyperdual),intent(in)::x(1:N)
  type(Hyperdual),intent(out)::f
 
  f = exp(x(1)*x(2))*tan(x(3))
 
  return
end subroutine func

subroutine Hessian(N,x,func,nabla,Hesse)
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  double precision,intent(in)::x(1:N)
  double precision,intent(out)::nabla(1:N)
  double precision,intent(out)::Hesse(1:N,1:N)
  external::func

  integer::i,j,k
  type(Hyperdual),allocatable::xH(:)
  type(Hyperdual)::rH
 
  allocate(xH(1:N))
  do i=1,N
     xH(i)%x0 = 0d0
     xH(i)%x1 = 0d0
     xH(i)%x2 = 0d0
     xH(i)%x3 = 0d0
  enddo
 
  do i=1,N
     xH(i)%x0 = x(i)
     xH(i)%x3 = 0d0
  enddo

  do i=1,N
     do k=1,N
        xH(k)%x1 = 0d0
        xH(k)%x2 = 0d0
     enddo
     xH(i)%x1 = 1d0
     xH(i)%x2 = 1d0
     
     call func(N,xH,rH)
     nabla(i) = rH%x1
     Hesse(i,i) = rH%x3
  enddo

  do i=1,N
     do j=i+1,N
        do k=1,N
           xH(k)%x1 = 0d0
           xH(k)%x2 = 0d0
        enddo
        xH(i)%x1 = 1d0
        xH(i)%x2 = 0d0
        xH(j)%x1 = 0d0
        xH(j)%x2 = 1d0
        call func(N,xH,rH)
        Hesse(i,j) = rH%x3
        Hesse(j,i) = Hesse(i,j)
     enddo
  enddo  
 
  return
end subroutine Hessian

Hyperdual.f90にモジュールを入れ、メインプログラムをmain.f90に入れたとすると、以下の結果を得ます。

> gfortran Hyperdual.f90 main.f90  
> ./a.out
  0.1158128336233602E-01
  0.1158128336233602E-01
 -0.7720855574890680E-02
 -0.7720855574890680E-02
  0.8491012233306163E-02
  0.8491012233306162E-02
  0.3474385008700806E-01 -0.1930213893722670E-01  0.2547303669991849E-01
  0.3474385008700806E-01 -0.1930213893722670E-01  0.2547303669991849E-01
 -0.1930213893722670E-01  0.1544171114978136E-01 -0.1698202446661233E-01
 -0.1930213893722670E-01  0.1544171114978136E-01 -0.1698202446661233E-01
  0.2547303669991849E-01 -0.1698202446661233E-01  0.2644793608458059E-01
  0.2547303669991849E-01 -0.1698202446661233E-01  0.2644793608458058E-01

実行結果の奇数行目はHyper-dual Numberによる計算結果、偶数行目は解析解を表します。
また、6行目までは一階微分、7行目以降はヘッセ行列を表します。

参考文献


[1]関数
\(\displaystyle
f(x,y)=\frac{\text{ln}(xy^2)e^x}{\sqrt{\sin^3{x}+\sin^3{y}}}
\)

の偏微分\(\frac{\partial^2 f}{\partial x\partial y}\)
の計算
https://www.wolframalpha.com/input/?i=D%5Bln(x*y%5E2)*e%5E(x)%2Fsqrt(sin%5E3(x)%2Bcos%5E3(y)),+y,x%5D

[2]Jeffrey A. Fike and Juan J. Alonso,~”The Development of Hyper-Dual Numbers for Exact Second-Derivative Calculations”, 49th AIAA Aerospace Sciences Meeting(2011)
http://adl.stanford.edu/hyperdual/fike_aiaa-2011-886_slides.pdf,
J. A. Fike and J. J. Alonso. The Development of Hyper-Dual Numbers for Exact Second Derivative Calculations. AIAA paper 2011-886, 49th AIAA Aerospace Sciences Meeting, January 4-7, 2011.http://adl.stanford.edu/hyperdual/Fike_AIAA-2011-886.pdf

[3]JeffreyA.Fike,~”Derivative Calculations Using Hyper-Dual Numbers”, Sandia National Laboratories (2016)
https://www.osti.gov/servlets/purl/1368722

[4]Jeffrey Fike, Aerospace Design Lab, http://adl.stanford.edu/hyperdual/

[5]松本佳彦, 新しい数をつくる, (2018) http://ymatz.net/assets/docs/20180629-jtpa-slide-mod

束縛条件下の運動 – 自由度がうまく落とせない運動

このページは
束縛条件下の運動 – ホロノミックな束縛と非保存力
の続きです。

拘束がある場合の2次元の運動


ラグランジュの方程式を解く際に適切な座標変換が見付からず自由度が落とせない場合を考えます。その時の運動方程式は

の形で止まります。式(75)において、未知の関数は\(x(t), y(t), \lambda(t)\)の3つであり、条件式は3つなので解くことは出来るはずです。
方針は、式(75)から\(\lambda(t)\)を消去すれば良いのです。

束縛条件\(f(x,y,t)=0\)が成立する点\(x=x_0, y=y_0\)周りで、微小時間\(\Delta t\)の変化を考えます。
点\((x_0, y_0)\)まわりでテーラー展開すれば、

を得ます。よって

が満たされる変化でなければなりません。すなわち、\(x,y\)の時間変化は式(77)をいつも満たしていなければならないのです。両辺を\(\Delta t\)で割って極限をとれば、

が成立することが分かります。式(75)を書きかえれば、

を得ます。式(79c)もまた時間変化しても成立し続けなければなりません。よって、左辺の時間微分を取れば、

を得ます。整理して

を得ます。ここで関係式

を用いました。式(79a), (79b)から未定乗数\(\lambda(t)\)を消去すれば

を得ます。よって、式(81)と式(84)から連立方程式

を立てることが出来ます。これをあらわに解けば、運動方程式

を得ます。ここで、式(86), (87)の右辺の\(-\frac{\partial V}{\partial x}, -\frac{\partial V}{\partial y}\)は\(x, y\)方向の保存力であり、右辺の残りの項は束縛条件によって決まる束縛力であることが分かります。

拘束がある場合の3次元の運動


スカラー関数\(f\)で表される拘束がある条件のもとで, 運動方程式

またはベクトル表記で

について考えます。ここで表記を略すために

の略記を用いました(\(y, z\)方向についても同様)。
また、\(\nabla\)はナブラ演算子で、

です。2次元の場合と同様に、式(88d)の時間微分から

が導けます。ここで、\(\hat{H}f\)は関数\(f\)のヘッセ行列であり、

を表します。式(92)は、式(81)を3次元に拡張した形になっています。

式(88a),(88b),(88c)から未定乗数を消去すると従属な3つの関係式

を得ることが出来ます。この3式で独立なのは2つです。
独立な式として式(94a)と式(94b)を選ぶことにし、式(92)と共に書けば、連立方程式

を得ます。この式をあらわに解けば、運動方程式

を得ることが出来ます。ここで、\(n_x, n_y, n_z\)は壁の法線ベクトル\(\mathbf{n}=\frac{\nabla f}{|\nabla f|}\)の\(x,y,z\)成分です。偏微分ではないことに注意してください。

ベクトル表記であれば

と書くことが出来ます。

単振り子の例


ここでは単振り子を例に挙げ、式(86)に代入した時本当に振り子を表す方程式になっているのか、また束縛条件の選び方に依らず、同じ運動を記述しているのかを例に挙げます。
保存力は重力で

を考えます。

考える束縛条件は
\(
f(x,y,t)=x^2+y^2-l^2=0
\)


\(
f(x,y,t)=\sqrt{x^2+y^2}-l=0
\)

です。

束縛条件1


束縛条件として

を考えます。この束縛条件に対する偏微分を計算して、

を得ます。これを式(86), (87)に代入すれば、

であり、整理すれば

を得ます。

束縛条件2


束縛条件として

を考えます。この束縛条件に対する偏微分を計算して、

を得ます。これを式(86), (87)に代入すれば、

であり、整理すれば

を得ます。

考察


同じ束縛条件なので、式(105)と(109)の示す運動は同じになるはずです。
それを示すには、2式の差異である

を示せればよい、ということです。そのために拘束力に対して座標変換

を考えます。右辺と左辺をそれぞれ計算すれば、

となり、同じ運動を記述していることが分かります。

Fortran90によるプログラム


Fortran90で2次元平面を動く質点の束縛運動のプログラムを作ります。
そのプログラムは
https://slpr.sakura.ne.jp/qp/supplement_data/constraint.tar.gz
に置いてあります。

一応補足しておきますが、初期状態の位置は\(f(x, y, t)=0\)を満たすx,yでなければなりませんし、初期速度も\(f(x, y, t)=0\)の傾きと一致していなければなりません。
上の条件を満たさずとも何も警告は出さず、計算自体は行われますが、その計算は束縛運動ではなくなります。

上のプログラムの中で、変更すると思われるのは以下の2つです。
1. fp2d
2. fw2d

1. fp2dは、ポテンシャル\(V(x,y)\)を記述しており、デフォルトのプログラムでは、重力による位置エネルギー
\(
V(x,y)=mgy
\)
に設定しています。

2. fw2dは、束縛条件\(f(x,y, t)=0\)を記述しており、デフォルトのプログラムでは、
\(
f(x,y,t)=y – [2(cos(1.5*x)-1) + 0.1x^2]=0
\)
に設定しています。

もしも非保存力を入れたい場合は、grkの最後の方を変更してください。

数値計算では、
時間発展は刻み幅制御の陽的ルンゲクッタ、
gradient, Hesse行列、束縛条件に関する時間の二階微分は、Hyper-dual numbersを利用してほぼ厳密な値を得ています。

program main
  ! Author : sikino
  ! Date   : 2019/10/12 (yyyy/mm/dd)
  ! URL    : http://slpr.sakura.ne.jp/qp/constraint-condition3/
  use Parameters
  use Hyperdualmod
  implicit none
  integer::i,N,info,Nt
  double precision::t,h,tol,ta,tb
  double precision,allocatable::x(:),tgrid(:)
  external::grk
  double precision::E,Ek,Ep

  N=4
  allocate(x(1:N))

  ! Time range [ta, tb]
  ta = 0d0
  tb = 30d0
  ! Divide time range equally among Nt
  Nt = 201
  ! Tolerance of the RK method
  tol=1d-8

  ! Initial condition
  x(1)=0d0 ! x (ta)
  x(2)=4d0 ! x'(ta)
  x(3)=0d0 ! y (ta)
  x(4)=0d0 ! y'(ta)

  ! Time grid
  allocate(tgrid(1:Nt))
  do i=1,Nt
     tgrid(i) = (i-1)*(tb-ta)/dble(Nt-1)+ta
  enddo
 
  do i=2,Nt
     info = 0
     t = tgrid(i-1)
     h = (tgrid(i) - t)*0.1d0
     do while(info.le.0)
        call drkf45(grk,t,h,N,x,tgrid(i),info,tol)
     enddo
     
     ! Kinetic energy
     Ek = 0.5d0*mass*(x(2)**2+x(4)**2)
     ! Potential energy
     Ep = mass*g*x(3)
     ! Total energy
     E  = Ek + Ep
     write(10,'(8e25.10e3)')t,x(1),x(2),x(3),x(4),E,Ek,Ep
  enddo
 
  stop
end program main

subroutine fp2d(N,xH,fH)
  use Parameters, only:mass, g
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  type(Hyperdual),intent(in)::xH(1:N)
  type(Hyperdual),intent(out)::fH

  ! Potential

  type(Hyperdual)::x,y,t
  x = xH(1)
  y = xH(2)
  t = xH(3)
  ! If conservative force,
  !  fH doesn't depend on the time.
  !-----------------------------
 
  fH = mass * g * y
 
  return
end subroutine fp2d

subroutine fw2d(N,xH,fH)
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  type(Hyperdual),intent(in)::xH(1:N)
  type(Hyperdual),intent(out)::fH

  ! Wall of the function

  type(Hyperdual)::x,y,t
  x = xH(1)
  y = xH(2)
  t = xH(3)
  !-----------------

  fH = y - (2d0*(cos(1.5d0*x)-1d0) + 0.1d0*x*x)
!  fH = (y-5d0)**2 + x**2 - 25d0
 
  return
end subroutine fw2d

デフォルトのまま計算すると、fort.10に以下の出力が現れます。アニメーションは、同じファイルに入っているanime.pltで実行できます。

エネルギーは計算の範囲ではほとんどありません。全力学的エネルギーは、計算時間内でおおよそ5桁目が変化していました。

参考文献


6. 振り子 -物理学の見つけ方
9. 自由な座標 -物理学の見つけ方

最速のクイックソート(Fortran)

クイックソートが速いのは分かっていますが、コーディングの仕方によって速度は変わります。
本稿では、ネット上で公開されているどのクイックソートが早いのか調べていきます。

最も早いクイックソートはNUMPACのプログラムでした。

比較プログラム


比較するプログラムはネット上で公開されているクイックソート+αです。
対象は、倍精度実数をソートするプログラム、です。

1. 再帰を用いたクイックソート(f90)

t-nissieの日記: 【電脳】Fortranで書いたクイックソート

2. 再帰を用いないクイックソート(f90)

[Fortran]再帰を使わないquicksortその2 -fortran66の日記

※ただし、私が割と変えていますので、オリジナルそのものではないということを注記しておきます。この変更によって、実行速度はほぼほぼ変わらないことは確認しています。

3. Netlibのクイックソートdsort.f

https://www.netlib.org/slatec/src/dsort.f

4. NUMPACのクイックソート(sortdk.f)

SORTPACK(SORTxK,SORTxy,SRTVxz) (スカラー又はベクトルデータの内部ソーティング)

5. ヒープソート

fortran90によるヒープソートとバブルソート -シキノート

以上の5つのプログラムを比較していきます。
メインプログラムはこちら↓

コンパイルは

gfortran -O3 (ソートのプログラム達) main.f90

で行いました。

結果


結果を示します。
横軸にデータ数、縦軸にソートに掛かった時間を示しました。
”ソートに掛かった時間”とは、同じデータ数でソートを100回繰り返した時の1回当たりの時間です。


それぞれ、
赤線:再帰有りクイックソート
青線:再帰無しクイックソート
緑線:Netlibのクイックソート
紫線:NUMPACのクイックソート
黒線:ヒープソート
です。
最も早いのがNUMPAC, 最も遅いのがヒープソートだと分かりました。

ヒープソートもクイックソートも大体\(O(n\log n)\)ですので、グラフの傾きは両者でほとんど変わりません。再帰有りの方が遅い、という結果は面白いです。
先入観で”再帰は遅い”と考えていたのですが、そんなことはありませんでした。

続いて、同じデータを対数で見てみると以下の通りです。

走査したデータサイズの範囲では変化は有りません。更に大規模になれば、違いが見えてくるかもしれません。