github: https://github.com/devernay/cminpack
主页: http://devernay.github.io/cminpack/
使用手册: http://devernay.github.io/cminpack/man.html
lmdif
,lmdif1_
- 最小化非线性函数平方和
include <minpack.h>
void lmdif1_(void (*fcn)(int *m, int *n, double *x, double *fvec, int *iflag),
int *m, int *n, double *x, double *fvec,
double *tol, int *info, int *iwa, double *wa, int *lwa);
void lmdif_(void (*fcn)(int *m, int *n, double *x, double *fvec, int *iflag),
int *m, int *n, double *x, double *fvec,
double *ftol, double *xtol, double *gtol, int *maxfev, double *epsfcn, double *diag,
int *mode, double *factor, int *nprint, int *info, int *nfev, double *fjac,
int *ldfjac, int *ipvt, double *qtf,
double *wa1, double *wa2, double *wa3, double *wa4 );
lmdif_
的目的是最小化m个n元非线性方程的平方和,使用的方法是LM算法的改进。用户需要提供计算方程的子程序。Jacobian矩阵会通过一个前向差分(forward-difference)近似计算得到。
lmdif1_
是相同的目的,但是调用方法更简单一些。
这些函数是通过FORTRAN写的,如果从C调用,需要记住以下几点:
void fcn(int m, int n, double *x, double *fvec, int *iflag) {
/* 计算函数在x点的值,通过fvec返回。*/
}
iflag的值不能被fcn所修改,除非用户想要终止lmdif
/lmdif1_
。在这个例子中iflag设置为负整数。m:函数个数;
n:变量个数(n<=m)
x:长度为n的数组,设置为初始的估计解向量。输出的时候x内容为最终估计的解向量。
fvec:输出长度为m的数组,内容为最终输出x计算得到的函数解。
tol:作为输入,非负数。用于函数终止的条件判断:
info:作为输出。如果用户终止了函数的执行,info将被设置为iflag的值(负数)(详细见fcn的描述),否则,info的值如下几种情况:
iwa:长度n的工作数组;
wa:长度lwa的工作数组;
lwa:作为输入,整数,不能小于mn+5n+m;
[NOTE] 这三个输入我也不知道作用,从样例来看不需要初始化。
暂时不用这部分,跳过。
/* driver for lmdif1 example. */
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <minpack.h>
#define real __minpack_real__
// 用户自定义的函数f()
// real -> __cminpack_real__ -> 浮点数(double)
void fcn(const int *m, const int *n, const real *x, real *fvec, int *iflag);
int main()
{
int j, m, n, info, lwa, iwa[3], one=1;
real tol, fnorm, x[3], fvec[15], wa[75];
// 函数个数15; 变量数3
m = 15;
n = 3;
/* the following starting values provide a rough fit. */
// 初始位置
// 1.e0 = 1.0e0 = 1.0
x[0] = 1.e0;
x[1] = 1.e0;
x[2] = 1.e0;
// 为什么要设置75?
lwa = 75;
/* set tol to the square root of the machine precision. unless high
precision solutions are required, this is the recommended
setting. */
// (建议打印一下看值是多少)
tol = sqrt(__minpack_func__(dpmpar)(&one));
// 需要注意指针
__minpack_func__(lmdif1)(&fcn, &m, &n, x, fvec, &tol, &info, iwa, wa, &lwa);
// 最终的2范数(即平方和开根号)
fnorm = __minpack_func__(enorm)(&m, fvec);
printf(" final l2 norm of the residuals%15.7g\n\n", (double)fnorm);
printf(" exit parameter %10i\n\n", info);
printf(" final approximate solution\n");
for (j=1; j<=n; j++) {
printf("%s%15.7g", j%3==1?"\n ":"", (double)x[j-1]);
}
printf("\n");
return 0;
}
// The problem is to determine the values of x(1), x(2), and x(3)
// which provide the best fit (in the least squares sense) of
// x(1) + u(i)/(v(i)*x(2) + w(i)*x(3)), i = 1, 15
// to the data
// y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39,
// 0.37,0.58,0.73,0.96,1.34,2.10,4.39),
// where u(i) = i, v(i) = 16 - i, and w(i) = min(u(i),v(i)). The
// i-th component of FVEC is thus defined by
// y(i) - (x(1) + u(i)/(v(i)*x(2) + w(i)*x(3))).
void fcn(const int *m, const int *n, const real *x, real *fvec, int *iflag)
{
/* function fcn for lmdif1 example */
int i;
real tmp1,tmp2,tmp3;
// 实际的y值
real y[15]={1.4e-1,1.8e-1,2.2e-1,2.5e-1,2.9e-1,3.2e-1,3.5e-1,3.9e-1,
3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0};
assert(*m == 15 && *n == 3);
if (*iflag == 0) {
/* insert print statements here when nprint is positive. */
/* if the nprint parameter to lmder is positive, the function is
called every nprint iterations with iflag=0, so that the
function may perform special operations, such as printing
residuals. */
// 这段没有很看懂,在??情况下打印信息
return;
}
/* compute residuals */
for (i=0; i<15; i++) {
tmp1 = i+1;
tmp2 = 15 - i;
tmp3 = tmp1;
if (i >= 8) {
tmp3 = tmp2;
}
fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
}
}
[NOTE] 其他内容有待更新
原文:https://www.cnblogs.com/bemfoo/p/11203993.html