第 1 回 減衰方程式
支配方程式: cx dt
dx =− (1)
初期条件: x0
x = (2)
c 定 数 c 値 系 振 舞 い 変 わ し し パ ー タ
(parameter) 呼
練習 1 講義時間中 行う
初期条件 (2) , (1) 解析解を求 横軸 t, 縦軸 x し 解析解
フを け 手書 い
数値解法
今回 最 簡単 オイ ー法 オイ ー前進差分 を使う (1)を差分近似す
,
) ) (
( )
( cx t
t t x t t
x =−
∆
−
∆
+ . (2)
(2) 分母を払う t t cx t x t t
x( +∆ )= ()− ( )∆ (3)
(3)を初期条件
0 0)
(t t x
x = = (4)
解く
手順
ア ゴ 算法 呼 こ あ
[1] パ ータ c, 計算開始時刻t0，計算終了時刻tE, 時間刻 ∆t 値を設定す
[2] 初期条件(4) 値を設定す
[3] 計算終了時刻 ∆t 値 , 計算回数 n を求 [4] 初期条件 値, t0 x0を出力す
[5] t +∆t 値を求
[6] (3)を使 , x( t +∆t) 値を求 [7] t +∆t x( t +∆t) 値を出力す
[8] x( t +∆t) 値を x( t ) 値 置 換え アップ ー 呼 こ あ [9] n 回 [5]-[6]を繰 返す
練習 2 講義時間中 行う
t=t0 t+2∆t 上述 算法 従い手計算 数値解を 電卓等使用可
演習１ 講義終了後 行う
(1) 末尾 添付したプロ を入力し, 実行用シェ プ を使 実
行す 実行例
$ damp0d.run.sh
-rw-rw-r-- 1 am am 850 9 月 28 09:57 damp0d.f90
Compiling damp0d.f90 ... Done Compile.
-rwxrwxr-x 1 am am 786K 9 月 28 10:06 damp0d.exe
damp0d.exe is running ...
# t0= 0.00000
# te= 10.00000
# c = 0.10000
# dt= 0.05000
# n = 200
Output file : damp0d_result.txt
Done damp0d.exe.
(2) 横軸 t, 縦軸 x し 結果を図 す
(3) 練習 1 求 た解析解 (1) 数値解 差を フ せ
(4) dt 値を dt=0.5 変更し , (1) (3)を繰 返せ 結果 い 説明せ
(5)プロ 各部 内容 何を行 い を口頭 説明せ 次回発
表せ
プロ 例：damp0d.f90
program damp0d
character(len=5000)::outfile real::c,t0, te, x0, dt
namelist /para/ c,t0, te, x0, dt, outfile
read(*,para)
n=(te-t0)/dt
open(20,file=outfile)
write(20,'(A,f10.5)')'# t0=',t0 write(20,'(A,f10.5)')'# te=',te write(20,'(A,f10.5)')'# c =',c write(20,'(A,f10.5)')'# dt=',dt write(20,'(A,i10 )')'# n =',n write(20,'(A )')'# t, x'
t=t0 x=x0
write(20,'(f10.5,e14.6)')t,x
do i=1,N
t=dt*float(i)
xnew=x-c*x*dt !Time integration
x=xnew !Update
write(20,'(f10.5,e14.6)')t,x !Output result
enddo !i
print *
print '(A)','# Parameters:' print '(A,f10.5)', '# t0=',t0 print '(A,f10.5)', '# te=',te print '(A,f10.5)', '# c =', c print '(A,f10.5)', '# dt=',dt print '(A,i10 )', '# n =',n
print '(A)', 'Output file : ',trim(adjustl(outfile)) print *
end program damp0d
実行用 プ ：damp0d.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 c=0.1
t0=0.0 te=10.0 x0=1.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
¶ c=${c}, t0=${t0}, te=${te}, x0=${x0}, dt=${dt},
outfile="${outfile}",
&end EOF
# COMPILE SOURCE FILE echo Compiling $src ... ifort $opt $src -o $exe if [ $? -ne 0 ]; then echo
echo Compile error! echo
echo Terminated. echo
exit 1 fi
echo "Done Compile." echo
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 version 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
range=0/10/0/1.5 size=X4/4
xanot=a5f1 yanot=a0.5f0.1
anot=${xanot}:"t":/${yanot}:"x":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
# PRINT HEADER export LANG=C
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" 実行例
$ plot.sh damp0d_result.txt Bash script ./plot.sh starts.
pstext: Record 6 is incomplete (skipped)
Input : damp0d_result.txt Output : damp0d_result.ps
Done ./plot.sh
