Add secant method to optimize

src/optimize.f90

@@ -6,7 +6,7 @@ module optimize
 use utils, only: stop_error
 implicit none
 private
-public bisect
+public bisect,secant
 
 interface
     real(dp) function func(x)
@@ -44,4 +44,53 @@ real(dp) function bisect(f, a, b, tol) result(c)
 c = (a_ + b_)/2
 end function
 
+real(dp) function secant(f, a,b,tol, maxiter) result (c)
+!Solves f(x) = 0 on using the a and b as starting values
+procedure(func) :: f
+real(dp),intent(in) :: a,b,tol
+integer,optional :: maxiter
+
+integer :: maxiter_
+real(dp) :: xstart, xnext
+real(dp) :: fstart, fnext
+integer ::  i
+
+if(present(maxiter)) then
+    maxiter_ = maxiter
+else
+    maxiter_ = 100
+endif
+
+
+xstart = a
+xnext = b
+fstart = f(xstart)
+if(abs(fstart ) < tiny(1.0_dp)) then
+    c = xstart
+    return
+endif
+fnext = f(xnext)
+if(abs(fnext ) < tiny(1.0_dp)) then
+    c = xnext
+    return
+endif
+
+do i = 1, maxiter_
+
+    if( abs(fnext - fstart) < tiny(1.0_dp)) then
+       call stop_error("secant: division by zero")
+    endif
+    c = (xstart*fnext - xnext*fstart)/(fnext - fstart)
+    if (abs(c - xnext) < tol) return ! root found 
+    !update variables
+    xstart = xnext
+    fstart = fnext
+    xnext = c
+    fnext = f(c)
+enddo
+!max iterations number reached
+call stop_error("secant: method did not converge")
+end function secant
+
+
 end module

tests/optimize/CMakeLists.txt

@@ -3,6 +3,8 @@ include_directories(${PROJECT_BINARY_DIR}/src)
 project(optimize)
 
 add_executable(test_bisect test_bisect.f90)
+add_executable(test_secant test_secant.f90)
 target_link_libraries(test_bisect fortran_utils)
-
+target_link_libraries(test_secant fortran_utils)
 add_test(test_bisect ${PROJECT_BINARY_DIR}/test_bisect)
+add_test(test_secant ${PROJECT_BINARY_DIR}/test_secant)

tests/optimize/test_secant.f90

@@ -0,0 +1,51 @@
+program test_secant
+use types, only: dp
+use utils, only: assert
+use optimize, only: secant
+implicit none
+
+real(dp) :: x
+
+x = secant(f1, 0.0_dp, 0.1_dp, 1e-12_dp)
+ print *, x
+ call assert(abs(x - 0.73908513_dp) < 1e-8_dp)
+
+x = secant(f2, -4._dp, 4._dp, 1e-12_dp)
+print *, x
+call assert(abs(x - 0.86547403_dp) < 1e-8_dp)
+
+x = secant(f3, 0._dp, 3._dp, 1e-12_dp)
+print *, x
+call assert(abs(x - 2.40482556_dp) < 1e-8_dp)
+
+x = secant(f3, 3._dp, 6._dp, 1e-12_dp)
+print *, x
+call assert(abs(x - 5.52007811_dp) < 1e-8_dp)
+
+x = secant(f4, -4._dp, 4._dp, 1e-12_dp)
+print *, x
+call assert(abs(x) < 1e-8_dp)
+
+contains
+
+real(dp) function f1(x)
+real(dp), intent(in) :: x
+f1 = cos(x) - x
+end function
+
+real(dp) function f2(x)
+real(dp), intent(in) :: x
+f2 = cos(x) - x**3
+end function
+
+real(dp) function f3(x)
+real(dp), intent(in) :: x
+f3 = bessel_jn(0, x)
+end function
+
+real(dp) function f4(x)
+real(dp), intent(in) :: x
+f4 = x
+end function
+
+end program