4
$\begingroup$

I need to implement a non-linear fitting algorithm in Fortran and chose to use MINPACK's flavor of the Levenberg-Marquardt algorithm as a basis for the least-squares stuff. However, I seem to misunderstand something very basic about how these routines are to be used. My understanding is that I have to call the routine lmdif1 and define a function fcn to be called by the MINPACK routines that calculates the square of the difference between the data to be fitted and the model function; lmdif1 would then return the best solution for the free fitting parameters in an array x.

I have written a small test program (given below) to try this out. The program produces a set of "data" by calculating values from the model function (exp(ax) in this case) with a given parameter a (=-0.5 in this case) and changing them by a small random "error". This dataset is then fitted with the same model function, with a being the free parameter. But it doesn't really work, the fits are usually very bad, unless I almost give the "correct" solution as a starting guess.

Considering that Levenberg-Marquardt is said to be quite robust and seeing how well it works e.g. in Gnuplot, I wonder what I'm doing wrong here. Any ideas? The code below should work when linking to the MINPACK routines.

Tom

! test lmdif
program testanf
implicit none
integer, parameter :: m=10,n=1,lwa=m*n+5*n+m
integer :: info,iwa(n)
double precision :: x(n),fvec(m),wa(lwa),tol=1d-3
external :: fcn
x(1)=-0.5d0
call lmdif1(fcn,m,n,x,fvec,tol,info,iwa,wa,lwa)
call translate_info(info,n,tol)
write(*,*) 'a=',x(1)
end program testanf

subroutine fcn(m,n,x,fvec,iflag)
implicit none
! m: number of data points
! n: number of fitting parameters
integer, intent(in) :: m,n,iflag
integer :: i
real, allocatable, save :: t(:),fdata(:)
real :: tmin,tmax,dt,rand
double precision, intent(inout) :: x(n)
double precision, intent(out) :: fvec(m)
character(len=1) :: c
logical, save :: init=.true.
if (init) then
!  create "measured data"
   tmin=0.
   tmax=10.
   dt=(tmax-tmin)/m
   call random_seed(put=(/1,2,3,4,5,6,7,8/))
   allocate(fdata(m),t(m))
   do i=1,m
      t(i)=i*dt
      fdata(i)=fmodel(-0.5d0,t(i))
   end do
   tmax=0.1*maxval(fdata)
   do i=1,m
      call random_number(rand)
      fdata(i)=fdata(i)+(rand-0.5)*tmax ! add measurement errors
      write(10,*) t(i),fdata(i)
   end do
   init=.false.
end if
! minimize the residuals
do i=1,m
   fvec(i)=(fdata(i)-fmodel(x(1),t(i)))**2
end do
if (iflag < 0) deallocate(t,fdata)
return
contains
! model function
  real function fmodel(a,t)
  implicit none
  real, intent(in) :: t
  double precision, intent(in) :: a
  fmodel=exp(a*t)
  return
  end function fmodel
end subroutine fcn
$\endgroup$
  • 1
    $\begingroup$ Here is how to call minpack from modern Fortran: github.com/certik/minpack/blob/master/examples/… $\endgroup$ – Ondřej Čertík Aug 5 '13 at 2:25
  • $\begingroup$ Thanks for pointing me to this resource. The main mistake in my program seems to be that I write the square of the difference into fvec, as opposed to just the difference. However, in terms of the quality of the result, that didn't help much - the crucial improvement in this respect seems to be to use double precision throughout, as in your example. If further tests I'll do confirm that, I'd find it quite stunning that it would make such a difference. $\endgroup$ – TomR Aug 5 '13 at 12:09
  • $\begingroup$ Yeah, we would have to debug it. I know that my example works, as I tested it pretty well, so you can try to modify it for your function and see. $\endgroup$ – Ondřej Čertík Aug 5 '13 at 21:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.