1 /*
2 * Library: lmfit (Levenberg-Marquardt least squares fitting)
3 *
4 * File: demo/curve1.c
5 *
6 * Contents: Example for the solution of 2 nonlinear equations in 2 variables.
7 * Find the intersection of a circle and a parabola.
8 *
9 * Note: Any modification of this example should be copied to the wiki.
10 *
11 * Author: Joachim Wuttke <j.wuttke@fz-juelich.de> 2013
12 *
13 * Licence: see ../COPYING (FreeBSD)
14 *
15 * Homepage: apps.jcns.fz-juelich.de/lmfit
16 */
17
18 #include "lmmin.h"
19 #include <stdio.h>
20 #include <stdlib.h>
21
evaluate_nonlin1(const double * p,int n,const void * data,double * f,int * info)22 void evaluate_nonlin1(
23 const double *p, int n, const void *data, double *f, int *info )
24 {
25 f[0] = p[0]*p[0] + p[1]*p[1] - 1; /* unit circle x^2+y^2=1 */
26 f[1] = p[1] - p[0]*p[0]; /* standard parabola y=x^2 */
27 }
28
29
main(int argc,char ** argv)30 int main( int argc, char **argv )
31 {
32 int n = 2; /* dimension of the problem */
33 double p[2]; /* parameter vector p=(x,y) */
34
35 /* auxiliary parameters */
36 lm_control_struct control = lm_control_double;
37 lm_status_struct status;
38 control.verbosity = 31;
39
40 /* get start values from command line */
41 if( argc!=3 ){
42 fprintf( stderr, "usage: nonlin1 x_start y_start\n" );
43 exit(-1);
44 }
45 p[0] = atof( argv[1] );
46 p[1] = atof( argv[2] );
47
48 /* the minimization */
49 printf( "Minimization:\n" );
50 lmmin( n, p, n, NULL, evaluate_nonlin1, &control, &status );
51
52 /* print results */
53 printf( "\n" );
54 printf( "lmmin status after %d function evaluations:\n %s\n",
55 status.nfev, lm_infmsg[status.outcome] );
56
57 printf( "\n" );
58 printf("Solution:\n");
59 printf(" x = %19.11f\n", p[0]);
60 printf(" y = %19.11f\n", p[1]);
61 printf(" d = %19.11f => ", status.fnorm);
62
63 /* convergence of lmfit is not enough to ensure validity of the solution */
64 if( status.fnorm >= control.ftol )
65 printf( "not a valid solution, try other starting values\n" );
66 else
67 printf( "valid, though not the only solution: "
68 "try other starting values\n" );
69
70 return 0;
71 }
72