exp(x)

Fortranプログラミング入門
-課題の解説-
関根 晃太
問題
nを整数型，xを倍精度実数型の変数を宣言せよ．nとx
の値をそれぞれread文で読み込む．そのとき，指数関数
のテイラー展開
を計算せよ．さらに，数学関数exp(x)の値と計算した値
をそれぞれ表示し，比較せよ．ファイル名はexptay.f90と
せよ．
P-1
解答
program lowexp
implicit none
real(8) :: x, y, ans
integer :: i, j, n
実行:
write(*,*) 'Please input number n:'
n=100, x=1
read(*,*) n
exp(x) = 2.712818284590451
write(*,*) 'Please input number x'
ans
= *Infinity
read(*,*) x
ans = 1d0
y=1
j=1
do i=1, n
j=j*i
階乗の計算のため，
y=y*x
n>=13でオーバーフロー…
ans = ans + y/j
end do
write(*,*) exp(x), ans
stop
P-2
end program lowexp
解答
program lowexp2
implicit none
real(8) :: x, y, ans
integer :: i, j, n
実行:
write(*,*) 'Please input number n:'
n=100, x=1
read(*,*) n
exp(x) = 2.712818284590451
write(*,*) 'Please input number x'
ans
= 2.712818284590455
read(*,*) x
ans = 1d0
y=1
j=1
do i=1, n
y = y * x/i
階乗の計算を回避!
ans = ans + y
end do
write(*,*) exp(x), ans
stop
end program lowexp2
P-3
問題
課題4
nを整数型，xを倍精度実数型の変数を宣言せよ．nとxの値をそれぞれread文
で読み込む．そのとき，指数関数のテイラー展開
を計算せよ．さらに，数学関数exp(x)の値と計算した値をそれぞれ表示し，比較
せよ．
にてxを負の値にすると精度が悪くなる．xの値が負のと
き精度が良くなるように改造せよ．
ファイル名をmexptay.f90とせよ.
P-4
解答
program lowexp2
implicit none
real(8) :: x, y, ans
integer :: i, j, n
write(*,*) 'Please input number n:'
read(*,*) n
write(*,*) 'Please input number x'
read(*,*) x
ans = 1d0
y=1
j=1
do i=1, n
y = y * x/i
ans = ans + y
end do
write(*,*) exp(x), ans
stop
end program lowexp2
実行:
n=100, x=-20
exp(x) = 2.06115362243855786E-009
ans
= 5.62188447213041761E-009
1桁もあっていない…
P-5
解答
program lowexp2
implicit none
real(8) :: x, y, ans
integer :: i, j, n
write(*,*) 'Please input number n:'
read(*,*) n
write(*,*) 'Please input number x'
read(*,*) x
ans = 1d0
y=1
j=1
do i=1, n
y = y * x/i
ans = ans + y
end do
write(*,*) exp(x), ans
stop
end program lowexp2
yとansを出力してみると…
y
-20
200
-1333.33
6666.667
-26666.7
88888.89
-253968
634920.6
-1410935
2821869
ans
1
-19
181
-1152.33
5514.333
-21152.3
67736.56
-186232
448688.9
-962246
1859624
近い大きな数の引き算は誤差が大きくなる
P-6
解答
program lowexp2
implicit none
real(8) :: x, y, ans
integer :: i, j, n
write(*,*) 'Please input number n:'
read(*,*) n
write(*,*) 'Please input number x'
read(*,*) x
ans = 1d0
y=1
j=1
do i=1, n
y = y * x/i
ans = ans + y
end do
write(*,*) exp(x), ans
stop
end program lowexp2
引き算を発生させないために…
のとき
とし，計算すればよい!!
P-7
解答
program highexp
implicit none
real(8) :: x, y, ans, ax
integer :: i, j, n
write(*,*) 'Please input number n:'
read(*,*) n
write(*,*) 'Please input number x'
read(*,*) x
ans = 1d0
y=1
ax = abs(x)
do i=1, n
y = y * ax/i
ans = ans + y
end do
if (x < 0) then
ans = 1/ans
end if
write(*,*) exp(x), ans
stop
end program highexp
P-8