Devil Maker: Tokyo_はじめに

たまには, 息抜きにゲームでも.

日本のソシャゲは, どれも最終的にはガチャで良いものを引かないと強くなれない飽和状態によくなる.
また, ガチャを引くアイテムも, ゲームが進行するごとに手に入り難く, 課金を遠回しに要求している.

お金を使わないと強くなれないというのは, ゲームとしてどうなのだろうか.

今回は, 以前からプレイしているゲームを紹介する.

DevilMaker: Tokyoは, 神話などをモチーフにしたカードゲームである.
海外のソーシャルゲームであり. 韓国版と英語版のものがある.
日本語版も以前はあったが, 現在は運営されてない.

おすすめする理由

1, ガチャが重要視されてない.
このゲームは, ガチャで☆3~5のカードがドロップするが, 最高レア度は☆7である.
課金したからといって, 最高レア度を入手することはできず,強化して行けばどれも強いカードになる.
また, ガチャを引くためのアイテムも, 購入せずともゲーム内でどんどん手に入るので困ることは無い.
(リセマラは重要でない)

2, 課金は主にプレイ時間の短縮を目的としている.
課金の主な内容が, ゲーム内アイテムの販売である. これは普通にプレイしていれば入手できるものばかりである.
他にも経験値増加やドロップ率増加アイテムも販売している. これらも, ランダムでダンジョンで同じ効果が出現する.
アイテム収集する時間, ゲームをずっとプレイする時間はないけど, ゲームは楽しみたい人用に, 課金アイテムが存在する.
つまり, このゲームは, 課金はプレイ時間を買うことに相当する.
たまに有料のイベントガチャもあるが, 最強カードがドロップなんてことは, 無い.
(イベントで上位ランクに入れば, ゲーム内課金通貨も報酬で手に入る. )

3, 好きなカードがほぼ確実に手に入る.
たとえランクの低いカードを手に入れても強化・進化をして行けば, どれも最高ランクまで強化可能である.
また, 毎月好きな☆5カードが手に入るシステムもあるので, 物欲センサーがあろうとなかろうと関係ない.

4, 絵が綺麗
正直, 適当な絵だったら, プレイしていない.

ゲームが外国語というのはやや抵抗があるだろうが, ゲーム自体は難しくないので(艦これみたいに半自動進行)すぐになれるだろう.

これから, このゲームについて, 少しずつだが, 説明していく.

二重振り子

座標の取り方は下図のように取ります。棒の伸び縮みは無いものとします。
二重振り子座標
どういう解き方でもいいですが、ここでは

  1. デカルト座標\(L(x,y)\)でラグランジアンを記述
  2. デカルト座標から座標変換し、\((r,\theta)\)でラグランジアンを記述
  3. 新たな座標系で運動方程式を導く

という順で解いていきます。


スポンサーリンク


1, デカルト座標でのラグランジアンLは(運動エネルギーK)-(位置エネルギーU)と書けるため、
\(
L(x_1,\dot{x}_1,y_1,\dot{y}_1,x_2,\dot{x}_2,y_2,\dot{y}_2)\\
\displaystyle =\frac{1}{2}m_1(\dot{x}_1^2+\dot{y}_1^2)+\frac{1}{2}m_1(\dot{x}_1^2+\dot{y}_2^2)-(-m_1gy_1-m_2gy_2)
\)
と書けます。

2, デカルト座標から座標変換
式を簡単にするために座標変換を行います。新しい座標\((r_1,\theta_1,r_2,\theta_2)\)とデカルト座標\((x_1,y_1,x_2,y_2)\)の関係式は
\(
\begin{align}
x_1&=r_1\sin{\theta_1}\\
y_1&=-r_1\cos{\theta_1}\\
x_2&=r_1\sin{\theta_1}+r_2\sin{\theta_2}\\
y_2&=-r_1\cos{\theta_1}-r_2\cos{\theta_2}
\end{align}
\)
という関係があります。各々を時間で微分すれば、
\(
\begin{align}
\dot{x}_1&=\dot{r}_1\sin{\theta_1}+r_1\dot{\theta}_1\cos{\theta_1}\\
\dot{y}_1&=-\dot{r}_1\cos{\theta_1}+r_1\dot{\theta}_1\sin{\theta_1}\\
\dot{x}_2&=\dot{r}_1\sin{\theta_1}+r_1\dot{\theta}_1\cos{\theta_1}+\dot{r}_2\sin{\theta_2}+r_2\dot{\theta}_2\cos{\theta_2}\\
\dot{y}_2&=-\dot{r}_1\cos{\theta_1}+r_1\dot{\theta}_1\sin{\theta_1}-\dot{r}_2\cos{\theta_2}+r_2\dot{\theta}_2\sin{\theta_2}
\end{align}
\)
これらをラグランジアン\(L(x_1,\dot{x}_1,y_1,\dot{y}_1,x_2,\dot{x}_2,y_2,\dot{y}_2)\)に代入します。すると、新たな座標系でのラグランジアン\(L(r_1,\dot{r}_1,\theta_1,\dot{\theta}_1,r_2,\dot{r}_2,\theta_2,\dot{\theta}_2)\)が得られます。
\(
\begin{align}
L(r_1,\dot{r}_1,\theta_1,\dot{\theta}_1,& r_2,\dot{r}_2,\theta_2,\dot{\theta}_2) \\
=&\frac{1}{2}m_1(\dot{r}_1^2+r_1^2\dot{\theta}_1^2)+\frac{1}{2}m_2\left[\dot{r}_1^2+\dot{r}_2^2+r_1^2\dot{\theta}_1^2+r_2^2\dot{\theta}_2^2 \right. \\
&\left.+2(\dot{r}_1 r_2 \dot{\theta}_2-r_1\dot{r}_2 \dot{\theta}_1)\sin{(\theta_1-\theta_2)}+2(\dot{r}_1 \dot{r}_2 +r_1 r_2 \dot{\theta}_1 \dot{\theta}_2)+\cos{(\theta_1-\theta_2)}\right] \\
&+m_1gr_1\cos{\theta_1}+m_2g(r_1\cos{\theta_1}+r_2\cos{\theta_2})
\end{align}
\)

僕は先ほど式を簡単にするために座標変換をする、といいました。しかし、新しい座標系におけるラグランジアンはどう見ても元のデカルト座標系のラグランジアンに比べて複雑です。この理由は物理的な意味から来ています。
振り子をつないでいる棒が伸び縮みしないとすると系の自由度は角度\(\theta_1,\theta_2\)の2つです。
となると運動方程式は最高で2本の独立した方程式になるはずです。
しかし、デカルト座標の場合うまく自由度を落とすことができず、運動方程式は4つになってしまいます。
そこで棒が伸び縮みを簡単に表すことができる座標系に移ることでうまく方程式の数を減らせます。

新しい座標系でのラグランジアンで棒の伸び縮みがないという条件を表すには
\(
\begin{align}
r_1&=l_1\ (l_1\mbox{は定数}) \\
r_2&=l_2\ (l_2\mbox{は定数})
\end{align}
\)
と書けるわけで、また、
\(
\begin{align}
\dot{r}_1&=0 \\
\dot{r}_2&=0
\end{align}
\)
となるわけです。

\(m_1=m_2=m,\ l_1=l_2=l\)という場合を特に考えると、ラグランジアンは
\(
\displaystyle L(\theta_1,\dot{\theta}_1,\theta_2,\dot{\theta}_2)=ml^2\left[\dot{\theta}_1^2+\frac{1}{2}\dot{\theta}_2^2+\dot{\theta}_1\dot{\theta}_2\cos{(\theta_1-\theta_2)}\right]+mgl(2\cos{\theta_1}+\theta_2)
\)
と書けます。あとはラグランジュの運動方程式を当てはめて計算します。

3, 新たな座標系で運動方程式を導く
保存力場中でのラグランジュの運動方程式は
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\theta}_1}\right)-\frac{\partial L}{\partial\theta_1}&=0 \\
\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\theta}_2}\right)-\frac{\partial L}{\partial\theta_2}&=0
\end{aligned}
\right.
\end{eqnarray}
\)
なので、代入し、\(\ddot{\theta}_1,\ddot{\theta}_2\)に関する運動方程式にすれば
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\ddot{\theta}_1&=\frac{1}{2}\left\{-\ddot{\theta}_2\cos{(\theta_1-\theta_2)}-\dot{\theta}_2^2\sin{(\theta_1-\theta_2)}-2\frac{g}{l}\sin{\theta_1} \right\} \\
\ddot{\theta}_2&=-\ddot{\theta}_1\cos{(\theta_1-\theta_2)}+\dot{\theta}_1^2\sin{(\theta_1-\theta_2)}-\frac{g}{l}\sin{\theta_2}
\end{aligned}
\right.
\end{eqnarray}
\)
となります。これは非線形の2階連立微分方程式です。カオスです。解けません。
数値計算で解きます。用いるのは4次ルンゲ・クッタ法です。

実際に解いてみると、こうなります。

2重振り子の実験もされています。

スポンサーリンク

カオスが現れる場合、初期値鋭敏性という性質があります。
初期値がほんのちょっと変わるだけでその後の時間発展の様子が大きく変わる性質です。
ちょっと値を変えると上の動画とはまるで違う動きになります。
この動画は上の動画よりも中心に近いほうの粒子の初期速度が4%違います

機械精度の違いでさえも十分な時間発展の後には大きな違いが現れます。複雑であると同時に面白い現象です。

WordPressでgnuplotを使えるようにする

これは便利!
プラグインを入れることでWordPressでgnuplotが使えます。

使うプラグインは『GNUPlot』というものです。
GNUPlot wordpress plugin v1.1に詳しい説明があります。

プラグイン→新規追加でGNUPlotを導入し、有効化します。
その後サーバーへアクセスし、“wp-content”ディレクトリの中に”cache”という名前のディレクトリを作ります(wp-content/cache)。
これでできるはず。
間違って
wp-content/plugins/gnuplot-wordpress-plugin/

wp-content/plugins/
にディレクトリを作らないように注意しましょう。

使う際は
[gplot]ここにコマンド[/gplot]
として書いてあげればokです。
どうやらデータをプロットすることはできないようです。
Responses to “GNUPlot wordpress plugin v1.1”でのやり取りを見ると、

4. aneubau Says:
February 13th, 2009 at 8:41 am

Is it possible to load data from a file on the server with this plugin ?
5. Mohamed Ibrahim Says:
February 13th, 2009 at 11:05 pm

No its not. I run an edited version of GNUPlot. If you found a way, let me know 🙂

というやり取りがありました。5.で回答しているのがこのプラグインの製作者です。
製作者自体はデータを読み込むことを想定して作っているわけではないが、もしかしたら読み込ませる方法があるのかも。という回答をしています。

また、生成されたgnuplotの画像のデータはディレクトリcache内に残るようなので定期的に消すべきでしょう。
間違って、試しに生成された画像も残るようです。

僕がサーバーの容量を圧迫し始めたら消すので、このページを編集する人は失敗とか気にしないでねミ
また、サイズはset size 0.8,0.5などを用いて小さく変えるほうがいいと思います。
以下は
[gplot]
set size 0.5,0.3
set xr[-3:3]
set yr[-1:1]
plot sin(x)
[/gplot]

とした時の表示です↓

[gplot]
set size 0.5,0.3
set xr[-3:3]
set yr[-1:1]
plot sin(x)
[/gplot]

リーマン面だってこの通り!!
[gplot]
set view 60,240,1,1
set size 0.8,0.5
i={0,1}
set parametric
set isosamples 70
set xlabel “x axis”
set ylabel “y axis”
set zlabel “z axis”
unset sur
set palette defined(0″#00008b”,1″#2ca9e1″,2″#38b48b”,3″#ffff00″,4″#eb6101″,5″#c9171e”)
set ticslevel 0
set pm3d
set pm3d depthorder
set zr[-2.5:2.5]
splot u*cos(v),u*sin(v),real(sqrt(u)*exp(i*0.5*v))
[/gplot]

set view 60,240,1,1
set size 0.8,0.5
i={0,1}             #虚数単位iの定義
set parametric      #媒介変数を使いますという宣言
set isosamples 70   #描写するときに70点使用します
set xlabel "x axis" #x軸の名前
set ylabel "y axis" #y軸の名前
set zlabel "z axis" #z軸の名前
unset sur           #表面の線を見えなくする
set palette defined(0"#00008b",1"#2ca9e1",2"#38b48b",3"#ffff00",4"#eb6101",5"#c9171e")
                    #カラーマップの変更
set ticslevel 0     #z軸の下をxy平面に一致させる
set pm3d            #pm3dを使います
set pm3d depthorder #3dプロット時の陰面処理
set view equal xy   #x軸とy軸のメモリ間隔を一致させる
set zr[-2.5:2.5]
splot u*cos(v),u*sin(v),real(sqrt(u)*exp(i*0.5*v))

gnuplotで画像出力

gnuplotで画像を出力したいことがあります。
その方法を書きたいと思います。

想定する環境

  • gnuplot: Version 4.6 patchlevel 4
  • Imagemagick(convertコマンドを使用)

がインストールされている状況を考えます。特に、gnuplotではdofor構文を用いるため少なくともVersion4.6以上出ないと本稿に沿った画像出力はできませんのでご了承ください。

画像出力について

画像出力するための方法として、epsファイルを出力してからjpgやpng形式に変換します(僕がそれしか知らないためです)。

基本的には
ターミナル(x11とかwxtとかaquaとかpostscriptとかgifとか…)の指定→出力ファイル名の指定→グラフをファイルに出力
という流れになります。
手入力で画像出力する場合はgnuplotで

set terminal postscript eps enhanced color
set output "aaa.eps"
replot

とすれば今出力されているグラフの画像が”aaa.eps”というファイル名で得られます。
もしもjpgやpng形式に変換したければImagemagickのconvertコマンドを用いて、

convert -density 300x300 aaa.eps aaa.jpg

とすればjpg形式が得られます。jpgをpngに書き換えればpngが得られます。

しかし、何枚も画像がほしい場合には不向きです。スクリプトを書きましょう。
“out.plt”というファイルを作り、その中にこう書きます。

fname="$0"
set terminal postscript eps solid enhanced color
set output fname.".eps"
replot
set out
set terminal wxt enhanced

# 1. PNG with background
cv=sprintf("convert -density 300x300 %s.eps %s.png",fname,fname)

# 2. PNG with background
#cv=sprintf("convert -density 150x150 %s.eps %s.png",fname,fname)

# 3. PNG with background white
#cv=sprintf("convert -density 150x150 -background white -flatten -alpha off %s.eps %s.png",fname,fname)

system cv

このファイルをgnuplotのcallコマンドを用いて、

call "out.plt" "bbb"

として呼び出します。callは引数を指定できるコマンドです。
すると”bbb.eps”と”bbb.png”ファイルが作成されます。便利だね!

縮約

縮約はアインシュタインが一般相対性理論を説明するために最初に導入した数式のお約束です。
これが考え出されたのはアインシュタインが

和記号を書くのに飽きた

からであると予想されます。
このお約束なんですが,大学の物理科ではあまり触れられない癖に教授たちは常識のように語るので,
身につけておいて損はないですね。

さて,あるテンソル量\(A_{i}\)と\(B_{i}\)があったとしましょう。
この二つの縮約をとるということは,同じテンソル同士の和を取る,という意味になります。
アインシュタインの記法に則れば次のように書きます。

\(
A_{i}B^{i}=\sum_{i}{A_{i}B_{i}}
\)

ようするに同じ添え字が右下と右上に来たら和をとりましょうという約束です。
上の式は内積を表していることがわかります。
これを使うと外積の\(i\)成分も次のようになります。

\(
(A \times B)_i= \epsilon_{i,j,k}A^{j}B^{k}
\)
\(
\epsilon_{i,j,k}= 1 ~~ {\rm at} ~~ i,j,k=(1,2,3),(2,3,1),(3,1,2)
\)
\(
\epsilon_{i,j,k}= -1 ~~ {\rm at} ~~ i,j,k=(1,3,2),(3,2,1),(2,1,3)
\)
\(
\epsilon_{i,j,k}= 0 ~~ {\rm at} ~~ i,j,k=(otherwise)
\)

この\(\epsilon_{i,j,k}\)をレヴィチビタの完全反対称テンソルといいます。

右サイドバーの縮小とメインエリアの拡張

WordPressのテーマ『twentyfourteen』において、style.cssの11.0 Media Queriesにある↓の部分を変更することで右側のサイドバーの領域を狭めることができます。中央のメインコンテンツエリアを広げることができます。

右側のサイドバーを狭めるのは、下の.content-sidebar{…}だけを変えればいいです。
コンテンツ領域拡張

これを変えるだけではもともとのサイドバー領域の余白部分が増え、サイドバーのコンテンツ領域が狭まるだけです。
中央のメインコンテンツエリアのコンテンツ領域を増やすためには下の.site-content{…}の変更をしたうえで、
サイトをカスタマイズするためのプラグイン『Fourteen Extended』を使用して、
テーマのカスタマイズ→TwentyFourteen Content Options→で474となっている値を874にすればOKです。

@media screen and (min-width: 1008px) {
    .search-box-wrapper {
        padding-left: 182px;
    }

    .main-content {
        float: left;
    }

    .site-content {
        /*margin-right: 29.04761904%; 変更前*/
        margin-right: 15.0%;    /*変更後*/
        margin-left: 182px;
    }

    .site-content .entry-header {
        margin-top: 0;
    }

    .site-content .has-post-thumbnail .entry-header {
        margin-top: 0;
    }

    .content-sidebar {
        /*margin-left: -29.04761904%;   変更前*/
        /*width: 29.04761904%; 変更前*/
        margin-left: -20.00%;   /*変更後*/
        width: 20.00%;  /*変更後*/
    }

モンテカルロ法

モンテカルロ法は,乱数を使って系の平衡状態にせまる一つの方法です。

ここでは二次元正方格子のイジング模型について紹介します。
イジング模型は隣り合う格子のスピン同士の相関エネルギーのみを考慮した模型です。
ハミルトニアンは以下のようになります。

\(
\displaystyle
H=-J \Sigma_{< i , j >}{s_{i}s_{j}}
\)
 
あるスピン状態を持つ系の状態Aについて,そのハミルトニアンを\(H_A\)とすると,
ボルツマン因子と分配関数\(Z\)を用いて状態の実現確率を計算することが出来ます。

\(
P(A)=\frac{\exp{(-\frac{H_A}{k_B T})}}{Z}
\)
 
Aからランダムに一つの格子におけるスピンを一つだけフリップした状態を状態Bとして,
状態Aと状態Bの実現確率比を乱数と比較します。
\(P(B)/P(A)>{\rm random}\)ならば状態Bを採用し,そうでなければ状態Aを採用します。
ここまでの流れを1mcs(モンテカルロステップ)といい,これを何度も繰り返すことで平衡状態に近づきます。
以下の二つのgifは,相転移温度より高い温度と低い温度で,
モンテカルロステップにより二次元正方格子のスピンがどのように振る舞うかをみているアニメーションです。
黒がダウンスピン,白がアップスピンです。


温度T=10eVのスピンの挙動
温度T=10eVのスピンの挙動

こちらは温度T=10eVで,自発磁化を持たない事が分かります。


温度T=0.001eVのスピンの挙動
温度T=0.001eVのスピンの挙動

こちらは温度T=0.001eVで,スピンが揃って自発磁化を持つ事が分かります。

画像をクリックするとmcsスタート!!

スピンのふるまいを追うと,自発磁化を持つ温度と持たない温度でスピンのふるまいが違う事が分かります。
二次元のイジング模型はオンサーガーにより厳密に解かれている(相転移温度T=2.26eV at J=1)ので,
それを踏まえて下のプログラムで遊んでみるといいですよ。

これのプログラムソースコード(保証は一切しません!!)
*注意 このプログラムを実行する前にising_anというディレクトリを作ってください。
*追記 2015 2/13 磁化とエネルギーを計算出来るようにしました。
計算方法は重み付き選択法というのを使ってます。

\(
\langle A \rangle = \frac{1}{T} \sum_{t=t_0+1}^{t_0 + T}A_{\alpha_{t}}
\)
 
十分平衡状態になっている\(t_0\)以降のモンテカルロステップで,物理量の統計平均を取る方法です。
平衡以前の物理量を計算しないmcsをからまわしと言います。

格子点51×51で,\(J=2\)で計算しました。
スピンの大きさを1として,相互作用をダブルカウントしているので注意しましょう。
即ち\(J=1\)の結果が欲しければプログラム上ではintJ=0.125にすればよいかも。

**********parameter**********
i_max::x方向の格子点の数
j_max::y方向の格子点の数
mcs_max::最大モンテカルロステップ数
mcs2step::mcs2step以降のmcsを平衡状態と見なし,物理量を計算します
intj::\(J/4\)
temp::温度
mov_max::mcs_max/movmax毎のスピン状態をgifアニメーションで見れるように区切ります。
**********parameter**********

単位格子あたりのエネルギーの温度依存性グラフを添付します。

ene

相転移温度付近でエネルギーの変化が大きい事が分かります。

単位格子あたりの磁化の絶対値の温度依存性グラフを添付します。

mag

相転移温度付近で自発磁化が発生し始めている事が分かります。

gnuplotを用いて,格子点上のスピンがどのように平衡状態へ近づくかをダイナミックに確認する事が出来ます。

gnuplot anime.gnuplot
 
を実行するとアニメが始まります。

gnuplot animegif.gnuplot
 
でgifアニメが作成出来ます。

!————————————————
! title = ising2d.f90
! developer = Amesyabody
! released = 2015 2/13
! programming language = fortran
!
!explanation
!2 dimentional square lattice Ising model's dynamical magnetism
!calculated by Monte Carlo method.
!————————————————
   
program main
    implicit none
    integer(4)::i,j,i_max,j_max,mcs,mcs_max,mov,mov_max,mcs2step
    parameter(i_max=50,j_max=50,mcs_max=1000000,&
    mcs2step=500000,mov_max=100)
    real(8)::intj,temp
    parameter(intj=0.25d0,temp=0.1d0)
    real(8)::spin(0:i_max,0:j_max)
    real(8)::ratio,ham1,ham2,hamres,dum,ene
    integer(4)::px,mx,py,my,rndx,rndy
    character(10)::movc,movc2

    real(8)::mag,mag2,magres
   
    !-----make random-----
    integer(4),allocatable::seed(:)
    integer cnt,seedsize
    real(8)::rnd,rndi,rndj,rnds
   
    call random_seed(size=seedsize)
    allocate(seed(seedsize))
    write(*,*)seedsize

    dum=rndmf(0)

    call system_clock(count=cnt)
    seed=cnt!+idnint(dum*1000000)
!            dum=rndmf(seed(0))
!            seed=idnint(seed*dum)
    call random_seed(put=seed)

    !------------------

mag2=0.0d0
ene=0.0d0

do mov=0,mov_max

write(movc,'(i3)')mov

open(100+mov,file="ising_an/spin"//trim(adjustl(movc))//".dat"&
,status="unknown")

end do

open(5000,file="ising_an/anime.gnuplot",status="unknown")
open(5001,file="ising_an/animegif.gnuplot",status="unknown")

    !-----primitive spin set-----
   
    do i=0,i_max
        do j=0,j_max
            call random_number(rnd)
!           write(*,*)rnd
            if(rnd>0.5d0)then
                spin(i,j)=1.0d0
            else
                spin(i,j)=-1.0d0
            end if
!           write(*,*)spin(i,j)!,cnt,seed(1)

        write(100,'(i3,i3,f5.1)')i,j,spin(i,j)

        end do

        write(100,*)" "

    end do
   
    !----------------------------
   
    do mcs=1,mcs_max
   
        do i=0,i_max
       
            do j=0,j_max
           
                px=i+1
                if(i==i_max)px=px-i_max-1
                py=j+1
                if(j==j_max)py=py-j_max-1
                mx=i-1
                if(i==0)mx=mx+i_max+1
                my=j-1
                if(j==0)my=my+j_max+1

                ham1=ham1-intj*spin(i,j)*spin(px,j)&
                -intj*spin(i,j)*spin(i,py)&
                -intj*spin(mx,j)*spin(i,j)&
                -intj*spin(i,my)*spin(i,j)

                do mov=1,mov_max-1

                    if(mcs==mov*mcs_max/mov_max)&
                    write(100+mov,'(i3,i3,f5.1)')i,j,spin(i,j)

                end do

                if(mcs>mcs2step)then

                    mag2=mag2+spin(i,j)

                end if
               
            end do

            do mov=1,mov_max-1

                if(mcs==mov*mcs_max/mov_max)write(100+mov,*)" "

            end do

        end do

        if(mcs>mcs2step)then

            magres=magres+abs(mag2)
            mag2=0.0d0
            ene=ene+ham1

        end if

        call random_number(rndi)
        call random_number(rndj)
        rndx=idnint(rndi*i_max)
        rndy=idnint(rndj*j_max)

!       write(*,*)rnd
!write(*,*)spin(rndx,rndy)
        spin(rndx,rndy)=-spin(rndx,rndy)
!write(*,*)spin(rndx,rndy)
        do i=0,i_max
       
            do j=0,j_max
           
                px=i+1
                if(i==i_max)px=px-i_max-1
                py=j+1
                if(j==j_max)py=py-j_max-1

                mx=i-1
                if(i==0)mx=mx+i_max+1
                my=j-1
                if(j==0)my=my+j_max+1

                ham2=ham2-intj*spin(i,j)*spin(px,j)&
                -intj*spin(i,j)*spin(i,py)&
                -intj*spin(mx,j)*spin(i,j)&
                -intj*spin(i,my)*spin(i,j)
   
            end do
           
        end do
       
        ratio=exp((-ham2+ham1)/temp)
!       write(*,*)ham1,ham2,ratio
       
        ham1=0.0d0
        ham2=0.0d0
       
        call random_number(rnds)

        if(ham1 > ham2)then!metroporis
        goto 200
        end if
       
        if(rnds>ratio)then
        spin(rndx,rndy)=-spin(rndx,rndy)
        end if

200     continue

    end do
   
    do i=0,i_max
        do j=0,j_max
       
        mag=mag+spin(i,j)
       
        px=i+1
        if(i==i_max)px=px-i_max-1
        py=j+1
        if(j==j_max)py=py-j_max-1
        mx=i-1
        if(i==0)mx=mx+i_max+1
        my=j-1
        if(j==0)my=my+j_max+1

        hamres=hamres-intj*spin(i,j)*spin(px,j)&
        -intj*spin(i,j)*spin(i,py)&
        -intj*spin(mx,j)*spin(i,j)&
        -intj*spin(i,my)*spin(i,j)

        write(100+mov_max,'(i3,i3,f5.1)')i,j,spin(i,j)
       
        end do

        write(100+mov_max,*)" "

    end do
   
    mag=mag/(dble(i_max*j_max))
   
!   write(*,'(3f10.5,f15.5)')temp,intj,mag,hamres

    ene=ene/(dble((mcs_max-mcs2step)*(i_max+1)*(j_max+1)))
    magres=magres/(dble((mcs_max-mcs2step)*(i_max+1)*(j_max+1)))

    write(*,'(A,f10.5,A,f10.5)')"energy=",ene,"      magnet=",magres

    write(5000,*)"set size sq"
    write(5000,*)"set pm3d map"
    write(5000,*)"set palette rgbformulae 21,22,23"
    write(5000,*)"set xlabel 'x'"
    write(5000,*)"set ylabel 'y'"
    write(5000,*)"set zra[-1:1]"
    write(5000,*)"set zlabel 'spin'"
    write(5000,*)"splot 'spin0.dat'"
    write(5000,*)"pause 0.1"
    do mov=1,mov_max
    write(movc,'(i3)')mov
    write(5000,*)"splot 'spin"//trim(adjustl(movc))//".dat'"
    write(5000,*)"pause 0.1"

    end do

    write(5001,*)"set tics font 'Times,15'"
    write(5001,*)"set xlabel font 'Times,15'"
    write(5001,*)"set ylabel font 'Times,15'"
    write(5001,*)"set size sq"
    write(5001,*)"set pm3d map"
    write(5001,*)"set palette rgbformulae 21,22,23"
    write(5001,*)"set xlabel 'x'"
    write(5001,*)"set ylabel 'y'"
    write(5001,*)"set zra[-1:1]"
    write(5001,*)"set cbrange [-1:1]"
    write(5001,*)"set zlabel 'spin'"
    write(5001,*)"set term gif animate optimize"
    write(5001,*)"set output 'is_t10.gif'"
    write(5001,*)"do for[i=0:100] {"
    write(5001,*)'file = sprintf("spin%01d.dat", i)'
    write(movc2,'(i10)')mcs_max/mov_max
    write(5001,*)'time = sprintf("s=%d[mcs]",i*'//trim(adjustl(movc2))//')'
    write(5001,*)"set title time"
    write(5001,*)"splot file"
    write(5001,*)"}"

    contains

real(8) function rndmf(seeds)
implicit none

integer(4)::a,m,q,p,n,ndiv,j,k,seeds
real(8)::rm,rmax
parameter(a=16807,m=2147483647,rm=1.0/m)
parameter(q=127773,p=2836,n=32,ndiv=1+(m-1)/n)
parameter(rmax=1.0-1.2e-7)
integer(4)r(n),r0,r1

if (seeds .ne. 0)then
r1=abs(seeds)
do j=n+8,1,-1
k=r1/q
r1=a*(r1-k*q)-p*k
if(r1 .lt. 0)r1=r1+m
if(j .le. n)r(j)=r1
end do
r0=r(1)
end if


k=r1/q
r1=a*(r1-k*q) -p*k
if(r1 .lt. 0)r1=r1+m
j=1+r0/ndiv
r0=r(j)
r(j)=r1
rndmf=min(rm*r0,rmax)

end function
   
end program

2分法、false position法、ニュートン法

ゼロ点を探す、または方程式の解、根を求める方法です。
ここではその方法の説明と、fortran90でのプログラムを載せておきます。

  1. 2分法
  2. false position法
  3. (有限差分による)ニュートン法
  4. 二分法とfalse position法、ニュートン法のプログラム

方程式の解を見つけるとは、その関数のゼロ点を見つける事に帰着されます。
具体的に
\(cos(x)=0\)
という方程式を満たすxを探すことを考えます。
これは、\(f(x)=cos(x)\)という関数を考え、\(f(x)=0\)になることを考えればいいのです。
\(f(x)=cos(x)\)のグラフはこのようになります。
cosx_1
もしもxの範囲を0~10に限れば、解は
\((\frac{\pi}{2}, \frac{3\pi}{2}, \frac{5\pi}{2})\)
の3つです。このくらいの関数では簡単にわかります。

しかし、方程式
\(cos(2\sqrt{x}\sin(x))=0\)
を満たすxを見つけなさい、だったら容易には分かりません。グラフはこうなります。
cos2Jxsinx_1

これを求める方法として有名なのが、2分法とニュートン法(微分は有限差分)です。
有名かどうかは分かりませんが、2分法よりも早い方法としてfalse position法というものがあり、これも紹介します。
ここで紹介するのは、上記の3つですが導関数が分からない場合、2分法の確実性と高次の方法の早さを合わせたBrent法(Van Wijngaarden-Dekker-Brent法, William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery著, 丹慶勝市ら訳(1993) 『Numerical Recipes in C [日本語版]C言語による数値計算のレシピ』 技術評論社 261pp.)が推奨されています。

収束の確実性と収束の早さの関係はおおよそこんな感じです。

確実性:
2分法 >> false position法 > 割線法(Secant法) >> ニュートン法
収束の早さ:
ニュートン法 >> 割線法(Secant法) > false position法 >> 2分法

Brent法はミックスした方法なので外しています。

2分法


これは解を含むxの範囲に対し解が存在する範囲を半分ずつ挟み込み、狭めていく方法です。
グラフで概念を表すとこういう手順になります。
二分法

収束速度について

2分法では、1回の計算あたり範囲が\(1/2\)になるため、漸化式
\(\varepsilon_n=\frac{1}{2}\varepsilon_{n-1}\)
が成立するため、
\(\varepsilon_n=\left(\frac{1}{2}\right)^n \varepsilon_{0}\)
が導けます。ここで\(\varepsilon_{0}\)は初期のxの範囲(\(\varepsilon_0=|b-a|\))を表します。
収束精度と計算回数の関係を決めるためにはnについて解いて、
\(
\begin{align}
\varepsilon_n &=\left(\frac{1}{2}\right)^n \varepsilon_{0} \\
2^n &=\frac{\varepsilon_0}{\varepsilon_n} \\
n &=\log_2\left(\frac{\varepsilon_0}{\varepsilon_n}\right) \\
\mbox{または}& \\
n &=\frac{\ln\left({\frac{\varepsilon_0}{\varepsilon_n}}\right)}{\ln{2}}
\end{align}
\)
より、例えば初期の解の範囲(\(\varepsilon_0\))を1, 解の収束を\(10^{-14}\)まで収束させようとすれば、大体\(n=47\)となり、47回、2分法を繰り返す必要があります。

false position法


false position法というものを紹介します。(2016/11/08 これは嘘です。falsepositionに関しては参考しないでください。)
こんな感じで収束させます。
falseposition_compressed

これは、片方の点を固定して常に符号の違う点同士で挟みます)、一次方程式で区間を近似、そのゼロ点を元に解に近づいていく方法です。
2分法より早いです。2分法ほど確実性はありませんが、収束は全然早いです。振る舞いがいい関数であれば、10回もやれば十分すぎます。
また、割線法(Secant法)という手法もあります。これは、false position法ほどの確実性はありませんが、収束はfalse position法よりも若干早い方法です。
(参考:『Numerical Recipes in C [日本語版]C言語による数値計算のレシピ』 技術評論社)

ニュートン法


これは接線を利用して収束させていく方法です。
グラフで概念を表すとこういう手順になります。
ニュートン法

収束精度について

※詳細はあとで書きます。

完全な微分が分かっている場合、ニュートン法は2次の収束を見せます。
これは、ニュートン法を1回行うごとに一致桁数が2倍になる、という意味です。
すなわち、1回目に本当の解と1桁目が一致していたら、2回目には2桁一致、3回目には4桁一致、4回目には8桁一致することになります。
なので、ニュートン法において10回計算して収束しないとかはあり得ないのです。
ただし、微分を有限差分
\(\frac{dy}{dx}\sim \frac{f(b)-f(a)}{b-a}\)
で求めていた場合は収束は遅くなるようです。それでも10回位計算すれば機械精度まで収束するでしょう…たぶん。

プログラム


\(cos(2\sqrt{x}\sin(x))=0\)
の解を求める2分法とfalse position法、それに有限差分によるニュートン法のプログラムです。
特異点があってもそれは根ではない、と判断して排除してくれます。

出力結果はこうなります。
zeropoint

変数の説明です。
ad:ゼロ点を探すx軸の範囲
bd:ゼロ点を探すx軸の範囲
rd:見つかったゼロ点の値を入れる配列
count:見つかったゼロ点の総数
is: 1 → 2分法で見つける
2 → ニュートン法で見つける
Nr: 取得したいゼロ点の数
count=5の時、Nr=3にした場合

  1 2 3 Nr 4 5 count
本来のゼロ点 0.94 2.64 3.57 5.96 6.59
配列rdに入るもの 0.94 2.64 3.57 0 0

count=5の時, Nr=7にした時

  1 2 3 4 5 count 6 7 Nr
本来のゼロ点 0.94 2.64 3.57 5.96 6.59    
配列rdに入るもの 0.94 2.64 3.57 5.96 6.59 0 0
program main
  implicit none
  integer::i,count,Nr
  double precision::ad,bd
  double precision,allocatable::rd(:)
  double precision::f
  external::f

  Nr=20
 
  allocate(rd(1:Nr))
  rd(1:Nr)=0.d0
 
  ad=0.d0
  bd=10.d0
 
  call rootfinding(f,ad,bd,Nr,rd,count,"bisection")
  do i=1,count
     write(6,*)rd(i)
  enddo
 
  call rootfinding(f,ad,bd,Nr,rd,count,"newton_raphson")
  do i=1,count
     write(6,*)rd(i)
  enddo

  call rootfinding(f,ad,bd,Nr,rd,count,"false_position")
  do i=1,count
     write(6,*)rd(i)
  enddo

  stop
end program main

function f(x)
  implicit none
  double precision::f,x
 
  f=cos(2.d0*sqrt(x)*sin(x))
 
  return
end function f

subroutine rootfinding(func,ad,bd,Nr,rd,count,method)
  !Date : 2015/07/16
  !Developer : sikino
  implicit none
  interface
     function func(ix)
       implicit none
       double precision,intent(in)::ix
       double precision::func
     end function func
  end interface
  integer,intent(in)::Nr
  double precision,intent(in)::ad,bd
  integer,intent(out)::count
  double precision,intent(out)::rd(1:Nr)
  character(*),intent(in)::method

  integer::N,i,j
  double precision,parameter::hd=1d-2
  !                           ^resolution of root
  double precision,parameter::eps=1.d-14
  !                           ^convergence epsilon
               
  double precision::h,x0,x1,x2,y0,y1,y2,Dy0,tx1,ty1

  !ad < bd
  if(ad.ge.bd)then
     write(6,'(A,2e15.6e3)')"must be ad < bd, your ad,bd --> ",ad,bd
     write(6,'(A)')"program stop at rootfinding"
     stop
  end if

  !Announsment
  write(6,'(A,A)')"----",trim(method)
 
  N=int(abs(ad-bd)/hd)
  count=0

  x0=ad
  y0=func(x0)
  do i=1,N
     x1=x0+hd
     y1=func(x1)
     if(y0.lt.0.d0.neqv.y1.lt.0.d0)then
        tx1=x1
        ty1=y1
        if(trim(method).eq."bisection")then
           !bisection rule
           do j=1,60
              x2=0.5d0*(x0+x1)
              y2=func(x2)
              if(y2.lt.0.d0.eqv.y0.lt.0.d0)then
                 x0=x2
              else
                 x1=x2
              endif
              if(abs((x1-x0)/x2).lt.eps)then
                 !If y2 is large, y2 is singular point.  
                 if(y2.le.1d0)then
                    count=count+1
                    rd(count)=x2
                 endif
                 exit
              endif
           enddo
           if(j.eq.60)write(6,*)"+--cannot convergence at root-finding method bs--+"
        elseif(trim(method).eq."newton_raphson")then
           !Newton-Raphson method
           do j=1,10
              h=1.d-4
              y0=func(x0)
              Dy0=(func(x0+h)-y0)/h
              x2=x0-y0/Dy0
              if(abs((x0-x2)/x2).lt.eps)then
                 if(y0.le.1d0)then
                    count=count+1
                    rd(count)=x2
                 endif
                 exit
              endif
              x0=x2
           enddo
           if(j.eq.10)write(6,*)"+--cannot convergence at root-finding method nr--+"
        elseif(trim(method).eq."false_position")then
           !false position method
           if(abs(y1).gt.abs(y0))then
              do j=1,30
                 x2=x0-y0*(x1-x0)/(y1-y0)
                 if(abs((x2-x0)/x2).lt.eps)then
                    if(y0.le.1d0)then
                       count=count+1
                       rd(count)=x2
                    end if
                    exit
                 endif
                 x0=x2
                 y0=func(x0)
              enddo
           else
              do j=1,30
                 x2=x0-y0*(x1-x0)/(y1-y0)
                 if(abs((x2-x1)/x2).lt.eps)then
                    if(y1.le.1d0)then
                       count=count+1
                       rd(count)=x2
                    endif
                    exit
                 endif
                 x1=x2
                 y1=func(x1)
              enddo
           endif
           if(j.eq.30)write(6,*)"+--cannot convergence at root-finding method fp--+"
        else
           write(6,'(A,A)')"unknown type of method, your method--> ",trim(method)
           write(6,'(A,A)')"program stop"
           stop
        end if
        if(count.ge.Nr)exit
        x1=tx1
        y1=ty1
     endif
     x0=x1
     y0=y1
  enddo

  return
end subroutine rootfinding