時間発展演算子の導出

量子力学における時間発展演算子\(\displaystyle \hat{U}_{(t)}=e^{-i\frac{\hat{H}}{\hbar}t}\)の導出です。
これは、ハミルトニアンは時間依存しない場合に限ります。時間依存するときは時間順序積という概念が出てきます。


シュレディンガー方程式
\(
\begin{align}
i\hbar\frac{\partial}{\partial t}\psi_{(x,t)}=\hat{H}\psi_{(x,t)}
\end{align}
\)
は既知のものとして進めます。

形式的な導出でよく説明されるのは以下の手順です。

シュレディンガー方程式の両辺を\(i\hbar\)で割って、更に波動関数\(\psi_{(x,t)}\)で割り両辺を時刻\(0\sim t\)まで積分します。
すなわち、
\(
\begin{align}
\int_{\psi_{(0)}}^{\psi_{(t)}}\frac{d\psi}{\psi}&=\int_0^t\frac{\hat{H}}{i\hbar}dt^{\prime} \\
\ln\frac{\psi_{(t)}}{\psi_{(0)}}&=\frac{-i}{\hbar}\hat{H}t \\
\end{align}
\)

なので
\(
\displaystyle \psi_{(x,t)}=e^{-i\frac{\hat{H}}{\hbar}t}\psi_{(x,0)}
\)

証明終了,,

(;’;゚;ж;゚;`;)ブホォッ
(; ^ω^)は…

何で波動関数で割ってくれちゃってるんですか?ハミルトニアンは演算子ですよ。形式的にしてもひどくない?

ということで、ちゃんと求めましょう。結果は↑のものと同じになります。不思議ですね。


ハミルトニアンが時間依存しないという条件の下求めます。
まず、時間依存するシュレディンガー方程式より、
\(
\displaystyle i\hbar\frac{\partial}{\partial t}\psi_{(x,t)}=\hat{H}\psi_{(x,t)}
\)

となります。無限小時間\(\Delta t\)に対して、シュレディンガー方程式は、
\(
\displaystyle i\hbar\lim_{\Delta t\rightarrow 0}\frac{\psi_{(x,t+\Delta t)}-\psi_{(x,t)}}{\Delta t}=\hat{H}\psi_{(x,t)}
\)

です。以降\(\lim\)は省略します。

式を変形して、
\(
\displaystyle \psi_{(x,t+\Delta t)}=\left(1-i\frac{\hat{H}}{\hbar}\Delta t\right)\psi_{(x,t)}
\)

また、もう\(\Delta t\)だけ時間を進めると、
\(
\displaystyle \psi_{(x,t+2\Delta t)}=\left(1-i\frac{\hat{H}}{\hbar}\Delta t\right)^2\psi_{(x,t)}
\)

という結果が得られます。

ここで、\(\Delta t\)をN回作用させることを考えます。
すなわち、\(t^{\prime}=N\Delta t\)とおき、\(\Delta t=\frac{t^{\prime}}{N}\)として考えれば、任意の時間\(t^{\prime}\)に対して、
\(
\displaystyle \psi_{(x,t+t^{\prime})}=\left(1-i\frac{\hat{H}}{\hbar}t^{\prime}\frac{1}{N}\right)^N\psi_{(x,t)}
\)

という式が成立します。

無限回、微小区間\(\Delta t\)を動かす操作を考えればいいので、\(N\rightarrow \infty\)として考えれば、
\(
\displaystyle \psi_{(x,t+t^{\prime})}=\lim_{N\rightarrow \infty}\left(1-i\frac{\hat{H}}{\hbar}t^{\prime}\frac{1}{N}\right)^N\psi_{(x,t)}
\)

公式
\(
\displaystyle \lim_{x\rightarrow \pm\infty}\left(1+\frac{a}{x}\right)^x=e^a
\)

より、
\(
\displaystyle \psi_{(x,t+t^{\prime})}=e^{-i\frac{\hat{H}}{\hbar}t^{\prime}}\psi_{(x,t)}
\)

\(t=0\)と置き、\(t^{\prime}=t\)と文字を置き換えれば、
\(
\displaystyle \psi_{(x,t)}=e^{-i\frac{\hat{H}}{\hbar}t}\psi_{(x,0)}
\)

となり、ちゃんと時間発展演算子が導けました。めでたしめでたし。


ハミルトニアンが時間依存する場合、
\(
\displaystyle i\hbar\frac{\partial}{\partial t}\psi_{(x,t)}=\hat{H}_{(t)}\psi_{(x,t)}
\)

となります。無限小時間\(\Delta t\)に対して、シュレディンガー方程式は、
\(
\displaystyle i\hbar\lim_{\Delta t\rightarrow 0}\frac{\psi_{(x,t+\Delta t)}-\psi_{(x,t)}}{\Delta t}=\hat{H}_{(t)}\psi_{(x,t)}
\)

次の時刻では
\(
\displaystyle \psi_{(x,t+2\Delta t)}=
\left(1-i\frac{\hat{H}_{(t+\Delta t)}}{\hbar}\Delta t\right)\left(1-i\frac{\hat{H}_{(t)}}{\hbar}\Delta t\right)\psi_{(x,t)}
\)

となります。2乗にはなってくれません。
と、どんどん作用させていくわけです。時間依存しない時と明らかに違ったものになりますが、詳しくは次の機会のお話で。

ルンゲクッタ法の説明と刻み幅制御

ルンゲクッタ法(Runge-Kutta method、RK法)とは?
僕の知る限りの知識で紹介します。

微分方程式
\(
\displaystyle \frac{d^2}{dt^2}x=t
\)

を倍精度実数、刻み幅制御で解く場合のコードはこれです。複素数の場合は本文参照。

もくじ

  1. ルンゲクッタ法の系統
  2. 陽的ルンゲクッタ法の段数と次数について
  3. 誤差について
  4. Butcher tableによるルンゲクッタ法の記述
  5. 埋め込まれた陽的ルンゲクッタ法
  6. ルンゲ-クッタ-フェールベルグ法による刻み幅の自動制御(適応刻み幅制御)
  7. ルンゲ-クッタ-フェールベルグ法による刻み幅の自動制御のプログラム(実数/複素数 で変数を 自動/指定 で解く 1/2 階微分方程式のプログラム)
  8. 4倍精度ルーチン
  9. 不連続な点を含む場合
  10. 刻み幅制御のベンチマーク(振り子)
  11. 陽的ルンゲクッタ法の導出
  12. 参考文献

理論はいいから4次ルンゲクッタ法の数値計算ではどうやるかだけ教えてくれ!という方は4次ルンゲクッタ法へどうぞ。

ルンゲ-クッタ法の系統


ルンゲ-クッタ法は微分方程式の数値計算解を得るための手法のことを指します。

通常の参考書で微分方程式を解くために良く紹介されているのは、オイラー法と中点法、4次ルンゲクッタ法でしょうか。
オイラー法も中点法も4次ルンゲクッタ法も、”陽的ルンゲクッタ法” と呼ばれる枠組みの1つです。

オイラー法は正確には “陽的1段1次ルンゲクッタ法” と呼ばれ、
中点法は “陽的2段2次ルンゲクッタ法”、
4次ルンゲクッタ法(RK4)は、”陽的4段4次ルンゲクッタ法” と呼ばれています。

“段”と”次”とはなんなんでしょう?それは、
計算の大変さ(段)と、計算の正確さ(次)
です。”“の値が小さければ小さいほど計算時間が少なくて済みますし、”“の値が高ければ高いほど計算が正確です。

オイラー法は1という計算コストで正確さ1が得られますし、RK4は4という計算コストで正確さ4が得られます。

4次ルンゲクッタ法が使われる理由

理由は実装が簡単でそれなりの精度を持つから。です。

陽的ルンゲクッタ法において、pという計算コスト(p段)で、pより大きな正確さq(q次)を得ることはできません
Derivation of Runge–Kutta methodsによれば、
\(q\)次の正確さ(q次のオーダー)を得たい場合、最低限必要な段数\(p_{\mbox{min}}(q)\)は

という関係にあります。

ここで注目するべきは4次の時までは計算コストに比例して計算精度が上がっていきます。
しかし、4次以上では、計算コストの増加と計算の正確さが見合わなくなっていきます。故に
計算効率が良いのは4次だろう、
と予想できます。
また、重要な理由として、4次ルンゲクッタ法に現れる係数が5次以降と比べて圧倒的にシンプルであることが挙げられます。4次では\(0,1/2,1/6\)程度の係数だけが使われ、プログラム作成時の入力ミスがほぼ生まれません。しかし、5次では\(28561/56430, -7200/2197\)といった係数が数多く出てきます。
これらの理由から,4次ルンゲクッタ法(RK4)が数値計算科学の世界でよく使われるのです。

陽的ルンゲクッタ法に限って言えばプログラムの実装が非常に簡単であることが挙げられます。陰的ルンゲクッタ法と呼ばれるアルゴリズムもあり、これは陽的ルンゲクッタよりも優れていますが、計算量が多くなり、若干複雑なアルゴリズムになります。

Q. オイラー法もものすごく細かい分点を取れば、その計算効率はRK4と同じなんじゃないの?
A. 刻み幅の乗数で効いてくるのでそうではありません。高次を使っても計算が信頼できるのであれば、大きなステップで進めるほうが早いです。例え、目標精度への計算時間が同じだとしても、計算機の有効桁数によって否定されてしまいます。
RK4で典型的にとられる時間ステップの間隔は、おおよそ\(10^{-2}\sim 10^{-4}\)程度であり、RK4のエラーのオーダーは\(O(h^5)\)です。
そして、科学計算で使う際の有効桁数は倍精度型で16桁です。
1ステップ当りの誤差は\(h\)の減少に伴い、解が\(h^4\)の早さで収束していく、と言えます。
だから16桁の計算では\(h=10^{-1}\to h=10^{-4}\)に変化させる時、誤差は\(O(h^5)=c 10^{-5}\to c 10^{-20}\)と変化します。
おおよそ\(c\approx 1\)と見積もれば、(有効桁数16桁を多少超えてしまいますが、)有効桁数いっぱいまで正しい値が出るであろうと期待できます。

これに対し、オイラー法で同じことをするには\(h\)を\(10^{-16}\)にしなくてはいけません。
\(t\)の値が\(10^{-16}\)変わった時に、桁落ちの問題を回避できるほど関数\(x\)の値に変化が生まれるか?
が問題になり、まぁそんな急激な変化は生まれないでしょう、と予想できます。これでは桁落ちの問題を回避するほどの変化は到底望めません。

よって計算の効率と有効桁数の限界から、RK4なのです。

また、あまりに高次の方法を使うとルンゲ現象に代表される不安定性といったことが起こるかもしれません。
高次は高精度という意味とイコールではないことに注意しましょう。この高次≠高精度については等間隔の分点における数値積分の時に書いたので気になる方はどうぞ。。

誤差について


4次ルンゲクッタ法の、1ステップ当りの誤差は\(h^5\)に比例,すなわち\(O(h^5)\)です。
しかし、通常は細かいルンゲクッタ法を何度も繰り返して計算します。
区間\([a,b]\)を刻み幅\(h\)の4次ルンゲクッタ法で\(N\)回のステップを繰り返し計算する場合、誤差は
\(
\displaystyle N\cdot O(h^5)=\frac{b-a}{h}\cdot O(h^5)=O(h^4)\)

となります。よって、\(N\)回繰り返すような計算では、オーダーが1つ落ちる事に注意しましょう。


スポンサーリンク


陽的ルンゲクッタ法の段数と次数について


さて、ここまで“段”は計算コスト、で“次”は計算の正確さ、という曖昧な表現でしたが、その表現をちゃんと知りましょう。
段と次を知るためにはルンゲクッタ法の計算方法を知る必要があります。
具体例を載せます。
\(
\displaystyle \frac{dx}{dt}=f(t,x)
\)

の、\(t_{n}\rightarrow t_{n}+h\ \ (=t_{n+1})\)における初期値問題に対する、
陽的1段1次ルンゲクッタ法(オイラー法)の計算スキームは、
\(
x_{n+1}=x_{n}+h\cdot f(t_{n},x_{n})
\)
です。

陽的4段4次ルンゲクッタ法(RK4)の計算スキームは、
\(
\begin{align}
k_1&=f(t_n, x_n) \\
k_2&=f(t_n+h/2, x_n+h k_1/2) \\
k_3&=f(t_n+h/2, x_n+h k_2/2) \\
k_4&=f(t_n+h, x_n+h k_3) \\
x_{n+1}&=x_{n}+{(k_1+2k_2+2k_3+k_4)}h/6
\end{align}
\)
として与えられます。

一般的に、陽的s段のルンゲクッタ法とは
\(
\begin{align}
g_i&=x_n+h\sum_{j=1}a_{i,j}k_j\ \ \ (j\lt i, \ i=1,2,…,s) \\
k_i&=f(t_n+c_ih, g_i) \\
x_{n+1}&=x_n+h\sum_{i=1}^s b_ik_i
\end{align}
\)
として書けます。
ここで行列形式で与えられる係数\(a_{i,j}, b_{i},c_{i}\)によって、そのs段ルンゲクッタ法が持つ次数が決められます。段数はここから由来します。

点\((t_n, x(t_n))\)周りで関数をテーラー展開し、その関数が点\((t_n+h\ \ (=t_{n+1}), x(t_{n+1}))\)で作る点を近似解とするのがルンゲクッタ法です。
故に、\(x(t_{n+1})\)は、
\(
\begin{align}
x(t_{n+1})=x(t_n)+\left.\frac{h}{1!}\frac{dx}{dt}\right|_{t=t_n}+\left.\frac{h^2}{2!}\frac{d^2x}{dt^2}\right|_{t=t_n}+\left.\frac{h^3}{3!}\frac{d^3x}{dt^3}\right|_{t=t_n}+\left.\frac{h^4}{4!}\frac{d^4x}{dt^4}\right|_{t=t_n}+…
\end{align}
\)
と書けます。
ここで、テーラー展開としてどの程度一致させて\(x(t_n+h)\)を決定するか?を表すのが次数に当たります。

言葉で書くなら、

局所的な精度は,高次項が解のテイラー(Taylor)展開とどの程度マッチしているかによって測定される.オイラー法は一次精度なので,誤差は一次高く\(h^2\)というベキ乗から始まる. ~NDSolveの”ExplicitRungeKutta”メソッドより

ともあります。

Butcher tableによるルンゲ-クッタ法の記述


行列形式で与えられるルンゲ-クッタ法での係数\(a_{i,j}, b_{i},c_{i}\)は何なのか?
具体的に記述してみましょう。
オイラー法(1段1次)はもっとも単純で、係数は
\(
\begin{align}
a_{1,1}&=0 \\
b_{1}&=1 \\
c_{1}&=0
\end{align}
\)
です。これを一般的な表記法の式に当てはめれば、
\(
\begin{align}
g_1&=x_n+h a_{1,1}k_1 \\
k_1&=f(t_n+c_1h, g_1) \\
x_{n+1}&=x_n+h b_1k_1
\end{align}
\)
となります。

中点法は、
\(
\begin{align}
a_{1,1}&=0 \\
a_{1,2}&=0 \\
a_{2,1}&=1/2 \\
a_{2,2}&=0 \\
b_{1}&=0 \\
b_{2}&=1 \\
c_{1}&=0 \\
c_{2}&=1/2
\end{align}
\)
という組で与えられます。

この係数行列の組はまとめてButcher tableと呼ばれる表記をするのが便利です。

これは、\(a_{i,j}, b_{i},c_{i}\)を
Buchertable
としてまとめて書く表記法です。

再び、オイラー法はButcher tableで書くと
euler_butcher
とまとめて書くことができます。
中点法は
midpoint_butcher
RK4は
rk4_butcher
です。

高次のルンゲクッタ法(10,12,14次)


4次、5次…とずっとあるわけです。
こんなページがありました[3]。
High-Order Explicit Runge-Kutta Methods
この上のページには

  • 17段10次(8次が埋め込まれてる)
  • 25段12次(10次が埋め込まれてる)
  • 35段14次(12次が埋め込まれてる)

といったButcher tableにおける係数の値が書かれています。埋め込まれてる、の意味は次の節で説明します。
ただし、上のページのbutcher tableは
highorderbutcher2
となっているので注意が必要です。

埋め込まれた陽的ルンゲクッタ法


埋め込まれた“という表現が出てきたのでその説明を行いましょう。
日本語では『埋め込まれた陽的ルンゲクッタ法』、英語では『embedded explicit runge-kutta method』と呼ばれるものがあります。
これは、p段q次陽的ルンゲクッタ法を作ったら、別の次数の陽的ルンゲクッタ法も、係数行列\(a_{i,j}, c_{i}\)を使って作れるじゃありませんか!
というものです。

Butcher tableは、この場合extended Butcher tableと呼ばれ、こういう形式で書かれます。
embedded_rk
この埋め込まれたルンゲクッタ法のいいところは、

  1. 計算誤差の評価ができる
  2. 刻み幅を自動的に制御できる、適応刻み幅制御。(応用として。)

という点です。ルンゲクッタ法によって得られた解が真の解とどのくらい違っているのか?が評価できるんです。

例えば、4次のルンゲクッタ法を使って得られた解\(x^{(4)}(t)\)と5次のルンゲクッタ法を使って得られた解\(x^{(5)}(t)\)があったとします。
もしも、\(x^{(4)}(t)\)と解\(x^{(5)}(t)\)の解の差を調べ、その差が無かったらその数値計算解は真の解に限りなく近い、と判断することができ、差が大きかったらその解は真の解から離れていて、数値計算の精度が足らない、と判断することができます。どちらも1つだけの解では出来なかったことです。

精度が足らない場合、刻み幅を小さくすれば精度が上がります。また逆に、精度が十分に足りている場合、刻み幅を大きくし、計算時間を減らすことができます。
これが適応刻み幅制御なのです。

違った次数のルンゲクッタ法を、まるで別々に計算してもいいのですが、そうすると計算時間が単純に考えておおよそ2倍になります。
解を評価するために2倍の計算時間が必要というのは良くない計算効率です。
そこで考えられたのが埋め込まれたルンゲクッタ法なのです。

具体例を挙げましょう。
一番簡単な埋め込まれたルンゲクッタ法は、ホイン法と呼ばれています。
butcher_heun
1行目は2次のオーダーを持ち、2行目は1次のオーダーを持ちます。

また、4次と5次を持つ埋め込まれたルンゲクッタ法は、ルンゲ-クッタ-フェールベルグ(Runge-Kutta-Fehlberg)法と呼ばれています。
その埋め込まれたルンゲクッタ法は、
butcher_rkf45
と書かれます。1行目は5次のオーダー、2行目は4次のオーダーを持ちます。

ルンゲ-クッタ-フェールベルグ法による刻み幅の自動制御(適応刻み幅制御)


さて、次数の違う2つのルンゲクッタ法を用いて、適応刻み幅制御を行いたいと考えます。
刻み幅を制御するにあたって、適当に精度良かったから2倍にしてもまだ大丈夫だろ、とか差が大きすぎるから刻み幅半分にしよう、ということをやってはいけません。
適当にやったら計算時間が余計にかかり、精度が良くない変な結果が得られます。

[5~9]によれば、ルンゲ-クッタ-フェールベルグ法において区間\(i\)での最適な刻み幅\(h’\)は区間\(i\)の誤差評価の結果を使って、
\(
\displaystyle h’=\left(\frac{\varepsilon h}{2|x^{(5)}_{i}-x^{(4)}_{i}|}\right)^{1/4}
\)

と予想できます。もちろん、この\(h’\)は区間\(i\)の最適な幅ですが、関数に劇的な変化は無いだろうとして、次の区間の計算の一番初めに用いる刻み幅を推定するのです。
なので、\(i+1\)番目の計算区間では、計算するときはこの\(h’\)の値を使えばいいんです。
(ちなみに、m次ルンゲクッタ法の場合では
\(
\displaystyle h’=\left(\frac{\varepsilon h}{2|x^{(m+1)}_{i}-x^{(m)}_{i}|}\right)^{1/m}
\)

と推測されます。)

詳しい理由は分かりませんが、5次オーダーではなく、4次です。5次のオーダーは誤差評価のためだけに用いられているようです。
ルンゲ-クッタ-フェールベルグ法の計算スキームは[7]に詳しく書かれています。
日本語訳して、その計算スキームを書けば下のようになります。
刻み幅制御rkf45_4order

ルンゲ-クッタ-フェールベルグ法による刻み幅の自動制御(プログラム)


講義のレポート等の宿題で使うのは僕の意向と異なるので使用はご控えください。
研究目的、趣味、確かめの場合はミスがあるかもしれないことを念頭に置いたうえならば使用と改変をご自由にしてください。
このプログラムを使用して生じた責任は取りません。

fortran90によるプログラムです。ほぼ上の説明をそのままプログラミングしたものです。

実例を8つ載せます。

  • 実数で、時間の刻みを自動で決める1階微分方程式のプログラム
  • 実数で、時間の刻みを指定する1階微分方程式のプログラム
  • 複素数で、時間の刻みを自動で決める1階微分方程式のプログラム
  • 複素数で、時間の刻みを指定する1階微分方程式のプログラム
  • 実数で、時間の刻みを自動で決める2階微分方程式のプログラム
  • 実数で、時間の刻みを指定する2階微分方程式のプログラム
  • 複素数で、時間の刻みを自動で決める2階微分方程式のプログラム
  • 複素数で、時間の刻みを指定する2階微分方程式のプログラム

です。
時間の指定する、というのは、指定した時刻の間は自動で刻み幅を制御しながらやりますが、出力するのは指定した時間のみ、という意味です。

下のモジュールRKmodと共にコンパイルし、実行してください。

  • 実数で、時間の刻みを自動で決める1階微分方程式のプログラム

    下のプログラムは一階微分方程式
    \(
    \displaystyle \frac{dy}{dx}=y\cos(x)
    \)
    を初期条件
    \(
    \left. y(x)\right|_{x=0}=1
    \)
    の条件下で解くものです。
    解析解は\(y(x)=\exp(\sin(x))\)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(y(x)\), 刻み幅, 積分\(\int_{x_0}^{x}y'(x)dx\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    expsin1
    が出力されます。紫色の線は解析解\(\exp(\sin(x))\)を表しています。関数の変化が激しいところで刻み幅が自動的に小さくなっているのが分かるでしょう。

  • 実数で、時間の刻みを指定する1階微分方程式のプログラム

    下のプログラムは一階微分方程式
    \(
    \displaystyle \frac{dy}{dx}=y\cos(x)
    \)
    を初期条件
    \(
    \left. y(x)\right|_{x=0}=1
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(\sin(x))\)となります。
    ただし、出力の間隔を指定します。

    上のプログラムを実行すると、
    “fort.11″に
    位置, 関数\(y(x)\), 積分\(\int_{x_0}^{x}y'(x)dx\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    expsin2
    が出力されます。紫色の線は解析解\(\exp(\sin(x))\)を表しています。
    ここで、時間の間隔は等間隔に荒く指定しました。ですが、点と点の間は自動で刻み幅制御が行われているので、精度は指定したものとなっています。等間隔でなくても構いません。

  • 複素数で、時間の刻みを自動で決める1階微分方程式のプログラム

    下のプログラムは一階微分方程式
    \(
    \displaystyle \frac{dy}{dx}=\left(\frac{1}{2x}+i\right)y
    \)
    を初期条件
    \(
    \left. y(x)\right|_{x=1}=\exp(i)
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(ix)\sqrt{x}\)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(Re\{y(x)\}\), 関数の虚部\(Im\{y(x)\}\), 刻み幅、積分の実部\(Re\{\int_{x_0}^{x}y'(x)dx\}\), 積分の虚部\(Im\{\int_{x_0}^{x}y'(x)dx\}\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    expsqrt1
    が出力されます。水色、オレンジの点は数値計算で解いた解の実部、虚部を表し、紫、緑の実線は解析解の実部、虚部を表します。

  • 複素数で、時間の刻みを指定する1階微分方程式のプログラム

    下のプログラムは一階微分方程式
    \(
    \displaystyle \frac{dy}{dx}=\left(\frac{1}{2x}+i\right)y
    \)
    を初期条件
    \(
    \left. y(x)\right|_{x=1}=\exp(i)
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(ix)\sqrt{x}\)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(Re\{y(x)\}\), 関数の虚部\(Im\{y(x)\}\), 積分の実部\(Re\{\int_{x_0}^{x}y'(x)dx\}\), 積分の虚部\(Im\{\int_{x_0}^{x}y'(x)dx\}\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    expsqrt2
    が出力されます。水色、オレンジの点は数値計算で解いた解の実部、虚部を表し、紫、緑の実線は解析解の実部、虚部を表します。

  • 実数で、時間の刻みを自動に決める2階微分方程式のプログラム

    下のプログラムは二階微分方程式
    \(
    \displaystyle \frac{d^2y}{dx^2}= \frac{dy}{dx}\cos{x}-y\sin(x)
    \)
    を初期条件
    \(
    \begin{align}
    \left. y(x)\right|_{x=0}=1 \\
    \left. \frac{dy(x)}{dx}\right|_{x=0}=1
    \end{align}
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(\sin{x}),\frac{dy(x)}{dx}=\cos{x}\exp(\sin{x})\)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(y(x)\), 導関数\(\frac{dy(x)}{dx}\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    expsinx_d2
    が出力されます。水色、オレンジの点は数値計算で解いた,解\(y(x)\)、導関数\(\frac{dy(x)}{dx}\)を表し、紫、緑の実線は解析解、とその導関数を表します。

  • 実数で、時間の刻みを指定する2階微分方程式のプログラム

    下のプログラムは二階微分方程式
    \(
    \displaystyle \frac{d^2y}{dx^2}= \frac{dy}{dx}\cos{x}-y\sin(x)
    \)
    を初期条件
    \(
    \begin{align}
    \left. y(x)\right|_{x=0}=1 \\
    \left. \frac{dy(x)}{dx}\right|_{x=0}=1
    \end{align}
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(\sin{x}),\frac{dy(x)}{dx}=\cos{x}\exp(\sin{x})\)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(y(x)\), 導関数\(\frac{dy(x)}{dx}\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    expsinx_d2_d
    が出力されます。水色、オレンジの点は数値計算で解いた,解\(y(x)\)、導関数\(\frac{dy(x)}{dx}\)を表し、紫、緑の実線は解析解、とその導関数を表します。

  • 複素数で、時間の刻みを自動で決める2階微分方程式のプログラム

    下のプログラムは二階微分方程式
    \(
    \displaystyle \frac{d^2y}{dx^2}=\left(-\frac{1}{4x^2}+\frac{i}{x}-1\right)y
    \)
    を初期条件
    \(
    \begin{align}
    \left. y(x)\right|_{x=1}=\exp(i) \\
    \left. \frac{dy(x)}{dx}\right|_{x=1}=\left(\frac{1}{2}+i\right)\exp(i)
    \end{align}
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(ix)\sqrt{x},\frac{dy(x)}{dx}=\left(\frac{1}{2x}+i\right)\exp(ix)\sqrt{x} \)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(Re\{y(x)\}\), 関数の虚部\(Im\{y(x)\}\), 微分の実部\(Re\{\frac{dy(x)}{dx}\}\), 微分の虚部\(Im\{\frac{dy(x)}{dx}\}\), 刻み幅
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    sqrtxexpix_d2
    が出力されます。水色、オレンジの点は数値計算で解いた,解\(y(x)\)の実部、解\(Re\{\frac{dy(x)}{dx}\}\)の実部を表し、紫、緑の実線は解析解の実部、とその導関数の実部を表します。

  • 複素数で、時間の刻みを指定する2階微分方程式のプログラム

    下のプログラムは2階微分方程式
    \(
    \displaystyle \frac{d^2y}{dx^2}=\left(-\frac{1}{4x^2}+\frac{i}{x}-1\right)y
    \)
    を初期条件
    \(
    \begin{align}
    \left. y(x)\right|_{x=1}=\exp(i) \\
    \left. \frac{dy(x)}{dx}\right|_{x=1}=\left(\frac{1}{2}+i\right)\exp(i)
    \end{align}
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(ix)\sqrt{x},\frac{dy(x)}{dx}=\left(\frac{1}{2x}+i\right)\exp(ix)\sqrt{x} \)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(Re\{y(x)\}\), 関数の虚部\(Im\{y(x)\}\), 微分の実部\(Re\{\frac{dy(x)}{dx}\}\), 積分の虚部\(Im\{\frac{dy(x)}{dx}\}\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    sqrtxexpix_d2e
    が出力されます。水色、オレンジの点は数値計算で解いた,解\(y(x)\)の実部、解\(Re\{\frac{dy(x)}{dx}\}\)の実部を表し、紫、緑の実線は解析解の実部、とその導関数の実部を表します。

4倍精度ルーチン

置いておきます。
メインプログラム

program main
  use qRKmod
  implicit none
  integer::i,N,Nx,info
  real*16::dx,tol,xbound,x0,tx
  real*16,allocatable::x(:)
  complex*32,allocatable::y(:),ss(:)
  complex*32::rkfc
  external::rkfc
 
  N=2
  allocate(y(1:N))
 
  x0=1q0; dx=1000q0; Nx=2
  allocate(x(0:Nx-1))
  do i=0,Nx-1
     x(i)=x0+dble(i)*dx
  enddo

  y(1)=exp(cmplx(0q0,1q0,kind=16))
  y(2)=cmplx(0.5q0,1q0,kind=16)*y(1)

  call qrk_preparation("rkf45")

  tol=1q-18
  write(10,'(3e30.17e4)')x(0),real(y(1)),imag(y(1))
  do i=1,Nx-1
     info=0; tx=x(i-1); xbound=x(i)
     call qrkf45_e(rkfc,tx,y,xbound,info,tol)
     if(info.eq.-1)write(6,'(A,f10.5,A,f10.5)')"strange point between ",x(i-1)," and ",x(i)
 
     write(10,'(3e30.17e4)')x(i),real(y(1)),imag(y(1))
  enddo
  write(6,*)x(0),y(1)
 
  call qrk_deallocation("rkf45")
 
  stop
end program main
!----------------
function rkfc(N,x,y,s)
  implicit none
  integer,intent(in)::N,s
  real*16,intent(in)::x
  complex*32::y(1:N)
  complex*32::rkfc
 
  rkfc=0q0
  if(s.eq.1)then
     rkfc=y(2)
  elseif(s.eq.2)then
     rkfc=(-1q0/(4q0*x*x)+cmplx(0q0,1q0,kind=16)/x-1q0)*y(1)
  else
     write(6,*)"unknown s, program stop"
     stop
  endif
 
  return
end function rkfc
スポンサーリンク

不連続な点を含む場合


不連続な点を含む場合、境界条件を指定しないと解くことはできません。

さて、ここで微分方程式
\(
\begin{eqnarray}
\frac{dy}{dx}=
\left\{
\begin{aligned}
0\;\;(x\le 0)\\
1\;\;(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
\)

を初期条件\(y(-0.095)=0\)の下で考えます(意図的に境界条件は考えません)。
これを4次ルンゲクッタと適応刻み幅ルンゲクッタで解いてみましょう。
プログラム上ではそのまま解くことが出来ます。
実際に解かせてみますと、
rk4_rkf451_disc
となり、2つの結果(赤:4次ルンゲクッタ、緑:刻み幅制御ルンゲクッタ)は異なってしまいます。これは1階微分の不連続性のため発生します。
不連続点\(x=0\)で関数\(y(x)\)に境界条件を指定しない限り、どちらも正しい解なのです。

さて、なぜこんなことが発生するのでしょうか?以下のように問題を表すことにします。

不連続な点を含む1階の微分方程式を考えます。
ここで不連続、という意味は関数\(y(x)\)の一階微分が、点\(x’\)で
\(
\displaystyle \left. \frac{dy}{dx}\right|^{x=x’+0}_{x=x’-0}=a,\;\;\;(a\ne 0)
\)

であるような点を指しているとします。

上記の例題を考えてみましょう。
上記の例題では\(a=1\)です。微分方程式を解析的に解いてみますと、
\(
\begin{eqnarray}
y(x)=
\left\{
\begin{aligned}
C_0\;\;(x\le 0)\\
x+C_1\;\;(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
\)

となります。ここで\(C_0, C_1\)は定数です。
\(C_0, C_1\)は\(y(x)\)が解きたい問題の境界条件によって決まります。

例えば、\(y(x)\)は全領域に対して繋がっている、という条件を課しましょう。この場合、不連続点\(x’\)で
\(
\displaystyle \left. y(x)\right|^{x=x’+0}_{x=x’-0}=0
\)

という境界条件を満たさなければなりません。この条件を課すと、\(C_1=C_0\)となり、初めて関数\(y(x)\)を一意に決めることが出来ます。

1階微分方程式を解く場合、適応刻み幅制御では関数\(y(x)\)は計算領域内で繋がっている事が課されています。しかし、4次ルンゲクッタではその条件は課されません。\(C_1\)の値は初期条件に依存し、一意に関数が決まりません。
どちらが悪いという話ではありません。

通常は適応刻み幅でも、4次ルンゲクッタでも\(y(x)\)にどこか連続ではない変な点がある場合、その点で区間を別々に分けて解きます。その後、境界条件に従って値を調節して全体の関数を構成します。

ベンチマーク用


微分方程式の解法がどれくらい正しそうかのベンチマーク問題として振り子(角度が大きい時)を考えましょう(振り子の詳しい解説はこちら)。
以下の\(\omega=1\)としたときの運動方程式
\(
\displaystyle \frac{d^2\theta}{dt^2}=-\sin\theta
\)


\(
t=0,\; \theta=0; \frac{d\theta}{dt}=1.9\cdots (k=0.95)
\)
の初期条件の下解いた場合、
1周期\(T\)は
\(
T=4K(0.95)=10.360044923498004876778\cdots
\)
となります。
この値はwolfram alphaから求めました。
4*EllipticK[0.95*0.95] wolfram alpha

刻み幅制御を行い、45000周期目の値を考えます。45000周期目は時刻
\(
T_{45000}=466202.0215574102\cdots
\)
です。刻み幅制御による精度を\(10^{-12}\)に設定し、数値計算を行わせます。
プログラムはこちら

すると実行結果として”fort.10″に

0.4662020113E+006  -0.2103356901E-001   0.1899883579E+001   0.4224109363E-002
0.4662020155E+006  -0.1300808922E-001   0.1899955473E+001   0.4223843994E-002
0.4662020198E+006  -0.4982881533E-002   0.1899993468E+001   0.4223658015E-002
0.4662020240E+006   0.3042061693E-002   0.1899997567E+001   0.4223520921E-002

というデータが出力されます。
1列目が時刻\(t\)、2列目が\(\theta(t)\),3列目が\(\frac{d\theta(t)}{dt}\),4列目が刻み幅\(h\)です。
1回のステップでの要求精度12桁に対し、最終的な結果は8桁まで正しい値を出しています。
この時、計算回数は77,852,488回\(~10^{8}\)回行われているので、最終的な結果として4桁は少なくとも正しいと考えられます。
4桁合っていればいい状況で8桁もあっているのは、被積分関数が変な振る舞いをしないからでしょう。

また、60000周期で確認してみると(\(T_{60000}=621602.695409880292\cdots\))

0.6216026888E+006  -0.1531918479E-001   0.1899938246E+001   0.4223959920E-002
0.6216026930E+006  -0.7293808996E-002   0.1899986003E+001   0.4223717084E-002
0.6216026973E+006   0.7312355417E-003   0.1899999862E+001   0.4223582630E-002
0.6216027015E+006   0.8756011575E-002   0.1899979827E+001   0.4223462029E-002

です。
1回のステップでの要求精度12桁に対し、最終的な結果は8桁まで正しい値を出しています。
この時、計算回数は103,803,513回\(~10^{8}\)回行われているので、最終的な結果として4桁は少なくとも正しいと考えられます。
・・・まだまだ大丈夫そうですね。

少し特殊な初期条件(\(k=1\))でやってみましょう。
この\(k=1\)は、振り子の質点がちょうど真上に来て静止する非常に不安定な状態です。
何秒間静止していられるか試してみましょう。刻み幅の制御等は上記条件と同じです。
横軸に時間\(t\)、縦軸に\(\theta(t)\)を取った時のグラフです。
ellk1
すぐに破綻しました。正しい値は\(t=20\)位までですかね。これは、機械的な誤差があることによって不安定な平衡点からほんの少しだけ外れてしまったからです。だからカオスとかを考えるときとかは注意しなければなりません。

陽的ルンゲクッタ法の導出


ルンゲクッタ法の導出は煩雑です。単に複雑なだけです。
導出過程について詳しく述べられているページは、早川様が記述された以下のpdfを参照すると良いと思います。
Runge-Kutta法についてのノート(早川尚男)
計算過程を含め記述されているので分かりやすいです。

参考文献


[1]Derivation of Runge–Kutta methods
[2]NDSolveの”ExplicitRungeKutta”メソッド
[3]High-Order Explicit Runge-Kutta Methods
[4]List of Runge–Kutta methods
[5]Runge-Kutta-Fehlberg Method (RKF45)
[6]Runge-Kutta-Fehlberg method
[7]Lecture:13Runge-Kutta-Fehlberg Method
[8]GPU acceleration of Runge Kutta-Fehlberg and its comparison with Dormand-Prince method
[9]William H. Pressら著『ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ– 』(1993)

↑この本が一番有益だと思います。Fortran版もありますが、英語しかありません。ちなみに、英語で,若干古いバージョンでいいのならば
公式のホームページ
Numerical Recipes in C (1992)
Numerical Recipes in Fortran 77 and Fortran 90 (1992, 1996)
で無料で公開されています。

運動量演算子の導出

量子力学の運動量演算子\(\hat{p}\)が
\(
\displaystyle \hat{p}=-i\hbar \frac{d}{d x}
\)

であらわされることを導出します。
ここでの”導出”は運動量演算子として妥当なものを導出する、と言った方がいいかもしれません。

平行移動演算子\(\hat{T}\)を出発点とします。
この平行移動演算子はある波動関数\(\psi_{(x)}\)を平行移動させる演算子で、これを作用させると
\(
\hat{T}_{(\Delta x)}\psi_{(x)}=\psi_{(x+\Delta x)}
\)

となる演算子です。

よく説明される、シュレディンガー方程式の平面波解からの出発はしません

この\(\hat{T}\)の持つであろう性質を考え、その特性から\(\hat{T}\)の具体的な形を推定し、運動量演算子を導きます。

平行移動演算子\(\hat{T}\)が満たすべき性質


当たり前と思われる性質から、平行移動演算子に関して以下の4つが言えます。

  1. 波動関数\(\psi_{(x)}\)が1に規格化済みであれば、平行移動した波動関数\(\psi_{(x+\Delta x)}\)もまた規格化されているはずである。
    すなわち、
    \(
    \begin{align}
    \langle\psi_{(x)}|\psi_{(x)}\rangle &= \langle\psi_{(x)}|\hat{T}^{\dagger}_{(\Delta x)}\hat{T}_{(\Delta x)}|\psi_{(x)}\rangle \\
    &= \langle\psi_{(x+\Delta x)}|\psi_{(x+\Delta x)}\rangle \\
    &\rightarrow \hat{T}^{\dagger}_{(\Delta x)}\hat{T}_{(\Delta x)}=\mathbf{1}
    \end{align}
    \)
    となり、3本目の式から\(\hat{T}\)はユニタリー演算子であることがわかります。
  2. \(\Delta x\)進めた後、\(\Delta x^{\prime}\)を進める操作は\(\Delta x + \Delta x^{\prime}\)進める操作に等しいはずである。
    すなわち、
    \(
    \hat{T}_{(\Delta x^{\prime})}\hat{T}_{(\Delta x)}=\hat{T}_{(\Delta x+\Delta x^{\prime})}
    \)
  3. 負の方向への平行移動\(\hat{T}_{(-\Delta x)}\)は正の方向への平行移動\(\hat{T}_{(\Delta x)}\)の逆のはずである。
    \(
    \hat{T}_{(-\Delta x)}=\hat{T}^{-1}_{(\Delta x)}
    \)
  4. \(\Delta x\)が0のとき、恒等操作なはずである。
    \(
    \begin{align}
    \hat{T}_{(0)}\psi_{(x)}&=\psi_{(x)} \\
    &\rightarrow \hat{T}_{(0)}=\mathbf{1}
    \end{align}
    \)

平行移動演算子の具体的な形


1,より、\(\hat{T}_{(\Delta x)}\)はユニタリー演算子である事が分かりました。ユニタリー演算子を作る1つの方法として、エルミート演算子\(\hat{A}_{(\Delta x)}\)を用いて、
\(
\hat{T}_{(\Delta x)}=e^{i\hat{A}_{(\Delta x)}}
\)

と置けば、\(\hat{T}_{(\Delta x)}\)はユニタリー演算子となります(エルミート演算子とユニタリー演算子の性質)。
この仮定の下、
\(
\hat{T}_{(\Delta x)}\psi_{(x)}=\psi_{(x+\Delta x)}
\)

に代入して、
\(
e^{i\hat{A}_{(\Delta x)}}\psi_{(x)}=\psi_{(x+\Delta x)}
\)

無限小移動を考えて\(\Delta_{(x)}\rightarrow 0\)と約束すれば、
\(
(1+i\hat{A}_{(\Delta x)})\psi_{(x)}=\psi_{(x+\Delta x)}
\)

変形して、
\(
\psi_{(x+\Delta x)}-\psi_{(x)}=i\hat{A}_{(\Delta x)}\psi_{(x)}
\)

辺々を\(\Delta x\)で割れば、
\(
\begin{align}
\displaystyle \lim_{\Delta x\rightarrow 0}\frac{\psi_{(x+\Delta x)}-\psi_{(x)}}{\Delta x} = \lim_{\Delta x\rightarrow 0} i\frac{\hat{A}_{(\Delta x)}}{\Delta x}\psi_{(x)}
\end{align}…(i)
\)

ここで右辺を見ます。\(\hat{A}_{(\Delta x)}\)は\(\Delta x\)のなんらかの関数であるはずですが、\(\hat{A}_{(\Delta x)}\)が分母の\(\Delta x\)と打ち消す形でなければ0、もしくは発散してしまう事になります。
具体的に、もしも\(\hat{A}_{(\Delta x)}=(\Delta x)^2\)だったら、
\(
\displaystyle \lim_{\Delta x\rightarrow 0} \frac{(\Delta x)^2}{\Delta x}=\lim_{\Delta x\rightarrow 0} \Delta x =0
\)

であるし、もしも\(\hat{A}_{(\Delta x)}=(\Delta x)^{1/2}\)だったら
\(
\displaystyle \lim_{\Delta x\rightarrow 0} \frac{(\Delta x)^{1/2}}{\Delta x}=\lim_{\Delta x\rightarrow 0} \frac{1}{\sqrt{\Delta x}} =\infty
\)

となります。
すなわち、言いたいことは、\(\hat{A}_{(\Delta x)}\)は\(\Delta x\)に関して1次であるべきで、
\(
\hat{A}_{(\Delta x)}=\hat{B}\cdot \Delta x
\)…(ii)
と書けるはずです。ここで\(\hat{B}\)は\(\Delta x\)に依らない演算子です(下に補足説明あり)。

よって、式(i)は\(\Delta x\rightarrow 0\)のとき、
\(
\begin{align}
\displaystyle \lim_{\Delta x\rightarrow 0}\frac{\psi_{(x+\Delta x)}-\psi_{(x)}}{\Delta x} &= \lim_{\Delta x\rightarrow 0} i\frac{\hat{A}_{(\Delta x)}}{\Delta x}\psi_{(x)} \\
\frac{d \psi_{(x)}}{dx}=i\hat{B}\psi_{(x)}
\end{align}…(iii)
\)

と書けるため、
\(
\psi_{(x)}=e^{i\hat{B}\cdot x}\psi_{(0)}
\)

となります。(この式で右辺の\(\psi_{(0)}\)が0である必要はなく、定数であればいいんです。)

平行移動演算子と運動量演算子との関係


次元について考えましょう!
今、\(\hat{T}_{(\Delta x)}=e^{i\hat{B}\cdot \Delta x}\)であり、\(\Delta x\)の次元は\(\mathrm{[m]}\)であるため、\(\hat{B}\)の次元は\(\mathrm{[m^{-1}]}\)でなければなりません。よって、\(\hat{B}\)は波数の次元を持つ演算子であると考えられるため、\(\hat{B}\)を\(\hat{k}\)と記述することにします。
式(iii)より、
\(
\displaystyle \frac{d \psi_{(x)}}{dx}=i\hat{k}\psi_{(x)}
\)

辺々に\(-i\hbar\)を掛けると、
\(
\displaystyle -i\hbar \frac{d}{dx} \psi_{(x)}=\hbar \hat{k}\psi_{(x)}
\)
…(iv)

!!!

\(p=\hbar k\)は運動量を表すものでした。よって,\(\hat{p}=\hbar \hat{k}\)と書けば、式(iv)より
\(
\begin{align}
\hat{p}\psi_{(x)}&=-i\hbar\frac{d}{d x} \psi_{(x)} \\
&\rightarrow \hat{p}=-i\hbar\frac{d}{d x}
\end{align}
\)
となり、運動量演算子が導けました。
各成分について言えるので、多次元の場合では
\(
\hat{p}_x=-i\hbar\frac{\partial}{\partial x}
\)
としても同じことです。
申し訳ないですが、数学的に厳密か?は保証しません。

補足説明


果たして運動量演算子として考えた\(e^{i\hat{B}\cdot \Delta x}\)は\(\hat{T}\)の特性1.~4.を満たすでしょうか?

  1. \(
    \left(e^{i\hat{B}\cdot \Delta x}\right)^{\dagger} \left(e^{i\hat{B}\cdot \Delta x}\right)=e^{-i\hat{B}^{\dagger}\cdot \Delta x}e^{i\hat{B}\cdot \Delta x}=\mathbf{1}
    \)
    ここで\(\hat{B}\)はエルミート演算子を考えているので、\(\hat{B}^{\dagger}=\hat{B}\)です。
    \(
    e^{-i\hat{B}^{\dagger}\cdot \Delta x}e^{i\hat{B}\cdot \Delta x}=\mathbf{1}
    \)
    ok!
  2. \(
    \hat{T}_{(\Delta x^{\prime})}\hat{T}_{(\Delta x)}=\hat{T}_{(\Delta x+\Delta x^{\prime})}
    \)
    か?
    \(
    (\mbox{左辺})=e^{i\hat{B}\cdot \Delta x^{\prime}} e^{i\hat{B}\cdot \Delta x}=e^{i\hat{B}\cdot (\Delta x^{\prime}+\Delta x)}
    \)
    ok!
  3. \(
    \hat{T}_{(-\Delta x)}=\hat{T}^{-1}_{(\Delta x)}
    \)
    か?
    両辺の左から\(\hat{T}_(\Delta x)\)を掛けて、
    \(
    \begin{align}
    &\hat{T}_{(\Delta x)}\hat{T}_{(-\Delta x)}=\mathbf{1} \\
    &= e^{i\hat{B}\cdot \Delta x}e^{i\hat{B}\cdot (-\Delta x)}=e^{i\hat{B}\cdot 0}=\mathbf{1}
    \end{align}
    \)

    ok!
  4. \(e^{i\hat{B}\cdot 0}=1\)
    ok!

よって、4つすべてを満たすので、\(\hat{T}_{(\Delta x)}=e^{i\hat{B}\cdot \Delta x}\)は平行移動演算子として適当であると考えることができるのです。

[PSO2]シェーダー品質の簡易と標準の違い

最近PSO2のスクリーンショットを撮ることにはまりましてね。
あまりこつとか知らないもので、『数打ちゃ当たる』でとっているわけです。

ただ、綺麗にとるためには何をすればいいか?ということで検索してみたらあったんです。
【PSO2】秘法!豚が教える『萌カワ目線でスクリーンショット』テクニック!
記事自体は古いので、カメラの視点を指定できることとかは書いていませんが、視点変更や、設定に関して参考になりました。

その中の重要なポイントに、
シェーダー品質を標準にすること!!
という文があったわけなんですね。
当然自分のは標準になっているだろうと思ったらなってなかったのです。

でもそんな効果ないやろ~って思っていたわけですが、明らかに以前と違いが判るんです。
SSを撮る場合だけではなく、プレイしているだけでも完全に分かるんです。

シェーダー品質による実際の違い

一例として紹介しますね。
TMGのPAデッドアプローチの効果の違いです。

これがシェーダー品質を簡易にした時で
pso20150411_005211_005_c

これがシェーダー品質を標準にした時です。
pso20150411_161048_000_c

全然違うじゃない!綺麗すぎるよ!感動!!!

この分動作が重くなったかと言われればそれは感じなかったですね。
環境は、
CPU: Intel® corei7 2600K

ビデオカード: AMD Radeon HD6800series
です。


スポンサーリンク


レジスト系によるダメージ減少量

レジスト系のダメージ減少量はどのくらい?
乙りたくない一心で作り上げた僕の装備は本当に効果があるのか?

装備はこの2つで試しました。レジスト系がある場合と無い場合です。

レジスト系なし

レジスト無_c2

レジスト系あり

レジスト系あり_c2
打撃系のレジストは合計で19%です。

これでSHでの、
ウーダン攻撃2
ウーダンのこの攻撃で比較、試しました。

結果はこれです。
レジストによるダメージ量減少_c

おおよそ150/184≒0.815なので、おおよそ18-19%のダメージカットが成功しています。数値通りですね。
中央値は大きく変わっていますが、ダメージのばらつき具合はあんまり変わらないんですね。

試行回数は少ないですが、まぁ数学的な精密性を求めてるわけじゃないので・・・。

一つ言えることは、死ぬことは少なくなりました。HPも上げたしね。
オールレジストⅢと付けたい人、レジスト系は本当に効果があるか試したい方、参考にしてみてください。

じゃあ実際Guでダメージどのくらい与えられるの?

ガンナーwwオワコンwwwww

うっさいわ!ダメージこんくらいあたえられるんや!これ見い!!

PAサテライトエイム使用時, GuRa, ブロッカブレッド使用です。
対クラパーダ
pso20150409_220845_012_a_c

対ゴルドラーダ
pso20150409_220747_010_a_c

対ウォルガーダ
pso20150409_220802_011_a_c

ちょうどいざないが始まったので撮りました。
弱点当たった時のおおよその与ダメージです。弱点当てられる確率はだいたい5分5分といったところでしょうかね…
(まぁ弱点以外当たったら目も当てられないですが…アンジャドゥリリは嫌いです。↓こいつ、(pso2 wikiより)
アンジャドゥリリ_0_c
サテライトエイムが車輪に吸収されちゃうので…)

あ、あとずっと滞空できるのでやってて楽しいし!

数値積分(等間隔)

数値計算での積分方法、特に等間隔の分点の場合であるニュートン・コーツ積分(とロンバーグ積分)に関する理論とのプログラムを載せます。

  1. ニュートンコーツ型の公式
  2. ルンゲ現象
  3. ロンバーグ積分
  4. サブルーチン”integral”のコードと使い方
  5. サブルーチン”romberg”のコードと使い方

等間隔の分点で積分を行う場合、よく使われる方法は台形則、もしくはシンプソン積分です。

台形則

台形則は、分点間を台形近似して面積を求める方法であり、下図のようなイメージで積分を近似する方法です。
外形_台形_c
積分\(\displaystyle \int_a^b f(x) dx\)を、
\(
\displaystyle \int_a^b f(x) dx \sim \sum_{i=0}^{N-1} \frac{f(x_{i+1})+f(x_i)}{2}h
\)

として近似します。この時、誤差のオーダーは\(O(h^3)\)となります。
言い換えると、
台形積分とは、隣り合う分点間を1次関数で近似して求める積分
言えます。

シンプソン積分

シンプソン積分は、隣り合う分点間を2次関数で近似して求める積分です。イメージは
外形_シンプソン_c
な感じです。数式では
\(
\displaystyle \int_a^b f(x) dx \sim \sum_{i=0}^{N/2} \frac{f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2})}{3}h
\)

となります。

高次≠高精度

これ以降もあります。分点間を3次、4次、…としていけば高次の積分をすることが可能になります。
この分点間の次数をどんどん上げていって積分をする方法をまとめてニュートンコーツの公式、と呼びます。
例えばの台形則はニュートンコーツの公式の1次に相当し、シンプソン則はニュートンコーツの公式の2次に相当しています。

じゃあ100次の公式作れば精度が凄い公式ができるんじゃないか?となりますが、これは違います。
高次であることと精度が高いことは違うのです。だから使うな、というわけではなくて知っておいてほしいことなんです。
ではなぜか?これはwikipediaのニュートンコーツの公式のページからとってきたものですが、この式に理由があります。
ニュートンコーツ型式_c

問題となるのが誤差項です。例えば台形近似を見ますと、誤差項は\(h^3 f^{(2)}(\xi)\)とあります。hは刻み幅なので、この項は刻み幅を小さくすれば小さくするほど精度がよくなることを言っています。
が、問題は\(f^{(2)}(\xi)\)の方です。この項は本当の関数が分からないと評価できません。
多項式だったら微分するごとに次数が1つづつ減少するのでこの項は高次になれば小さくなります。じゃあもしも微分するごとに値が増加していったら…?この場合、hを小さくしてもそれを上回るくらい\(f^{(2)}(\xi)\)が大きくなることがあります。これをルンゲ現象と呼びます。

ルンゲ現象

例えば関数
\(
\displaystyle f(x)=\frac{1}{1+25x^2}
\)

を考えます。wikipediaのページを参考にすれば、その導関数と値は、
ルンゲ現象_c
となります。高次になればなるほど値が増加する関数なのです。
実際には\(h^n f^{(m)}(\xi)\)の積の値が小さいか大きいか?なので、一概には高次は悪い!とは言えません。高次は危ないくらいでしょう。

これを可能な限り回避するために、低次の積分を組み合わせて使う方法がとられます。
こうすることで高次の値を増加を抑えるのです。

そのほかの手法もいろいろあります。分点の位置に縛られなければガウス求積法が最高でしょう。素晴らしい方法です。ぜひ調べてみてください。


スポンサーリンク


ロンバーグ積分

さて、高次を追い求めるのもいいですが、一線を駕す…とまではいきませんが、ロンバーグ(Romberg)積分というものがあります。
このロンバーグ積分とは、台形則による計算を基本とし、無限小区間での積分結果を補外によって求めよう、というものです。
ここでいう補外の概念は以下のような感じです。

刻み幅1で台形則の結果から積分値S1が分かる

刻み幅0.5で台形則による結果から積分値S2がわかる

刻み幅0.25で台形則による結果から積分値S3がわかる

刻み幅0.175で台形則による結果から積分値S4がわかる

とある程度計算していきます。そうすると刻み幅hを変化させるとそれに応じて積分値がS1→S2→S3→S4…と推移していくわけです。
この推移の仕方を計算すれば、刻み幅無限小の場合の積分値Sが分かる、こういう仕組みになっています。収束加速法と言ったり、リチャードソン補外、なんて言ったりします。
この補外による方法はかなり優れています。先ほどの高次がどうこう、といった問題が発生しないのです。台形則しか使ってないんですから。有限桁計算に適した方法である、とかwikipediaに載ってたりします。

ロンバーグ積分の詳しい計算方法や理論は
ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ


を見ると詳しく書かれています。

日本語では第6章 数値積分と数値微分にロンバーグ積分の記述があります。
英語ではRomberg integralで調べると良いでしょう。参考までに、英語では
Romberg integration exampleや、Romberg Integrationなにかがいいと思います。

2016/03/01
上のRomberg Integrationを参考にして書いたプログラムはこちらです。これは、解析関数fをRomberg積分するプログラムです。
解析関数ではなく、サイズ\(2^{n}+1\)の配列に格納されている場合はこのページの下の方にあるプログラムを参照してください。

call romberg(a,b,s,pre)

で呼び出し、区間a~bの積分\(s=\int_a^b f(x) dx\)を精度preで求めるプログラムです。

サブルーチン”integral”のコードと使い方

下のサブルーチン”integral”は1,2,3次元の短冊近似、台形則、シンプソン積分、シンプソン3/8積分、ブール則(ニュートン・コーツ型積分)をカバーします。
ただし、分点の個数に注意してください。分点の個数を気にせず使える近似方法は、短冊近似、台形近似のみです。

ロンバーグ積分では引数の関係上、ニュートン・コーツ型積分との間をうまく処理できなかったため、別のサブルーチンとして書きます。もしもロンバーグ積分を使いたい場合はもう少し先に進んでください。

まずはニュートン・コーツ型積分です。
下のサブルーチンを使ってください。
使い方は、

 call integral(配列の大きさ,配列,刻み幅,積分結果,積分手法)

という感じです。
実際のプログラムでは、以下のように使用してください。配列yは複素数解列(ただし実軸上の積分)でもokです。

1次元

サブルーチンintegralの使い方
積分則 呼び方
短冊近似(長方形近似)
 call integral(size(w,1),w,h,s,"box")
配列yの大きさに指定はない。
台形近似
 call integral(size(y,1),y,h,s,"trapezoid")
配列yの大きさに指定はない。
シンプソン則
 call integral(size(y,1),y,h,s,"simpson")
配列yの大きさが3,5,7,…個、\(2n+1\)個でないといけない
シンプソン3/8則
 call integral(size(y,1),y,h,s,"simpson38")
配列yの大きさが4,7,10,…個、\(3n+1\)個でないといけない
ブール則
 call integral(size(y,1),y,h,s,"boole")
配列yの大きさが5,9,13,…個、\(4n+1\)個でないといけない

2次元配列の場合は以下のように指定してください。

 call integral(size(w,1),size(w,2),w,hx,hy,s,"simpson")

3次元の場合はこう指定してください。

 call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"simpson")

積分のためのモジュール

module integral_mod
  !developer --> sikino
  !date --> 2015/04/07
  implicit none
  interface integral
     module procedure &
          dintegral, &
          dintegral2d, &
          dintegral3d, &
          cintegral, &
          cintegral2d, &
          cintegral3d
  end interface integral
contains
  subroutine dintegral(N,y,h,s,method)
    integer,intent(in)::N
    double precision,intent(in)::h,y(1:N)
    character(*),intent(in)::method
    double precision,intent(out)::s
    integer::i
    double precision::y0,y1,y2,y3
       
    s=0.d0; y0=0.d0; y1=0.d0; y2=0.d0; y3=0.d0
   
    if(trim(method).eq."box")then
       s=h*sum(y(1:N-1))
    elseif(trim(method).eq."trapezoid")then
       y1=y(1)+y(N)
       do i=2,N-1
          y2=y2+y(i)
       enddo
       s=(y1+2.d0*y2)*h*0.5d0
    elseif(trim(method).eq."simpson")then
       if(mod(N,2).ne.1)then
          write(6,*)"=====cannot calculation with simpson"
          write(6,*)"=====program stop"
          stop
       endif

       y1=y(1)+y(N)
       do i=2,N-1,2
          y2=y2+y(i)
       enddo
       do i=3,N-2,2
          y3=y3+y(i)
       enddo

       s=(y1+4.d0*y2+2.d0*y3)*h/3.d0

    elseif(trim(method).eq."simpson38")then

       if(mod(N,3).ne.1)then
          write(6,*)"=====cannot calculation with simpson38"
          write(6,*)"=====program stop"
          stop
       endif

       y0=y(1)+y(N)
       do i=2,N-2,3
          y1=y1+y(i)
       enddo
       do i=3,N-1,3
          y2=y2+y(i)
       enddo
       do i=4,N-3,3
          y3=y3+y(i)
       enddo
       s=(y0+3.d0*(y1+y2)+2.d0*y3)*3.d0*h/8.d0        

    elseif(trim(method).eq."boole")then
       if(mod(N,4).ne.1)then
          write(6,*)"=====cannot calculation with boole"
          write(6,*)"=====program stop"
          stop
       endif

       y0=y(1)+y(N)
       do i=5,N-4,4
          y1=y1+y(i)
       enddo
       do i=2,N-1,2
          y2=y2+y(i)
       enddo
       do i=3,N-2,4
          y3=y3+y(i)
       enddo

       s=(14.d0*y0+28.d0*y1+64.d0*y2+24.d0*y3)*h/45.d0
    else
       write(6,*)"=====cannot calculation in integral"
       write(6,*)"=====program stop"
       stop
    end if

    return
  end subroutine dintegral
 
  subroutine dintegral2d(Nx,Ny,z,hx,hy,s,method)
    implicit none
    integer,intent(in)::Nx,Ny
    double precision,intent(in)::hx,hy,z(1:Nx,1:Ny)
    character(*),intent(in)::method
    double precision,intent(out)::s
    integer::i
    double precision::ty(1:Ny),r(1:Nx)

    s=0.d0
    ty(1:Ny)=0.d0
    r(1:Nx)=0.d0
    do i=1,Nx
       ty(1:Ny)=z(i,1:Ny)
       call integral(Ny,ty,hy,s,method)
       r(i)=s
    enddo
    call integral(Nx,r,hx,s,method)
       
    return
  end subroutine dintegral2d
 
  subroutine dintegral3d(Nx,Ny,Nz,w,hx,hy,hz,s,method)
    implicit none
    integer,intent(in)::Nx,Ny,Nz
    double precision,intent(in)::hx,hy,hz,w(1:Nx,1:Ny,1:Nz)
    character(*),intent(in)::method
    double precision,intent(out)::s
    integer::i
    double precision::tyz(1:Ny,1:Nz),r(1:Nx)
       
    s=0.d0
    tyz(1:Ny,1:Nz)=0.d0
    r(1:Nx)=0.d0
    do i=1,Nx
       tyz(1:Ny,1:Nz)=w(i,1:Ny,1:Nz)
       call integral(Ny,Nz,tyz,hy,hz,s,method)
       r(i)=s
    enddo
    call integral(Nx,r,hx,s,method)    
       
    return
  end subroutine dintegral3d
 
  subroutine cintegral(N,y,h,s,method)
    integer,intent(in)::N
    complex(kind(0d0)),intent(in)::y(1:N)
    double precision,intent(in)::h
    character(*),intent(in)::method
    complex(kind(0d0)),intent(out)::s

    double precision::res,ims

    s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

    call integral(N,dble(y),h,res,trim(method))
    call integral(N,dimag(y),h,ims,trim(method))
   
    s=dcmplx(res,ims)
   
    return
  end subroutine cintegral

  subroutine cintegral2d(Nx,Ny,z,hx,hy,s,method)
    integer,intent(in)::Nx,Ny
    complex(kind(0d0)),intent(in)::z(1:Nx,1:Ny)
    double precision,intent(in)::hx,hy
    character(*),intent(in)::method
    complex(kind(0d0)),intent(out)::s

    double precision::res,ims

    s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

    call integral(Nx,Ny,dble(z),hx,hy,res,trim(method))
    call integral(Nx,Ny,dimag(z),hx,hy,ims,trim(method))
   
    s=dcmplx(res,ims)
   
    return
  end subroutine cintegral2d

  subroutine cintegral3d(Nx,Ny,Nz,w,hx,hy,hz,s,method)
    integer,intent(in)::Nx,Ny,Nz
    complex(kind(0d0)),intent(in)::w(1:Nx,1:Ny,1:Nz)
    double precision,intent(in)::hx,hy,hz
    character(*),intent(in)::method
    complex(kind(0d0)),intent(out)::s

    double precision::res,ims

    s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

    call integral(Nx,Ny,Nz,dble(w),hx,hy,hz,res,trim(method))
    call integral(Nx,Ny,Nz,dimag(w),hx,hy,hz,ims,trim(method))
   
    s=dcmplx(res,ims)
   
    return
  end subroutine cintegral3d
end module integral_mod

サブルーチン”integral”を用いた例題

実際の使い方。参考にどうぞ。
下のコードは2次元の積分
\(
\displaystyle \int\int xe^{-x^2}e^{-y}dxdy=\frac{e^{-x^2}}{2}e^{-y}+\mbox{const}
\)

を数値積分します。積分範囲を\(x=0\sim 3, \; y=0 \sim 5\)にした場合、その解析解は
\(
\begin{align}
\int_0^3 dx\int_0^5 dy xe^{-x^2}e^{-y}&=\frac{1}{2}(1-e^{-5}-e^{-9}+e^{-14})\\
&=0.4965697373627734784608751894356320535936864348993604\cdots
\end{align}
\)

となります。実行すると

./a.out
simpson   :    0.496569737366759E+000
romberg,3 :    0.496569737362773E+000

となります。コードは、

program main
  use integral_mod
  implicit none
  integer::i,j,n,Nx,Ny
  double precision::hx,hy,ans
  double precision::xmin,xmax,ymin,ymax
  double precision,allocatable::x(:),y(:)
  double precision::s
  double precision,allocatable::w(:,:)
  double precision::f
  external::f
 
  n=10

  xmin=0d0
  xmax=3d0

  ymin=0d0
  ymax=5d0

  Nx=2**(n)
  Ny=2**(n)
  allocate(x(0:Nx),y(0:Ny),w(0:Nx,0:Ny))
  x=0d0; y=0d0; w=dcmplx(0d0,0d0)

  hx=(xmax-xmin)/dble(Nx)
  hy=(ymax-ymin)/dble(Ny)
 
  do i=0,Nx
     x(i)=xmin+hx*dble(i)
  enddo
  do i=0,Ny
     y(i)=ymin+hy*dble(i)
  enddo

  do i=0,Nx
     do j=0,Ny
        w(i,j)=f(x(i),y(j))
     end do
  end do

  call integral(size(w,1),size(w,2),w,hx,hy,s,"simpson")
  write(6,'(A,e25.15e3)')"simpson   : ",s
 
  call integral(size(w,1),size(w,2),w,hx,hy,s,"romberg",3)
  write(6,'(A,e25.15e3)')"romberg,3 : ",s
 
  return
end program main

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

3次元の計算の例題です。
下のコードは
\(
\displaystyle \int\int\int x^4y^4z^4 dxdydz=\left(\frac{1}{5}\{0.7^5-(-1)^5\}\right)^3
\)

を数値積分します。積分範囲を\(x,y,z=-1\sim 0.7\)にした場合、その値は
\(
\displaystyle \int_{-1}^{0.7}\int_{-1}^{0.7}\int_{-1}^{0.7} x^4y^4z^4 dxdydz = 0.127496010896795\cdots
\)

となります。実行すると、

$ gfortran integralmod.f90 main.f90
$ ./a.out
analysis  :   0.127496010896795E-01
simpson   :   0.10507613E-006  0.10507613E-006
romberg,3 :  -0.17347235E-017 -0.17347235E-017

と得られます。中身はこうです。分点の数はx,y,z軸各々で違っていて構いません。

program main
  use integral_mod
  implicit none
  integer::i,j,k,n,Nx,Ny,Nz
  double precision::hx,hy,hz,ans
  double precision,allocatable::x(:),y(:),z(:)
  complex(kind(0d0))::s
  complex(kind(0d0)),allocatable::w(:,:,:)
 
  n=6
  Nx=2**(n)
  Ny=2**(n-1)
  Nz=2**(n+1)
 
  allocate(x(0:Nx),y(0:Ny),z(0:Nz),w(0:Nx,0:Ny,0:Nz))
  x=0d0; y=0d0; z=0d0; w=dcmplx(0d0,0d0)
 
  hx=1.7d0/dble(size(w,1)-1)
  hy=1.7d0/dble(size(w,2)-1)
  hz=1.7d0/dble(size(w,3)-1)
 
  do i=0,Nx
     x(i)=hx*dble(i)-1d0
  enddo
  do i=0,Ny
     y(i)=hy*dble(i)-1d0
  enddo
  do i=0,Nz
     z(i)=hz*dble(i)-1d0
  enddo

  do i=0,Nx
     do j=0,Ny
        do k=0,Nz
           w(i,j,k)=dcmplx((x(i)*y(j)*z(k))**4,(x(i)*y(j)*z(k))**4)
        end do
     end do
  end do

  ans=(0.7d0**5.d0+1.d0)/5.d0
  ans=ans**3.d0
  write(6,'(A,e23.15e2)')"analysis  : ",ans
 
  call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"simpson")
  write(6,'(A,2e17.8e3)')"simpson   : ",dble(s)-ans,dimag(s)-ans
 
  call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"romberg",3)
  write(6,'(A,2e17.8e3)')"romberg,3 : ",dble(s)-ans,dimag(s)-ans
 
  return
end program main

ロンバーグ積分のプログラム


ロンバーグ積分
 call romberg(jx,x,y,s,pre)
配列yの大きさは2,3,5,10,…,\(2^{jx}+1\)個でないといけない ロンバーグ積分では収束精度”pre”を指定すること。もしも収束精度に達しなかった場合、警告と共に、与えられた値での収束限界の積分結果を返す。積分精度を高めたければ刻み幅を小さくする操作(例えば分点数を増やす等)をしてください。”pre”は、積分結果の絶対値が1より大きくなる場合は相対誤差を、小さくなる場合は絶対誤差を取ります。この理由は以前の結果の補正を加え、収束をさせるためであり、非常に小さい誤差の場合いつまでたっても機械誤差のため収束しないからです。
module romberg_mod
  !developer --> sikino
  !date --> 2016/03/01
  !         2016/03/03
  !
  ! 1D case :
  !  romberg(jx,x,y,s,pre)
  !          |  | | | +- (in)  precision (e.g. 1d-8)
  !          |  | | +--- (out) integration result
  !          |  | +----- (in)  y(1:2**jx+1) f(x)
  !          |  +------- (in)  x(1:2**jx+1) x
  !          +---------- (in)  related to array size
  !
  ! 2D case :
  !  romberg(jx,jy,x,y,z,s,pre)
  !          |  |  | | | | +-- (in)  precision (e.g. 1d-8)
  !          |  |  | | | +---- (out) integration result
  !          |  |  | | +------ (in)  z(1:2**jx+1,1:2**jy+1) f(x,y)
  !          |  |  | +-------- (in)  y(1:2**jy+1) y
  !          |  |  +---------- (in)  x(1:2**jx+1) x
  !          |  +------------- (in)  related to array size y
  !          +---------------- (in)  related to array size x
  !
  ! 3D case :
  !  romberg(jx,jy,jz,x,y,z,w,s,pre)
  !          |  |  |  | | | | | +-- (in)  precision (e.g. 1d-8)
  !          |  |  |  | | | | +---- (out) integration result
  !          |  |  |  | | | +------ (in)  w(1:2**jx+1,1:2**jy+1,1:2**jz+1) f(x,y,z)
  !          |  |  |  | | +-------- (in)  z(1:2**jz+1) z
  !          |  |  |  | +---------- (in)  y(1:2**jy+1) y
  !          |  |  |  +------------ (in)  x(1:2**jx+1) x
  !          |  |  +--------------- (in)  related to array size z
  !          |  +------------------ (in)  related to array size y
  !          +--------------------- (in)  related to array size x
  !
  implicit none
  interface romberg
     module procedure &
          dromberg, &
          dromberg2d, &
          dromberg3d
  end interface romberg
contains
  subroutine dromberg(jx,x,y,s,pre)
    implicit none
    integer,intent(in)::jx
    double precision,intent(in)::x(1:2**jx+1),y(1:2**jx+1),pre
    double precision,intent(out)::s

    integer,parameter::jm=6 !--> precision: O(h^(2*jm))
    integer,parameter::nm=2**jm
    double precision,allocatable::tx(:),ty(:)
    double precision::ts
    integer::k

    s=0d0
    if(jx.ge.jm)then
       allocate(tx(1:nm+1),ty(1:nm+1)); tx=0d0; ty=0d0
       do k=0,2**(jx-jm)-1
          ts=0d0
          tx(1:nm+1)=x(k*nm+1:(k+1)*nm+1)
          ty(1:nm+1)=y(k*nm+1:(k+1)*nm+1)

          call romberg_sub(jm,tx,ty,ts,pre)
          s=s+ts
       enddo
       deallocate(tx,ty)    
    else
       call romberg_sub(jx,x,y,s,pre)
    endif

    return
  end subroutine dromberg

  subroutine romberg_sub(jx,x,y,s,pre)
    ! reference "http://www.math.usm.edu/lambers/mat460/fall09/lecture29.pdf"
    implicit none
    integer,intent(in)::jx
    double precision,intent(in)::x(1:2**jx+1),y(1:2**jx+1),pre
    double precision,intent(out)::s

    integer::i,j,k,n,dn
    double precision::h,ps,tmp
    double precision::T(1:jx+1,1:jx+1)

    n=2**jx+1

    h=x(n)-x(1)
    dn=(n-1)/2

    T(1,1)=0.5d0*h*(y(1)+y(n))
    s=T(1,1)
    ps=s
    h=0.5d0*h
    do j=2,jx+1

       ! trapezoidal rule
       tmp=0d0
       do i=1,2**(j-2)
          tmp=tmp+y(1+(2*i-1)*(dn))
       enddo
       T(j,1)=0.5d0*T(j-1,1)+h*tmp

       do k=2,j
          ! Richardson extrapolation
          T(j,k)=T(j,k-1)+(T(j,k-1)-T(j-1,k-1))/(dble(4**(k-1))-1d0)
       enddo
       s=T(j,j)

       ! precision check
       if(abs(s).ge.1d0)then
          if(abs((ps-s)/s).le.pre)exit
       else
          if(abs((ps-s)).le.pre)exit
       endif
       ps=s
       h=0.5d0*h
       dn=dn/2
    enddo

    if(j-1.eq.jx)then
       write(6,'(A)')" -+-+- didn't converge at romberg integral -+-+- "
       write(6,'(A)')"       Please change stepsize h of array x  "
    endif

    return
  end subroutine romberg_sub

  subroutine dromberg2d(jx,jy,x,y,z,s,pre)
    implicit none
    integer,intent(in)::jx,jy
    double precision,intent(in)::pre,x(1:2**jx+1),y(1:2**jy+1)
    double precision,intent(in)::z(1:2**jx+1,1:2**jy+1)
    double precision,intent(out)::s

    integer::i,nx,ny
    double precision::ty(1:2**jy+1),r(1:2**jx+1)
    nx=2**jx+1
    ny=2**jy+1

    s=0.d0
    ty(1:ny)=0.d0
    r(1:nx)=0.d0
    do i=1,nx
       ty(1:ny)=z(i,1:ny)
       call romberg(jy,y,ty,s,pre)
       r(i)=s
    enddo
    call romberg(jx,x,r,s,pre)
       
    return
  end subroutine dromberg2d
 
  subroutine dromberg3d(jx,jy,jz,x,y,z,w,s,pre)
    implicit none
    integer,intent(in)::jx,jy,jz
    double precision,intent(in)::pre,x(1:2**jx+1),y(1:2**jy+1),z(1:2**jz+1)
    double precision,intent(in)::w(1:2**jx+1,1:2**jy+1,1:2**jz+1)
    double precision,intent(out)::s
    integer::i,nx,ny,nz
    double precision::tyz(1:2**jy+1,1:2**jz+1),r(1:2**jx+1)
       
    nx=2**jx+1; ny=2**jy+1; nz=2**jz+1
    s=0.d0
    tyz(1:ny,1:nz)=0.d0
    r(1:nx)=0.d0
    do i=1,nx
       tyz(1:ny,1:nz)=w(i,1:ny,1:nz)
       call romberg(jy,jz,y,z,tyz,s,pre)
       r(i)=s
    enddo
    call romberg(jx,x,r,s,pre)  
       
    return
  end subroutine dromberg3d
end module romberg_mod

ロンバーグ積分(romberg_mod)の例


必要なのは、上で紹介したモジュール “romberg_mod” と下のメインプログラムです。
以下のプログラムは1次元の定積分
\(
\int_1^10 \frac{1}{x^2} dx = 0.9
\)
を分点数\(2^8+1(jx=8)\)個でロンバーグ積分するものです。
精度は\(O(h^{2\cdot min{jm,jx}})\)です。
ここで\(jm\)はモジュールromberg_modの中にパラメータとして宣言されています。

コンパイルは

ifort romberg_mod.f90 main.f90

とでもすればいいでしょう。

program main
  use romberg_mod
  implicit none
  integer::i,jx,n
  double precision::h,a,b,s
  double precision,allocatable::x(:),y(:)
  double precision::f
  external::f
 
  a=1d0
  b=10d0
  jx=8
  n=2**jx+1
  allocate(x(1:n),y(1:n))
  h=(b-a)/dble(n-1)
  do i=1,n
     x(i)=a+h*dble(i-1)
     y(i)=f(x(i))
  enddo
 
  call romberg(jx,x,y,s,1d-8)
  write(6,*)s,"romberg"
 
  stop
end program main

function f(x)
  implicit none
  double precision::f,x
 
  f=1d0/(x*x)
 
  return
end function f

実行例(要求精度は\(10^{-8}\))

>./a.out
  0.900000000062642  
>

2次元ではこう。

program main
  use romberg_mod
  implicit none
  integer::i,j,Nx,Ny,jx,jy
  double precision::hx,hy
  double precision::xmin,xmax,ymin,ymax
  double precision,allocatable::x(:),y(:)
  double precision::s
  double precision,allocatable::w(:,:)
  double precision::f
  external::f
 
  jx=10
  jy=8
  xmin=0d0
  xmax=3d0

  ymin=0d0
  ymax=5d0

  Nx=2**(jx)
  Ny=2**(jy)
  allocate(x(0:Nx),y(0:Ny),w(0:Nx,0:Ny))
  x=0d0; y=0d0; w=dcmplx(0d0,0d0)

  hx=(xmax-xmin)/dble(Nx)
  hy=(ymax-ymin)/dble(Ny)
 
  do i=0,Nx
     x(i)=xmin+hx*dble(i)
  enddo
  do i=0,Ny
     y(i)=ymin+hy*dble(i)
  enddo

  do i=0,Nx
     do j=0,Ny
        w(i,j)=f(x(i),y(j))
     end do
  end do

  call romberg(jx,jy,x,y,w,s,1d-8)
  write(6,*)s
 
  stop
end program main

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

実行例

>./a.out
  0.496569737374163
>
スポンサーリンク

サバゲーフィールドの情報

サバゲーに行った時にメモした内容と知りたいことです。
日付に注意してください。
そして、刻々と変わるはずなので、参考程度にとどめてください。

  • 場所: OPS http://svg-ops.jp/index.html

  • 住所 東京都稲城市坂浜734、(新ゆり天望の丘墓苑の少し先)
    2015/03/31(火) に参加した時の情報

    最寄り駅から到着まで

    新百合ヶ丘駅からはOPSのHPに記載があります。
    京王よみうりランド駅からフィールドまでタクシーで1520円
    ※八幸タクシー使用(0120086850)。電話番号はタクシー乗り場に記載有。

    平日の火曜日、午前8時20分頃に電話しても大体『出払っている』と言われます。(2回中2回、ほかのタクシー会社も同様)
    15~20分程度待っていればタクシーは来ます。
    8時35分頃に駅でタクシーがつかまり、8時50分くらいにはフィールドに到着したと思います。

    到着時

    到着時に参加費を払う。平日だったので場所代は2500円。
    お昼ご飯4種類の中から1つ選ぶ。任意で500円。

    参加人数とゲームの種類

    参加人数 60人位(かなり少なめ。いつもだいたい100~200人はいるようです。)
    定例会ゲーム数 10~14戦 程度
    チーム分けは朝礼時に真ん中ここら辺、として右半分、左半分で分けられました。

    ナイフキルなし、フリーズコールなし、フレンドリーファイアは謝ったうえで当てられた側、当てた側どちらもヒット判定
    ショートフラッグ戦 10分
    フラッグ戦 15分
    ※この日、やったゲームはこの2種類でしたが、10種類以上のゲームパターンが用意されていました。おそらくこの日は人数が少なかったためだと思われます。

    終わりに

    ※ちなみに、帰りに京王よみうりランドまで乗合したいと思い参加者に尋ねましたが、誰もいませんでした。みなさん新百合ヶ丘駅へ行った模様です。
    京王よみうりランド駅近くの飲食店は絶望的に少ない

    費用

    銃を借りないで1人で行く場合
    タクシー 行き1430~1520円、帰り1430~1520+300(送迎)円
    場所代 2500円
    お昼代 500円
    計 6340円程度


アイキャッチ画像を横幅いっぱいに広げる

アイキャッチ画像を画面いっぱいに広げたい!

以前右サイドバーの縮小とメインエリアの拡張によって中央のコンテンツエリアを広げました。
そうしたら弊害が生まれたんです。
エリアを広げることに成功しましたが、アイキャッチ画像がコンテンツエリアの領域いっぱいに広がらなかったんです。

上のアイキャッチ画像中の青線はコンテンツエリアの領域を示しています。

これがモヤモヤしていたので調べて直してみました。
WPのテーマはTwenty fourteenです。

変更するファイルはfunction.phpです。
style.cssではないです。

やること

function.phpファイルの中に、

    // Enable support for Post Thumbnails, and declare two sizes.
    add_theme_support( 'post-thumbnails' );
    set_post_thumbnail_size( 672, 372, true );
    add_image_size( 'twentyfourteen-full-width', 1038, 576, true );

という部分があると思います。ここのset_post_thumbnail_size()を変更し、

    // Enable support for Post Thumbnails, and declare two sizes.
    add_theme_support( 'post-thumbnails' );
    //set_post_thumbnail_size( 672, 372, true );
    set_post_thumbnail_size( 1038, 519, false );
    add_image_size( 'twentyfourteen-full-width', 1038, 576, true );

と記述します。この記述をしたうえで、アイキャッチ画像を再びアップロードすれば横幅いっぱいに広がるようになります。

ここからは補足と関数set_post_thumbnail_sizeの説明です。

補足と関数set_post_thumbnail_sizeの説明


上で何をしたのか、説明をします。
デフォルトの

set_post_thumbnail_size( 672, 372, true );

の意味は、アイキャッチ画像のサイズを
横 672px
縦 372px
の大きさで画像をトリミングする(true)
という意味になっています(テンプレートタグ/the post thumbnail)。

これのせいで強制的に横幅672pxに変更され、幅が開いてしまっていたのです。しかも勝手にトリミングされるので用意した画像通りの出力がされないのです。
これを解決するために変更します。

set_post_thumbnail_size( 1038, 519, false )
としました。これの意味は
アイキャッチ画像のサイズを
横 1038px
縦 519px
の大きさになるまで画像をリサイズする(false)
という意味です。
リサイズといっても縦横比は固定されたままのリサイズになります。
例えば横2000px,縦500pxの画像を用意てアップロードしたら、設定したサイズまで縦、横を満たすよう変更されます。
この場合は縦幅はすでに満たされているため、横幅が1038pxになるまで縮小されるという意味です。

上の条件で横幅いっぱいに表示されるためには、
横:縦 = 2 : 1
よりも横長で、横幅1038px以上の画像ならばOKということです。

参考


複数サイズのサムネイル(アイキャッチ画像)を表示させる方法