1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
| !****************************************************
! 非線形方程式の解法（二分法）
! * 方程式： y = x + cos(2 * x) - sin(4 * x)
!
! date name version
! 2018.10.11 mk-mode.com 1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
program nonlinear_equation_bisection
implicit none
! SP: 単精度(4), DP: 倍精度(8)
integer, parameter :: SP = kind(1.0)
integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
integer(SP), parameter :: NMAX = 50
real(DP), parameter :: EPS = 1.0e-6
logical :: stat
integer(SP) :: n
real(DP) :: x, y, x_0, x_1, sgn
write (*, "(a)", advance="no") "x_0, x_1 : "
read (*,*) x_0, x_1
! x_0 > x_1 であれば、交換
if (x_0 > x_1) then
x = x_0
x_0 = x_1
x_1 = x
end if
! f(x_0) と f(x_1) の符号が逆なら、終了
if (f(x_0) * f(x_1) >= 0.0) then
write (*,*) "f(x_0) * f(x_1) >= 0.0"
stop
end if
sgn = sign(1.0_8, f(x_1) - f(x_0))
stat = .false.
do n = 1, NMAX
x = (x_0 + x_1) * 0.5_8
y = f(x)
! 収束判定
if (abs(x_1 - x_0) < EPS) then
stat = .true.
write(*, '("収束 [", i4, "]")') n
exit
else
write(*, '("誤差 [", i4, "] = ", e20.8)') n, abs(y)
end if
! 次の値を推定
if (y * sgn < 0.0) then
x_0 = x
else
x_1 = x
end if
end do
! 結果出力
if (.not. stat) then
write (*, *) "近似不可！"
end if
write (*, '("近似値 = ", e20.8)') x
write (*, '("誤差 = ", e20.8)') abs(y)
stop
contains
! 方程式
! * f = x + cos(2 * x) - sin(4 * x)
!
! :param real(8) x
! :return real(8) f
real(DP) function f(x)
implicit none
real(DP), intent(in) :: x
f = x + cos(2 * x) - sin(4 * x)
end function f
end program nonlinear_equation_bisection
|
