物理科学科 2年次 2007年 後期(水5-6) シミュレーション物理(2単位 選択)
担当者(所属)：高 橋 信 介
(理工学部): taka(_#_)stu.hirosaki-u.ac.jp (_#_)を @ に置換えて送信してください
Last modified 09:45, Feb.26, 2008
目次
おしらせ：レポート2の不完全な部分は急ぎ再提出すること。(1/22)
おしらせ：研究室にて期末試験を返却しています。(2/11)
おしらせ：レポート3のヒントを出しました。(2/11)
レポート提出状況
第01回 OCt 3, 2007 / (Last modified Oct.3,2007)
第02回 Oct.10, 2007 / (Last modified Oct.10,2007)
- 講義
- 運動方程式の解法
2階の微分方程式を，連立1階の微分方程式にする
時間間隔��tごとに計算をする。
- Euler法
��tの左端のvを使って右端のyを求める
��tの左端のaを使って右端のvを求める
- Euler-Cromer法
��tの右端のvを使って右端のyを求める
��tの左端のaを使って右端のvを求める
- 演習
- 演習そのまえに(file整理のためのUnixコマンド)
- ls sample*
sampleから始まる類似のファイル名をすべて表示する
- mkdir Prog
Progというディレクトリ(フォルダ)を作成する
- mv sample* Prog
前期の授業で作ったファイルをProgというディレクトリに移動する
- cd Prog
ディレクトリProgに移る
- ls
- pwd
現在のディレクトリをしめす
- cd
ホームディレクトリに戻る
- mkdir Simulation
Simulationというディレクトリを作成する
この授業で作るファイルはSimulationというディレクトリ内で作業する
- ls
- emacs ex01.f &
- fortranプログラムの編集・コンパイル・実行
暗黙の型宣言は使わないことにする
implicit noneを使う
- Euler法
ex01.f
program ex01
implicit none
real*8 t0, y0, v0, t, y, v, a, g, dt
integer i, imax
data t0, y0, v0, g, dt/0d0, 1d1, 0d0, 9.8d0, 1d0/
data imax/100/
t = t0
y = y0
v = v0
a = -g
i = 0
write(*,*) i, t, y, v, a
do i=1,imax
a = -g
y = y + v*dt
v = v * a*dt
t = t + dt
write(*,*) i, t, y, v, a
enddo
stop
end
- Euler-Cromer法
ex01-2.f
do i=1,imax
a = -g
v = v * a*dt
y = y + v*dt
t = t + dt
write(*,*) i,t,y,v,a
enddo
- gnuplotによる描画
予め
f77 ex01.f
./a.out > data01.dat
f77 ex01-2.f
./a.out > data01-2.dat
で描画データをファイルに作成する
plot "data01.dat" using 2:3, "data01-2.dat" using 2:3
- 予習のために・連絡
- 次回はオイラー・リチャードソン法を扱います
- サブルーチにする方法を復習してくること
- レポートなし
第03回 Oct 17, 2007 / (Last modified Oct.25,2007)
- 講義
運動方程式の解法・その２
- Euler Richardson法について
区間��tの中間点tmidにおける速度vmid、位置ymidを求めて、中間点での加速度amidを計算する
a{n} = F( y{n}, v{n}, t{n} )/m
vmid = v{n} + a{n}*��t/2
ymid = y{n} + v{n}*��t/2
amid = F( ymid, vmid, t+��t/2 )/m
v{n+1} = v{n} + amid*��t
y{n+1} = y{n} + vmid*��t
t{n+1} = t{n} + ��t
- 一次元・空気抵抗なしの場合
m*a = - m*g
- 一次元・空気抵抗あり（速度の１乗）の場合
m*a = - m*g - k*v
- 演習
- 一次元・空気抵抗なし
ex02.f
do i=1,imax
a = -g
vmid = v + a*dt/2
ymid = y + v*dt/2
amid = -g
v = v * amid*dt
y = y + vmid*dt
t = t + dt
write(*,*) i,t,y,v,a
enddo
- 一次元・空気抵抗あり（速度の１乗）
ex02-2.f
do i=1,imax
a = -g - k/m*v
vmid = v + a*dt/2
ymid = y + v*dt/2
amid = -g - k/m*vmid
v = v * amid*dt
y = y + vmid*dt
t = t + dt
write(*,*) i,t,y,v,a
enddo
- レポート課題 なし
- 連絡・注意
第04回 Oct 24, 2007 / (Last modified Oct.25,2007)
- 講義
- 速さの1乗に比例する空気抵抗を考慮した場合の運動方程式の追加説明
上向きをプラスの方向とすると、
上昇中、v>0、空気抵抗下向き、-k*v<0
降下中、v<0、空気抵抗上向き、-k*v>0
上昇・降下いずれも m*a = -m*g -k*v
- 解析解の説明。
具体的な解きかた。
およびその関数の形・特徴について
- 空気抵抗を-m*g*(v/vt)で表した場合について。vtは最終の速さ。
たとえば、vt = 3 m/sec
- 速さの2乗に比例する空気抵抗を考慮した場合の運動方程式の説明
空気抵抗を-m*g*(v/vt)**2で表現する。
上昇中 m*a = -m*g -m*g(v/vt)**2
降下中 m*a = -m*g +m*g(v/vt)**2
- 演習
- 1次元・速度1乗の空気抵抗あり（mg(v/vt)）
ex02-2.f
do i=1,imax
a = -g - g*(v/vt)
vmid = v + a*dt/2
ymid = y + v*dt/2
amid = -g - g*(vmid/vt)
v = v * amid*dt
y = y + vmid*dt
t = t + dt
write(*,*) i,t,y,v,a
enddo
- 1次元・速度2乗の空気抵抗あり（m*g*(v/vt)**2）
ex03.f
do i=1,imax
a = -g - g*abs(v)*v/vt**2
vmid = v + a*dt/2
ymid = y + v*dt/2
amid = -g - g*abs(vmid)*vmid/vt**2
v = v * amid*dt
y = y + vmid*dt
t = t + dt
write(*,*) i,t,y,v,a
enddo
abs(v)はvの絶対値を求める関数。
- コンパイル・実行・結果をファイルへ
f77 ex03.f
./a.out > data03.dat
- gnuplotのコマンドをファイルで編集・保存・loadして表示。
> emacs ex03.plt &
set grid
set xlabel 'time(sec)'
set ylabel 'y(m)'
set xrange[0:10]
set yrange[-10:50]
plot 'data03.dat' using 2:4, 'data02-2.dat' using 2:4
上書き保存をする。
> gnuplot
gnuplot> load "ex03.plt"
gnuplot> exit
- レポート課題あり
レポート番号：rep01
プログラムのファイル名：ex04.f
速さの２乗に比例する空気抵抗を考慮して、２次元の軌跡を求めよ。
下記の砲丸投げを想定した場合の、投げる角度と着地距離について論じよ
運動方程式
ｘ方向： m*ax = - Fd*cosθ
ｙ方向： m*ay = -m*g - Fd*sinθ ただし Fd = k2*v**2
加速度は
ｘ方向： ax = -k2/m * v * vx
ｙ方向： ay = -g - k2/m * v * vy
砲丸投げを想定した初期条件：
m=7(kg), k2=0.01, y0=h=2(m), v0=30(m/s), θ0=30°, 45°, 60°
��t=0.01
注意： ｘ方向とｙ方向の運動方程式は独立では
ありません。速度ベクトルの逆向きに空気抵抗（ブレーキ）が働きます。その
向きは時間とともに変化していきます。これに対し、空気抵抗が無い場合は、
ｘ方向とｙ方向の運動方程式は独立になります。このことを良く検討してから、
プログラムを作ることが、重要です。あわせて、空気抵抗がない場合に想定される
軌跡と比較してみると良いでしょう。プログラムのどこを修正すればよいか？
たとえば、こんな感じ！
t = t0
x = x0
y = y0
vx = v0*cosθ0 ! 角度の単位はラジアン
vy = v0*sinθ0
do i=1,imax
v = sqrt( vx**2 + vy**2 )
ax = - k2/m * v * vx
ay = -g - k2/m * v * vy
vxmid = vx + ax*dt/2
vymid = vy + ay*dt/2
xmid = x + vx*dt/2
ymid = y + vy*dt/2
vmid = sqrt( vxmid**2 + vymid**2 )
axmid = ...vmid...
aymid = ...vmid...
vx = ...
vy = ...
x = ...
y = ...
t = t + dt
write(*,'(....)') t, x, y, vx, vy, ax, ay
enddo
グラフ作成：
放出角θをいろいろ変えた軌跡を描く。
角度θと到達距離の関係を示すグラフ(横軸：θ、縦軸：到達距離)
- 連絡・注意
- メールを使ったレポ−ト送信について
件名(Subject:)に必ず学籍番号とレポート番号を記入する
件名の例： 06s20xx rep01
これに従わないとレポートが行方不明になるかもしれないので注意
送信先のアドレスは taka(_#_)stu.hirosaki-u.ac.jp へ。但し、(_#_) を @ に修正して送信。
レポートの送信前に計算結果をよく検討することが必要です。
計算結果に明らかな誤りを含んだままで、なんの説明の無いレポートを送信することは慎むべきです。
- 提出するもの
- プログラムリスト。これはメールの本体へコピペする。
- 計算結果のグラフ。epsファイルとして作成し添付する。
- 結果についての考察。これもメールの本体へ書き込む。
第05回 Oct. 31, 2007 / (Last modified Nov.01,2007 17:00)
- 講義
- 単振動の運動方程式を数値的に解く
- 抵抗力が無い場合
f = - k x
ma = - k x
x(t) = A cos(ω0t+δ) （解析解）
ただし ω0^2 = k/m
a = - ω0^2x としてEuler-Richardson法にかける （数値解）
全エネルギー E = (1/2)mv^2 + (1/2)kx^2 = (1/2)kA^2 は時間とともに変化しない。
- 速さに比例する抵抗力が働く場合
f = - kx - mγv
ma = - kx - mγv
a = - ω0^2 x - γv としてEuler-Richardson法にかける （数値解）
- 速さに比例する抵抗力および外力F(t)による駆動力が働く場合
f = - kx - mγv + (1/m)F(t)
ma = - kx - mγv + (1/m)F(t)
(1/m)F(t) = A0 cos ωt とすれば
a = - ω0^2 x - γv + A0 cos ωt
- 演習
ex05.f
- 抵抗力が無い場合
t = t0 <-- 0.0
x = x0 <-- 1.0
v = v0 <-- 0.0
w2 = w02 <-- 9.0
dt = 0.01
imax = 600
Em0 = 0.5 * v**2 + 0.5 * w2 * x**2 <-- 単位質量あたりの全エネルギー
write(*,'(1p6e15.7)') t, x, v, a, Em, 1.0
do i=1, imax
a = -w2 * x
vmid = v + a*dt/2
xmid = x + v*dt/2
amid = - w2 * xmid
v = v + amid*dt
x = x + vmid*dt
t = t + dt
Em = 0.5*v**2 + 0.5*w2*x**2
dE = Em - Em0
ratio = dE / Em0
write(*,'(1p6e15.7)') t, x, v, a, Em, ratio
enddo
t-x 曲線を描く
x-v 曲線を描く
t-ratio 曲線を描くbr>
dt を 0.2, 0.1, 0.05, 0.01, 0.0005 などに変えて、上記曲線がどう変化するか観察し、
考察しておくこと。
運動開始 6秒までに全エネルギーの値の精度が 0.1% 以下に保つた
めには、時間幅 dt はどれだけが適切か？
解析解との比較を試みる
- 速さに比例する抵抗力が働く場合
a = -w2 * x - gamma*v
.
amid = - w2 * xmid - gamma*vmid
.
.
t-x 曲線を描く
x-v 曲線を描く
t-ratio 曲線を描く
速さに比例する抵抗力および外力F(t)による駆動力が働く場合
レポート課題なし
連絡・注意
- 有限時間幅 dt ごとの数値計算ですので、必ず真の解とのずれがあります。
ここでは、抵抗力が無い場合は全エネルギーが一定になっていることを使って、確かめられます。
- 振動現象について良く復習しておくこと
- この減衰振動には、解析解があります。これも復習！
- 次回は強制振動およびLRC回路について、扱う予定ですので良く予習をしておくこと
第06回 Nov. 07, 2007 / (Last modified Nov.09,2007 17:00)
- 講義
- 外力F(t)により駆動された減衰調和振動子
a = - ω0^2 x - γv + F(t)/m
F(t)/m = A0 cos (ω・ t)
A0, ω を変化させて、 t-x-v曲線を描き、その形の変化を考察する
- RC回路（直列）
キルヒホッフの法則を復習
R・(dQ(t)/dt) = R・I(t) = Vs(t) - Q(t)/C Q(t)についての一階の微分方程式
Vs(t) = cos (ω・t)
- 演習
- 外力F(t)により駆動された減衰調和振動子の場合
ex05.f に追加修正する
.
gamma = 0.5
omega = 2.0 <--- ωのこと
amplitude = 1.0 <--- A0のこと
.
a = -w2 * x - gamma*v + amplitude * cos ( omega * t )
.
amid = - w2 * xmid - gamma*vmid + amplitude * cos ( omega * t )
.
.
t-x 曲線を描く
t-v 曲線を描く
x-v 曲線を描く
t-x-v 曲線を描く
- レポート課題なし
- 連絡・注意
- 設定したパラメータと表示された曲線の形状との関係を、よく考察すること。
第07回 Nov. 14, 2007 / (Last modified Nov.16,2007 18:00)
- 講義
- rep01の解説（次回にプログラム例を配布します）
砲丸の落下点の求め方
- RC回路について(low pass filter, high pass filter, 積分回路, 微分回路)
- 演習
- 1階の微分方程式
- 外部からの入力電圧と取り出す(観察する)端子間の出力電圧の時間的変化を観察し考察する
- レポート課題なし
- 連絡・注意
- 次回こそLRC回路(これは2階の微分方程式になります)
第08回 Nov. 21, 2007 / (Last modified Nov.22,2007 14:00)
- 講義
- rep01(ex04.f)のプログラム例を配布しました。自作したものと比較し、試してみること。
サブルーチンの利用法を良く復習しておくこと
- LRC回路の説明 (変数および係数の対応関係に注意)
L・(d^2Q/dt^2) + R・(dQ/dt) + Q/C = Vs(t)
プログラムを作る際は、
ex05.f で扱った
m(d^2 x/dt^2) = - kx - Γ(dx/dt) + F(t)
との比較をしてみる。
また
ex06.f で扱った
R・(dQ(t)/dt) = R・I(t) = Vs(t) - Q(t)/C
とも比較する。
- 演習
- ex07.f とする
ex05.f 及び ex06.f を参考として、オイラーで作成する
たとえば、R=100Ω、C=3.0μF、L=2.0mH
Vs(t) = V0 cos (ω・t)、V0=1.0V
周波数(又は角振動数ω)をいろいろ変えて、電流値I(t)を観察する
横軸：時間、縦軸：電流の図を作成する
- レポート課題なし（ex07.fを完成できた場合は提出すること。提出は全員必須ではありません）
- 連絡
- １２月５日（水）のこの授業は休講です
（弘前大学ドリーム講座のため県内の高校へ出張するためです）
第09回 Nov. 28, 2007 / (Last modified Dec.01,2007 12:00)
- 講義
- ex07.fを自由提出としたところ4名の提出だけでした。
できないのか？ できていても出さないのか？
- Euler Richardson法について再度の説明
- 2階の微分方程式の解法
- 変数の対応関係
| 運動方程式 | LRC回路のキルヒホッフの法則 |
| m(d^2x/dt^2) = - kx - Γ(dx/dt) + F(t) |
L(d^2Q/dt^2) + R(dQ/dt) + Q/C = Vs(t) |
| x | Q |
| v | I |
| 加速度 | d^2Q/dt^2 |
- 演習
- R=100Ω C=3.0μF L=2.0mH Vs(t)=V0cosωt V0=1.0V として
Euler Richardson法を用いて数値解を求めて考察する
- 横軸：t 縦軸：I
- 横軸：t 縦軸：電圧【 Vs(t), L, R, C それぞれにかかる電位差 】
- 横軸：周波数 縦軸：定常状態のImax
- レポート課題
レポート番号：rep02
Program file name : ex07.f
- 連絡・注意
Dec. 05, 2007 休講 / (Last modified Nov.22,2007 14:00)
第10回 Dec. 12, 2007 / (Last modified Dec.13,2007 15:00)
- 講義
- 確率論的シミュレーション
一様乱数の発生法（の復習）
- 乱数の度数分布 n(i )
- 確率密度の分布 p(i )
分割数、サンプリング数に依存しないグラフ
p(i ) = n(i ) / dx / Nsample
- 演習
- 組み込み関数
iseed = 123457
call srand(iseed)
do i=1,100
enddo
end
- 合同法
iy = 1234567
do i=1,100
call ransu(iy,r)
write(*,*) r
enddo
end
subroutine ransu(iy,r)
integer iy
real*8 r
iy = iy * 48828125
if(iy.lt.0) iy = ( iy + 2147483647 ) + 1
r = dfloat( iy ) * 4.6566128e-10
return
end
- 期待される頻度が平均値の周りにどれくらいばらついているかを観察する。
分割数 = 100, Nsample = 100万 の場合
分割数 = 500, Nsample = 100万 の場合
分割数 =1000, Nsample = 100万 の場合,,,,
- レポート課題 無し
- 連絡・注意
- 次回は各頻度値と平均値の差のばらつき度合いを数値化します（-->分布を見る）
この分布を正規分布にあてはめてみますので、正規分布（ガウス分布）を予習してくること
- gnuplot について
plot "data.dat" u 1:2 with histeps
この場合、x座標の値には各区間の中央の値を記入する
第11回 Dec. 19, 2007 / (Last modified Jan.09,2008 12:00)
- 講義
- レポート２の解説、プログラム例の配布
- 次のことを予測して計算を実行することが大事。
(1) 繰り返し計算が何回おこなわれているか？
(2) 計算時間はどれだけになるか？
(3) ファイルに書き出した容量はどれだけになっているか？
- 一様乱数の分布のばらつき度合いを調べる。
度数の平均値と度数の分布
- 演習
- レポート課題なし
- 連絡・注意
第12回 Jan. 09, 2008 / (Last modified Jan.17,2008 20:15)
- 講義
- 一様乱数の分布のばらつき度合いを調べる。
--> 正規分布をしている
- データの解析
--> 時間差の分布をみる --> 指数分布になっている
--> 日付変更となる時刻データの扱いに注意
- 演習
- ex08.f ex08.plt
- 070109.DAT
open(1,file='070109.DAT')
do i=1,20000
read(1,'(3i2)') hh,mm,ss
enddo
close(1)
- 連絡・注意
第13回 Jan. 16, 2008 / (Last modified Jan.17,2008 15:15)
- 講義
- 一様乱数の分布のばらつき度合いを調べる。
--> 正規分布曲線を重ねて描く
--> 縦軸が分割数Nやサンプリング数Imaxに依存しないようにする
--> 分割数 N，サンプリング数 Imax を変化させ、分布の形を観察する。
--> 平均値と標準偏差値との関係を観察する。
- データ解析について(次回のための解説)
時間差分布 --> 平均時間差を求めておく
1時間ごとの飛来数の分布 --> バラツキがあっても平均約１５０位。
指数分布の曲線を重ねて描く
この指数分布に従う乱数を発生して、飛来時間のシミュレーションをする
- 演習
- 連絡・注意
- 次回は通常授業をします。
期末試験は1月３０日に延期します。
- 次回にレポート３の課題を出します。
第14回 Jan. 23, 2008 / (Last modified Jan.17,2008 15:15)
期末試験の予定を、次回1月３０日に延期します。
第15回 Jan. 30, 2008 / (Last modified Jan.23,2008 15:15)
期末試験を実施します。
レポート提出状況 ( ○印：レポ−トメール受信， △印：内容不十分， _：未着信 )
(Last modified 09:45, Feb.26, 2008)
レポート2の不完全な部分は急ぎ再提出すること。未提出の者も必ず提出すること。(1/22)
修正提出および未提出レポートの提出締め切りは２月２２日(金)とします。(2/11)
| 出題日 |
10/24 |
11/28 |
01/23 |
_/__ |
_/__ |
|---|
| 締切日 |
11/09 |
12/14 (last check 2/22 18:15) |
02/22 (last check 2/26 09:45) |
_/__ |
_/__ |
|---|
| レポート番号 |
rep01 |
rep02 |
rep03 |
_____ |
_____ |
|---|
| Program file name |
ex04.f |
ex07.f |
ex09.f |
_____ |
_____ |
|---|
| 配点(暫定満点) |
10点 |
10点 |
15点 |
____ |
____ |
|---|
| program / figure |
p | f | p | f | p | f | p | f | p | f |
| 03 | ○ | ○ | _ | _ | _ | _ | _ | _ | _ | _ |
| 04 | ○ | ○ | ○ | ○ |
○ | △(fig0901, fig0903, fig0904-2,以外は開かないよ？) |
_ | _ | _ | _ |
| 05 | ○ | ○ | ○ |
△(図3まだ変だ！横軸を対数目盛りで描くと、左右対称になる)2/19 |
○ | ○ |
_ | _ | _ | _ |
| 06 | ○ | ○ | △(もう少し考察欲しいね) |
図3もうちょっとですね |
○ | ○ |
_ | _ | _ | _ |
| 07 | ○ | ○ | ○ |
図3もうちょっとですね |
○ | ○ |
_ | _ | _ | _ |
| 12 | ○ | ○ | △(もっと考察必要) |
図1,2,3どれも開かないぞ
今回も，ex07.eps, ex07-1.eps, ex07-2.epsどれも開かないよ？(2/22,13:00) |
△(考察足りない) | ○ |
_ | _ | _ | _ |
| 16 | ○ | ○ | ○ | ○ |
△ | △(2番目の図は縦軸を対数表示で) |
_ | _ | _ | _ |
| 17 | ○ | ○ | ○ | ○ |
○ | ○(もう少し考察が欲しかったね) |
_ | _ | _ | _ |
| 19 | ○ | ○ | ○ |
図3が無い |
△ |
△(もう少しがんばりませんか) |
_ | _ | _ | _ |
| 20 | ○ | ○ | △(考察必要) |
図1,2,3どれも開かないぞ |
_ | _ | _ | _ | _ | _ |
| 25 | ○ | ○ | ○ |
図1,図3開かないぞ |
△ |
△(もう少しがんばりませんか) |
_ | _ | _ | _ |
| 31 | ○ | ○ | ○ |
○(図3 今度はok!) |
○ | ○(３番目の図ちょっとへんかな！) |
_ | _ | _ | _ |
| 32 | ○ | ○ | ○ | ○ |
△ |
△(考察でも取り上げているのだが，report09-3-1.eps, report09-3-2.eps で赤線と緑色の線が離れているのはへんだぞ。いづれの場合でも一致するはず。) |
_ | _ | _ | _ |
| 33 | ○ | ○ | _ | _ | _ | _ | _ | _ | _ | _ |
| 35 | ○ | ○ | ○ | ○ |
△(考察不十分) | ○(縦軸・横軸のlabelくらい書きなさい) |
_ | _ | _ | _ |
| 38 | ○ | ○ | ○ | ○ |
○ | ○ |
_ | _ | _ | _ |
| 48 | ○ | ○ | ○ |
△(図3の左半分に不要な点が多数ありますね)2/19 |
○ | ○ |
_ | _ | _ | _ |
| program / figure |
p | f | p | f | p | f | p | f | p | f |
| レポート番号 |
rep01 |
rep02 |
rep03 |
_____ |
_____ |
|---|
taka(_#_)stu.hirosaki-u.ac.jp