2．振動方程式
2.1 出
今回 振動現象を表現 簡単 形 方程式を考え 。振動論 関 知識 ，
波動現象を理解 基礎 ，今回 内容 ，大気 海洋 け 波，光 =
電磁波 性質を理解 うえ 非常 役 立 。
い ，あ 物体 バネ い 状況を考え 。 々 経験的 ，バネを引 張
際 ，バネ 伸び 大 け 大 いほ バネ 復元力 大 く を知 い 。こ
を数学的 表現 フッ 法則 あ 。フッ 法則 ，バネ 伸びをx
，
バネ 復元力=−kx, (2.1)
表さ 。右辺 負号 ，バネ 伸び 方向 逆向 バネを引 戻 う 力 働
くこ を意味 い 。ここ ，k バネ定数 呼 定数 あ 。今，物体 居
位置 x け動い 。こ こ を 物体 変位 x あ 言う。今，バネ 物
体 い ，物体 変位 x あ バネ x け伸び こ 。
物体 度 ，物体 変位x 時間微分 =微 時間内 x 変化 表さ 。
わ ， dt dx
v= / , (2.2)
あ 。 ，物体 加 度 a 度 時間微分 =微 時間内 v 変化 ，
dt dv
a= / , (2.3) あ 。式(2.4) (2.5) ，
2 2
dt x d dt dx dt
a d =
= , (2.4)
。今，物体 質量をm 。物体 働く力 バネ 復元力 け 仮定 ，式
(2.2) 示 ニュ トン 運動 第 2 法則 ， kx
dt x md =−
2 2
, (2.5)
いう方程式 得 。式(2.5)を，ここ 便宜的 振動方程式 称 。
練習 1 m= k =1, x=1.0，v=0.0 (at t=0) (2.5) 解析解(x, v)を求 。 練習 2 (2.7)式 表さ 振動 周期を m, k, πを使 表せ。
2.2 数値解
式(2.7)を数値的 解く 若干 工夫 必要 あ 。式(2.5) 時間t 2 階微分 表
さ ，こ を直接数値積分 や い あ 。ここ ，式(2.5)を う
2 式 分け 考え 。 v
dt
dx/ = , (2.6) kx
dt
dv/ =− , (2.7)
式(2.12),(2.13)を前進差分 近似 , , )
) ( ( )
( v t
t t x t t
x =
∆
−
∆
+ , (2.8)
) ) (
( )
( x t
m k t
t v t t
v =−
∆
−
∆
+ , (2.9)
。今回 式(2.8) (2.9)を数値的 解い 。
練習 3 m=1, k=0.1, x(t=0)=1, v(t=0)=0, ∆t=0.1 , t=0 0.3 x vを手計算 計算
せ 。電卓 使用可。
演習 1
末尾 添付 プロ ム(osc0d.f90) (2.8), (2.9) 相当 部分を作 せ 。
演習 2
完 させ プロ ム 内容 各部 う 処理を行 い ，口頭 説明せ 次
回行う 。
演習 3
初期条件 練習 1 同 。時間刻 を∆t=0.05 プロ ム osc0d.f90 を動 ,
以 示 図を作図せ 。
i)横軸 t, 縦軸 x, v を フ 時系列 呼ぶ 。
ii)横軸 x，縦軸 vを フ 相図 呼 こ あ 。
iii)横軸 t, 縦軸 x 解析解 x 数値解 差を フ
計算回数(n) 少 く 振動 5 周期以 う 設定せ 。
実行用シェ ス プト(osc0d.run.sh) 作図用ス プト サンプ (plot.sh)を末尾 添付
。
osc0d.run.sh 実行例
$ osc0d.run.sh
-rw-rw-r-- 1 am am 1.1K 9 月 29 10:12 osc0d.f90
Done Compile.
-rwxrwxr-x 1 am am 788K 9 月 29 10:21 osc0d.exe
osc0d.exe is running ...
# Parameters:
# m = 1.00000 # k = 1.00000 # t0= 0.00000 # te= 30.00000 # x0= 1.00000 # v0= 0.00000 # dt= 0.05000 # n = 600 Output file : osc0d_result.txt
Done osc0d.exe.
プロ ム例(osc0d.f90)
(2.8), (2.9) 相当 部分 自分 追加 program osc0d
character(len=5000)::outfile real::m, k, t0, te, x0, v0, dt
namelist /para/m, k,t0, te, x0, v0, dt, outfile
read(*,para)
pi = atan(1.)*4.
n=(te-t0)/dt
open(20,file=outfile)
write(20,'(A,f10.5)')'# m =',m write(20,'(A,f10.5)')'# k =',k write(20,'(A,f10.5)')'# t0=',t0 write(20,'(A,f10.5)')'# te=',te write(20,'(A,f10.5)')'# x0=',x0 write(20,'(A,f10.5)')'# v0=',v0 write(20,'(A,f10.5)')'# dt=',dt write(20,'(A,i10 )')'# n =',n write(20,'(A )')'# t, x, v'
t=t0 x=x0 v=v0
write(20,'(f10.5,2e14.6)')t,x,v
do i=1,N
t=t0+dt*float(i)
!ADD YOUR OWN CODE HERE
x=xnew v=vnew
write(20,'(f10.5,2e14.6)')t, x, v
enddo !i
print '(A)','# Parameters:' print '(A,f10.5)',' # m =',m print '(A,f10.5)',' # k =',k print '(A,f10.5)',' # t0=',t0 print '(A,f10.5)',' # te=',te print '(A,f10.5)',' # x0=',x0 print '(A,f10.5)',' # v0=',v0 print '(A,f10.5)',' # dt=',dt print '(A,i10 )',' # n =',n print *
print '(A)', 'Output file : ',trim(adjustl(outfile)) print *
end program osc0d
実行用シェ ス プト (osc0d.run.sh)
#!/bin/bash
#
# Description:
#
src=$(basename $0 .run.sh).f90 exe=$(basename $0 .run.sh).exe nml=$(basename $0 .run.sh).nml
f90=ifort
opt="-CB -traceback -fpe0"
# SET PARAMETERS m=1.0
k=1.0 t0=0.0 te=30.0 x0=1.0 v0=0.0 dt=0.05
outfile=$(basename $0 .run.sh)_result.txt
# CHECK IF SOURCE FILE SURELY EXISTS. if [ ! -f $src ]; then
echo
echo ERROR in $0 : No such file, $src echo
exit 1 fi
echo
ls -lh $src echo
# CREATE NAMELIST FILE cat <<EOF>$nml
¶ m=${m}, k=${k}, t0=${t0}, te=${te}, x0=${x0}, v0=${v0}, dt=${dt},
outfile="${outfile}",
&end EOF
# COMPILE SOURCE FILE echo Compiling $src ... ifort $opt $src -o $exe if [ $? -ne 0 ]; then echo
echo
echo Terminated. echo
exit 1 fi
echo "Done Compile." echo
ls -lh $exe echo
# RUN THE PROGRAM echo
echo $exe is running ... echo
$exe < $nml
if [ $? -ne 0 ]; then echo
echo ERROR in $exe: Runtime error! echo
echo Terminated. echo
exit 1 fi echo
echo "Done ${exe}." echo
作図用ス プト サンプ (plot.sh) Generic Mapping Tools 4.4.0 使用
#!/bin/bash
# Description:
#
# Author: am
#
# Host: localhost
# Directory: /work2/am/TEACHING/Numerical.Modelling/01.damp0d
#
# Revision history:
# This file is created by /usr/local/mybin/ngmt.sh at 09:58 on 09-28-2017.
echo "Bash script $0 starts."
gmtset MEASURE_UNIT INCH gmtset LABEL_FONT_SIZE 18 gmtset HEADER_FONT_SIZE 20 gmtset ANOT_FONT_SIZE 18 gmtset TICK_PEN 4
dash="t15_5:15" sp="¥040"
range=0/30/-3.0/3.0 size=X4/4
xanot=a5f1 yanot=a1.0f0.5
anot=${xanot}:"t":/${yanot}:"x${sp}or${sp}v":WSne
in=$1
if [ ! -f $in ]; then
echo Error in $0 : No such file, $in exit 1
fi
out=${figdir}$(basename $in .txt).ps
awk '{if ($1!="#") print $1,$2}' $in |
psxy -R$range -J$size -B$anot -W3 -X1.5 -Y3 -P -K >$out
awk '{if ($1!="#") print $1,$3}' $in |
psxy <<EOF -R -J -W3 -O -K >>$out 20 -1.2
22 -1.2 EOF
pstext <<EOF -R -J -O -K >>$out 22.5 -1.2 12 0 0 CM x
EOF
psxy <<EOF -R -J -W3${dash} -O -K >>$out 20 -1.4
22 -1.4 EOF
pstext <<EOF -R -J -O -K >>$out 22.5 -1.4 12 0 0 CM v
EOF
# PRINT HEADER export LANG=C xoffset=
yoffset=5 comment=
currentdir=`pwd`
if [ ${#currentdir} -gt 90 ]; then curdir1=${currentdir:1:90}
curdir2=${currentdir:91} else
curdir1=$currentdir curdir2="¥ "
fi
now=`date`
host=`hostname`
#comment=" "
time1=$(ls -l $in | awk '{print $6, $7, $8}')
pstext -JX6/1.2 -R0/1/0/1.2 -N -X${xoffset:-0} -Y${yoffset:-0} << EOF -O >> $out
0 1.1 9 0 1 LM $0 $@ 0 0.95 9 0 1 LM ${now} 0 0.80 9 0 1 LM ${host} 0 0.65 9 0 1 LM ${curdir1} 0 0.50 9 0 1 LM ${curdir2}
0 0.35 9 0 1 LM Input: ${in} (${time1}) 0 0.1 9 0 1 LM ${comment:-""}
EOF
# Output: ${out}
echo
echo Input : $in echo Output : $out echo
echo "Done $0"
