下記の例では1/xという関数をx=1からx=2まで積分しています。答えは、log(2)-log(1)=log(2/1)=0.69314718055994530941
です。
bcコマンドを使った計算結果のチェック
$ bc -l
bc 1.06
Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'.
l(2)-l(1)
.69314718055994530941
シェルスクリプト：simspon.sh
内部でFortranプログラムとnamelistファイルを生成しています。
#!/bin/bash
# Description:
#
# Author: am
#
# Host: aofd165.fish.nagasaki-u.ac.jp
# Directory: /work1/am/TEACHING/TOOLS/Simpson
#
# Revision history:
# This file is created by /usr/local/mybin/nbscr.sh at 18:12 on 03-14-2015.
# 参考文献
# 戸川隼人, 数値計算リテラシ, サイエンス社. p.95
echo "Shell script, $(basename $0) starts."
echo
exe=integ_simpson
namelist=${exe}.namelist.txt
src=${exe}.f90
fort=ifort
opt="-traceback -check"
cat <<EOF>$namelist
¶
x0=1.0
x1=2.0
n=10
&end
EOF
cat <<EOF>$src
program $exe
real,allocatable::x(:),y(:)
integer nx
real x0,x1,integ_fx
namelist /para/x0,x1,n
read(*,nml=para)
allocate (x(0:n),y(0:n))
dx=(x1-x0)/float(n)
print *,'n=',n
print *,'dx=',dx
do i=0,n
x(i)=x0+dx*float(i)
y(i)=f(x(i))
print *,x(i),y(i)
enddo
call simpson(y,n,dx,integ_fx)
print *,'integ_fx=',integ_fx
stop
end program
subroutine simpson(y,n,h,s)
integer,intent(in)::n
real,intent(in)::y(0:n),h
real,intent(inout)::s
real s1,s2
s1=0.0
do i=1,n-1,2
s1=s1+y(i)
enddo
s2=0.0
do i=2,n-2,2
s2=s2+y(i)
enddo
s=(y(0)+4.0*s1+2.0*s2+y(n))*h/3.0
end subroutine
real function f(x)
real,intent(in)::x
f=1.0/x
!f=sin(x)
end function
EOF
$fort $opt $src -o $exe
if [ $? -eq 0 ]; then
${exe} < $namelist
else
echo Compile Error.
exit 1
fi
echo "Done $0"
echo
実行例
$ ./simpson.sh
Shell script, simpson.sh starts.
n= 10
dx= 0.1000000
1.000000 1.000000
1.100000 0.9090909
1.200000 0.8333333
1.300000 0.7692308
1.400000 0.7142857
1.500000 0.6666667
1.600000 0.6250000
1.700000 0.5882353
1.800000 0.5555556
1.900000 0.5263157
2.000000 0.5000000
integ_fx= 0.6931502
Done ./simpson.sh