計算物理 2018
物理学を計算させて理解する
C言語を使って計算させる
Cygwin版
竹下徹

1

計算機内部の表現 ;8bit
計算機なかでは数字も文字もbitの集合
0

1

２進法の一桁を表す：ビット bit

２進法の8桁を表す：バイト byte

２の8乗＝２５６通りの組み合わせ＝２５６種の文字
ASCII: アスキー ： American Standard Code forInformation Interchange0011 0000 : "0"0011 0001 : "1"0011 1001 : "9".....
0100 0001 : "A"

0110 0001 : "a"0100 0010 : "B"

0110 0010 : "b"...
0101 1010 : "Z"

0111 0010 : "z"2

計算機内部の表現 ;UNICODE
日本語：１６ビット：２バイト；全角
２の１６乗＝65536通りの組み合わせ：ASCIIを含む
ISO : International Organization forStandardization . 国際標準化機構. ISO-10646
ASCII

0000 0000 0011 0000 : "0"0000 0000 0011 0001 : "1"0000 0000 0011 1001 : "9".....
0000 0000 0100 0001 : "A"
0000 0000 0110 0001 : "a"0110 0001 0001 1011 : "愛" 611B

3

計算機内部の表現 ;数字 :IEEE754
３２ビット：４バイト
２の３２乗通りの組み合わせ整数
0000 0000 0000 0000 0000 0000 0000 0000 : 00000 0000 0000 0000 0000 0000 0000 0001 : 10000 0000 0000 0000 0000 0000 0000 1001 : 9.....
1111 1111 1111 1111 1111 1111 1111 1111 : -1
実数
s

e

a

1bit

8bit

23bit

実数 X

(10)

s
e -127=(-1) [0.a](2) x 2 (10)実数 6.0=0 10000001 10000000000000000000000(10)

実数 -0.375(10)=1 01111101 10000000000000000000000
4Cygwin• UNIX on windows• オープンソース(GNU GPL)
でPCと対話する• command命令• 大文字小文字は区別する• optionは - の後につける
windows:コマンドプロンプト
5unix command•ls ：listの略
folder内の今いるdirectoryあるいはファイルやdirectoryを表示する!

•ls -l ：表示の仕方が細かい!
•pwd ：今いるdirectory を表示 / で上下関係を区切る

左が上側、右が下位のdirectory !

•cd XXX ： change directoryの略
directory 名

directory

XXXへ移動!

•cd ..

：一つ上のdirectoryへ上がる
6/home/user/XXX

C言語• 自由度の高い低級言語：関数からなる• 古い(1972)- 信頼性高い＜＞面倒• ユーザー多数• 汎用性高い：パソコンソフト（OSを含む）、ゲームソフト、機械制御、自動車、家電、システム管理…

C言語の書き方.hはヘッダーを表す• #include <stdio.h> で始まる
• main() {..文...} が全体を囲う• 文は変数や関数から成る、main()も関数• 変数は、整数intか、実数ﬂoatか、文字charか• 計算は、右辺で行い、結果を左辺に入れる
stdio.hはprintf()などの関数を含む

文は；で終わる：フリーフォーマット

!
a*b=x-c/d; は意味なし
a*bという変数が出来ている

x=a*b+c/d;
8

作業手順• エディターでテキストファイルを編集保存
xxx.c という名保存はUSBかcloud• コンパイルは gcc xxx.c• 結果の出力は a.exe• 実行は ./a
• 出力を画面でなく、ファイルへ行うときはテキストファイルを計算機が理解できる形式に変換する>を使う./a > yyy

yyyはファイル名
yyyの中身は9

more yyyで見える

抜けるには CNTL C

C事初め• 作法がある• プログラムは上から下へ一行 #include <stdio.h>
スペース

づつ実行

•

int
main(){#include <stdio.h> 最低限
printf(“文字”);
main { で始まり、}で終わる}
画面への出力は printf(“ここに文字”);••• 各行は; セミコロンで終わる
10

C事初め

• 計算
• 変数は整数型 int か、実数型

•

•
•
•

#include <stdio.h>
int main() {スペース
ﬂoat か定義する
int i,j,k; 整数型変数名i,j,k
画面への出力は printf(“%?” ﬂoat a,b,c;
実数型変数名a,b,c変数名);? = d整数、f実数
i=1;％d 整数型
a=0.1;％f 実数型
printf(“ %d %f\n“,i,a);print文の\nは改行
}
11

変数

•
• 実数 ﬂoat y;• 計算文 y=x+1;
整数

int x;

スペース

演算

記号

例

和

+

z=x+y;

差-

z=x-y;

積*

z=x*y;

商/

z=x/y;

あまり%

z=x%y;
12xという名の整
数変数に整数１を代入するx=1;xという名の整数変数に入った数値とさら
に１を足して、
結果をyに代入する
y＝x＋1;

関数文<stdio.h>• printf() : 端末画面に印字する<math.h>sin(x):
サイン
x•<math.h>• pow(x,n): power : x ベキ乘
<math.h>sqrt(x):square
root
:
xの平方根
•
<stdio.h>• scanf(): キーボードから入力
n

13

Cの計算
• 左辺Aへ右辺Bの計算結果を代入 A=B
計算
加算
減算
乗算
除算
あまり
１加算
１減算
和の代入
差の代入
積の代入
商の代入
あまりの代入

記号
+
*
/
%
＋＋
ーー
+=
-=
*=
/=
%=

例
c=a+b
c=a-b
c=a*b
c=a/b
c=a%b
k++
k-k+=5
k-=5
a*=5
a/=5
a%=5
14

c=a-(a/b)*b [int a,b,c]
k=k+1
k=k-1
k=k+5
k=k-5
a=a*5
a=a/5
a=a%5

練習問題
• int 型の変数 x に数値 11 を代入し、xの値を “x=11”のように表示するプログ
ラムを作成しなさい。

• 1割る3 を計算させよ、変数は整数や
実数で試せ

• 宿題は

tohrutt@shinshu-u.ac.jpへ送る
15くりかえす• for ( k=1; k<3;k++) { printf();}• printf()をkが１から３を超えない範囲• k++ は k=k+1
1• printf(“%d /n”,k)なら216• 繰り返しfor文
整数 k ,実数 x,y•sin(x)
は数学関数

結果
k=1
x=1/10.0=0.1#include <math.h> が必要
y=sin(0.1)
k=2
for ( k=0; k<10;k++)x=2/10.0=0.2{
y=sin(0.2)
x=k/10.0;
k=3
y=sin(x);printf(“ %d %f %f \n“, k,x,y); x=3/10.0=0.3
y=sin(0.3)}./a > aplot.dat> redirection

....

17
出力ファイル名aplot.datはどんな名前でもよい、しかし.cだけは使わないgnuplot• グラフツール• plot sin(x) でサイン関数を赤線でプロット• plot “yyy” でデータファイル名yyyから数値を読み込みプロット• plot [xmin:xmax][ymin:ymax] “yyy” :グラフの最大値と最小値の設定• plot “yyy” using 1:3 で１列目vs３列目でプロット
18

カオス シミュレーションの• 実数変数xはある集団の個体数• 実数変数aは子供の増え方を表す増倍率であるととても奇妙な振る舞い• この系を（１で規格化されている）

for文で計算させるをする

x aが大きい場合飢餓と繁殖を繰り返す

for ( k=1; k<10;k++){x=a*x*(1.0-x);
printf(“ %d %f \n“, k,x);
}
19

横軸世代

k

世代と個体数x• x=a*x*(1-x)の系は面白い• a が１より小さいと滅亡 x->0
• aが１〜？の間では一定値に収束する• aが？〜??付近以上になるとカオスである（前ページ）• さらに大きくなるとxが１を超えるxa20

for 文その２• 整数 k,s
int k,s;• 和をとる:初期化s=0 s=0;for(
k=1;k<5;k++)• s=s+k はs+=k{• と書いてもよい
s=s+k;
printf(“ %d %d \n“,k,s) ;
}

21

和の計算
N

s
1
int k,N=??;
s=/ k
k=1
ﬂoat x,s;
ここでsは実数s=0.0;for (k=1;k<N+1;k++){x=1.0/k;
注意：1/kは整数なので０s=s+x;
printf(“%d %f %f ¥n“,k,x,s);}N
N-> ∞ で発散(log)
2

1
1
1
π
次を示せ； 1+ 2 + 2 + 2 ... == ς (2)234622

無限級数和∞1
Z(s) = ∑ s
k
k=1
問2:

問1: s=2,4を計算で示せ∞

1€
Z(s) = ∑ s
k
k=1

問3:

" 1% 2"1% 4
Z(2) = $ 'π ,Z(4) = $ 'π
# 6&
# 90 &

が収束するsの条件を見つけよ

1 1 1
π
1− + − + .... =3 5 7
4を計算で示せ

2
3
問4:
x x
x
x
+
+ ....
マクローリン 展開を示せ e = 1+ +1! 2! 3!
23

複数のxについて等号を計算して確かめる

階乗計算

long k,s;s=1;
for( k=1;k<19,k++){s=s*k;
printf(“ %d %d\n “,k,s) ;}• n! 計算• 破綻• 倍精度計算nを実数にしてみよ• 整数 long 印刷 %ld
やはりおかしいはず• 実数 double 印刷 %lf• それでも計算がおかしいのはどこからか？
• おまけ C =n!/r!/(n-r)!を計算せよ
s=s*k
で前回までの積sにkを掛けている
数学的帰納法にように、
k=1,2,3と進むと
1!, 2!=1*2,3!=2!*3 ,,,

n

r

24

2

12

２重for文xxxe = 1 + 1! + 2! + ... + 12!
!

logy

for(j=0;j<50;j++) for-1
xの定義{ x=j/10.0;
このままでは正しく動きません！
s=1;
waの初期化をしてください。
f=1.0;
for-2
for( k=1;k<13;k++){ s=s*k;f=f*x;
wa=wa+f/s;}

k!
xk
12
k
x
/ k!
k=1
set logscale yxprintf(“ %f %f\n “,x,wa) ;}25

分布：rの関数• 二項分布 確率=p (0<p<1の実数)
p =pow(p,r)• P(n,p;r) = C p (1-p)• C =n!/r!/(n-r)!を計算せよ• n=20 & p=0.7 で分布のグラフを描け• Poisson分布P(r)=λ exp(-r)/r!をプロット
n

n

r

r

(n-r)

r

r

せよ

r

２重for文の練習問題• xが0,1,2……のとき、下記のy=f(x)をマクローリン展開した式で計算し、plot
せよ,関数も重ね書きせよ。

•• (2) y=sin(x)=自分で書き下せ• (3) y=log(1+x)=自分で書き下せ(1) y=1/(1-x)=1+x+x2+x3+…

• カオス問題をaを変化させながら、p18のplotを作れ
27

計算失敗

x=10.0• 扱えない大きさに注意• ベキ乗計算:実数• pow(x,y) =x 計算関数
• 階乗計算• 整数はint (32bit) では11桁(10進)• 整数はlong(64bit) では19桁(10進)2210000000000000000000000.000000
23 99999999999999991611392.000000
24 999999999999999983222784.000000

308
10000000000000000109790636294404554174049230967731184633681068290315758540491149153716332897849468889906124966972117251561159028374314008832
83070091981460460312716645029330271856974896995885590433383844661650011784268976262129451776280911957867074581227839701717844151052918028932
07873272974885715430223118336.000000
309 infy-2,147,483,648 〜 2,147,483,647‒ 9,223,372,036,854,775,808 〜 9,223,372,036,854,775,807

フーリエ級数
任意の周期関数は三角関数の和で表される

f(x)

∞

4 sin{(2k −1)x}f (x) = ∑
π k=1
2k −1
矩形関数に近づくxplot [-pi:pi][] 4/pi*(sin(x)/1),4/pi*(sin(x)+sin(3*x)/3),4/pi*(sin(x)/1+sin(3*x)/3+sin(5*x)/
5)∞

問題： f (x) = 1+

∑ 2cos(kx)を計算して図示せよ

k= 0
29

•

プログラム例∞f (x) = 1+ ∑ 2cos(kx)k= 0

どんな関数でしょう？
kmax=25

kmax=5030

for (j=1;j<jm;j++){x=-5.0 + j /10.0;
for(k=0;k<km;k++){f=f+2*cos(k*x);}}
kmax=200フーリエ級数の問題
次の関数(1)と(4)をフーリエ変換して、図を描け
大きさ

1)
2)
3)
4)

右図1を作る関数を見つけよnについて和を取って図示せよ 0 ω
大きさ
右図２を作る関数を見つけよ
a

図１
x
4ω
図２

y~-x-20 1 2 ３周波数

キーボードから入力• 数値xを読みこむ• scanf()が入力関数<> printf()に対応
ﬂoat x;printf(" input a real number : ");scanf("%f",&x);実数を読む &:アンパサンドとxcomment out 先頭に // :行だけ除外•• あるいは /* …実行文（複数行OK)… */ で囲う :実行文から除外される
32

差分：数値微分と書かれる•• limをとって有限のhで書いたものを差分と微分はf (x + h) - f (x) df (x)
lim
=
h"0
h
dx

h=0.01;for(k=0;k<n;k++){hが十分小さいと差分〜微分
x=a+h*k;小さい実数hで計算する
df=((x+h)-(x))/h;
printf(“%f%f”);x,df
右は f(x)＝x, dfが差分
呼ぶ

•
•
•

}
33

main(){
ﬂoat f(),x;....f(x).......
f(); 名前は自由}
main 文の外に関数文を置く ﬂoat f(x)関数文,名前はfこともできるﬂoat x;複雑な関数なら間違いがない {ﬂoat y;
変数のつながりに注意
y=sin(x)*exp(-x);
多変数関数：f(x,y)
return(y);

関数文

•
•
•
•
•

定義された変数は{}内で有効、外では無効
34}

微分/差分の問題y• 次を数値的に示せ• dsin(x)/dx=cos(x)• de^x/dx=e^x :^はベキ乗
x
x+h x
x-h
• dlog(x)/dx=1/x• darctan(x)/dx = 1/(x*x+1)• df(x,y)/dx: f(x,y)=x*y/(x*x+y*y)• 中心差分：fʼ(x)={f(x+h)-f(x-h)}/2h =[{f(x)f(x-h)}+{f(x+h)-f(x)}]/2h 3点の差分の平均

（数値）積分• 積分＝面積計算• 長方形近似by=f(x)

S = ∫ f (x)dxadx=(b-a)/n;一定for(k=0;k<n;k++){x=a+dx*k;s=s+f(x)*dx;}

s

f(x)

s s s
s s

s

x
b

dx

a
台形近似

y=sqrt(a*a-x*x)が半円の面積 s=s+(f(x)+f(x+dx))*dx/2;
y=-x*xの面積
36

課題• Plank分布B(m)：黒体からの光（波長m)の分布
2hc
1
2

!

B( λ) =

5λehc / λkT−1

Tは温度• B (m)を波長の全域で積分すると温度Tのエネルギーを計算できる：次の二つを示せ€この値は温度Tの４乗に比例する：シュテファン＝•ボルツマンの法則•# B (m) dm

積分せよ

B(m)の頂点のλはT-1に比例する：ヴィーンの変位則
dB (m) 微分せよ
h=c=k=1の単位で計算せよ
dm37

Plank 分布の積分と微分
case : h=k=c=1
Stefan-Boltzmann則
gnuplot: set logscale y
積分結果

Wien変位則
微分結果
T-1線

T4線

T

38

T

debug：デバグ• bugバグ＝虫：プログラム中の間違いだがcompileで来てしまう• de は 「反対」の意味de-bugでbugをなくす事を言います• 最も単純なdebug法は• 実行文中に printf(" %f %f %f ...",s,x,y..)；を入れて使っている変数の値を打ち出
してください。
なにを打ち出して、なにになるべきかを考えて眺めてください。!
なにが間違いがわかります39

よく起きるbug• float a,b,c,x,y; a=b*c/x*y;//これは• yの位置に注意• float a,b,c,x,y; a=b*c/x/y;• a=b*c/(x*y);

//これは

と同じ

• int a=1,x=3; ﬂoat y; y=a/x; y=0.0です• 自分で試して確認してください。

40

b#c#y
a=
x
b#c
a = x#y•if
文
もし なんとか なら 式
１を実行、そうでないなら
式２を実行• if (なんとか) {• 式1;}• else {<<< そうでないなら• 式2;} 省略可能
• z=|x-y| を計算

scanf(“%f”,&x);
scanf(“%f”,&y);if(x>y) {z=x-y;}else{
z=y-x;}
printf(“%f\n”,z);
差の絶対値を取る

if(w>x>y) は使えない if(w>x&& x<y) を使う
41

if (式)a>bif(a>b)a>=bif(a>=b)a<=bif(a<=b)整数

a=bif(a==b)整数

not equal

if(a!=b)(式1) AND ( 式2)if( (式1) &&(式2) )(式1) OR ( 式2)if( (式1) || (式2) )要注意：実数同士が＝＝になることはまずない
42

if文の課題• [1] if(a=b)はどういう動きをするか述べよif(a=b)は数値式a＝bは、aという変数にbの中身を代入! if(数値) 数値が０であれば偽として扱われ、０以外であれが真として動く• 本当はif(a==b)と書きたかった• いろんなa,bの組み合わせで試せ• [2] 整数を入力して、偶数か奇数か答えるプロ
グラムを作れ：できたら授業中にメールしてくださいif文の使用例

極値探し• 関数を微分• 有限hで近似• 極値を探す場合
for(k=0;k<?; k++) {• 微分係数gが正の間印刷 x=a+h*k;g=(f(x+h)-f(x))/h;
印刷が止まると極値に近い•
if(g>0)• f=x -2x-1 の極値を計算 { printf(“%f %f %f3

微分係数の符号が正の間中出力する原始的なプログラム：止まった
点で符号が変化
賢い方法:breakの利用\n”,x,f(x),g);}44}

if() break;• break 文• forのループから抜ける時• 右のプログラムを動かせ• 課題 関数が極値を持つ点を見つけるプログラム
を作り、３次関数でそれを示せ

int x,y=?,ix=??;for(i=0;i<ix;i++){x=i;if(x>y) break;}printf(“%f\n”,x)
条件が満たされたらfor文の外に飛び出す：break課題• Plank分布B(m)：黒体からの光（波長m)の分布
2hc
1
2

!

•

B( λ) =5λh=c=k=1の単位で計算せよehc / λkT−1

Tは温度B(m)の頂点のλはT-1に比例する事を示せ：
ヴィーンの変位則
Wien変位則dB (m)€ ２重for文
dm
温度Tのループと波長λのループが必要 微分して微分係数の符号が変化したら
if文でbreakする
その点を印刷せよ

T-1

T
46

微分方程式の解法(1)
右辺はtの関数(2) 左辺はt+Δtの関数

時間をdtづつ動かす

右辺はtでの値
(3)

(4)
(5)

自由落下

dv
m
= mv "(t) = mg
dt

dx
=vʼ
v=
dt

v(t + Δt) = v(t) + Δt × g
fʼ=vʼ=g€(6) x(t + Δt) = x(t) + Δt × v
fʼ=xʼ=v

47

for(k=1;k<km;k++){v=v+dt*g; 加速度=gx=x+dt*v;
t=t+dt;}空気抵抗• 速さに比例する場合
dv! m= mv "(t) = mg − α × v(t)
dt

m=1.0

for(k=0;k<?; k++){t=t+dt;v=v+dt*(g-a*v);
x=x+dt*v;}• 宿題１：終端速度はmg/a• 速さの２乗に比例する場合• 符号に注意 v**2>0• 投げ上げ問題 でif文で分ける• 宿題２：速度の２乗の抵抗で投げ問題で次の
式を示せ

1
1
1
= 2+ 22v g v i v48f

vg:地上に戻った速さ
vi:初速
vf:終端速度

v

空気抵抗
v

• 速度に比例：終端速度~mg/a• 速度の２乗に比例：終端速度~sqrt(mg/a)m=1.0 v
a=1.0

t

49

t

単振動と抵抗• 単振動の微分方程式
dv
m
!

dt

for(k=0;k<?; k++){t=t+dt;v=v+dt*(-b*x);x=x+dt*v;}= mv "(t) = −b × x(t)• 単振動を示せ• 空気抵抗がある場合

dv
m!
= mv "(t) = −b × x(t) − α × v(t) for(k=0;k<?; k++){dt•

減衰振動を示せ

t=t+dt;v=v+dt*(-b*x-a*v);

a=0.1

t

振幅
50

x=x+dt*v;}

振り子• 振れ角θが小さいときは• 振れ角θが小さいときは• これは簡単には手で解けない• 計算機はできる• 動かしてみてください• 楕円関数ですか？
• sinθ〜θでなくてもOK?
• 初速度を大きくすると回るか？

θ

x

ℓ

x
mgsinθ
mg

１次元問題•• x〜０での挙動はおかしい• 正しくない：差分の限界？
力F=-k/x/x として

v,x

x
v

!

ｔ

ｖ

x

t
x
x<0,v<0x<0t

平面上の運動• 惑星の楕円運動• 連立微分方程式 を２次元で作る
r=sqrt(x*x+y*y);• g=1,m=1のとき

安定軌道の初期値x=1, y=0, vx=0, vy =の最小値は？

y

楕円軌道

a=g*m/r/r;
c=x/r;
s=y/r;
vx=vx-a*c*dt;
vy=vy-a*s*dt;

実はyが長いx=x+vx*dt;
初期値x=1, y=0, vx=0, vy =?

x53

y=y+vy*dt;

平面上の運動• 課題１ 前ページのように万有引力での楕円運動を計算せよ。

• 課題２ 課題１で初期条件を変えて、円運動する条件を探せ、またこれより初速度が小さいとどうなるか？

• 進んだ課題（解ける人のみ）ポテンシャ
ルをバネや一定重力に変えるとどんな運動をするか？
54

乱数• 一様乱数:どの数値も同じ確率• 一様乱数発生関数 rand()• 乱数初期化 srand(int)• rand()の最大値をさがせ• 最大値で割ると0~１の実数一整数

様乱数ができる• 自然現象の裏にあるので• シミュレーションに使える
55

int r,rand();srand(整数);
for(k=0;k<?; k++){r=rand();printf(“ %d %d\n”,k,r);}rm=231-1=2147483647

さいころ• ランダムな１­６の整数を作る道具• int A%6 はAを６で割ったあまり• １の目の出る割合を調べよ• 6の目の出る割合を調べよ• どちらも同じ：一様56

int r,rand();srand(??);
for(k=0;k<?; k++){r=rand()%6;
printf(“ %d %d\n”,k,r);
}

面積計算

乱数の応用

•
• 円の内側の点の個数nを数える• n/M は rr =3.14159/4

0-1の一様乱数の組みをM作る
1
4

2

srand(int)
M個の乱数を生成
for(k=0;k<M; k++){x=rand()/rmax;
y=rand()/rmax;
r=sqrt(x*x+y*y);if (r<1.) n=n+1;}

正方形内の乱数M個

M個

1.0

n個
r<1.0の乱数の数

0.0.,x,y,r;ﬂoat rm=2147483647

1#1

y

0<x<1.0, 0<y<1.0

1.0 x

57

面積計算問題
モンテカルロ法とも言う

b-a
課題：楕円の面積を計算せよ
2
2

y
x

• x y
-b
+
=
1
の面積はrab• a2 b2• 一様乱数 z=(ﬂoat)rand()/rm は 0.0~1.0なので• -a~aまでの一様乱数 xを作るには a>0として• x=a*(2.0*z-1.0)とする
実数

rmは2147483647

58

a

面積計算
• 関数の面積＝積分
y=f(x)
c• 関数の積分領域よりs• 広い領域内に一様乱数• を生成するab• そのうち関数より小さい物の数を数えて• 割り合いを計算する• 精度は悪いがどんな形にも適用可f(x)微分できない関数も積分できる：例えば
59

f(x)=|x|x乱数を使う面積計算• 問題：原点を中心として半径２の円のxの区間 a<= x <=bでの面積Sをモンテyカルロ法で計算せよ。
2

S
-2

a

0

-2

b

2

x

よっぱらいの歩み srand(int)• 乱数の応用：１次元xよっぱらい for(k=0;k<?; k• 実数 乱数r 0<= r <=1++)• r<0.5 なら右へ:x=x+1 ランダム r=rand()/rm;if(r<0.5)
r>=0.5なら左へ:x=x-1•
x=x+1;• 複数のよっぱらいのある歩数で elseの到着地点の平均値と標準偏差
位置x
左へひっくり返りそうなよっぱらい

x=x-1;
printf(“ %d %d\n”,k,x);}歩数k61

2次元でのよっぱらいy• 座標を２つ用意 x,y• ２次元平面内で４通りの選択肢x• 乱数x 0<x<1• 0.<=x<0.25 なら右へ、0.25<=x<0.5なら左へ• 0.5<=x<0.75 なら上へ、0.75<=x<1.なら下へ• N歩で到着する点は中心からの距離R、Rの分布と平均値は？

62

２次元よっぱらい• 原点から出発、１歩を４通りからランダムに選択• １００歩あるいたら1000人の到着点(100歩後)分布
中心が多い

63

配列array a[3]• 3成分ベクトルaをa[3]とかく• 2つのベクトルの内積xを計算• 同様に外積も計算できる
ﬂoat a[3]={ 0.,0.,1.} もOK
for(i=0;i<3;i++) a[i]= 0. ;もOK
64

ﬂoat a[3],b[3];a[0]=0.;a[1]=0.;a[2]=1.;b[0]= ;
b[1]= ;
b[2]= ;
x=0.0;
for(i=0;i<3; i++){x=x+a[i]*b[i];}配列の応用• 配列内の総和を求めよ• 配列内の最小値を求めよ• ２つの配列a[3],b[3]の• 和のベクトルを計算せよ• 外積ベクトルを計算せよ• 1次元よっぱらいの到着点で-10~＋10のすべての点を埋めるまで何人必要か？(歩数が偶数の時は、偶数のみ、奇数なら奇数点のみのすべて）
65histogramヒストグラム

頻度(回数）• 頻度分布• 一様分布 x実数：0<=x<=1.0• 横軸は有限幅の「ビン」に分ける 0• 例えば幅＝0.1、0<=x<0.1,
0.1<=x<0.2,,,, 10個のビン

1.0

頻度(回数）• 縦軸は何回一様乱数xがビンに入っ
たか。
0.0
66

x

...

x
1.0

一様乱数のテスト• x=rand()/rm は一様乱数• その意味は、xの分布が０.〜１.で一様• h[10] でh[0]は 0.=<x<0.1となった回数
1000個の乱数
h[0]=116h[1]= 85h[2]= 91.......
h[8]=119h[9]=102100h[]with step
最後のビンは手で入れたx67

histogram• ヒストグラムの整数配列 h(10)• まず中身を０にする• 実乱数x, int j=x*10;• h[j]=h[j]+1 :１個加える
bin数10int h[10] ; としただけでは0になっていない

j番目のビンに

jは横方向の整数でビン番号h[j]の中身は何回１がたされたか、
縦軸を確率にしたかったら、kmで割れ

for(k=0;k<10; k++){printf(“ %d %d \n”,k,h[k]);}printf(“ %d %d \n”,k,h[k-1]);
最後のビンは手で入れた!

68

int h[10];for(k=0;k<10;k++){h[k]=0;}srand(?)
for(k=0;k<km;k++){x=rand()/rm;
int j=x*10;
h[j]=h[j]+1;}

中心極限定理• N個の一様乱数の平均値の分布(histgram)は gauss分布（分布によらないが。。。）• 分布の平均値は一様分布の平均値• 分布の標準偏差は• 一様分布の標準偏差の• 1/sqrt(N)になる• これらを示せ

gnuplot “ ” with step

0.0
69

<x>
0.5

1.0

一様乱数を使う例• ２つの一様乱数の和のヒストグラムはガウス関数となる• ２つの一様乱数の差のヒストグラム？
• ２つの一様乱数の積のヒストグラム？
２つの一様乱数の商のヒストグラム？

非一様分布関数の作り方!

• 非一様乱数の例
１0<f(x)<=1-2x:0<=x<0.5•• 0<f(x)<=2x-1: 0.5<=x<=1.0
x
０
0.5
１
• に従う乱数を作れ
１• (0,0)-(1,1)の正方形内に y• 一様乱数(x,y)を作り
x０採用
y<f(x)を満たすxのみ使う！
•
不採用
１• このxは分布関数を満たし、非一様乱数71

中心極限定理II•

分布関数f(x)<=4(x-0.5)2をに従う非一様乱数xを作れ• この分布に従う２つのX=(x1＋x2)の分布は１• ガウス分布である• 平均値<X>=<x1>+<x2>=1
０
0.5• 標準偏差=sqrt(s1*s1+s2*s2)=72

１xコンデンサーの電場• Laplace eq. ２次元 51x51のメッシュ(点）
0-50• 座標交点のポテンシャル V[int x][int y]２次元配列
境界条件：中に１点
V=V+•
V[51][51]• 境界条件：箱の外周はGRD• V[0][any]=0., V[any][0]=0. yV+• V[50][any]=0., V[any][50]=0.
50
.
!
!
.
.
2
1
0

01,2,...
73

20

x 30

50

２次元Laplace equation
式１
式２
式３

差分（偏微分）

式４

二回の差分
式５
式６

式７

V(x,y)=(周りの４つのVの和)/4=(周りの４つのVの平均)
74

x+δ
x
x-δ
y-δ y y+ δ

初期化：全空間のポテンシャルを０

for( i=0;i<51;i++) /* clear array */{ for (j=0;j<51;j++){ v[i][j] = 0.0;}}

境界条件：周囲も０!

ラプラス差分方程式:
実はある点の周りのポテンシャルの平均を取る事

2
3

4
1

y-1 y y+1

x+1
x
x-1

vd=v[x-1][y]; 1
vu=v[x+1][y]; 2
vr=v[x ][y+1]; 3
vl=v[x ][y-1]; 4
v[x][y]=(vl+vr+vu+vd)/4.0 ;すべてのx,yについて繰り返す
75

ただし、電位を決めた点はのぞくfor (n=0;n<100;n++){安定するまで繰り返す
for(ix=1;ix<50;ix++){
すべてのx,yについて
ix,iy=0~50あるしかし
for(iy=1;iy<50;iy++){ix,iy=1~49
if(ix==25 && iy==25){
だけ計算
v[ix][iy]=10.0 ;境界の電位は常に０
break; }if文から抜ける
中心点の電位を10V[25,25]
V[][] splot “ﬁle-name” with lines条件•

for(ix=0;ix<51;ix++) {for(iy=0;iy<51;iy++){printf(“ %d %d %f”,ix,iy,v[ix][iy]):
}}

ix
76

iy

1本棒の電位（ポテンシャル:V)
V+=10 : ix=20~30, iy=25
V(x,y)

50
.
!
!
.
.
2
1
0

V+

y
x

y

y

splot ファイルの中は、x y V(x,y)
77

01,2,...

20x 30

50