1 /* multiroots/test_funcs.c
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 #include <gsl/gsl_vector.h>
23 #include <gsl/gsl_matrix.h>
24 #include <gsl/gsl_multiroots.h>
26 #include "test_funcs.h"
28 /* For information on testing see the following paper,
30 J.J More, B.S. Garbow, K.E. Hillstrom, "Testing Unconstrained
31 Optimization Software", ACM Transactions on Mathematical Software,
32 Vol 7, No 1, (1981) p 17-41
36 /* Rosenbrock Function */
38 gsl_multiroot_function_fdf rosenbrock =
45 rosenbrock_initpt (gsl_vector * x)
47 gsl_vector_set (x, 0, -1.2);
48 gsl_vector_set (x, 1, 1.0);
52 rosenbrock_f (const gsl_vector * x, void *params, gsl_vector * f)
54 double x0 = gsl_vector_get (x, 0);
55 double x1 = gsl_vector_get (x, 1);
58 double y1 = 10 * (x1 - x0 * x0);
60 gsl_vector_set (f, 0, y0);
61 gsl_vector_set (f, 1, y1);
63 params = 0; /* avoid warning about unused parameters */
69 rosenbrock_df (const gsl_vector * x, void *params, gsl_matrix * df)
71 double x0 = gsl_vector_get (x, 0);
75 double df10 = -20 * x0;
78 gsl_matrix_set (df, 0, 0, df00);
79 gsl_matrix_set (df, 0, 1, df01);
80 gsl_matrix_set (df, 1, 0, df10);
81 gsl_matrix_set (df, 1, 1, df11);
83 params = 0; /* avoid warning about unused parameters */
89 rosenbrock_fdf (const gsl_vector * x, void *params,
90 gsl_vector * f, gsl_matrix * df)
92 rosenbrock_f (x, params, f);
93 rosenbrock_df (x, params, df);
99 /* Freudenstein and Roth function */
101 gsl_multiroot_function_fdf roth =
108 roth_initpt (gsl_vector * x)
110 gsl_vector_set (x, 0, 4.5); /* changed from the value in the paper */
111 gsl_vector_set (x, 1, 3.5); /* otherwise the problem is too hard */
115 roth_f (const gsl_vector * x, void *params, gsl_vector * f)
117 double x0 = gsl_vector_get (x, 0);
118 double x1 = gsl_vector_get (x, 1);
120 double y0 = -13.0 + x0 + ((5.0 - x1)*x1 - 2.0)*x1;
121 double y1 = -29.0 + x0 + ((x1 + 1.0)*x1 - 14.0)*x1;
123 gsl_vector_set (f, 0, y0);
124 gsl_vector_set (f, 1, y1);
126 params = 0; /* avoid warning about unused parameters */
132 roth_df (const gsl_vector * x, void *params, gsl_matrix * df)
134 double x1 = gsl_vector_get (x, 1);
137 double df01 = -3 * x1 * x1 + 10 * x1 - 2;
139 double df11 = 3 * x1 * x1 + 2 * x1 - 14;
141 gsl_matrix_set (df, 0, 0, df00);
142 gsl_matrix_set (df, 0, 1, df01);
143 gsl_matrix_set (df, 1, 0, df10);
144 gsl_matrix_set (df, 1, 1, df11);
146 params = 0; /* avoid warning about unused parameters */
152 roth_fdf (const gsl_vector * x, void *params,
153 gsl_vector * f, gsl_matrix * df)
155 roth_f (x, params, f);
156 roth_df (x, params, df);
163 /* Powell badly scaled function */
165 gsl_multiroot_function_fdf powellscal =
172 powellscal_initpt (gsl_vector * x)
174 gsl_vector_set (x, 0, 0.0);
175 gsl_vector_set (x, 1, 1.0);
179 powellscal_f (const gsl_vector * x, void *params, gsl_vector * f)
181 double x0 = gsl_vector_get (x, 0);
182 double x1 = gsl_vector_get (x, 1);
184 double y0 = 10000.0 * x0 * x1 - 1.0;
185 double y1 = exp (-x0) + exp (-x1) - 1.0001;
187 gsl_vector_set (f, 0, y0);
188 gsl_vector_set (f, 1, y1);
190 params = 0; /* avoid warning about unused parameters */
196 powellscal_df (const gsl_vector * x, void *params, gsl_matrix * df)
198 double x0 = gsl_vector_get (x, 0);
199 double x1 = gsl_vector_get (x, 1);
201 double df00 = 10000.0 * x1, df01 = 10000.0 * x0;
202 double df10 = -exp (-x0), df11 = -exp (-x1);
204 gsl_matrix_set (df, 0, 0, df00);
205 gsl_matrix_set (df, 0, 1, df01);
206 gsl_matrix_set (df, 1, 0, df10);
207 gsl_matrix_set (df, 1, 1, df11);
209 params = 0; /* avoid warning about unused parameters */
215 powellscal_fdf (const gsl_vector * x, void *params,
216 gsl_vector * f, gsl_matrix * df)
218 powellscal_f (x, params, f);
219 powellscal_df (x, params, df);
225 /* Brown badly scaled function */
227 gsl_multiroot_function_fdf brownscal =
234 brownscal_initpt (gsl_vector * x)
236 gsl_vector_set (x, 0, 1.0);
237 gsl_vector_set (x, 1, 1.0);
241 brownscal_f (const gsl_vector * x, void *params, gsl_vector * f)
243 double x0 = gsl_vector_get (x, 0);
244 double x1 = gsl_vector_get (x, 1);
246 double y0 = x0 - 1e6;
247 double y1 = x0 * x1 - 2;
249 gsl_vector_set (f, 0, y0);
250 gsl_vector_set (f, 1, y1);
252 params = 0; /* avoid warning about unused parameters */
258 brownscal_df (const gsl_vector * x, void *params, gsl_matrix * df)
260 double x0 = gsl_vector_get (x, 0);
261 double x1 = gsl_vector_get (x, 1);
263 double df00 = 1.0, df01 = 0.0;
264 double df10 = x1, df11 = x0;
266 gsl_matrix_set (df, 0, 0, df00);
267 gsl_matrix_set (df, 0, 1, df01);
268 gsl_matrix_set (df, 1, 0, df10);
269 gsl_matrix_set (df, 1, 1, df11);
271 params = 0; /* avoid warning about unused parameters */
277 brownscal_fdf (const gsl_vector * x, void *params,
278 gsl_vector * f, gsl_matrix * df)
280 brownscal_f (x, params, f);
281 brownscal_df (x, params, df);
287 /* Powell Singular Function */
289 gsl_multiroot_function_fdf powellsing =
296 powellsing_initpt (gsl_vector * x)
298 gsl_vector_set (x, 0, 3.0);
299 gsl_vector_set (x, 1, -1.0);
300 gsl_vector_set (x, 2, 0.0);
301 gsl_vector_set (x, 3, 1.0);
305 powellsing_f (const gsl_vector * x, void *params, gsl_vector * f)
307 double x0 = gsl_vector_get (x, 0);
308 double x1 = gsl_vector_get (x, 1);
309 double x2 = gsl_vector_get (x, 2);
310 double x3 = gsl_vector_get (x, 3);
312 double y0 = x0 + 10 * x1;
313 double y1 = sqrt (5.0) * (x2 - x3);
314 double y2 = pow (x1 - 2 * x2, 2.0);
315 double y3 = sqrt (10.0) * pow (x0 - x3, 2.0);
317 gsl_vector_set (f, 0, y0);
318 gsl_vector_set (f, 1, y1);
319 gsl_vector_set (f, 2, y2);
320 gsl_vector_set (f, 3, y3);
322 params = 0; /* avoid warning about unused parameters */
328 powellsing_df (const gsl_vector * x, void *params, gsl_matrix * df)
330 double x0 = gsl_vector_get (x, 0);
331 double x1 = gsl_vector_get (x, 1);
332 double x2 = gsl_vector_get (x, 2);
333 double x3 = gsl_vector_get (x, 3);
335 double df00 = 1, df01 = 10, df02 = 0, df03 = 0;
336 double df10 = 0, df11 = 0, df12 = sqrt (5.0), df13 = -df12;
337 double df20 = 0, df21 = 2 * (x1 - 2 * x2), df22 = -2 * df21, df23 = 0;
338 double df30 = 2 * sqrt (10.0) * (x0 - x3), df31 = 0, df32 = 0, df33 = -df30;
340 gsl_matrix_set (df, 0, 0, df00);
341 gsl_matrix_set (df, 0, 1, df01);
342 gsl_matrix_set (df, 0, 2, df02);
343 gsl_matrix_set (df, 0, 3, df03);
345 gsl_matrix_set (df, 1, 0, df10);
346 gsl_matrix_set (df, 1, 1, df11);
347 gsl_matrix_set (df, 1, 2, df12);
348 gsl_matrix_set (df, 1, 3, df13);
350 gsl_matrix_set (df, 2, 0, df20);
351 gsl_matrix_set (df, 2, 1, df21);
352 gsl_matrix_set (df, 2, 2, df22);
353 gsl_matrix_set (df, 2, 3, df23);
355 gsl_matrix_set (df, 3, 0, df30);
356 gsl_matrix_set (df, 3, 1, df31);
357 gsl_matrix_set (df, 3, 2, df32);
358 gsl_matrix_set (df, 3, 3, df33);
360 params = 0; /* avoid warning about unused parameters */
366 powellsing_fdf (const gsl_vector * x, void *params,
367 gsl_vector * f, gsl_matrix * df)
369 powellsing_f (x, params, f);
370 powellsing_df (x, params, df);
378 gsl_multiroot_function_fdf wood =
385 wood_initpt (gsl_vector * x)
387 gsl_vector_set (x, 0, -3.0);
388 gsl_vector_set (x, 1, -1.0);
389 gsl_vector_set (x, 2, -3.0);
390 gsl_vector_set (x, 3, -1.0);
394 wood_f (const gsl_vector * x, void *params, gsl_vector * f)
396 double x0 = gsl_vector_get (x, 0);
397 double x1 = gsl_vector_get (x, 1);
398 double x2 = gsl_vector_get (x, 2);
399 double x3 = gsl_vector_get (x, 3);
401 double t1 = x1 - x0 * x0;
402 double t2 = x3 - x2 * x2;
404 double y0 = -200.0 * x0 * t1 - (1 - x0);
405 double y1 = 200.0 * t1 + 20.2 * (x1 - 1) + 19.8 * (x3 - 1);
406 double y2 = -180.0 * x2 * t2 - (1 - x2);
407 double y3 = 180.0 * t2 + 20.2 * (x3 - 1) + 19.8 * (x1 - 1);
409 gsl_vector_set (f, 0, y0);
410 gsl_vector_set (f, 1, y1);
411 gsl_vector_set (f, 2, y2);
412 gsl_vector_set (f, 3, y3);
414 params = 0; /* avoid warning about unused parameters */
420 wood_df (const gsl_vector * x, void *params, gsl_matrix * df)
422 double x0 = gsl_vector_get (x, 0);
423 double x1 = gsl_vector_get (x, 1);
424 double x2 = gsl_vector_get (x, 2);
425 double x3 = gsl_vector_get (x, 3);
427 double t1 = x1 - 3 * x0 * x0;
428 double t2 = x3 - 3 * x2 * x2;
430 double df00 = -200.0 * t1 + 1, df01 = -200.0 * x0, df02 = 0, df03 = 0;
431 double df10 = -400.0*x0, df11 = 200.0 + 20.2, df12 = 0, df13 = 19.8;
432 double df20 = 0, df21 = 0, df22 = -180.0 * t2 + 1, df23 = -180.0 * x2;
433 double df30 = 0, df31 = 19.8, df32 = -2 * 180 * x2, df33 = 180.0 + 20.2;
435 gsl_matrix_set (df, 0, 0, df00);
436 gsl_matrix_set (df, 0, 1, df01);
437 gsl_matrix_set (df, 0, 2, df02);
438 gsl_matrix_set (df, 0, 3, df03);
440 gsl_matrix_set (df, 1, 0, df10);
441 gsl_matrix_set (df, 1, 1, df11);
442 gsl_matrix_set (df, 1, 2, df12);
443 gsl_matrix_set (df, 1, 3, df13);
445 gsl_matrix_set (df, 2, 0, df20);
446 gsl_matrix_set (df, 2, 1, df21);
447 gsl_matrix_set (df, 2, 2, df22);
448 gsl_matrix_set (df, 2, 3, df23);
450 gsl_matrix_set (df, 3, 0, df30);
451 gsl_matrix_set (df, 3, 1, df31);
452 gsl_matrix_set (df, 3, 2, df32);
453 gsl_matrix_set (df, 3, 3, df33);
455 params = 0; /* avoid warning about unused parameters */
461 wood_fdf (const gsl_vector * x, void *params,
462 gsl_vector * f, gsl_matrix * df)
464 wood_f (x, params, f);
465 wood_df (x, params, df);
471 /* Helical Valley Function */
473 gsl_multiroot_function_fdf helical =
480 helical_initpt (gsl_vector * x)
482 gsl_vector_set (x, 0, -1.0);
483 gsl_vector_set (x, 1, 0.0);
484 gsl_vector_set (x, 2, 0.0);
488 helical_f (const gsl_vector * x, void *params, gsl_vector * f)
490 double x0 = gsl_vector_get (x, 0);
491 double x1 = gsl_vector_get (x, 1);
492 double x2 = gsl_vector_get (x, 2);
499 t1 = atan(x1/x0) / (2.0 * M_PI);
503 t1 = 0.5 + atan(x1/x0) / (2.0 * M_PI);
507 t1 = 0.25 * (x1 > 0 ? +1 : -1);
510 t2 = sqrt(x0*x0 + x1*x1) ;
512 y0 = 10 * (x2 - 10 * t1);
516 gsl_vector_set (f, 0, y0);
517 gsl_vector_set (f, 1, y1);
518 gsl_vector_set (f, 2, y2);
520 params = 0; /* avoid warning about unused parameters */
526 helical_df (const gsl_vector * x, void *params, gsl_matrix * df)
528 double x0 = gsl_vector_get (x, 0);
529 double x1 = gsl_vector_get (x, 1);
531 double t = x0 * x0 + x1 * x1 ;
532 double t1 = 2 * M_PI * t ;
533 double t2 = sqrt(t) ;
535 double df00 = 100*x1/t1, df01 = -100.0 * x0/t1, df02 = 10.0;
536 double df10 = 10*x0/t2, df11 = 10*x1/t2, df12 = 0;
537 double df20 = 0, df21 = 0, df22 = 1.0;
539 gsl_matrix_set (df, 0, 0, df00);
540 gsl_matrix_set (df, 0, 1, df01);
541 gsl_matrix_set (df, 0, 2, df02);
543 gsl_matrix_set (df, 1, 0, df10);
544 gsl_matrix_set (df, 1, 1, df11);
545 gsl_matrix_set (df, 1, 2, df12);
547 gsl_matrix_set (df, 2, 0, df20);
548 gsl_matrix_set (df, 2, 1, df21);
549 gsl_matrix_set (df, 2, 2, df22);
551 params = 0; /* avoid warning about unused parameters */
557 helical_fdf (const gsl_vector * x, void *params,
558 gsl_vector * f, gsl_matrix * df)
560 helical_f (x, params, f);
561 helical_df (x, params, df);
567 /* Discrete Boundary Value Function */
571 gsl_multiroot_function_fdf dbv =
578 dbv_initpt (gsl_vector * x)
581 double h = 1.0 / (N + 1.0);
583 for (i = 0; i < N; i++)
585 double t = (i + 1) * h;
586 double z = t * (t - 1);
587 gsl_vector_set (x, i, z);
592 dbv_f (const gsl_vector * x, void *params, gsl_vector * f)
596 double h = 1.0 / (N + 1.0);
598 for (i = 0; i < N; i++)
600 double z, ti = (i + 1) * h;
601 double xi = 0, xim1 = 0, xip1 = 0;
603 xi = gsl_vector_get (x, i);
606 xim1 = gsl_vector_get (x, i - 1);
609 xip1 = gsl_vector_get (x, i + 1);
611 z = 2 * xi - xim1 - xip1 + h * h * pow(xi + ti + 1, 3.0) / 2.0;
613 gsl_vector_set (f, i, z);
617 params = 0; /* avoid warning about unused parameters */
623 dbv_df (const gsl_vector * x, void *params, gsl_matrix * df)
627 double h = 1.0 / (N + 1.0);
629 for (i = 0; i < N; i++)
630 for (j = 0; j < N; j++)
631 gsl_matrix_set (df, i, j, 0.0);
633 for (i = 0; i < N; i++)
635 double dz_dxi, ti = (i + 1) * h;
637 double xi = gsl_vector_get (x, i);
639 dz_dxi = 2.0 + (3.0 / 2.0) * h * h * pow(xi + ti + 1, 2.0) ;
641 gsl_matrix_set (df, i, i, dz_dxi);
644 gsl_matrix_set (df, i, i-1, -1.0);
647 gsl_matrix_set (df, i, i+1, -1.0);
651 params = 0; /* avoid warning about unused parameters */
657 dbv_fdf (const gsl_vector * x, void *params,
658 gsl_vector * f, gsl_matrix * df)
660 dbv_f (x, params, f);
661 dbv_df (x, params, df);
666 /* Trigonometric Function */
668 gsl_multiroot_function_fdf trig =
675 trig_initpt (gsl_vector * x)
679 for (i = 0; i < N; i++) /* choose an initial point which converges */
681 gsl_vector_set (x, i, 0.05);
686 trig_f (const gsl_vector * x, void *params, gsl_vector * f)
691 for (i = 0; i < N; i++)
693 sum += cos(gsl_vector_get(x,i));
696 for (i = 0; i < N; i++)
698 double xi = gsl_vector_get (x,i);
699 double z = N - sum + (i + 1) * (1 - cos(xi)) - sin(xi);
701 gsl_vector_set (f, i, z);
704 params = 0; /* avoid warning about unused parameters */
710 trig_df (const gsl_vector * x, void *params, gsl_matrix * df)
714 for (i = 0; i < N; i++)
716 for (j = 0; j < N; j++)
719 double xi = gsl_vector_get(x, i);
720 double xj = gsl_vector_get(x, j);
723 dz = sin(xi) + (i + 1) * sin(xi) - cos(xi);
727 gsl_matrix_set(df, i, j, dz);
731 params = 0; /* avoid warning about unused parameters */
737 trig_fdf (const gsl_vector * x, void *params,
738 gsl_vector * f, gsl_matrix * df)
740 trig_f (x, params, f);
741 trig_df (x, params, df);