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.
! 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