我需要在 Fortran 中实现一个非线性拟合算法,并选择使用 MinPACK 的 Levenberg-Marquardt 算法作为最小二乘法的基础。但是,我似乎误解了有关如何使用这些例程的一些非常基本的东西。我的理解是我必须调用例程 lmdif1 并定义一个函数 fcn 由MINPACK例程调用,该函数计算要拟合的数据与模型函数之间的差异的平方;然后 lmdif1 将返回数组 x 中自由拟合参数的最佳解。
我编写了一个小测试程序(如下所示)来尝试一下。该程序通过使用给定参数a(在这种情况下=-0.5)计算模型函数(在这种情况下为exp(ax))的值并通过一个小的随机“错误”来改变它们来生成一组“数据”。然后用相同的模型函数对该数据集进行拟合,其中 a 是自由参数。但它并没有真正起作用,拟合通常非常糟糕,除非我几乎给出“正确”的解决方案作为开始猜测。
考虑到据说 Levenberg-Marquardt 非常健壮并且看到它在 Gnuplot 中的工作情况如何,我想知道我在这里做错了什么。有任何想法吗?当链接到MINPACK 例程时,下面的代码应该可以工作。
汤姆
! 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