J. Plasma Fusion Res. Vol.96, No.8 (2020)４
０５‐４１７

講座 シミュレーションのための乱数入門

５．乱数の変換と利用例
5. Transformation of Random Numbers and its Applications
佐 竹 真 介，菅 野 龍 太 郎，本 多充１）

SATAKE Shinsuke, KANNO Ryutaro and HONDA Mitsuru１）
自然科学研究機構核融合科学研究所，１）量子科学技術研究開発機構（原稿受付：２
０
２
０年５月２
２日）

この章では，プラズマ・核融合分野，特に，磁場によって閉じ込められた核融合プラズマの理工学におけるシミュレーションにおいて，乱数がどのように利用されているのかを紹介します．シミュレーションにおいては，一様乱数だけでなく，正規乱数など様々な乱数が目的に応じて利用されています．そこで，まずは，乱数を用いたシミュレーションを行う際に必須である，様々な乱数を生成する「乱数の変換」について説明します．次に，シミュレーションにおいて，具体的にどのように乱数が用いられているのかを解説します．乱数を用いるシミュ
レーション手法は，一般に，モンテカルロ法と呼ばれます．ここでは，モンテカルロ法に基づくシミュレーショ
ンについて，基礎的な例を挙げて説明します．また，近年，様々な科学技術分野において注目を集めている遺伝的アルゴリズムと，そのアルゴリズムを用いた大域的最適解を求めるシミュレーションについても説明します．
Keywords:
random number transformation method, Monte-Carlo method, Dirichlet-Poisson solver, stochastic differential equation,
genetic algorithm, global optimizationに持つ単調増加関数になります（ 'が不連続な場合につい５．
１ 乱数の変換ては後述）．よって，,は !の逆関数 !!"を使って

本講座の第２章から４章で紹介してきた疑似乱数や物理乱数は基本的に +!
!
".の実数値一様乱数でした（利用((-+*,#!!")

するサブルーチンによって )
，!
!
!
!
"*
!
"*ないし )!
".の場
合もありますが，ここでは厳密な区別はしません）．しかのように表すことができます．さて，今 #が区間 !
!
".

し，様々なシミュレーションに乱数を利用する際には，一に一様分布する確率変数とすると，#$+となる確率は様乱数以外の様々な確率密度関数に従う乱数が必要になり#$+,,*で す．こ の 時，こ れ に 対 応 し て，
$&
+
#+#!)

ます．本節では，そのような任意の確率密度関数に従う疑,*で$#!!")#*+*と な る 確 率 も 同 じ く !)$,#!!")

似乱数を発生させる方法を紹介します．

す．すなわち，

任意の確率密度関数に従う乱数を生成する方法として

!!")
#*
#"!)
,*
,*
$&
+
$,,
#$&
+
,
#!)

は，大きく分けて２つの代表的な方法が知られています．
まず１つ目は，逆変換法（直接法）と呼ばれる方法で，確が成り立ちます．上式より，区間 !
+!
".の一様乱数 +(,率密度関数 'に対し累積分布関数!"
を用いて ,,+(,は累積分布関数(#! )(*として与えた +,!),*に従う分布をとる，つまり確率密度関数 'に従う分布%になることがわかります．

+#!)
,*
*
%*
#! '
)
*
を考えることで，一様乱数 ++,から !を通じて 'に従う乱(

ところで，!!")+*は常に解析的に表せるわけではなく，数+,(,を生成するという方法です．以下，具体的に見てい不連続な関数である場合も考えられます．そのような場合

きましょう．

は，生成したい乱数+の値域,%!
&.を十分細かく分割し，(,数値積分を使うなど し て +
（)#!
!
!
!
"），
"
#
/!
,
) #!)
)*

ここでは簡単のため，確率密度関数の定義域が %!
&.の

+
，+
，の離散データのテーブルを用意
!#!
" #"

連続な関数であるとします．もし下限，上限の両方あるい

ただしはいずれかが !'，"' となる場合は，*"%，*#&で十分します．そして，一様乱数ルーチンで生成した +
(の値が

に'
*
)
*
&!とみなせるような下限，上限を設定しておきま+) $+)""に入る時，例えば線形補間を用いて，("+しょう．さて，確率密度関数の性質として，'
*
)
*
%!である

,
)""!,
)
,
+
)
)"
)*
(#,
(!+
+
)""!+
)

こと，定義域全体で積分すると "になることの２つがあります．したがって，累積分布関数 !),*は，値域を !
!
".
National Institute for Fusion Science, Toki, GIFU 509-5292, Japancorresponding author’s e-mail:satake@nifs.ac.jp

４０５!2020 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.96, No.8 August 2020!"を生のようにして 4の近似値を得ることで+
4
3
(,
($! )
(*

/

#
0
!/

#
!/

1$!)
/
*
$% '/0/0& # $"!& #成することが可能です．

!

なお，+
4
(,が離散値しか取らない場合も，上と同様の考えか ら 逆 変 換 法 を 用 い て 生 成 で き ま す．確 率 変 数 4がの値域が !
!
".なので先ほどの逆変換法が使えそうです．

+ 通りの離散値 4!
!
+）を取り，それぞれの発
#
1!
*（*$"

しかしこのままでは逆関数が簡単に記述できないので，生確率が -* であるとします（ただし # +
）．この
* $"
*$"-

10$"!1と置き換え，10$& # とすると，その逆関数は/$&!#)
10*で あ る こ と が わ か り ま す．と こ ろ で
*)
/
!
'*
3
!
4*の逆変換は)/)#!/

時，区間 !
!
".一様乱数 3(に対して，
!"
,
#,
* "3
*
(%# *$!*$!-

3$/&
+,'，

を満たす,を求め，この,に対して 4,を採用すること
($4

4$/,
(
*'

で与えられるので，一様乱数 +1(!
2(,を用いて，で，+4(,を生成します．ただし，上式の総和を取る時に便

，
3
1(*&
&2(*
*)
+,)
#
!#)
($&
，
4
1(*,
&2(*
*)
(
*)
#
!#)
($&

宜的に -!$!とおきます．
もう１つの代表的な方法は棄却法と呼ばれるもので，一様乱数を使った試行錯誤によって得たい分布に従う乱数をと与えると，+
は２次元空間の正規分布 '
3
3
!
4*に従う
4
)
(!
(,

発生させるものです．まず，先ほどと同様に，生成したい分布を持つ乱数となります．

，その確率密度関数を '
とし，さら
乱数 +
#!
$.
.(,の範囲が-他にも様々な関数形の確率分布を一様乱数から生成するに'0"%であることが既知であるとします．ここで，%)*方法が知られていますが，ここでは紹介しきれないので，は厳密な最大値である必要はなく，確実に '
0
"%となる
)
*

興味のある方は文献［１，２］などを参照してください．

ような適当な定数を設定すれば大丈夫です（近い方が乱数（佐竹）の生成効率はよいですが）．次に，一様乱数発生ルーチンを使って，二組の一様乱数 ()(#!
$.と%!
%.を用意し!
)(-５．
２ 乱数の利用例ます．これに対して，'()*と %))の大小関係を判定します．

ここからは，磁場によって閉じ込められた核融合プラズ

もし '
$%
()*
)
($(
)を採用し，そうでなければ)であれば 4

マの理工学における乱数を用いたシミュレーション手法をに対して同様の判定を行います．そし
破棄し，次の)()!
%)*

いくつか紹介します．シミュレーションにおいてどのようて，必要な個数の +4(,が採用されるまでこの判定を繰り返に乱数が用いられているのかを説明する前に，なぜシミュせば，確率密度関数を 'に従う乱数 +4(,が生成されます．

レーションで乱数を用いるのかについて考えてみます．

!0

この他にも，'
0)*$& や正規分布の場合など，特定の関まずは，シミュレーションによって扱われるプラズマの数の場合に使える乱数の変換法が色々と知られています．

諸現象が，どのような方程式により記述されているのか，

ここでは代表的なものとして，２次元空間における正規分

その概略を見ていきましょう．議論を単純化するために，

布'
3!
4*を与える，Box-Muller 法を紹介します．
)

完全電離プラズマが，それぞれ "##個の水素イオンと電子

Box-Muller 法では区間 !
!
".一様乱数を二組用意しまで構成されているとします．ここで，"&"!#"程度です．

す．それらを +1(!
2(,とします．この手法では 1(と 2(は無このとき，荷電粒子の運動は，Newton 方程式（および電磁相関で独立な一様乱数であることを前提としますが，第４場を解くために Maxwell 方程式）で記述することができま章で説明したような「品質の良い」疑似乱数であれば，実す［３‐５］．電場と磁場をそれぞれ "および !とし，)番目の際に独立な疑似乱数列を２系統使う必要はなく，１系統の荷電粒子の質量を +)，電荷を .)として，Newton 方程式か

疑似乱数列からその奇数番目を 1(，偶数番目を 2(，というら粒子の位置 $))00*および速度 #))*は，具合に取れば十分です．さて，平均が )，分散共分散行
!
!
!*

/
$
0
0
*
$#))
*
))

列が !" !"となる２次元正規分布関数は
! "

#

および

#

" !3 "4
')
3
!
4*
$ & #
&
#

/
#
0
*
$
))

で 与 え ら れ ま す が，こ れ を 極 座 標 系（/$&3#"4#，）で表現すると，
4#
3*
'$%
*!")

.)
")
0
!$))
0
0
0
!$))
0
+
*
*
"#))
*
#!)
*
*
,
+)

で与えられます．ここで，)$"!
!
!
" です．この New#$1!
ton 方程式と Maxwell 方程式を解くことができれば，原理

!/#

/
')
/
!
'*
$ &#
&
#

的には，プラズマの温度，密度，閉じ込め時間など，磁場閉じ込め核融合炉のプラズマ性能を評価するために必要なであり，'に関しては !
&*に一様に分布していることに!
#なります（ここで因子 /は )
3
!
4*
/
!
'*の座標変換のヤコ/)情報を得ることができるはずです．しかし，現実には，この方程式系を解くことは不可能です．なぜなら，全粒子間'ビアンから出てきます）．また，$//&!/##$"を満たす! 'の相互作用を考慮して，" 個すべての荷電粒子の Newton

ことは容易に確認できます．そこで，/方向に関する分布

方程式を解く必要があるからです．

#

を乱数から生成する方法を考えると，累積分布関数

そこで，プラズマの状態を近似的に記述する基礎方程式
４０６

Lecture Note

5. Transformation of Random Numbers and its Applications

S. Satake et al.

としては，水素イオンと電子のそれぞれに対して，２体相に近いと仮定しています．一方，プラズマの力学的平衡状

互作用のみを考慮した６次元位相空間内の１体粒子分布関態における磁場や，温度・密度分布そのものは，計算コス数%%!
$)の時間発展方程式である Boltzmann 方程式：)!
#(

トを考えて，通常，Boltzmann 方程式から導かれる３次元空間の流体方程式とMaxwell方程式を自己無撞着に解くこ$$ (# " )
$
*)
)
!
*)
%
(!
"$#!(
)
$ "%
! "$$ " (# %"#(
#)
$
)
$% &#
$$

とで求めます［１２］．この流体方程式に基づいて温度・密度分布を決定する輸送方程式が非線形偏微分方程式であるたが用いられます．ここで，#は粒子種を表し，"および !

めに，解を数値的安定に求めることが困難になる場合があは外場とプラズマの集団運動が作り出す電磁場とします．

ります［１３］．そのような場合に対処するため，５．
２．
２節で$個の Newton 方程式から Boltzmann 方程式を導出する議紹介する遺伝的アルゴリズムを用いた手法が開発されてい論については，文 献［３‐５］を 参 照 く だ さ い．上 式 の右辺ます．この場合，乱数は，物理現象のモデル化に用いるとは，２体相互作用（Coulomb 衝突）によって起こる速度空

いうより，非線形偏微分方程式を解くアルゴリズムにおけ間での拡散を表現する Landau 衝突項で，Fokker-Planck

るツールとして使用されているわけです．

型の式（＊１）として記述されます［３］．さらに，左辺にお以上のように，乱数は，運動論方程式を解く際の膨大ないて，荷電粒子のジャイロ運動に対して近似的な取り扱い

計算コストを減じるためのモデリングにおいて利用されたを行うことで，５次元位相空間の１体粒子分布関数の運動

り，流体方程式により記述される非線形な輸送方程式を解

論方程式（近似の取り方によって，ジャイロ運動論方程式，

くために利用されたりしています．

またはドリフ ト 運動 論 方 程 式）を 得 る こ と が で きます（菅野・本多）５．
２．
１ 拡散過程とその応用［５，６］．この運動論方程式によって，現実的にプラズマの諸現象を扱うことができるようになります．

本節では，乱数を用いた基礎的なモンテカルロシミュ

プラズマの振る舞いは，運動論方程式によって記述されレーションの代表的な例として，拡散方程式の解法とるわけですが，次に，シミュレーションにおいてなぜ乱数

Dirichlet-Poisson 混合問題への適用の２つを取り上げ，乱が用いられるのかを Landau 衝突項に注目して考えていき

数がどのように利用されているのか紹介します．紹介するます．例えば，プラズマ性能を決める主な要因の１つであ

シミュレーション手法の数学的基盤は，確率過程論です．

る微視的乱流輸送を扱うシミュレーションにおいては，運「確率過程」とは，時間と共に変化する確率変数（確率的に動論方程式の主に左辺の寄与により生じる１体粒子分布関値が定まる変数）です．この後に示すように，シミュレー

数の微細揺動に対する数値計算の高精度化のために，分布ションにおける乱数の役割は，数値計算上の単なるツール

関数そのもの を５次 元 の 連 続 体 と し て 扱 う こ と が多くではなく，確率過程論に基づいて拡散現象を粒子的にモデ［７］，Landau 衝 突 項 は，速 度 空 間 の 微 積 分 方 程 式 と しリングしたときに導かれる数学的性質の具現化を担っていて，差分スキームなどを使って数値的に解きます．このよます．

うな微視的乱流シミュレーションは，理化学研究所のスー

以下では，確率過程論における数学的結果（定理など）

パーコンピュータ「京」などの大規模計算機を利用して１００を用いた手法の説明に重点を置きますので，数学用語の厳

時間程度の計算時間を必要とします［８，９］．一方，最近の研密な定義や数学的結果の導出については，文献［１４‐１７］な究［１０］で，プラズマ性能を劣化させる不純物の輸送においどを参照ください．なお，本節で扱う手法では，簡単のたて，Coulomb 衝突の寄与による輸送現象（新古典輸送）のめ，空間を１次元としますが，'次元空間（'""）へ拡張

重要性が指摘されていますが，そのような現象に対するシできます．

ミュレーションにおいては，微視的乱流シミュレーション!拡散方程式の解法

ほどの速度空間分解能は要求されていないので，計算コス

プラズマ・核融合分野におけるモンテカルロシミュレー

トの低減を優先した粒子的なシミュレーション手法を用いションの代表的なものの１つは，拡散方程式の解法ではなることが多いです［１１］．その場合，Landau 衝突項によるいでしょうか．ここでは，荷電粒子の Coulomb 衝突や，熱拡散現象は，分布を構成するテスト粒子の速度空間におけ

拡散などのランダムネスが内在する物理現象に対する粒子

るランダムウォークによってモデル化されます．そのよう的なモデリングに基づいたシミュレーション手法を紹介しな Coulomb 衝突の効果を考慮したシミュレーションにおます．以下のような拡散方程式（Fokker-Planck 方程式）をいて，粒子のランダムウォークを具現化するため，乱数が解くことを考えることにします．

用いられるのです．このような乱数の利用に対する簡単な$%)$" $# #(
*)
!(
*)
%
)
!
*)
%
)
!
*)
(!
%! *
(
+
"
* *)
(
+
$
)
$
*
#$
*#
（１）

例を５．
２．
１節で紹介します．また，このシミュレーション手法を発展させると，科学技術分野で扱われることが多い
Dirichlet-Poisson 混合問題に適用することができ，そのよ

ここに，時間パラメータを )&!，空間を *とし，)!
*'#

うな応用についても紹介します．
上記のような運動論方程式を解く際には，多くの場合，で，#は 実 数 と し ま す．ま た，!(
，#(
*)
*)
'#お よ びプラズマを閉じ込めている磁場や，プラズマの温度および#(*)&!です．式（１）を解くため，分布 %を構成する粒子密度の分布は，既知である，もしくは，初期に与えたものの位置 *()))という量を導入し，*()の時間変化 の式を（＊１） Fokker-Planck 型とは，どのようなものかについては，５．
２．
１節で紹介します．

４０７

Journal of Plasma and Fusion Research Vol.96, No.8 August 2020

＋ノイズ項，と表現しましょう．ここ
+
#
-*
+
%-*
%+#!*
+
+
+

具体的な例で，拡散方程式を解くモンテカルロシミュ

で，この「ノイズ項」から生じるランダムウォークの歩幅レーションがどのようなものかを見ていきましょう．プラが拡散係数 #*-+を与えるように設定します．多数の粒子

ズマ中のテスト粒子（例えば，電子とします）が，背景プに対して，この式を解くと，時刻 +におけるそれぞれの粒

ラズマ（水素イオン と し ま す）と の Coulomb 衝 突 に よ子の位置 -*++が求まりますので，粒子の位置分布も得らり，速度空間をピッチ角散乱する場合を考えます．電子のれ，それが式（１）の解です．数学的に，もう少し厳密に言速さを ,，磁力線方向の速度を ,/として，*#,/#,で速度空えば，式（１）の解 &+!
-+を求めるために，式（１）に対する*間を表現することにします．*は，電子のピッチ角 &（速度

粒子描像の微分方程式：*+&を表します．ここベクトルと磁力線のなす角）の余弦 $

+
-*
+
-*
+
+
%-*
%+
%%*
+
#!*
+
+
")*
+
+
+

で，ピッチ角散乱によって速さ ,は変化しないと仮定して（２）います．*の取り得る値は，!"%*%"です．電子の速度空に注目し，この式を解いて得られる「確率過程 -*」の分++

間 *における分布 &を与える拡散方程式は，以下のように

布が解 &
+
!
-+であること（＊２）を利用します．ここで，*表現されます［１９］．

-+)*-+です．この確率過程 -*+#*-++は，拡散過程と#)*++&+" +# ( " # &
*+
*
&
*!
#! 0
!(
0 , !* -1 （６）
1
"
&'
+
*
+
+
#+
*# &'

も呼ばれます．式（２）の右辺第２項における%*++は，数学における意味での「Brown 運動」で，Wiener 過程とも呼ばれます［１４］．Wiener 過程 %*++は，①連続で，②離散時刻

簡単のため，ここでは，電子‐イオン衝突周波数 (
&'を定数

!#+
!+
#!
!"+
""+
#"4"+
(!""+
( を 適 当 に 取 っ て，%*と 仮 定 し ま す．式（６）に 対 応 す る 拡 散 過 程 **
，つ ま
+
+

として，任意の時刻 +++!%*( に対して %*(+(!"+が平均 !

り，テスト粒子の速度空間における位置の時間発展は，時

で分散 *
+
!
+
!
(!+
(!"+の正規分布 $*(!+(!"+に従い，③ま間ステップを %+として，式（２）のように与えられる確率

と%*
が互い
た，()
#'として%*
+
+
+
+
!%*
!%*
(+
(!"+
!"+
'+
'

微分方程式から自然に導かれる計算スキーム：に独立となるという性質を持ちます［１４］．式（２）の「ノイ

+
**
+
%+
%% （７）
**
+
+
"!*#*
+
+
"#(
"%+
+
#**
+
!(
&'
&',

ズ項」において %%*+#+%+のように書かない理由は，%*++

により与えることができます．ここで，%**
+
+
+
+
'**
"%+
+
!**
+

が至るところ微分不可能であるためです．
式（２）のような方程式は，確率微分方程式と呼ばれ，日としました．また，テスト粒子は，時刻 +#!で位置 **!+本の数学者，伊藤 清氏によって，１９４２年に数学的基礎が与から出発するとします．式（７）のような計算法は，Euler‐

えられました［１８］．式（２）は，確率過程 -*++を求めるた丸山スキームと呼ばれています［２０，２１］．もちろん，常微分

めの「積分」を表現していて，方程式の解を求める数値計算スキームと同様に，Runge-

+

+

!

!

Kutta スキーム（後で紹介します）など高次の近似法もあ

-*
*
-*
+
-*
*
*
%*
!+
%%*
+
+
"" )*
+
#-*
"" !*
+
+
+ （３）

ります［２０‐２２］．ここでは，Euler‐丸山スキームを利用することにします．

のように与えられます．式（３）の右辺第３項は，伊藤積分

式（７）のような計算スキームは，確率微分方程式の時間と呼ばれ，)*-**++における時刻 *の取り方に気を付けて，離散近似です．時間離散近似には，大別して，確率微分方+

程式の解の経路 ,**+++&!をよく近似する「強い近似」と解

)!"

%*
+
+
-*
+
+
,
!%*
" )*
-*
*
*
%%*
'
) ! )*
+
+
+
#(
(+
(""+
(+
%+3 !(#!
!
（４）の平均や分散などの分布特性を近似する「弱い近似」があり［２０］，%% に対して，強い近似では正規乱数 'を用いて%% #'%
2 +，弱い近似では，例えば，%% #$2%+のように取ります．ここで，弱い近似における $ は，コイン投げなのように定義されます．後で数値計算することを意識して，%+#+
%+，%*
#
)とし，+
+
+
#%*
"%% とし
( #(
(""+
(+

どで発生させた，値が ""または !"の二値乱数を用いて，

ま す．こ こ で，%% は，平 均０で 分 散１の 正 規 分 布

それぞれ確率 1/2 で与えられるとします．
（弱い近似にお$*!
!
"+に従う正規乱数'を用いて，%% #'%2 +と表現できます．確率微分方程式（２）の解の存在とその一意性のた

いては，%% に対して他の取り方もあります．
）どちらの

め，!*
-+および )*-+の満たすべき条件は，適当な定数切に判断すればよいでしょう．

近似を採用するかは，この後に示すように，問題ごとに適"!および "" に対して
!*
-+
)*
-+
-.
"".
.
.
".
.
%"!,
-

拡散方程式
（６）を数値シミュレーションで解く場合には，Euler‐丸山スキーム（７）における強い近似の採用は，（５a）あまり有効ではないことがわかっています．以下では，このことについて説明します．式（７）において強い近似を採

および
!*
-+
.+
)*
-+
.+
-!..
.
!!*
.
".
!)*
.
%"".

用すると，数値的に .
*.
$"となる不具合が生じます．すな（５b）

わち，数値計算において %+は有限な値であるので，１スです．ここで，１５］など!
.(!です．詳しくは，文献［１４，テップの計算の際に発生させた正規乱数の絶対値が，あるで論じられていますので，そちらを参照ください．

値'*.
$"が発生します．ここで，'!より大きくなると .
!

（＊２）証明については，例えば，文献
［１
４］
の問題８．
３を参照ください．この後で紹介する，式（１０）
を用いる計算法との関係も説明されています．

４０８

Lecture Note

5. Transformation of Random Numbers and its Applications

は，ある時刻 )における .の値を .&)"!('として，
'
#$&
)
!#

(",
"!,'
&
(&
#!('
#,

S. Satake et al.

は+
!
"!'
!""'
,の範囲に存在することになります．ここ
で，'#*
%)
#
#です［２３］．また，有限個のテスト粒子を用(*いた手法なので，得られた分布に統計的な誤差は生じま（８）

す．統計誤差は，テスト粒子数 #.に対し，"##- .に比例しで す．た だ し，!",#*%)""，!"("#で，ま た，(*て小さくなります．

.#$&"!('の 正 符 号 は .$"と な る ケ ー ス，負 符 号 は以上が，テスト粒子の位置分布によって拡散方程式の解."!"となるケースで選ぶこととします．.#"!(からを与えるモンテカルロシミュレーションの紹介ですが，こ出 発 す る テ ス ト 粒 子 の１ス テ ッ プ の 計 算 で 不 具 合れとは別の計算法もあります．関数 *&)!
+'を )$!で，（.$"）が発生する確率 $(,,は，$(,,#&
#
)
"
#
(,
)
&&
'
- 'で
!##

*&
)
!
+'
&&
+&
)
#%++
'
'
,

す．出発点を .#!""(とした場合も同じです．ここに，

（１０）

+'
+'で，(,+'は，誤差関数です．発生確率
(,
)
&&
)
)
#"!(,
&
&

であると定義すると，このとき，*&)!
+'に対して以下の関$(,,から，不具合は，*.*#"の近傍で比較的起こりやすいこ

係式が成り立つことを利用します［１４］．

とがわかります（＊３）．時間ステップ %)を小さくするこ/*#!*および *&
!
+'
+'
!
#&&
/
)

とで，１ス テ ッ プ 毎 の 発 生 確 率 $(,,は 小 さ く で き ま すが，ゼロにはならないので，根本的に解決したとは言えま（つまり，せん．数学的に %). !の極限を取れば，$(,,. !

（１１）

ここで，演算子!は，式（２）に対応して，式（７）を確率微分方程式として見たときには測度ゼロのイ#+'/"""&
+'/#
! #!&
/
+ #
/
+

ベント）となりますが，シミュレーションにおいてはゼロでない確率で発生するレアイベントとなります．そのた（１２）ありますし，どのような処理でも，拡散方程式（６）には含であり，%++&&+&)'',は，式（２）に従って時刻 (#!に位置+から出発した確率過程 (
+&
(
+&
)
'
)
'
'
(%!の分布による &&め，*.*$"となった場合の処理をどうするか考える必要がまれていない操作なので，*.*#"近傍のテスト粒子の分布の期待値で，この場合には，式（２）が与える時刻 )におけへの影響が心配です．また，拡散方程式（６）の基となるる +か ら ,へ の 遷 移 確 率 密 度 を '&)!
+!
,'と し て，Landau 衝突項のモデリングでは，小角度散乱を前提とし

!&&
,'
'&
)
!
+
!
,'
',と表現することもできます．最初に紹介ていますので［３，１９］，有限な%)の場合に，%% #)%- )におした計算法と式（１０）を用いる計算法は，Kolmogorov の前ける )の値によって散乱効果が大きくなりすぎること自進方程式と後退方程式の関係になっています［１４］．

体，レアイベントとは言え，モデリングの前提に反するの

先ほどのピッチ角散乱の問題に適用してみましょう．解ではないかという疑問もあります．

くべき方程式は以上の議論から，この問題では，Euler‐丸山スキームに

/#&
/
&
/&)
.'
. ""*
"!.#) #
&!
#!*
(
(*
/
. # (*/)/.

おける強い近似の採用は，あまり有効ではないことがわかりました．そこで，文献［２３］に従い，弱い近似で扱うことにし，%% #$-%)として，（１３）

で，対応する拡散過程は，式（７）（数値計算する際には，
!
.'
.'
!
&
#&&
#
式（９））で す．例 え ば，初 期 分 布 を &

.&
)
)
)
*
)
*
%)
%) （９）
"!.#&
"%)
'
#.&
'
!.&
'
'
)
$"(
(*
(*

+.'
#
#!*
+&
$で あ る と す る と，各 時 刻 )に お け る 解()

例えば，１ステップごとに !から "の範囲の一様乱数を

&
)
!
.'
&&
.&
)
&
#%.+
'
'
,は，図１のようになります．ここで，速度空間の位置 .における時刻 )の分布の値 &
)
!
.'は，式
&

のように計算することにします．ここで，$については，##
Tausworthe（トーズワース）法で発生させ，その値が "（９）により与えられる，時刻 )#!に .から出発したテス未満なら負符号，1/2 以上なら正符号とします．式（９）のト粒子の時刻 )での位置 .&).&)'における &&''の平均値でように計算すれば，*.*$"となる不具合は起こりません．

す．当然のことながら，図１の分布は，時刻 )#!でテスト

時刻 )#!で初期分布 &!
.'となるようにテスト粒子の初!
&粒子を初期分布 &&（９）によっ.'となるように配置し，式期位置 .&（９）によって与えられる時!'を設定すれば，式て時間発展させたテスト粒子の位置の分布と同じものにな

刻 )におけるテスト粒子の位置の分布が，拡散方程式（６）ります．最初に紹介した計算法の方が，シンプルでわかりの数値解 &)!
.'になります．時間ステップ%)を十分小さく&やすく感じるかもしれませんが，式（１０）を用いる計算法の取れば，分布の計算に関して，式（９）は，確率微分方程式

考え方は，より一般的な初期値境界値問題への応用に発展として見たときの式（７）の良い近似となり，シミュレー

させることができ，次に紹介する手法とも関連していまションにおいて，分布 &)!
.'は，初期分布 &
!
.'に依ら
!
&
&

す．詳しくは，例えば，文献［１５］で論じられていますので，!"（#*ず，１衝突時間 '以降に理論通り一様な分布に緩(*(*そちらを参照ください．
!Dirichlet-Poisson 混合問題への適用

和します．ただし，どのような初期分布であっても，.&)'

この節の後半では，科学技術分野で扱われることが多いの数値計算における最初の１ステップ以降は，テスト粒子（＊３）出発点を *（５b）を満たしていないため，解の一意性に問題が#"とした場合は，*.*.*$"となる不具合は発生しませんが，式あるように見えます．しかしながら，後で示す式（９）に従って，有限な時間ステップ %)で数値的に解く際には，右辺第２項がある
ために，そのような問題は生じません．

４０９

Journal of Plasma and Fusion Research Vol.96, No.8 August 2020

り，Runge-Kutta スキームを採用することで，これらを与える式（１６）に対する計算精度を向上させることができます．ここで，式（１６）に対する Runge-Kutta スキームは，以下の通りです［２１，
２２］．
" !
"
-'
+
+
)"")#(
%
"#+
(
$-'
(
" #+
" -#+'' ""!#(
"
#
#
（１８a）
!"$!'
-'
+
(
(
（１８b）

図１

ピッチ角散乱による電子分布 f (t, +)の時間発展を式（１０）に基づき計算．このシミュレーションにおいて，'ei#t ＝ 10−3と設定し，−1!+!1 の範囲の格子点（分割数５０）のそれぞれに対して，105 個のテスト粒子を用いて期待値 E+[ $(+(t))]
を求めました．実線は，初期分布 $(+) ＝ { 2 − sin ((+)}/4，記号＋は時刻 t/*eit ＝ 0.5 に お け る 分 布，記 号!は 時 刻
ei ＝ '
t/*ei ＝ 5.0 における分布です．

（１８c）

!#$!'
-'
+
%
(
"!"#+
")"-#+
#(
-'
+
%
)#$)'
(
"!"#+
")"-#+
#(

（１８d）
（１８e）

ただし，%
!
!
"(に従う互いに独立な正規乱数
"と %
#は，%'
です．これらの正規乱数は，例えば，Tausworthe 法で発生させた !から "の範囲の一様乱数を使って，
５．
１節で説明した Box-Muller 法で生成します．
以下の例題に適用してみましょう．
*#'*'-#""-" #!''%"-( !#'""-(!-$"'*# *-'"!-($&!#

Dirichlet-Poisson 混合問題に対して，モンテカルロ法をいかに適用するか紹介します．次の方程式を考えます．
!"&'
-(
,'
-(
-(
)
*
$#'

)"$)'
-'
+
(
(

（１４）

（１９）

!
），境 界 条 件 はた だ し，-&+!
!
",と し て（!$!
"$"
$'
!(
"(
$"お よ び $'$!と し ま す．式（１９）の 解 は，

ただし，演算子!および境界条件は，以下の通りとします．

'
-(
-"
-#"
#(
#(です．シミュレーションの結果を
'
$"!'
!'

#
,
!,$ ""'
-(,,
-(,
"!'
,
#
,
-#

図２に示します．ここで，図２では，得られた数値解の平（１５a）

滑化のため，最小二乗法を用いた多項式近似を行っていま

境界条件：-$-(-( （１５b）'$$'!!
" において，,す．多項式の項数の決定には，AIC（Akaike Information

ここで，空間 -については，-&+!!
",を仮定しています．

Criterion，赤池情報量規準）［２４］を利用しました．AIC は，式（１５a）から，これに対応する確率微分方程式は，次のよ

統計モデルの良さを評価し，良いモデルを選択するための

うになります［１４，１５］．

指標です．数値的に得た解 'を微分したい場合には，その

+
-'
+
+
-'
+
*-'
*&'
*+
(
$)'
(
(
(
"!'
(
(

ままでは数値解の統計誤差がノイズになるので，AIC など（１６）を利用しながら数値解を平滑化することが有効です．また だ し，"'-(-()'-(で す．方 程 式（１４）の 解 は，式$)'た，得られた数値解の妥当性は，異なる乱数列を使って複（１６）に従う確率過程 -'+(により数回，数値解を求めて確認します．例えば，２回計算して得られた数値解をそれぞれ '"，'#として，以下のような相

*

-

-'
*
-'
*
,'
-(
)'
+-,%' &'
**
(
(
(
(
$( #
&$
!

*

*

!

!

' #'
-'
*
-'
)
+-,%' &'
*)&**
(
(
(
(
!(-#
（１７）
$
と与えられます［１４，１５］．ここで，(-+#,は，式（１６）に従って，時刻 +$!に位置 -か ら 出 発 し た 確 率 過 程 )
-'
+
(
*
+%!
の分布による期待値を意味します．また，*は，式（１６）に従って動くテスト粒子（モンテカルロ法におけるサンプルの外へ最初に出た時粒子）それぞれに対して，区間 +!!
",刻によって与えられます．式（１７）により数値解,'-(を求めるには，すべてのテスト粒子それぞれに対する時刻 *の値が必要となるので，式（１６）の !'-(の値などから，計算時間がどれくらいか推定できても，乱数を用いてシミュレー

図２

ションしているために，計算時間を正確に定めることはできないという欠点があることに注意ください．
式（１７）を 用 い る 計 算 法 で は，各 テ ス ト 粒 子 の 経 路
-'
+
)
*
#!
$!
. と境界に最初に到達した時刻 *が重要であ
((
($"!
４１０

方程式（１
９）
の解析解と数値解．解析解を実線，数値解を破線で示しています［２５］．数値解は，式（１７）に基づき計算．
区間[0, 1]における格子点（分割数１００）のそれぞれに対して
103 個のテスト粒子を用いて，数値解を求めました．また，#t ＝ 10−4 と設定．得られた数値解の相対誤差は，-＜1/200
でした．

Lecture Note

5. Transformation of Random Numbers and its Applications

S. Satake et al.

敗し，!%
%&
%&かつ "#!とすれば解析解と良く一
#!###
%

対誤差
%
"

%
"

%
!

%
!

#
#
$#! '
%%! '
%%
"!#
#'
"'

致する数値解が得られました．もちろん，例題（１９）においても，左辺第３項を右辺に持っていき，ソース項として（２０）扱って問題無く解くことができます［２５］．
以上のように，本節では，確率過程論に基礎づけられたが十分小さいなら，数値解として妥当と言えます．
例題（１９）では解がわかっていましたが，もちろん，事前

モンテカルロシミュレーションを紹介しました．プラズに解がわからない問題を解くのが通常です．第４章で議論

マ・核融合分野では，ここで紹介したもの以外に，いろい

した，乱数の統計的品質が良いと判定される疑似乱数を利ろなモンテカルロ手法が開発・利用されています．例え用していれば，疑似乱数の性質が計算結果に影響すること

ば，文献［２６‐４４］などです．文献［２６‐３１］は Coulomb 衝突のは一般的には考えづらいです．それでも乱数の検定法で見モデリング，［３２‐３８］は運動論方程式の解法，［３９‐４２］は流抜けないような，疑似乱数の隠れた非ランダム性が計算結体方程式の解法，［４３］はプラズマ‐壁相互作用研究におけ

果に影響を与える可能性を排除したいのであれば，物理乱る動的モンテカルロ法，［４４］は中性子・光子輸送に関して数を利用することで，数値解の妥当性を判断する際に，そ書かれています．興味を持たれた方は，これらの文献を参のような影響を考えなくてよくなります．

照ください．
（菅野）

ところで，例題（１９）では，式（１７）における指数関数の引数に含まれる "が負の値で，計算の実行に問題は無かった

５．
２．
２ 輸送方程式求解への応用のですが，"が正の値を取る場合には，式
（１７）
のままでは

プラズマ輸送研究や統合モデリングにおいて，プラズマ

不具合（テスト粒子によっては，指数関数部分の数値計算の巨視的な時間発展や定常状態における分布を求めるためが破綻するなど）が発生することがあります．そのようなに輸送コードが広く用いられています．とりわけ，運転シ

不具合への対処の１つとして，式（１４）において "%
%&
$%
%&

ナリオ開発や実験解析においては定常状態を対象とするこを右辺に持っていき，ソース項 !%%&%&$%%&として処!"%とが多く，非定常輸送コードでプラズマの分布が変化しな理すると，数値解が落ち着くまで再帰的に計算を行う必要

くなるまで時間発展させて定常状態の分布を求めていまはありますが（つまり，ソース項における $%%&として，最す．輸送コードは放物型偏微分方程式である１次元の輸送

初は適当な近似解を与え，それにより得られた近似解をま方程式に基づいており，通常は有限差分法や有限要素法にたソース項の $%%&に代入して計算を繰り返すことで），う

よって時空間方向に離散化して数値的に解きます．温度分

まく克服できる場合があります［２５］．例として，以下のよ布を求めたい場合，圧力（内部エネルギー）を従属変数とうな方程式を考えます．

した輸送方程式を解きます．拡散係数（輸送係数とも言い

%## #
#!
"##
%%#

ます）がプラズマの内部状態に対して不変である場合，またはその応答が線形である場合は求解は容易です．しか（２１）

し，現実のプラズマの輸送は新古典輸送や乱流輸送によった だ し，%$(!
!
")と し て，境 界 条 件 は "%!&#"お よ びて支配され，それらは温度勾配などの熱力学的力に非線形

"%
%&
#%&
"&
&'%
#!"とします．式（２１）の解析解は， #%#$に依存しています．そのため，輸送方程式は非線形の偏微です．計算結果を図３に示します．ここで，#%&として最%分方程式となり，解析的に一般解を求めることは困難です

初に与えた近似解は，#%&%""です．シミュレーショ
%#!#し，数値的に解けたとしても数値的安定に求めることも困ンにおけるその他の設定は，図２と同じです．図３に示し

難になります．従来型解法の枠内で数値的安定に定常解をたように，式（１７）で !%%&#!かつ "###とすると計算に失求めるアルゴリズムが提案されてきましたが
［４５，
４６］，期
待した結果が得られないこともしばしばです．
本節では，輸送方程式を解いてプラズマの定常分布を求めるという問題を，定常状態を満たす温度と温度勾配の適切な組み合わせを求める大域的最適化問題へと転換し，その解法としてメタヒューリスティックなアルゴリズムである遺伝的アルゴリズム［４７，４８］を適用した例を紹介します．
ここまで読まれてきて，この話が一体どう乱数と関係するのか疑問に思われたかもしれませんが，遺伝的アルゴリズ
ムの中で乱数は重要な働きをします．以下では，まず遺伝的アルゴリズムの簡単な紹介をしたのちに，遺伝的アルゴ
リズムを用いた輸送方程式の求解法へと進んでいきます．
!遺伝的アルゴリズムとは

図３

遺伝的アルゴリズムは進化的アルゴリズムの中でも代表方程式（２１）の解析解と数値解．解析解を実線，式（１
７）
で
７）
で
F (x) ＝ 0 か つ "＝ #2 と し た 場 合 は 点 線，式（１F (x) ＝ − #2f (x)かつ "＝ 0 とした場合の数値解（得られた数値解の相対誤差 $＝ 1/100）を破線で示しています［２５］．

的なものであり，自然界における自然淘汰，遺伝，交叉，突然変異などの生命の進化過程を数値的アルゴリズムで模したものになります．遺伝的アルゴリズムでは必然的に進化
４１１

Journal of Plasma and Fusion Research Vol.96, No.8 August 2020

遺伝学の用語が多く使われますが，あくまで数値的な最適言います）の集団を生成します．扱う問題にも依りますが，

化アルゴリズムであり，生物学の文脈で使われる概念とは一般的に多くのテスト粒子数を必要とするモンテカルロ法

必ずしも厳密に一致しませんのでその点はご留意くださとは異なり，個体数は典型的に１００程度です．アルゴリズム

い．一口に遺伝的アルゴリズムと言っても実装は様々かとの流れは以下の通りです．
（１）パラメータ空間の中で個体思いますが，以下では PIKAIA［４９，５０］という Fortran 言語

をランダムに生成し，各個体の適応度を評価します．
（２）

における実装を基に，遺伝的アルゴリズムがどのようにし次に，適応度に応じて現在の集団からつがいを選択し交配て大域的最適解を求めるのかを見ていきましょう．

させ，次世代の集団を作ります．
（３）次世代の集団を現世最適解，すなわち解を求めるためには，当然ですが最適代と置き換え，（４）新しい集団の個体に対して適応度を評さの基準が必要になります．この品質の尺度を適応度と呼価します．集団の中で最適な（最大の適応度を持つ）個体

び，適応度を求めるための関数を適応度関数 #!*と言いま)が，設定した基準を満たすか設定した世代数に到達するます．大 域 的 最 適 化 問 題 を 解 く と は，パ ラ メ ー タ 空 間

で，（２）‐（４）の過程を繰り返します．最終世代の集団で最!*を最大化するパラメータ !を見つけるこ
!(.
"
!
"/で #
)
とを言います．なお，!は適当な規格化によって.
"
!
"/に収

適な個体の持つパラメータ &
!*
,
)!+&
-#
)
"が最適解とな
!""/
!(.
ります．
以上が遺伝的アルゴリズムのまるようにしています．

アルゴリズムになりますが，これだけだと乱数との関わりがはっきりしないかも

まず，これから生存競争を繰り広げる個体（表現型とも"#########################

乱

数

茶

話

５．
２．
１節で紹介したモンテカルロシミュレーションでです．

"#########################"#########################
は，拡散方程式を解く際に，対応する確率微分方程式に

試しに，拡散方程式（１３）で考えてみましょう．*)%*

従うテスト粒子の運動を利用しました．例えば，拡散方を−1/10 から＋1/10 の範囲の一様乱数として，次のよう

程式（１３）に対応した確率微分方程式（の計算スキーム）に式（９）に加えてみます．

として式（７）（数値計算する際には，式（９））を用いましたが，これ以外にも式（１３）に対応する確率微分方程式は

%
+
.)
%
%
.)
%
%
+
&%
&%
"!.#)
*
,
"&%
*
#.)
*
!+
*
"*)
*
,
$$+
(*
(*

あるのでしょうか？このことに関して，文献［１４］に興味

この式を用いて，%..
')
.)
%
*
*
/を計算すると，例えば，深い問いが載っています．それは，「いつ伊藤過程は拡%$!
$における分布は，図４のように図１と同じ結果(*#!

散過程となるか？」です．伊藤過程とは，次のように与になります．.%"とならない範囲で，*)%*-の最大値をえられる確率過程 ')%*を 指 し ま す．よ り 正 確な定義変えてみても，結果は変わりません．つまり，このようは，文献［１４］を参照ください．

な一様なノイズは，分布 #
%
"
.*
')
.)
%
)
#%..
*
*
/に影響し
ません．

%
%
%
%
'')
'%
'")
*
#()
*
"))
*
*

ここで紹介した文献［１４］の数学的結果は，微視的運動

ただし，()%%*お よ び ))*も 確 率 過 程 で，任 意の時刻

から巨視的スケールの分布が従う方程式へのモデリング%&!に対して確率１で以下を満たすとします．

にとって，重要な意味を持つように思えるのですが，い

%

%

!

!

かがでしょうか．
##()$$'$#' および # ))
*
*'$#'

（菅野）

文献［１４］による「問い」への答えですが，以下の拡散過
程 &)
%
*
%
&)
%
&)
%
%
'&)
'%
'")
*
#!)
*
*
",)
*
*
*

と ')
%
%
*との間に，()*に対する次の条件付き期待値につ
いて
!%/
()
%
'&)
%
%&.
*
#!)
*
*

かつ，))
%
*について
%
,)
'&)
%
))
%
))
%
'&)
*
*
*
*
*
*
#,)
図４

が，ほとんどすべての時刻 %，およびほとんどすべてのテスト粒子に対して成立する場合，&)
%
%
*と ')
*は，法則
の意味で一致する（つまり同じ分布を持つ）というのです．ここで，'&)%%*は %#!で位置 &から出発した ')*を意味し，!%は +; $%%‐加法族［１４］')$*,の生成する ,ピッチ角散乱によって時間発展する電子分布に対する一様なノイズの影響．実線は図１の時刻 t/ei ＝ 0.5 における分布，記号〇は一様 な ノ イ ズ を 加 え た 場 合 の 時 刻t/ei ＝ 0.5 における分布．一様乱数 *(t)は，Tausworthe法を利用して，−1/10 から＋1/10 の範囲で与えました．初期分布，各パラメータの設定値，およびテスト粒子の総数は，図１と同じです．

"#########################
４１２

Lecture Note

5. Transformation of Random Numbers and its Applications

S. Satake et al.

しれません．そのため，乱数との関連という観点から各段としましょう．座標値は有効数字４桁で取ったので，合わ階をもう少し詳細に見ていきたいと思います．

せ て８つ の 遺 伝 子 を 持 っ た 個 体 か ら 各 々 の 染 色 体

まず，個体とは何かについて説明します．個体とは遺伝

SP1: 12345678，SP2: 86427531 を作ります．このエンコード子の集合からなるものとして捉えられており，最適化した過程は自明かと思います．次に，一点交叉と呼ばれる演算

いパラメータの組の数だけ遺伝子の組があると思ってくだを行います．両親の染色体のつがいに対して乱数によってさい．例えば，２次元空間 &（ここで，%
）に
%
"
&'
"
&%(
"
!
")

決まった位置で染色体を切断し，付け替えます．この場合

おいて&
!
"
!
%
"
&'を求めるとします．念!
&!
&'に最も近い !$&ですと，一様整数乱数 "%(""))を生成し，例えば４が選ばのためですが，理想的な答えは&
!
"
!
!
&
!
&'なので自明です！

れたとすると，３番目と４番目の遺伝子座の間で組み替え

あくまで遺伝的アルゴリズムの概念 を わ か り や すく伝が 起 こ る た め，染 色 体 は そ れ ぞ れ SO1: 12327531，

え る た め に 設 け た 問 題 で す．適 応 度 関 数 の 設 定 方

SO2: 86445678 となります．PIKAIA では世代ごとに集団の法は色々あり ま す が，ユ ー ク リ ッ ド 距 離 の 二 乗 の逆数

数は増減しないため，２人の親から２人の子が生まれるこ

#

# !"

!'
$
%!!
!
&!!
!
&'"&
&') とすれば良いことがわかり&#(&とになります．つまり，これらは子の染色体となります．

ます．%， &は何桁かの数字（遺伝子）からなっており，今更に，ある変異率（通常は 0.005 程度です）で突然変異を起回の例では座標値になる訳ですが，その初期の座標値をこします．染色体を構成する遺伝子ごとに生成した一様乱"!
")の一様乱数から生成し，個体数分だけ繰り返します．
(

数 #%(
"
!
")が変異率を下回った場合に，その遺伝子を一

これが（１）で行っている作業です．オリジナルの PIKAIA

様整数乱数 "%("!
*)で置き換えます．これを一様な一点では線形合同法による乱数を用いていますが，筆者はそれ突然変異と呼びます．突然変異のさせ方は他にも色々あり

を Mersenne-Twister（メルセンヌ‐ツイスター）に換装しますが，詳細な説明は割愛します．ここでは SO1 に対してて使用しています．

は６番目の遺伝子が９に，SO2 に対しては２番目の遺伝子（２）の過程は遺伝的アルゴリズムの中でも最も重要なもが３になったとしましょう．これをデコードすると，子はのです．個体の集団の中から，１組のつがい，すなわち父それぞれ O1:&
，O2:&
!
"
!
!
"
!
!
"#$#
!
(*$"'
!
)$%%
!
&'()'となりまと母を選択します．父と母をランダムに選択するのではな

す．こうして生まれた子世代は一旦プールされ，親世代とく，適応度に応じて選ばれる確率を変える必要がありそう同数になるまで上記の過程が繰り返されます．なお，ランです．そのために，まずルーレット選択と呼ばれるサンプ

キング選択方式からわかるとおり，適応度の高い個体は１リング手法を用いると，どうなるか見てみましょう．そのつの世代で何度でも親に選ばれることがあります．

方法ではまず，全個体の適応度の和 !と，個体番号が若い

子世代に個体が出揃うと，（３）の世代交代を行います．

方から順に定義される適応度の累積分布を作成し準備して最も簡単で広く使われている方法は，集団を完全に子世代おきます．そして，一様乱数 #%("!)を生成すると，その!

の集団で置き換えてしまう，というものです．但し，エ乱数#は必ず離散的な累積分布のどこかの区間に対応する

リート選択という手法を使う場合は，親世代で最も適応度値を取り，その区間を示すインデックス値から特定の個体が高い個体は無条件で子世代へと引き渡されます．交叉や番号と一対一対応が付けられるために，乱数によって個体突然変異によってその世代での最適個体を壊してしまうのを選択できます．ルーレット盤は個体数で分割されておを避けることで，収束を悪くしないための工夫です．子世り，各個体の適応度の大きさによって区間の面積が決定さ

代からは１個体がランダムに排除されてしまいますが，仮れているイメージです（イメージしにく い 方 は［４８］のに選ばれた子世代の個体が親世代の最適個体よりも適応度

図２．
５をご覧ください）．これで問題は全くないように思いが高い場合は，親世代の最適個体は次世代に引き継がれまますが，実は適応度の大きさに応じて決めてしまうと，多せん．

くの問題が生じることが知られています．わかりやすい例

最後に，こうして作られた子世代の適応度を評価するのを挙げると，１つか２つの個体が極端に高い適応度を持っが（４）です．P1，P2，O1，O2 の適応度を評価すると，O2ておりその他がそうでない場合，このサンプリング方式でが最も大きな値を示していることがわかります．これははほぼ常に適応度の高い同じ個体が選ばれてしまうことにO2 が &!
"!
!
&!
&'に最も近い個体であることを示しており，なり，最適化過程にとって重要な集団の多様性が失われてこの例では，１世代を経てよりしまいます．それを避けるために，ランキング選択と呼ばたことになります．参考までに，このシミュレーションのれる方法が使われます．適応度の高い順に個体を並べるの結果を図５に示します．図５では各世代において最適個体

ですが，個体の選ばれやすさはランク（順位）の高さで決の持つ座標値が世代ごとに示されています．既に第１世代

まり，適応度の値そのものでは決まらない，という方式でから正解にかなり近く，１０世代程度で正解と見なして申しす．つまり個体間での選択確率の差が同じである，という

分ない精度に到達し，第２５世代で完全に正解に到達しま

わけです．ランクを横軸に取り，相対適応度を縦軸に取る

す．なお，適応度関数 $の定義から，&
!
"
!
&
!
&'における適!

と，右肩下がりの直線を引くことができますが，この勾配応度は無限大になりますが，数値的に工夫することで無限を制御することによって，淘汰圧の強さを調整することが大を避けるようにしています．ただし，遺伝的アルゴリズ

できます．

ムはメタヒューリスティックアルゴリズムであるため，正優秀な子どもが生まれこ の よ う に し て 選 ば れ た つ が い の 父 を解に必ず到達する保証はないことに注意してください．こ

P1: &，母を ,#+
%
"
&'
!
"
!
%
"
&'
!
"
!
!
"#$%
!
&'()'
!
)'%#
!
(&$"'
#&
&
#&

の例では &!
"!
!
&!
&'にたどり着きましたが，選択する乱数
４１３

Journal of Plasma and Fusion Research Vol.96, No.8 August 2020が算出する輸送係数に大きな影響を与えてしまいます．そこで微分に依らない解法である遺伝的アルゴリズムに着目し，大域的最適化手法を用いた定常輸送方程式の解法の開発に取り組みました［１８］．
以下では，簡単のため１粒子種のみを考えます．定常の輸送方程式は'# '$ &""%!#! ". '.!!))
+
&%+
*
' '
'
%
%

（２２）

と書けます．%は規格化小半径，' は体積，.は %微分，)は密度，$は拡散係数，&は温度，%は加熱パワーなどに対#応するソース項になります．)+&%+*は平衡計算で求められ

図５

各世代における最適個体の持つ遺伝子（x，y座標）．最初の１０世代でほぼ正解の(0.5, 0.5)に収斂し，２
５世代で完全に正解と一致します．

る メ ト リ ッ ク で あ り，こ こ で は 既 知 で あ る と し ま す．

シードによっては 0.4999 などで止まってしまうこともあり

分方程式となります．'")"%が既知である場合に，この方

)"
)."
$#$'
&"
&."
/(は一般に密度，温度やそれらの勾配などに非線形に依存しているため，式（２２）は非線形の偏微ます．しかし，それでも十分正解と見なせるでしょう．こ

程式を満たす &を求める，というのがここで考える問題でこまで簡単に遺伝的アルゴリズムの実装を見てきました

す．ここで，式（２２）の括弧内は熱流束 $ であり，が，数多の過程において一様乱数が用いられていることが

# '
$'
%(
$&
#!))
+
&%+
*
'
%

わかるかと思います．
!大域的最適化手法を用いた輸送方程式の求解

（２３）

を使って

遺伝的アルゴリズムは大域的最適化手法の一つとして位置付けられており，多峰性を持ったパラメータ空間におい

'
'.$(
!#! '
"%
'
'

て最適解（最大値とします）を見つけることを得意としています．多峰性の例として，*%,"!
"-における正規分布の重# #ね合わせからなる関数 +'
*(
*!*
&(+*,
###
!'
((#
(#"!()（２４）と書き換えられます．式（２４）を空間積分すると，を 考 え，そ の 最 大 値 を 見 つ け て み ま し ょ う．こ こ で，

%

'.$'
%(
%(
#$ '.%(%$#'

!
"
!
"
!
"
!
"
!
"
!
*
&""
!#"
*
&#!""
!
'
!
#
!
!#%
!
$
!
&
!
"%-と し ま
#,
,
""
#"

!

!
#を中心とする急峻な正規分布の頂す．最大値 0.9 は *#!

（２５）点に位置する一方，幅広い裾野を持った正規分布に極大値となります．ある磁気面を横切って外側に流れ出る熱流束あります．Newton 法などに代表される局所最適化手法では，その磁気面より内側に吸収されたパワーと等しいこと

は，極大値ではなく最大値を見つけられるかどうかはひとを表しています．式（２３）の温度 & を変化させることで式

えに初期値の選択に掛かってきますが，問題を解く前に最（２５）を満たすように$を変えることから，##'.を目標熱流

大値にたどり着く適切な初期値を知っていることは稀で束と呼びます．

す．局所勾配に依存する最適化手法では，*!0.282 に初期

では，どのように式（２５）を満たす &を見つければ良いの値を取ると極大値 0.3 へと向かってしまい真の最大値 0.9でしょうか．通常の手法ですと，空間方向にメッシュを作を見つけることはできません．言い換えると，ランダムに成し，有限差分法や有限要素法で空間離散化し，非線形項

初期値を選んだ場合，７割強の確率で最大値を求められなである $を線形化して & に対する連立一次方程式を作いことになります．大域的最適化手法は局所最適化手法とり，それを直接法ないし間接法で解き，得られた解（温度）比べ一般に計算コストが掛かる一方で，アルゴリズムにもを再び $に代入して連立一次方程式を解き直し，その作業依りますが局所最適化手法よりも広くパラメータ空間を俯を温度が収束するまで繰り返します［５１］．ここではそのよ瞰する特徴があります．極大値近傍に適応度の高い個体がうな一般的な方法を採らずに解&を求める方法を考えてみ!
#近傍に個集まっている状況下でも，突然変異などで *#!

ましょう．温度勾配を規格化した "
&#
'
%(
#
& で式
#
"& $'
!'

体が生じるとそちらの方が高い適応度を持つため，局所最（２３）を表記すると，

適解から抜け出して最大値を見つけ出せます．

#
$'
%(
$'
&"
#
"& (
& " # #.
"
#!))
+
&%+
*
"& '

遺伝的アルゴリズムは最大値を見つける際に微分値を必要としません．そのため，微分を高精度に評価することが（２６）となります．最後の等式は式（２５）から来ました．式（２６）で困難であったり，僅かな微分値の違いが結果を大きく左右してしまったりするような問題であっても，安定に解を見は $の&"#"& の依存性を陽に書きました．温度と温度勾配"つけることができます．そのような問題の例としては，分には物理的に現実的な定義域がありますから，そのパラ布の勾配に強い依存性を持つ輸送モデルを用いる場合が挙メータ領域内で各%において式（２６）を満たす'&"#"& (の組"

げられます．分布の勾配という微分の精度が，輸送モデル

み合わせを見つけることができれば，それが式（２２）の解と４１４

Lecture Note

5. Transformation of Random Numbers and its Applications#

S. Satake et al.
#

$'!$)
"!%,#
'
"， '
"
"#!
##!
$'
"

なります．この手法だと，直接温度の微分値である ""!$を探しているため，微分方程式の解を求めるのに微分が不（２８）を評価し，最終的に遺伝的アルゴリズムによって要であることがわかると思います．
解の組み合わせを求める手法ですが，極端なことを言え

!"
'
'
'
($
*%
#*
+
"!
#&

ば，そのパラメータ領域内を極度に細かく分割して全ての（２９）組み合わせに対して当てはめを行えば，どこかで解は見つを 最 大 化 す る 個 体 を 見 つ け ま す．そ の 個 体 の 持 つかるでしょう．現実的にはそこまで細かな分割はできませ$!
"!$ &が最終的な式（２２）の解となる訳です．'
"
%
#の意味は

んので，計算精度と掛かる計算時間を秤にかけて妥当な刻わかりにくいかもしれませんが， '#が小さいとはすなわ

み幅を決めざるを得ません．非線形性の強い輸送モデルをち，ある温度 $)を入力として用いて式（２７）を計算した時に用いる場合，僅かな勾配の差が大きな熱流束の差を生んで得られる$'が入力の$)とがほぼ同じである，ということをしまうため，先程まで極端に細かくないにしてもやはり十意味しています．そのような $ は式（２７）を 満 た し て い

分細かな分割が必要になります．規格化温度勾配の定義域

る，ということがわかるかと思います．

を*!
!
"!+として，刻み幅を 0.01 とかなり粗く取ったとして!計算例

も，１０００メッシュ必要になります．温度に対しても同じく

上記のアルゴリズムで温度分布を解いた計算例を一つ示

6１０００メッシュ用意すると，組み合わせは 10 にも上ります．

します．計算は，上記のアルゴリズムに基づく定常輸送

これを，空間分割数分繰り返さなくてはいけないため，膨

コード GOTRESS で行いました．加熱分布や平衡は JT-60U

大な試行が必要になります．さらに，この組み合わせの中のとある放電実験のデータから持ってきました．輸送モデには解として全くあり得そうに無い組み合わせも存在してルには CDBM［５２］を用い，径方向５０点で１粒子種（電子にいるため，無駄の多い手法であることがわかります．その相当します）の温度分布の計算を PC クラスタで実行した

ような無駄を省くために，この解の組み合わせを見つける

ところ，約 2.77 秒で計算を終えました．平均すると，各径アルゴリズムとして大域的最適化手法である遺伝的アルゴ

方向点当たり２１２世代で収束解が得られていることになりリズムを用いるのです．遺伝的アルゴリズムでは，最初こます．なお，世代ごとの個体数は１００であり，CDBM の評価

そモンテカルロ的にランダムに値を生成しますが，そこかを行った回数の総数は１
０６万６千回になります．図６が示らはより解に近い（適応度の高い）組み合わせを持った個

すとおり，各径方向位置において目標熱流束に対応する加体がより生き残るため，急速に最適解へと収斂していきま

熱パワー累積分布 "（実線）と輸送熱流束に対応するパす．

ワー累積分布 %,#（丸印）が一致していることがわかりま,

熱流束を求める時は， 正解 である目標熱流束 ""%がす．これは，式（２６）を満たしていることを示しています．

事前にわかっていたため，# を目標熱流束に合わせることその時の規格化温度勾配分布は図７に，温度分布は図８にができたのですが，正解がわかっていない温度 $はどのよ

示されています．温度分布や規格化温度勾配分布が共に滑うに求めたら良いでしょうか．一つの考え方は，""!$ の値らかに得られていることがわかるかと思います．この時のが定まったのだったら，それは温度に対する常微分方程式

規格化温度勾配分布は温度勾配を何らかの微分スキームでになるため，境界の温度さえわかれば後は全て積分で求め微分したものではなく，温度分布と独立に得られたものでられる，というものです．しかし，この手法は上手くいかあることを強調しておきます．また，図８には拡散係数分ないことがあることがわかっており［１８］，別の手法を考え布も示されています．GOTRESS 内部では新古典熱拡散係,る必要があります．式（２３）で % ##"であることに留意し数も計算しているため，CDBM で評価された乱流熱拡散係て外側境界から積分してみましょう．すると，数との和が示されています．なお，磁気軸においては規格

$&

$%
$&
$%&
#$%
"#

$

&$
"
,
#
('
#%
$!
"
!$ &%
"
)
$$)
(

化温度勾配，つまり輸送熱流束は０であることが自明であ（２７）という関係が得られます．$%$%&は計算領域の外側境界$#$%での固定された温度であり，Dirichlet 条件に相当するものです．式（２７）は非線形の積分方程式であり，被積分
関数中の #に陰に含まれる$（$)とします）と左辺の$（$'）
が同一となるような $を見つけなければなりません．本来
このような解を見つけることは難しいのですが，大域的最適化手法によって比較的容易に見つけ出すことができます．ところで，式（２７）は ""!$ にも陰に依存しているため，"!$ も適切に定まらなければ収束解にはたどり着きませ"ん．つまり，式（２６）も同時に満たさなくてはいけません．
（２６）
および（２７）を計算
各個体の持つ %$!
"!$ &に対して式"し，得られた # と $ に対して

図６

４１５

目標熱流束に対応する加熱パワー累積分布 P と輸送熱流束に対応するパワー累積分布 V´Q ．

Journal of Plasma and Fusion Research Vol.96, No.8 August 2020

参考文献

図７

図８

［１］津田孝夫：モンテカルロ法とシミュレーション（三訂版）（培風館，１９９５）．
［２］伏見正則，逆瀬川浩孝（監訳）：モンテカルロ法ハンド
ブック（朝倉書店，２
０１
４）．
［３］D.R. Nicholson, Introduction to Plasma Theory (Wiley, 1983);
小笠原正忠，加藤鞆一
［訳］：プラズマ物理の基礎（丸
善，１
９
８
６）
．
［４］D.G. Swanson, Plasma Kinetic Theory (Chapman and Hall/
CRC, 2008).
［５］プラズマ・核融合学会［編］：プラズマシミュレーショ
ン（京都大学学術出版会，２０
１
８）．
［６］渡邉智彦，洲鎌英雄：プラズマ・核融合学会誌 81, 534(2005).
［７］渡邉智彦：プラズマ・核融合学会誌 81, 686 (2005).
［８］前山伸也：プラズマ・核融合学会誌 91, 589 (2015).
［９］S. Maeyama et al., Phys. Rev. Lett. 114, 255002 (2015).
［１０］M. Nunami et al., Phys. Plasmas 27, 052501 (2020).
［１１］K. Fujita et al., Plasma Fusion Res. 14, 3403102 (2019).
［１２］J.P. Freidberg, Ideal Magnetohydrodynamics (Plenum, 1987).
［１３］M. Honda, Comput. Phys. Commun. 231, 94 (2018).
［１４］B. !ksendal, Stochastic Differential Equations: An Introduction with Applications (Springer, 2003); 谷口説男
［訳］
：確
率微分方程式―入門から応用まで（丸善，１
９
９
９）
．
［１
５］A. Friedman, Stochastic Differential Equations and Applications (Dover, 2004).
［１６］保江邦夫：確率論（数理物理学方法序説４）
（日本評論
社，２
０
０
１）
．
［１
７］保江邦夫：物理数学における微分方程式（数理物理学方法序説別巻）（日本評論社，２
０
０
２）
．
［１
８］高橋陽一郎［編］：伊藤清の数学（日本評論社，２
０
１
１）
．
［１
９］P. Helander and D.J. Sigmar, Collisional Transport in Magnetized Plasmas, Chapter 3 (Cambridge University Press,
2002).
［２０］P.E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations (Springer, 1999).
［２１］T. Misawa and H. Itakura, Phys. Rev. E 51, 254 (1995).
［２２］J.R. Klauder and W.P. Petersen, SIAM J. Numer. Anal. 22,1153 (1985).
［２３］A.H. Boozer and G. Kuo-Petravic, Phys. Fluids 24, 851(1981).
９８３）．
［２４］坂元慶行 他：情報量統計学，第４章（共立出版，１［２５］R. Kanno et al., Plasma Fusion Res. 6, 2403066 (2011).
［２６］T. Takizuka and H. Abe, J. Comput. Phys. 25, 205 (1977).
［２７］K. Nanbu and S. Yonemura, J. Comput. Phys. 145, 639(1998).
［２８］C. Wang et al., J. Comput. Phys. 227, 4308 (2008).
［２９］Y. Masaoka and S. Murakami, Plasma Fusion Res. 8,2403106 (2013).
［３０］G. Zhang and D. del-Castillo-Negrete, Phys. Plasmas 24,
092511 (2017).
［３１］S. Satake et al., Comp. Phys. Comm. 255, 107249 (2020).
［３２］X.Q. Xu and M.N. Rosenbluth, Phys. Fluids B 3, 627 (1991).
［３３］W.X. Wang et al., Plasma Phys. Control. Fusion 41, 1091
(1999).
［３４］S. Brunner et al., Phys. Plasmas 6, 4504 (1999).
［３５］S. Murakami et al., Nucl. Fusion 40, 693 (2000).
［３６］S. Murakami et al., Nucl. Fusion 46, S425 (2006).
［３７］S. Matsuoka et al., Phys. Plasmas 22, 072511 (2015).

規格化温度勾配 1/LT の定常解．

温度 T の定常解とその時の拡散係数 !．拡散係数は CDBM
で評価された乱流熱拡散係数に新古典熱拡散係数が重畳されたものです．磁気軸での拡散係数は０としています．

るため，磁気軸では計算しておらず，磁気軸での温度は温度勾配零の条件から算出しています．そのため磁気軸での拡散係数を算出する必要はなく，０と置いています．
本節では，一見およそ繋がりがなさそうな遺伝的アルゴ
リズムと偏微分方程式の求解が，解法の工夫によって結びつけられ，輸送方程式を解くという実際の問題へと適用できることをご 紹 介い た し ま し た．輸 送 コ ー ド と してのGOTRESS は複数粒子種を扱えたり並列計算が行えたりす
るなど多くの特徴を備えています．詳細は［５３］を参考にしてください．
（本多）

謝

辞

５．
２．
１節 に お け るシ ミ ュ レ ー シ ョ ン で 利 用 し た Tausworthe 法に基づく疑似乱数発生プログラムを提供いただ
いた，高丸尚教氏（中部大学 工学部 教授）に感謝いたします．また，５．
２．
１節に対して，松山顕之氏（量子科学技術
研究開発機構 六ヶ所核融合研究所 主幹研究員）より，貴重なコメントをいただきました．ここに感謝申し上げます．

４１６

Lecture Note

5. Transformation of Random Numbers and its Applications

S. Satake et al.

!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

しんすけ

佐 竹 真 介自然科学研究機構 核融合科学研究所 ヘリ
カル研究部 核融合理論シミュレーション
研究系 准教授，２００
３年総合研究大学院大学 博士（学術）．
モンテカルロ法を使った３次元磁場配位中の新古典輸送現象，新古典粘性のシミュレーションや最適化配位の研究が主なテーマ．乱数については深い思い入れがありますが，特にギャンブル好きというわけではありません．

!!!!!!!!!!!!!!

たけ!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!

ほん

だ

みつる

充

量子科学技術研究開発機構 那珂核融合研究所 先進プラズマ研究部 上席研究員．主にトカマク中の輸送現象や輸送シミュレー
ションを研究対象としていますが，機械学
習やベイズ推定，最適化問題にも興味を持っていて研究への応用を進めています．家族でカラオケボックスに行くことが多かったのですが，昨今の情勢でなかなかままなりま せん．９
５点以上を獲れる曲を増やすのが目標です．

!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

本 多

かん

の

りゅう

た

ろう

菅 野 龍 太 郎
自然科学研究機構 核融合科学研究所 ヘリ
カル研究部 核融合理論シミュレーション
研究系 准教授．乱数を用いた確率論的な計算手法であるモンテカルロ法全般に興味があります．最近の研究では，プラズマの衝突輸送現象に対するモンテカルロ法に基づいたドリフト運動論シミュレー
ションを行っています．

!!!!!!!!!!!!!!!!!!!!!!!!!!

さ［５２］M. Honda and A. Fukuyama, Nucl. Fusion 46, 580 (2006).
［５３］M. Honda and E. Narita, Phys. Plasmas 26, 102307 (2019).

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!

An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence, second ed. (MIT Press, Cambridge, 1992).
［４８］C.R. Reeves, J.E. Rowe, Genetic Algorithms: Principles andPerspectives (Kluwer Academic Publishers, New York,
2002).
［４９］P. Charbonneau, Astrophys. J. Suppl. Ser. 101, 309 (1995).
［５０］http://www.hao.ucar.edu/modeling/pikaia/pikaia.php［５１］J. Candy et al., Phys. Plasmas 16, 060704 (2009).

!!!!!!!!!!!!!!!!

［３８］H. Yamaguchi and S. Murakami, Nucl. Fusion 56, 026003(2016).
［３９］Y. Feng et al., J. Nucl. Mater. 266-269, 812 (1999).
［４０］A.M. Runov et al., Phys. Plasmas 8, 916 (2001).
［４１］M. Kobayashi et al., Contrib. Plasma Phys. 44, 25 (2004).
［４２］R. Tatsumi et al., Plasma Fusion Res. 15, 1403003 (2020).
［４３］プラズマ・核融合学会［編］：プラズマシミュレーショ
ン，第７．
４．
１節（京都大学学術出版会，２
０
１
８）
．
［４
４］桜井 淳：日本原子力学会和文論文誌 2, 555 (2003).
［４５］S.C. Jardin et al., J. Comput. Phys. 227, 8769 (2008).
［４６］G. V. Pereverzev and G. Corrigan, Comput. Phys. Commun. 179, 579 (2008).
［４７］J.H. Holland, Adaptation in Natural and Artificial Systems:!!!!!!!!!!!!!!!!!!!!!!!!!!

４１７