csimyu
おドン
csimyu
準備中
準備中
準備中

基本データ型

char文字か整数1文字-128から127
unsigned char符号なし0から255
int整数型-2147483648から2147483647, 20億
unsigned int整数型0から4294967295, 40億
double実数型1.7E-308から1.7E+308
long double実数型,修飾子1.7E-308から1.7E+308

型変換

異なるデータ型同士の数値の演算は最もサイズの大きな型を用いて計算が行われる.

多倍長整数(基本データ型では扱えないデータ処理)

典型的なC言語の処理系では長い整数であるlong型でも32bit. unsigned longでもlog10 2^32=9.6より10進数で概ね9桁以下の整数しか表現できない.
例)3^99を計算する方法
(1)unsigned long型の配列を使って多倍長整数を表現する. bignum[k]を用いて配列各要素に7桁の整数を割り当てる. 7(k+1)桁.
(2)bignum[k]を下位から3倍していく. 下位からの繰り上げを忘れずに計算.
3no99jou.cxx

誤差

桁落ち: ほぼ等しい引き算で生じる.
丸め: 実数を2進数で表現すると生じる.
情報落ち: 大きいものと小さいものを足すか引くかで生じる.
(1)桁落ち
(x+1)^0.5-x^0.5=1/((x+1)^0.5+x^0.5)
有理化して引き算をなくす
(2)丸め誤差
0.1を1000000回足す. 0.1は2進数で循環小数となるので100000.000001となってしまう.
このようなアルゴリズムは採用すべきでない. (3)情報落ち
絶対値の小さいもの同士から順に足していくと防げる.

2次方程式の解

桁落ち誤差防止のため有理化した解の公式を使って計算している.
nijihouteisiki_no_kai.cxx

常微分方程式

常微分方程式をある初期値のもとで数値的に徳とは初期値から初めてある刻み間隔で次の値を順々に求めていく作業のことを言う.
一回微分方程式を解く方法としてオイラー法(教育用にしか用いられない素朴な方法), ルンゲクッタ法などがある.

質点の落下(オイラー法)

(1)刻み幅, 速度の初期値v_0, 位置の初期値x_0をあらかじめ決めておく.
(2)オイラー法により次のステップの速度を求める.
  v=v_0+gh
(3)オイラー法により次のステップの位置xを求める.
  x=v+gh
(4) (2),(3)を繰り返して順に求めていく.
freefall.cxx

着陸船のシミュレーション

落下し始めてから数秒後に逆噴射をする. freefall_hunsya.cxx

ルンゲクッタの法

   準備中です.

ポテンシャルに基づく2次元運動のシミュレーション

電荷によって与えられたポテンシャルの中を運動する荷電粒子のシミュレーション.
オイラー方程式で解いた. denka_nijigen_undou.cxx


偏微分方程式

ラプラス方程式

△u(x, y)=0
例えば、一変数函数f"(x)=0の時f(x)は直線になる.
これと同じで, 二変数関数の場合は表面に膨らみや凹みのない枠に合わせてゴムの膜を張ったような形状となる.
ラプラス方程式の境界値問題とは関数を表す平面の端の値を決め, その内部をラプラス方程式を満足するように計算すること.
数値計算を行うためには, 問題を離散化しなければならない.
u(x, y)の定義される領域Dを格子状に離散化し, ラプラス方程式を差分方程式に直して, 各格子点の値を数値計算する.
u(x, y)の値を各格子点上でu_ijとして求める.ラプラス方程式を満足するように格子点u_ijの値を決めるには、隣接する上下左右の格子点の値に対して,u_ijの値が隣接点の値から離れることなく, 滑らかでなければならない.
このためにはu_ijの値が隣接する4点の平均値である必要がある.
  u_11=(0+0+u_21+u12)/4
  u_21=(0+u_11+u_31+u_22)/4
  ・
  ・
  ・
  u_33=(u_32+u_23+0.75+0.75)/4
このようなu_11からu_33を決定する9元連立方程式となる.

  

連立方程式を数値計算で解くガウスの消去法(吐き出し方)

gauss_no_syoukyohou.cxx

逐次近似

(1)u_ijの初期値を設定する.
(2)以下を適当な終了条件まで繰り返す.
  1. u_ijの周囲の4点の平均を計算し, u_ij_nextとする.
  2. u_ij_nextをu_ijにコピーする.
  3. これがu_ij_next=u_ijぐらいに収束するまで繰り返す.
laplace.cxx
laplace_syokiti.cxx
laplace_syokiti.txt

使い方: ./laplace < laplace_syokiti.txt
またはofsをcoutに書き換えて ./laplace_syokiti | ./laplace

領域Dが長方形以外の場合

図を切って、長方形の組み合わせとして考える.

セルオートマトンを使ったシミュレーション

セルオートマトンとは内部状態を保ったセルが他のセルとの相互作用により時間的に変化していくというモデル.
各セルは隣接するセルと情報を交換する.
セルオートマトンの世界には時間が存在する.セルの状態は時刻とともに変化する.
通常, セルオートマトンの世界に含まれるすべてのセルはあるタイミングで一斉に状態を更新する.
つまり、セルオートマトンの世界では, 時刻tは, t=0, 1, 2, ...と離散的に変化する.

一次元のセルオートマトン

セルの状態は0と1の2通りとする.
隣接する3つの状態は1時刻前のt_(k-1)の状態で決定されるが, 決定のためには, ルールが必要.
3つのセル全てが0の状態から, すべて1の状態まで2^3=8通りある.
111, 110, 101, 100, 011, 010, 001, 000
この8通りのそれぞれの場合に対して2通りの選択が可能なので, ルールの総数は2^8=256通りとなる.
] ルールの名称は, 並びを8桁の2進数として読んだ時の10進数ひょうげんで与える.
(00010010)_2= 18: ルール18となる.
cell_automaton.cxx
cell_automaton_65syokiti.cxx
cell.txt
使い方: ./hello ルール番号 < 初期状態ファイル.txt

ライフゲーム

ルール: (1)時刻_kにおいて, セルc_ijの周囲8セルの状態の総和S_ji_(t_k)が3ならば, 次の時刻t_(k+1)におけるセルc_ijの状態a_ij_(t_(k+1))は1.
(2)S_ji_(t_k)が2ならば, a_ij_(t_(k+1))は変化なし.
(3)上記以外の場合はa_ij_(t_(k+1))は0.

・セルの状態が1であることを, 生物が存在することに例えるならば, (1)と(2)は生物の誕生あるいは存続を意味する.
・セルの状態が2であれば、周囲の環境が適切である場合のみ生物が存続することをシミュレートしているとみなすことができる.
life_game_syokiti.cxx
life_game_syokiti.txt
life_game.cxx
使い方: ./life_game < life_game_syokiti.txt

擬似乱数

乱数と擬似乱数

擬似乱数は計算によって求めた一見ランダムに見える数の並び.
計算アルゴリズムが分かれば, 擬似乱数の値を予想することは容易.
しかし, シミュレーションの目的を達成できる範囲で数がランダムに並んでいれば, 擬似乱数を乱数(random number)として扱うことは可能である
線形合同法はC言語のライブラリ関数であるrand()乱数生成関数の多くの実装において生成の基礎となるアルゴリズムとして広く使われている.
線形合同法のアルゴリズムは非常に簡単. 乱数系列R_1, R_2, R_3, ...に置いて
  R_(i+1)=(aR_i+c)%m
  a, c, mは正の整数 によって, 次の値を順次計算する.
通常mは, R_iのビット幅に依存して決定される. 例えばR_iが32bitのunsigned int 型であれば, mを2^32とすることができる.
mを小さくすると乱数の周期が短くなり, 同じ並びが現れやすくなる.mを大きくとる方が有利.
aとcは乱数の性質に大きな影響を与える. 適切に選ぶと周期が長くランダムな乱数系列を得ることができる.
逆に不適切な選択をすると, 乱数としての性質が失われる. 「William H. Press 『ニューメリカルレシピ・イン・シー』技術評論社 1993」ではmを2^32とする場合について
  a=1664525
  c=1013904223
を適切な値と例示している.実際のrand関数の実装では, 線型結合法をそのまま用いるのではなく, 線型結合法の出力に適切な変換を加えている場合が多い.
乱数生成のアルゴリズムについては, Webで「M系列」, 「Mersenne Twister」を検索するとよい.
senkeigoudouhou.cxx
使い方: ./senkeigoudouhou (初期値)

積分

台形公式

数値積分の基本的な考え方を与える公式のひとつである.
台形公式では, 関数f(x)のある区間[x_i, x_(i+1)]を直線で近似する.すると区間[x_i, x_(i+1)]の積分値は台形部分の面積となる.
f(x)上の点x_0, x_1, x_2, ...,x_nが等間隔hで並んでいるとすると,

senkeigoudouhou.cxx

乱数を用いた数値積分

台形公式は単純な公式だが, f(x)が滑らかな関数であれば, 分割数を増やすことでそれなりの制度を持った計算ができる.
これに対し乱数を用いた数値積分は, 計算量が大きいにも関わらず, 制度は良くない. 乱数による数値積分はf(x)の性質が悪く他の計算が困難な場合に用いる特殊な方法である.
例)正方形の中にランダムに点をプロットしていく.打った点の総数に対する4分円の面積の近似値となる. 10個中8点が4分円の内部に含まれているとすると4分円の面積の近似値Iは
  I=8/10=0.8である.
C言語のライブラリ関数であるrand()関数を用いて乱数生成.strand()関数を用いて乱数の初期値を設定する.
記号定数RAND_MAXはrand()関数の返す値の最大値であり, stdlib.hで定義されている.
ransuu_sekibun.cxx

乱数と最適化

組み合わせ最適化問題の典型例としてナップサック問題を取り上げる.
制限重量のあるナップサックに価値と制限重量が決まっている複数の荷物を詰め込む問題.
ナップサックの制限重量を超えることなく, なるべく荷物の価値の合計値が大きくなるように荷物を詰めなければならない.
価値が最大となる荷物の組み合わせを最適解, 厳密解と呼ぶ.
ナップサック問題の最適解の解析的に求める方法は知られていないので基本的には, 最適解を求めるためには, 全ての組み合わせを調べなければならない.
荷物の個数が一つ増えるごとに組み合わせの数があ2倍になるので, 力ずくの総当たり法では, 荷物の数が数十個程度までしか求めることができない.
そこで探索の方法を工夫する手段として, 一般に動的計画法や分岐限定法などが用いられる.

乱数を用いて比較的優良な解を求める(最適)解を求める方法

乱数を使ってナップサックにランダムにつめ, その結果を評価する. これを繰り返すことで, 重量を越えないで価値が最大となる荷物の組み合わせを探す.
napsakku.cxx
napsakku.cxx
出力結果napsakku.cxx
使い方 ./napsakkku (荷物の個数) (制限重量)(試行回数の上限) < napsakku_data.txt

ランダムウォーク

酔歩と呼ばれる. 物理現象のシミュレションだけでなく経済学における経済現象のモデル化にも適応される.
2次元ランダムウォークでは単位時間ごとに乱数の値を座標値に加算する.
napsakku.cxx
使い方 ./randam_walk (試行回数n) 7ffffffe

エージェントシムレーション

エージェントは内部状態を持っていて, 外部とやり取りをすることのできるプログラムである.
複数のエージェントを同一環境内で動作させるマルチエージェントの枠組みを用いると大勢の人間の挙動や動物の群れの様相など, 物理的シミュレーションだけでは, 再現の難しい社会現象や集団挙動などをシミュレーションすることが可能である.

感染のシミュレーション

gnuplotを用いてグラフィカルに描写した.
infection.cxx

正方格子の古典イジングモデルのモンテカルロ法

マルコフ過程

確率変数が1ステップ前の状態にのみ依存し, 他の時刻間の相関を持たない, つまり隣り合ったステップ間の状態によってのみ決定されるような確率過程.
マスター方程式はマルコフ過程を前提として導出される.

マスター方程式

系列が時刻tにおいて離散的な状態iにある確率P(i,j) (i=1,2,3,...)
規格化

とする.
状態iが状態jへの単位時間当たりに遷移する確率w_(i→j)とするとマスター方程式は以下のようにかける.


マスター方程式の行列表示

但し, 確率行列Lは


性質

より

モンテカルロ法


Nが大きいので数値計算は無理である.そこでMこの状態を確率分布p_iにしたがって発生させ, 和をとる.

定常マルコフ過程

定常分布へ収束させる過程
Copyright (C) 2017 Odon Laboratry  All Rights Reserved.