第六章 数值计算方法举例
6-1二分法(Biction)
二分法是最简单的解法，该算法只有简单的几个步骤。
（1）先猜两个值a、b，使得f(a)*f(b)小于0，也就是f(a)、f(b)必须异号。这样才能保证在a与b间存在一个c值，使得f(c)=0。
（2）令c=(a+b)/2,如果f(c)=0就找到了一个解，工作完成。
（3）f(c)不为时，如果f(a)、f(b)异号，则以a、c为新的两个猜测值来重复步骤2；如果f(b)、f(c)异号，则以b、c为新的猜测值来重复步骤2。
编程举例：用二分法求解方程及方程的根。
1. program main
2. implicit none
3. real :: a, b !两个猜测值
4. real :: ans !算出的值
5. real, external :: bict, f1, f2
6. do while(.true.)
7. write(*,*)＂输入两个猜测值＂
8. read(*,*)a, b
9. !f(a)*f(b)<0的猜测值才是有效的猜测
10. if(f1(a)*f1(b)<0)exit
11. write(*,*)＂不正确的猜测＂
12. end do
13. !调用二分法求根的函数
14. ans=bict(a, b, f1)
15. ！显示结果
16. write(*,＂(＇x=＇,F6.3)＂) ans
17.
18. do while(.true.)
19. write(*,*)＂输入两个猜测值＂
20. read(*,*)a, b
21. !f(a)*f(b)<0的猜测值才是有效的猜测
22. if(f2(a)*f2(b)<0)exit
23. write(*,*)＂不正确的猜测＂
24. end do
25. !调用二分法求根的函数
26. ans=bict(a, b, f2)
27. ！显示结果
28. write(*,＂(＇x=＇,F6.3)＂) ans
29. stop
30. end
31. real function bict(a, b, func)
32. implicit none
33. real, parameter :: zero=0.00001
34. real :: a, b !输入的猜测
35. real :: c !用来计算(a+b)/2
36. real :: fa, fb, fc !用来记录f(a),f(b),f(c)
37. real, external :: func !所要求解的函数
38. ！先求出c, f(c)的值
39. c = (a+b)/2.0
40. fc=func(c)
41. !f(c)小于zero时，就视f(c)=0,结束循环
42. do while (abs(fc)>zero)
43. fa=func(a)
44. fb=func(b)
45. if (fa*fc<0) then
46. !f(a)*f(c)<0,以a, c值为新的区间
47. b=c
48. c=(a+b)/2.0
49. el
50. !不然就是以b, c为新的区间
51. a=c
52. c=(a+b)/2.0
53. end if
54. !求出新的f(c)值
55. fc=func(c)
56. end do
57. bict=c
58. return
59. end function
60. !求解用的函数f1(x)
61. real function f1(x)
62. implicit none
63. real :: x
64. f1=sin(x)-x
65. return
66. end function
67. !求解用的函数f2(x)
68. real function f2(x)
69. implicit none
70. real :: x
71. f2=x**2-x-20
72. return
73. end function
编程举例：用二分法计算在三维空间内由坐标原点发出的任意一条射线与椭球面
的交点（坐标），并指出迭代多少步。
program main
implicit none
！函数biscet返回值是数组，要用interface说明函数接口
interface
function bict(x)
implicit none
real ::x(3),bict(3)
end function
end interface
integer, parameter :: n=3
real :: ans(n),b(n)
real :: fb
real ,external :: f
write(*,*)'输入一个方向'
read(*,*)b
fb=f(b)
if (fb>0 ) goto 10
do while (.true.)
b=2*b
fb=f(b)
if (fb>0)exit
end do
10 ans = bict(b)
write(*,"('x(n)=',3(f6.3))")ans
stop
end
function bict(b)
implicit none
integer, parameter :: n=3
integer :: iter
real, parameter :: zero=1.0e-6
real :: a(n)=0.0
real :: b(n),c(n),bict(n)
real :: fa,fb,fc
real, external :: f
c=(a+b)/2
fc=f(c)
do while(abs(fc)>zero)
fa=f(a)
fb=f(b)
if(fa*fc<0.0)then
b=c
c=(a+b)/2
el
a=c
c=(a+b)/2
end if
fc=f(c)
iter=iter+1
end do
bict=c
write(*,"('iter=',i5)")iter
return
end function
real function f(x)
implicit none
integer, parameter :: n=3
real :: x(n)
f=(x(1)**2)/2+(x(2)**2)/3+(x(3)**2)/4-1.0
return
end function
6-2 牛顿法
牛顿法是利用线段来逼近的结果，计算过程如下：
（1） 先给一个猜测值a。
（2） 以为斜率，经过（a, f(a)）作一条直线，令这条直线与x轴的交点为b。检
查f(b)是否为0，如果是就找到一个解。
（3） f(b)不为0时，重新令b为新的猜测值a，回到步骤（2）来重复。
若初始猜测值a取的好，f(b)会越来越接近0，程序实验如下：
program main
implicit none
real :: a !初始猜测值
real :: ans !解
real, external :: newton, func, dfunc
write(*,*)＂输入初始猜测值＂
read(*,*)a
!输入初始猜测值及要求值的函数
ans=Newton(a, func, dfunc)
write(*,＂(＇x=＇,F8.4)＂)ans
stop
end program
！牛顿法的函数
real function newton(a, f, df)
implicit none
real, parameter :: zero=0.00001
real :: a
real, external :: f ！输入的求值函数
real, external :: df ! f(x)的导数
real :: b !逼近得到的解
real :: fb !记录f(b)
b=a-f(a)/df(a)
fb=f(b)
!在趋于0之前要一直做逼近工作
do while(abs(fb)>zero)
a=b
b=a-f(a)/df(a)
fb=f(b)
end do
newton=b
return
end function Newton
!要求值的函数
real function func(x)
implicit none
real :: x
func=x**3-2*x**2+5*x-6
return
end function func
real function dfunc(x)
implicit none
real :: x
dfunc=3*x**2-4*x+5
return
end function dfunc
