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

Fortran90で構造体を使う

Fortran90で自分で型を定義し、それを使用します。
本稿でいう”型”は、一般には構造体と呼ばれます。
Fortran使用者にとっては”型”という名称が分かりやすいでしょう。

整数型、倍精度型などと同列にあるもの、と考えて差し支えないです。

本稿では基本的なことしか書いていません。

二階偏微分の計算の時でも構造体の考えを使用していますので、もう少し例が欲しい方はこちらもどうぞ。
Hyper-dual numbersによる二階偏微分の計算 -シキノート

概要

簡単な実例



module CEOM_mod
  implicit none
  type Particle
     double precision::m,x,v
  end type Particle
end module CEOM_mod

program main
  use CEOM_mod
  implicit none
  type(Particle)::p1

  p1%m  = 1d0
  p1%x  = 2d0
  p1%v  = 3d0

  write(6,*)p1%m,p1%x,p1%v

  stop
end program main

実行結果

$ ./a.out
   1.0000000000000000        2.0000000000000000        3.0000000000000000
$

説明


上は、粒子p1が質量m, 位置x, 速度vを持つ、ということを想定しています。

そこで、Particleという新しい”型”(複素数型みたいに、実部と虚部という2つの独立な変数を1つの複素数型変数が持つ、みたいな変数の型)を定義しました。
Particle型を持つ変数は、すべてm, x, vを内包するというイメージです。

mainプログラムではparticle型が定義されているモジュールを呼び出すために

    use CEOM_mod

という文を入れ、実際にparticle型の変数p1を定義するために

  type(Particle)::p1

と書いています。変数p1が定義されると自動的にp1が持つm,x,vの固有な値も定義されます。固有の値をいじる場合は

  p1%m  = 1d0

と書いて、変数に%をつけてm,x,vのいずれかを書いて代入します。


配列



module CEOM_mod
  implicit none
  type Particle
     double precision::m,x,v
  end type Particle
end module CEOM_mod

program main
  use CEOM_mod
  implicit none
  integer::i,N
  type(Particle),allocatable::p(:)

  N=5
  allocate(p(1:N))
  do i=1,N
     p(i)%m  = 1d0
     p(i)%x  = dble(i)
     p(i)%v  = 0.1d0*i
  enddo

  do i=1,N
     write(6,*)i,p(i)%m,p(i)%x,p(i)%v
  enddo

  stop
end program main

実行結果

$ ./a.out
                    1   1.0000000000000000        1.0000000000000000       0.10000000000000001<br />
                    2   1.0000000000000000        2.0000000000000000       0.20000000000000001<br />
                    3   1.0000000000000000        3.0000000000000000       0.30000000000000004<br />
                    4   1.0000000000000000        4.0000000000000000       0.40000000000000002<br />
                    5   1.0000000000000000        5.0000000000000000       0.50000000000000000
$

サブルーチンの利用



module CEOM_mod
  implicit none
  type Particle
     double precision::m,x,v
  end type Particle
end module CEOM_mod

subroutine twicemass(N,p)
  use CEOM_mod
  implicit none
  integer,intent(in)::N
  type(Particle),intent(inout)::p(1:N)

  integer::i
 
  do i=1,N
     p(i)%m  = 2 * p(i)%m
  enddo
 
  return
end subroutine twicemass

program main
  use CEOM_mod
  implicit none
  integer::i,N
  type(Particle),allocatable::p(:)

  N=3
  allocate(p(1:N))
  do i=1,N
     p(i)%m  = 1d0
     p(i)%x  = dble(i)
     p(i)%v  = 0.1d0*i
  enddo

  call twicemass(N,p)
 
  do i=1,N
     write(6,*)i,p(i)%m,p(i)%x,p(i)%v
  enddo

  stop
end program main

実行結果

$ ./a.out
                    1   2.0000000000000000        1.0000000000000000       0.10000000000000001<br />
                    2   2.0000000000000000        2.0000000000000000       0.20000000000000001<br />
                    3   2.0000000000000000        3.0000000000000000       0.30000000000000004<br />
$

補足


オブジェクト指向になれている方は、構造体の中身を外部から直接アクセスできる点で褒められた書き方ではないでしょう。
Fortran90にオブジェクト指向の考え方を持ってくるのは、計算速度というFortran90の長所を潰してしまうので、このくらいで良いでしょう。

ガウス求積法の導出

ガウス求積法は、関数の数値積分を行う方法です。
N個の点で関数を参照すると、2N-1次以下の多項式で表される関数を厳密に数値積分することができる方法です。

ガウス求積法の凄さを少し述べます。
通常、3個の任意の点における値があれば、その3点を結ぶ2次関数を作ることができるので、うまい選び方をすればN個の点でN-1次の関数を厳密に表現できそうであり、それが厳密に積分できてもそこまで不思議ではありません。

ですが、ガウス求積法では分点の位置にも制限を加えることで積分の次数を倍にすることができます。
分点の位置を制限するにもかかわらず、高々二倍になるだけか…とも思われますが、これは同じ精度に達するまでに関数の評価点数を半分に減らせることと同じです。
評価点数を少なくできるというこの性質は、多次元の計算において顕著に現れます。なぜなら、1次元で半分になるので2次元では4分の1になります。

多少無理してもガウス求積法を使う価値があります。
科学技術計算ではガウス求積法を使うことをお勧めします。ガウス求積法を推奨します
ガウス求積法を使いましょう(圧)。

本稿のPDF版はこちらのリンクからどうぞ。

https://slpr.sakura.ne.jp/qp/wp-content/uploads/2021/03/Gauss_legendre_quadrature_derivation_20210904.pdf

導出手順


導出は以下の順番で行います。

  1. 点数と次数の関係
  2. 具体的な位置と重みの計算

『点数と次数の関係』では、M個の離散的な点を用いてN次多項式を厳密に積分できるとき、NとMの関係を導出します。
『具体的な位置と重みの計算』では、具体的に、M点の位置とその重みを導出します。

ガウス求積法を利用する方は、1. が分かればよいと思います。
ガウス求積法を研究する方は、2. まで分かればよいと思います。

また、本稿ではガウス求積法の中で最も基本的な
ガウス=ルジャンドル求積法
について記述します。他の方法も、ガウス=ルジャンドル求積法の導出を理解できれば、おのずと分かるでしょう。

先に結論をまとめます。
\(2N-1\)次以下の多項式\(f(x)\)の積分を、N個の点だけを用いて厳密に計算できます。
つまり、

を厳密に計算する点の位置\(x_i\)と\(w_i\)が存在し、N個の\(x_i\)とN個の\(\omega_i\)は

から求めることができます。ここで、\(P_n(x)\)はn次ルジャンドル多項式です。

点数と次数の関係


N次多項式を有限区間\([-1,1]\)で、M個の点を用いて数値的に厳密に積分することを考えます。つまり、関数\(f(x)\)がN個の既知の係数\(a_n\)を用いて

と書けている場合に、積分

を考えます。この積分を離散的な点だけを用いて計算します。つまり、適当なM個の点を用いて

となるような右辺の\(x_i, \omega_i\)の組み合わせを求めることが目標です。

式(3)をできるだけ計算していきましょう。

まず、式(3)の左辺を積分すると、

となります。

続いて、式(3)の右辺を計算すると、

と書けます。ここで、\(\omega_i\)は重みであり、適当な係数です。もしも任意の\(a_n\)について、式(3) の右辺と左辺が厳密に一致するならば、\(n=0, 1,\cdots, N\)について、

が成立します。この式は合計で\(N+1\)個の方程式の数があります。今、右辺の未知の係数は\(\omega_i, x_i, (i=1,\cdots, M)\)なので、右辺の自由度は\(2M\)個あります。
つまり、両辺を一致させるような係数の組み合わせがあるならば、

である必要があります。よって、もっとも小さくなるような\(M\)であるならば、上式の等号が成立します。できるだけ少ない点数で計算したいので、等しい場合を採用して

となります。つまり、M個の点を利用すれば\(2M-1\)次の多項式を厳密に計算できる点の位置\(x_i\)と\(w_i\)が存在します。

具体的な位置と重みの計算


具体的な分点の位置

積分

を求めるために、離散的なN個の点\(x=x_i\)における\(f(x)\)の値だけを利用します。つまり、

の点は既知とします。ここで、\(f_i\equiv f(x_i)\)と定義しました(この時点でまだ\(x_i\)は決まっていません)。
新たな関数\(f(x)\)を定義します。\(f(x)\)は、式(12)のN個の点を全て通る\(N-1\)次多項式であり、一意に決まります。具体的には、ラグランジュ多項式を用いて

と書けます。ここで、更に\(f(x)\)と\(f(x)\)の差の関数\(g(x)\)

を定義します。\(x\)の次数を確認しておくと、\(g(x), f(x)\)は\(2N-1\)次、\(f(x)\)は\(N-1\)次多項式です。
ここで、\(f(x_i)=F(x_i),~(\mbox{for all } i)\)なので、\(g(x)\)は

の性質を持ちます。積分を考えると

であり、\(f(x)\)を離散的に表すと、

です。今、\(f(x)\)は\(2N-1\)次であり、\(f(x)\)はN個の点で表現されているので、前節の議論から式(17)をゼロにするような\((x_i, \omega_i)\)の組み合わせが存在することが分かります。
今, \(g(x)\)の満たす性質として式(15)が分かっていますので、N個の\(x_i\)について、

と書き換えることができます。
ここで、\((x-x_1)(x-x_2)\cdots(x-x_N)\)はN次多項式で、\(G(x)\)は\(N-1\)次多項式です。つまり、

を満たすN個の\(x=x_i\)を探すことができれば嬉しいわけです。ここで、\(G(x)\)は\(N-1\)次多項式なので、\(N-1\)次のlegendre多項式を用いた基底関数で展開可能です。なので、

と書くことができます。また、

はN次なので、N次多項式の関数で書くことができます。

つまりN次多項式としてN次legendre多項式\(P_n(x)\)を選ぶことにします。
単にN次であるだけだと、N次以下の全てのlegendre多項式の和で書けなければなりませんが、点\(x_1, x_2, \cdots, x_N\)を\(P_N(x)=0\)を満たすN個の点を作為的に選ぶことにしましょう。

このように選ぶと、N次以下の和ではなく、N次のlegendre多項式だけで書くことができて、

とする定数\(c\)を見つけることができます。

これからの議論において定数\(c\)を具体的に決める必要はありませんが明確に決めたい場合

の関係から決まります。

式(19)に式(20), 式(22)を代入すると

となります。これはつまり、N次legendre多項式のゼロ点を求積法におけるゼロ点として選ぶと

を満たすように\(g(x)=0\)にできるということです。

ここまでの議論で、分点の位置\(x_i\)が求まりました。
残るは、まだ決定していない重み\(\omega_i\)を決めることです。

具体的な重み

上記までの議論から、\(2N-1\)次の多項式で記述される\(f(x)\)の積分は、N個の点を用いて

を厳密に積分する分点の位置と重み\(x_i, \omega_i\)の組があることを示し、分点の位置\(x_i\)は、N次legendre多項式のゼロ点、すなわち

を満たすN個のゼロ点であることが分かりました。次は重みを決めます。
式(16)より\(2N-1\)次の\(f(x)\)をN次のlegendre多項式のゼロ点を通る多項式の積分に置き換えできます。

つまり、式(13)を用いて

のように置き換えできます。つまり、重みは

ですので、これを求めることが問題となります。

積和は書き換えがやりにくいですので、意図的に消していきましょう。
式(22)から、

分かっています。意図的にk番目の式についてまとめると、

となります。また、式(35)の両辺を微分すると

となります。特に、\(x=x_k\)ならば、

式(36)と式(40)を組み合わせると、

なので、重み\(\omega_k\)は

となります。この積分を計算するために多項式の性質を用いて計算を進めましょう。
右辺の形は、クリストッフェル=ダルブーの定理から導けそうな形であることを当たりをつけて計算します(私のような凡人には思い付きません。天下り的に使います。導けるからそれで良いです)。

直交多項式(Legendre多項式やLaguerre多項式など)にはクリストッフェル=ダルブーの定理 (Christoffel–Darboux theorem) があります。それは、n次の直交多項式を\(P_n(x)\)と書くとき、\(P_n(x)\)には

という性質であります。ここで、\(k_n\)は、直交多項式\(P_n(x)\)の\(x^n\)の係数で、

で表されるときの係数です。この係数を具体的に求めてみましょう。そのために、
legendre多項式の性質の一つである三項漸化式

を用います。xの次数で両辺を比較すると

なので、次数を比較すれば最高次の係数について

の関係があります。よって、

を得ます。また、\(h_N\)は直交性から導かれる係数で、

を満たします。legendre多項式の場合、\(w(x)=1\)であり、\(h_n=\frac{2}{2n+1}\)です。

重みの計算に使いたいため、クリストッフェル=ダルブーの定理に対し、\(N\to N-1\)にして、legendre多項式の場合に適用、さらに\(y=x_k\)の場合を考えて、

より、

を得ます。ここで\(P_{N}(x_k)=0\)を利用しました。

式(43)に式(56)を代入して、

となります。\(P_n(x)\)の微分は次数が一つ下がるので、\(P_{n-1}(x)\)と同じ次数になります。つまり、両者に何か関係がありそうですので頑張って見つけてみます。

微分に関する三項漸化式

を用います。特に\(n=N,~x=x_k\)の場合、

の関係があります。代入すると重みは

と求めることができます。

まとめ


\(2N-1\)次以下の多項式\(f(x)\)の積分を、N個の点だけを用いて厳密に計算できます。
つまり、

を厳密に計算する点の位置\(x_i\)と\(w_i\)が存在し、\(x_i, \omega_i\)は

から求めることができます。ここで、\(P_n(x)\)はn次ルジャンドル多項式です。

参考文献など


恐らく、ガウスが発表した論文がこれではないかと思います。ドイツ語のwikipedia(https://de.wikipedia.org/wiki/Gau%C3%9F-Quadratur)にリンクがありました。
Methodus nova integralium valores per approximationem inveniendi. In: Comm. Soc. Sci. Göttingen Math. Band 3, 1815, S. 29–76

千葉豪著, Gauss-Legendre求積公式の導出, (2016)
http://roko.eng.hokudai.ac.jp/studentadm/chiba_data/PnSn/gl5.pdf

Abramowitz and Stegun, “Handbook of Mathematical Functions.”
https://www.math.ubc.ca/~cbm/aands/page_887.htm

L. Bos et al, “An Orthogonality Property of the Legendre Polynomials”, Constr Approx (2016)
https://math.indiana.edu/promotion/files/legendrepaperinprintCA.pdf

LECTURE 16GAUSS QUADRATURE
https://coast.nd.edu/jjwteach/www/www/30125/pdfnotes/lecture16_19v17.pdf

Identities and properties for associated Legendre functions
https://www.mat.univie.ac.at/~westra/associatedlegendrefunctions.pdf

Christoffel–Darboux formula -wikipedia
https://en.wikipedia.org/wiki/Christoffel%E2%80%93Darboux_formula

陰関数の等高線を数値計算で求める

まとめ


関数

を求める問題を考えます。
\(f(x_0, y_0)=0\)を満たす\((x_0, y_0)\)が与えられている場合、\((x_1\equiv x_0+\Delta x, y_1\equiv y_0+\Delta y)\)を求めることを考えていきます。
まず、

の変換を行い、初期値\((x_0, y_0)\)を\((x’_0, y’_0)\)に変換します。
ここで、\(\theta\)は定数で

です。\((x’_0, y’_0)\)を新たな初期値として微分方程式

を解き、\(x_0′\)から微小距離\(h\)だけ発展させた\(x_1’\equiv x_0′ + h\)とその値\(y_1’\equiv y_0′ + \Delta y\)を求めます。
その後、逆変換

を経ることで\(f(x_1,y_1)=0\)を満たす\((x_1,y_1)\)の値を得ることができます。

導出


陰関数

を満たす\((x,y)\)を求めます。ある点\((x_0, y_0)\)において\(f(x_0, y_0)=0\)を満たすことが分かっている場合、
陰関数の全微分より、\((x_0, y_0)\)の近傍で

が満たされています。\(f(x_0,y_0)=0\)であり、\(\Delta x\approx \Delta y \to 0\)ならば、\(O(\Delta^2)\)の項を無視して

が成立する変化のみが許されます。両辺を\(\Delta x\)で割って、\(\Delta x\to 0\)の極限をとれば、

を得ます。すなわち、式(9)を何らかの手段、例えば数値計算で解けば等高線を得ることができます。

発散への対処


さて数値計算によって微分方程式(9)を解いて等高線を描く場合、\(y\)が\(x\)の一価関数になることはほとんどありません。
そのため、式(9)の右辺の分母がゼロになることがあり、微分が発散するため計算を行うことができなくなることがあります。言い換えれば、発散しない範囲しか計算できなくなるということです。
これを回避するためには様々な方法があると思いますが、例えば変数変換を適宜行いながら積分を実行していくことでこの問題を回避することができます。

発散を回避するためのアイデアは以下の通りです。
\(y\)を\(x\)の関数と見たとき、微分が発散するということは、分母がゼロになることを意味しています。
しかし、その点において90°回転させて\(y\)を\(x\)の関数としてみれば微分はゼロに等しいため、問題なく計算を進めることができます。
図示すれば以下の通りになります。

左図は座標を変化せずにそのままの関数ですが、右図のようにある点を計算したいときにその点近くの傾きと等しいように回転した座標から見るようにします。こうすることで発散を回避することができます。

これから変数変換について考えていきます。新たな変数\(x’, y’\)を用いて\(x=x(x’, y’),~y=y(x’, y’)\)(またはその逆関数)と書けているならば、

と書いても良いです。すなわち式(9)を解く代わりに\(x’\)と\(y’\)に関する方程式

を解いても良いのです。具体的に変数変換


を考えます。これは座標を\(\theta\)回転させた座標系への変換です。もちろん、\(\theta=0\)のとき、\(x=x’, y=y’\)が成立し、恒等変換となります。微分方程式は

であり、\((x’, y’)\)の空間で解いていきます。

さて、回転角\(\theta\)は任意ですが、式(15)の右辺が発散しないような\(\theta\)を選ぶ必要があります。
また、数値計算を行う上では、傾きがほとんど0に近いほうが有利です。そのため、方針としては現在の傾きを計算し、その傾き分だけ逆方向に座標回転させることでこれを実現させます。
すなわち、

分だけ回転させます。現在の傾きと座標軸の傾きを一致させることで傾きがほとんどゼロの座標軸に回転させるのです。

プログラム

以下にプログラムを示します。
必要な情報は\(f(x,y)\)の\(x, y\)方向の偏微分\(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}\)の解析解です。関数名”fx”, ”fy”にそれぞれを入力してください。
計算は刻み幅制御の陽的ルンゲクッタ法を用いています。

具体例


具体的に関数\(f(x,y)=x^3-2xy+y^3\)の等高線を考えていきましょう。
関数の値がそれぞれ、\(0.01,1,2,5\)に等しいときの等高線を考えます。
具体的な初期値は、
\(\displaystyle
f(-2,1.3657)\approx 0.01,~f(-2,1.4646)\approx 1,~f(-2,1.5568)\approx 2,~f(-2,1.7977)\approx 5
\)
とします。
具体的に解いたのが左の図です。

”正解”としたgnuplotの等高線表示の結果と同じ結果が得られました。座標回転がうまく働いているようですね。微分の発散を抑えるだけでなく、数値計算量も抑えられているようです。

刻み幅\(h\)が正である場合と負である場合で進む方向が変わります。
また計算の途中で分岐がある場合、どうなるかは保証できません。

また、閉じた等高線が複数ある場合、その全てを計算することは本稿のプログラムだけではできません。初期値を求める別のプログラムを作成する必要があります。
更に、等高線がある値に等しい量を計算することはできません。初期値の値に等しい等高線を引きます。

追記


調べ物をしていたら

照井章, 佐々木建昭著『27.代数関数の陰関数描画について』
http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0941-27.pdf
という素晴らしい資料がありました。

この文献では2つの方法が紹介されており、
1.\(f(x,y_0)=0\)を解いて\(x=x_1,x_2,\cdots,x_n\)を解く
2.陰関数定理から求める
という方法です。資料では、1.の方法を採用していましたが、本稿の方法は2.の方法でした。
1.はもれなく等高線を見つけることができますので、考えられる結果を得たい場合によい方法です。
一方、2.は全てのゼロ点を見つけることができませんが、計算が早くできます。別の方法などで注目したい等高線が分かっている場合に良い方法でしょう。

効率的に解くためには、場合によって使い分けることが一番です。

秋間補間

秋間補間と呼ばれる補間方法があります。

こんな補間結果となります(akimaとmodified akimaの線が秋間補間の結果です)。

秋間補間の特徴は、

  • 与えられた点を必ず通る
  • オーバーシュートが発生しにくい
  • 補間したい点から遠い補間点の影響を受けない
  • 導関数の連続性は保証されない
  • 元のデータ点を大きく逸脱するような補間は苦手(線形的な補間に近い)
  • 極値におけるデータ点など、元の関数の特徴的な値が分かっているときに良い補間結果を与える

ということです。オーバーシュートとは、元のデータが急峻に変化した場合、その前後の補間結果が振動してしまい大きく外れてしまうことを意味します([5]にオーバーシュートの様子を見ることができます)。
ただし、3次スプライン補間のような微分の連続性が保証されません。

それでは計算方法を追っていきます。

導出


秋間補間は1970年に秋間浩によって提案された補間方法です[1]。地図の等高線を表示するために生成されたようです。その後修正された秋間補間が紹介されました[2]。

具体的な計算方法を見ていきます。
補間は補間したい点を中心に前後のデータ点を用いるので、両端では特別な扱いが必要です。それを区別するために領域を3つに分けて考えます。具体的には下の図のような感じです。

領域IIにおける補間


まずは領域IIを考えます。簡単に考えるために、以下の図のように補間に用いるための6点\((x_i, y_i), i=0,5\)があった時、中間の領域である\(x_2\le x\le x_3\)の中の値を補間して求めたいと考えます。

補間は3次関数で行います。すなわち\(x_2\le x\le x_3\)において

の形で補間を行います。係数\(a_0, a_1, a_2, a_3\)は補間に用いる6点のデータ点から計算できて、

の通り求められます。ここで、\(q_1, q_2, m_1,\cdots, m_5\)は


の通り求められます。\(w_1, w_2\)は使用する秋間補間の種類によって変わります。通常の秋間補間では、

の通り計算し[1,3,4]、修正秋間補間においては

で計算します[2,3,4]。

領域I, IIIにおける補間


領域I, IIIについては以下の関係を満たすように、本来は存在しない補間用のデータ点を予測します。
領域I:最も左端のデータ点を\(x_3\)、その右隣のデータ点を\(x_4, x_5\)と置きます。すると、\(x_3\)より左のデータ点\(x_2, x_1\)を関係式(7)より予測します。
領域III:最も右端のデータ点を\(x_3\)、その左隣のデータ点を\(x_2, x_1\)と置きます。すると、\(x_3\)より右のデータ点\(x_4, x_5\)を関係式(7)より予測します。

ここで、\(m_{ij}\)は

を意味します。すると、秋間補間で用いる(本来は存在しない)端のデータ点は、

領域I)

領域III)

と求められます。

プログラム例


上の式をプログラムすれば、このようになるかと思います。
注意してプログラムを作らないとゼロ割が発生します。\(x_i\ne x_j\)を仮定していれば、式(5),(6)の\(w_1, w_2\)を計算する際にこれが発生します。なので、この回避を取り入れてプログラムを作成しています。

計算例


参考文献[4]に秋間補間と修正秋間補間の振る舞いの例があります。それを本稿のプログラムを用いて計算すると同じ振る舞いを示すことが分かります。

それでは、3次スプライン補間と秋間補間、修正された秋間補間を比べて、秋間補間にどのような特徴があるかみましょう。

オーバーシュートの抑圧


本当の答えを定義することは難しいですが、変化する点(\(x=2.8\sim 3.2\))近傍で直線とすれば傾きは5になります。一階導関数はどの補間方法も連続(右から極限と左から極限が一致の意味で)ですが、秋間補間は二階導関数は連続ではありません。3次スプラインは二階導関数まで連続であるという仮定の下での補間なので、一応連続的な変化をしています。

秋間補間では3次スプライン補間よりはオーバーシュートが抑えられているものの、近傍 \((x=2.5, 3.5)\)において外れた値となっています。
修正された秋間補間では、その点が改善されオーバーシュートが発生していないことが分かります。

振動関数の補間

振動関数を補間する例として、元の関数が\(y=\sin(x)\)だったとします。

sin関数

この場合は元の関数が十分滑らかであり、高階微分も連続であることが分かっています。与えられたデータ点も誤差がないので、3次スプラインが有効でしょう。
秋間補間では本来極値をとるはずの\(x=n\pi/2,~(nは整数)\)あたりで値があまり大きくならないことから、かなり局所的な補間である特徴が分かります。すなわち、大きくオーバーシュートしない反面、予測は得意ではないということです。これを回避するためには、元の関数の極値の値など特徴的なデータが必要であることが分かります。

二次元平面の補間は本稿には載せていませんが、[4,6]にどのような違いが表れるか図があるので、興味がある方は参照してみてください。
確かに、本来の目的である地図の等高線作成には大きな効果が期待できそうなことが分かります。おそらく、取得するデータ点も標高など、地形の特徴的な点におけるデータを集めることが主となるので、秋間補間との相性が良いのでしょう。

参考文献


[1]Akima H. A New Method of Interpolation and Smooth Curve Fitting Based on Local
Procedures. Journal of the ACM, 17(4), 589 ― 602. (1970).
doi:10.1145/321607.321609

[2]Akima H. A method of univariate interpolation that has the accuracy of a third-degree
polynomial. ACM Transactions on Mathematical Software, 17(3), 341 ― 366. (1991),
doi:10.1145/114697.116810

[3]
makima 修正Akima区分的3次エルミート内挿 -MathWorks
interp1 1次データ内挿(テーブルルックアップ) -MathWorks

[4]Cleve Moler, Makima Piecewise Cubic Interpolation

[5]スプラインによる補間

[6]格子点に内挿して等値線をひく

連成振動

弦を考えます。弦は重さ\(M\), 長さ\(L\), ばね定数\(K\)を持ちます。

この弦を波動方程式に頼らないで定式化しましょう。
波動方程式は、2次元の波動のうち、1次元分の方向を\(y=y(x, t)\)として表現するため、x軸にとって多価関数となる弦を表現することが出来ません。物性物理学等で格子の振動を考える場合は、結晶を構成する原子が平衡位置から殆ど動かないのでこういった問題は起こりません(起こったとしても無視される位少ないので、大多数の振る舞いを考える場合には問題になりません)。
よって多次元の問題を考えるには波動方程式を元にして考えるのは不適切です。そのため、弦を単なる連成振動として考えて、完全に2次元の問題として扱いたいと思います。

図のような連成振動のモデルを考えます。

質量\(m\)を持つ\(i\)番目の質点が、位置ベクトル\(\mathbf{x}_i\)で指定され、隣り合う質点との間が、ばね定数\(k\)、自然長\(l\)を持つ理想的なばねで繋がれているとします\((i=1,2,\cdots, N)\)。この時、運動方程式を立てると
\begin{eqnarray}
m\frac{d^2 \mathbf{x}_i}{dt^2}&=&
+k(|\mathbf{x}_{i+1}-\mathbf{x}_{i}|-l) \frac{\mathbf{x}_{i+1}-\mathbf{x}_{i}}{|\mathbf{x}_{i+1}-\mathbf{x}_{i}|}
-k(|\mathbf{x}_{i}-\mathbf{x}_{i-1}|-l)\frac{\mathbf{x}_i-\mathbf{x}_{i-1}}{|\mathbf{x}_i-\mathbf{x}_{i-1}|}
\end{eqnarray}
となります。全く同じですが、見通しをよくするため式変形すれば、
\begin{eqnarray}
m\frac{d^2 \mathbf{x}_i}{dt^2}&=& k\Bigl[A_i \mathbf{x}_{i+1}-\bigl(A_i+A_{i-1}\bigr)\mathbf{x}_{i}+A_{i-1} \mathbf{x}_{i-1}\Bigr]
\end{eqnarray}
と表せます。ここで、
\begin{equation}
A_i = \frac{|\mathbf{x}_{i+1}-\mathbf{x}_{i}|-l}{|\mathbf{x}_{i+1}-\mathbf{x}_{i}|},~~(i=1,\cdots, N, A_N = 0)
\end{equation}
と置きました。

さて、こうしたモデルを考えた時に、連成振動として考えた\(N\)分割した1つ1つの小さいばね(”ミクロなばね”と名づけます)とその弦を1つにみなした時の大きなばね(”マクロなばね”と名づけます)の間の関係式を考えましょう。
仮定として、等間隔に分割し、弦の密度は一様であるとします。

この条件の下で、
マクロなばねの重さ\(M\), 自然長\(L\), ばね定数\(K\)

ミクロなばねの重さ\(m\), 自然長\(l\), ばね定数\(k\)
の間の関係式を導出します。

結論として、マクロとミクロなばねの間に成り立つ関係は
\begin{eqnarray}
L&=&N l \\
M&=&N m \\
K &=& k/N
\end{eqnarray}
となります。

自然長

自然長\(L\)の弦を\(N\)分割するので、
\begin{equation}
L=N l
\end{equation}
の関係があります。

重さ

重さ\(M\)の弦を\(N\)分割するので、
\begin{equation}
M=N m
\end{equation}
と書けます。

ばね定数

最後にばね定数の関係を考えましょう。この関係は、弾性エネルギーを考えることで導くことができます。
マクロなばねとして全体を見た場合に長さ\(D\)だけ縮んだ時に蓄えられる弾性エネルギーと
同じだけ縮んだ時にミクロな\(N\)個のばねに蓄えられる弾性エネルギーが同じであってほしいと要請します。すなわち、ミクロなばねの縮みの合計がマクロなばねの縮み\(D\)に等しいとすると、
\begin{equation}
D=N d
\end{equation}
と書けます。
実は、大胆な過程をしないとマクロなばねとミクロなばねの関係を求めることが出来ません。
なぜなら、マクロなばねはいつでも線形応答をするかのように現在考えていますが、ミクロなばねでは場所ごとに違っても良いという様に定式化しています。なのでもともと無理な話なのです。
ではどうやって妥当な関係式を導出すればよいかといえば、ばねが非常にゆっくり運動が行われるときを仮定するのです。ゆっくり動くときの振る舞いは全てのばねが同じ動きをするとして差し支えないでしょう。そうして定式化します。
よって、
\begin{eqnarray}
\frac{1}{2}K D^2 &=& N\cdot \frac{1}{2} k d^2 \nonumber \\
\to ~~ K &=& k/N
\end{eqnarray}
という関係式が導けます。

あとは数値計算を行うだけです。陽的ルンゲクッタ法で解いたプログラムを以下に示します。

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

実行は、

$ gfortran rkf45.f90 main.f90
$ ./a.out

でokだと思います。初期状態や、重さや自然長はinputファイルを見てください。

実行後、動画をgnuplot上で出力したければ、

$ gnuplot
$ call "movie.plt" 100

としてください。100は表示する時間ステップを表します。この引数の値の最大はinputファイルの中のNtまでです。

数値計算結果


実際に数値計算を行う条件として、自然長や重力があると波動方程式と若干異なる振る舞いになってしまうので、対応を見る上ではなくしています。すなわち、自然長\(L=0\)、重力\(g=0\)として計算します。

三角パルス

三角形の波がパルス状に存在するとします。
そして初期状態として、初速度がゼロであることを仮定しましょう。

この場合、波動方程式の考察から、
同じ関数であらわされる右向きの波\(f(x-vt)\)と左向きの波\(f(x+vt)\)の重ね合わせとして書かれることが分かっています。すなわち、解\(y=y(x,t)\)は
\begin{equation}
y(x,t)=f(x-vt)+f(x+vt)
\end{equation}
と書かれます。図で表せば以下のように時間発展していくことが期待されます。

これが本稿の連成振動モデルで実際にそうなっているのかを確かめるのは良い問いかけでしょう。
実際に連成振動を数値的に解くと以下の通りになります。

ゆっくりと時間発展を見ていくと

となります。明らかに右向きに進む波と左向きに進む波として分離します。分離した結果から、初期状態がそれらの波の重ね合わせで表現されていたことが分かるでしょう。

矩形パルス

矩形パルスを考えてみます。実際に数値計算をするとこんな感じです。

結局右向きと左向きの重ね合わせとしてあらわされることも確認できますが、決して単純なそれらの足し合わせで書かれているわけではないことがわかります。この振る舞いは、波動方程式に対応するように連成振動モデルを作っていないことから由来します。波動方程式としてあらわした二階微分の偏微分方程式の場合は多価関数をとることができません。なので縦に弦が伸びる場合、間に質点を置くことができないので質点の位置がいきなり飛ぶことになります。しかし、連成振動モデルではその制約をして初期状態を準備しているわけではありませんのでその違いが出て来ます。図で表せば、以下のようになります。

両者の違いがあるので、波動方程式の結論である
\begin{equation}
y(x,t)=f(x-vt)+f(x+vt)
\end{equation}
が成り立っていない、と解釈することができます。

少しずつ時間発展をみるとこの通りになります。

三角形の波

端まで三角形のパルスであった場合にどのように振る舞うかを見てみましょう。特に面白い振る舞いがあるわけではなく、上のパルスの、右向きと左向きの波の重ね合わせで説明ができるので、単にシミュレーション結果を載せるだけにとどめます。実際に計算してみると以下のようになります。

星形

初期状態が星型である場合も見てみましょう。振る舞いも同じです。

注意点としては、実際にこの振る舞いは起こりえないことに注意しておきましょう。なぜなら、本当の弦で行った場合、重なる場所があると弦の衝突が起こります。しかし、この連成振動モデルでは計算上、衝突は起こりません。

偏りがある三角形のパルス

偏りがある場合も波動方程式では記述することができません。こんな感じの初期状態です。

これを実際にシミュレートすると以下の通りになります。

3次元の場合

3次元の場合も載せておきます。ふるまい自体はそんなに変わるわけではないので、シミュレーション結果を載せるだけにとどめておきます。

弾道計算のコード(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.