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