1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
|
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
<html>
<!-- Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 The GSL Team.
Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.3 or
any later version published by the Free Software Foundation; with no
Invariant Sections and no cover texts. A copy of the license is
included in the section entitled "GNU Free Documentation License". -->
<!-- Created by GNU Texinfo 5.1, http://www.gnu.org/software/texinfo/ -->
<head>
<title>GNU Scientific Library – Reference Manual: Example programs for Multidimensional Root finding</title>
<meta name="description" content="GNU Scientific Library – Reference Manual: Example programs for Multidimensional Root finding">
<meta name="keywords" content="GNU Scientific Library – Reference Manual: Example programs for Multidimensional Root finding">
<meta name="resource-type" content="document">
<meta name="distribution" content="global">
<meta name="Generator" content="makeinfo">
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<link href="index.html#Top" rel="start" title="Top">
<link href="Function-Index.html#Function-Index" rel="index" title="Function Index">
<link href="Multidimensional-Root_002dFinding.html#Multidimensional-Root_002dFinding" rel="up" title="Multidimensional Root-Finding">
<link href="References-and-Further-Reading-for-Multidimensional-Root-Finding.html#References-and-Further-Reading-for-Multidimensional-Root-Finding" rel="next" title="References and Further Reading for Multidimensional Root Finding">
<link href="Algorithms-without-Derivatives.html#Algorithms-without-Derivatives" rel="previous" title="Algorithms without Derivatives">
<style type="text/css">
<!--
a.summary-letter {text-decoration: none}
blockquote.smallquotation {font-size: smaller}
div.display {margin-left: 3.2em}
div.example {margin-left: 3.2em}
div.indentedblock {margin-left: 3.2em}
div.lisp {margin-left: 3.2em}
div.smalldisplay {margin-left: 3.2em}
div.smallexample {margin-left: 3.2em}
div.smallindentedblock {margin-left: 3.2em; font-size: smaller}
div.smalllisp {margin-left: 3.2em}
kbd {font-style:oblique}
pre.display {font-family: inherit}
pre.format {font-family: inherit}
pre.menu-comment {font-family: serif}
pre.menu-preformatted {font-family: serif}
pre.smalldisplay {font-family: inherit; font-size: smaller}
pre.smallexample {font-size: smaller}
pre.smallformat {font-family: inherit; font-size: smaller}
pre.smalllisp {font-size: smaller}
span.nocodebreak {white-space:nowrap}
span.nolinebreak {white-space:nowrap}
span.roman {font-family:serif; font-weight:normal}
span.sansserif {font-family:sans-serif; font-weight:normal}
ul.no-bullet {list-style: none}
-->
</style>
</head>
<body lang="en" bgcolor="#FFFFFF" text="#000000" link="#0000FF" vlink="#800080" alink="#FF0000">
<a name="Example-programs-for-Multidimensional-Root-finding"></a>
<div class="header">
<p>
Next: <a href="References-and-Further-Reading-for-Multidimensional-Root-Finding.html#References-and-Further-Reading-for-Multidimensional-Root-Finding" accesskey="n" rel="next">References and Further Reading for Multidimensional Root Finding</a>, Previous: <a href="Algorithms-without-Derivatives.html#Algorithms-without-Derivatives" accesskey="p" rel="previous">Algorithms without Derivatives</a>, Up: <a href="Multidimensional-Root_002dFinding.html#Multidimensional-Root_002dFinding" accesskey="u" rel="up">Multidimensional Root-Finding</a> [<a href="Function-Index.html#Function-Index" title="Index" rel="index">Index</a>]</p>
</div>
<hr>
<a name="Examples-26"></a>
<h3 class="section">35.8 Examples</h3>
<p>The multidimensional solvers are used in a similar way to the
one-dimensional root finding algorithms. This first example
demonstrates the <code>hybrids</code> scaled-hybrid algorithm, which does not
require derivatives. The program solves the Rosenbrock system of equations,
with <em>a = 1, b = 10</em>. The solution of this system lies at
<em>(x,y) = (1,1)</em> in a narrow valley.
</p>
<p>The first stage of the program is to define the system of equations,
</p>
<div class="example">
<pre class="example">#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
struct rparams
{
double a;
double b;
};
int
rosenbrock_f (const gsl_vector * x, void *params,
gsl_vector * f)
{
double a = ((struct rparams *) params)->a;
double b = ((struct rparams *) params)->b;
const double x0 = gsl_vector_get (x, 0);
const double x1 = gsl_vector_get (x, 1);
const double y0 = a * (1 - x0);
const double y1 = b * (x1 - x0 * x0);
gsl_vector_set (f, 0, y0);
gsl_vector_set (f, 1, y1);
return GSL_SUCCESS;
}
</pre></div>
<p>The main program begins by creating the function object <code>f</code>, with
the arguments <code>(x,y)</code> and parameters <code>(a,b)</code>. The solver
<code>s</code> is initialized to use this function, with the <code>hybrids</code>
method.
</p>
<div class="example">
<pre class="example">int
main (void)
{
const gsl_multiroot_fsolver_type *T;
gsl_multiroot_fsolver *s;
int status;
size_t i, iter = 0;
const size_t n = 2;
struct rparams p = {1.0, 10.0};
gsl_multiroot_function f = {&rosenbrock_f, n, &p};
double x_init[2] = {-10.0, -5.0};
gsl_vector *x = gsl_vector_alloc (n);
gsl_vector_set (x, 0, x_init[0]);
gsl_vector_set (x, 1, x_init[1]);
T = gsl_multiroot_fsolver_hybrids;
s = gsl_multiroot_fsolver_alloc (T, 2);
gsl_multiroot_fsolver_set (s, &f, x);
print_state (iter, s);
do
{
iter++;
status = gsl_multiroot_fsolver_iterate (s);
print_state (iter, s);
if (status) /* check if solver is stuck */
break;
status =
gsl_multiroot_test_residual (s->f, 1e-7);
}
while (status == GSL_CONTINUE && iter < 1000);
printf ("status = %s\n", gsl_strerror (status));
gsl_multiroot_fsolver_free (s);
gsl_vector_free (x);
return 0;
}
</pre></div>
<p>Note that it is important to check the return status of each solver
step, in case the algorithm becomes stuck. If an error condition is
detected, indicating that the algorithm cannot proceed, then the error
can be reported to the user, a new starting point chosen or a different
algorithm used.
</p>
<p>The intermediate state of the solution is displayed by the following
function. The solver state contains the vector <code>s->x</code> which is the
current position, and the vector <code>s->f</code> with corresponding function
values.
</p>
<div class="example">
<pre class="example">int
print_state (size_t iter, gsl_multiroot_fsolver * s)
{
printf ("iter = %3u x = % .3f % .3f "
"f(x) = % .3e % .3e\n",
iter,
gsl_vector_get (s->x, 0),
gsl_vector_get (s->x, 1),
gsl_vector_get (s->f, 0),
gsl_vector_get (s->f, 1));
}
</pre></div>
<p>Here are the results of running the program. The algorithm is started at
<em>(-10,-5)</em> far from the solution. Since the solution is hidden in
a narrow valley the earliest steps follow the gradient of the function
downhill, in an attempt to reduce the large value of the residual. Once
the root has been approximately located, on iteration 8, the Newton
behavior takes over and convergence is very rapid.
</p>
<div class="smallexample">
<pre class="smallexample">iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03
iter = 1 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03
iter = 2 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01
iter = 3 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01
iter = 4 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01
iter = 5 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01
iter = 6 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01
iter = 7 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00
iter = 8 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00
iter = 9 x = 1.000 0.878 f(x) = 1.268e-10 -1.218e+00
iter = 10 x = 1.000 0.989 f(x) = 1.124e-11 -1.080e-01
iter = 11 x = 1.000 1.000 f(x) = 0.000e+00 0.000e+00
status = success
</pre></div>
<p>Note that the algorithm does not update the location on every
iteration. Some iterations are used to adjust the trust-region
parameter, after trying a step which was found to be divergent, or to
recompute the Jacobian, when poor convergence behavior is detected.
</p>
<p>The next example program adds derivative information, in order to
accelerate the solution. There are two derivative functions
<code>rosenbrock_df</code> and <code>rosenbrock_fdf</code>. The latter computes both
the function and its derivative simultaneously. This allows the
optimization of any common terms. For simplicity we substitute calls to
the separate <code>f</code> and <code>df</code> functions at this point in the code
below.
</p>
<div class="example">
<pre class="example">int
rosenbrock_df (const gsl_vector * x, void *params,
gsl_matrix * J)
{
const double a = ((struct rparams *) params)->a;
const double b = ((struct rparams *) params)->b;
const double x0 = gsl_vector_get (x, 0);
const double df00 = -a;
const double df01 = 0;
const double df10 = -2 * b * x0;
const double df11 = b;
gsl_matrix_set (J, 0, 0, df00);
gsl_matrix_set (J, 0, 1, df01);
gsl_matrix_set (J, 1, 0, df10);
gsl_matrix_set (J, 1, 1, df11);
return GSL_SUCCESS;
}
int
rosenbrock_fdf (const gsl_vector * x, void *params,
gsl_vector * f, gsl_matrix * J)
{
rosenbrock_f (x, params, f);
rosenbrock_df (x, params, J);
return GSL_SUCCESS;
}
</pre></div>
<p>The main program now makes calls to the corresponding <code>fdfsolver</code>
versions of the functions,
</p>
<div class="example">
<pre class="example">int
main (void)
{
const gsl_multiroot_fdfsolver_type *T;
gsl_multiroot_fdfsolver *s;
int status;
size_t i, iter = 0;
const size_t n = 2;
struct rparams p = {1.0, 10.0};
gsl_multiroot_function_fdf f = {&rosenbrock_f,
&rosenbrock_df,
&rosenbrock_fdf,
n, &p};
double x_init[2] = {-10.0, -5.0};
gsl_vector *x = gsl_vector_alloc (n);
gsl_vector_set (x, 0, x_init[0]);
gsl_vector_set (x, 1, x_init[1]);
T = gsl_multiroot_fdfsolver_gnewton;
s = gsl_multiroot_fdfsolver_alloc (T, n);
gsl_multiroot_fdfsolver_set (s, &f, x);
print_state (iter, s);
do
{
iter++;
status = gsl_multiroot_fdfsolver_iterate (s);
print_state (iter, s);
if (status)
break;
status = gsl_multiroot_test_residual (s->f, 1e-7);
}
while (status == GSL_CONTINUE && iter < 1000);
printf ("status = %s\n", gsl_strerror (status));
gsl_multiroot_fdfsolver_free (s);
gsl_vector_free (x);
return 0;
}
</pre></div>
<p>The addition of derivative information to the <code>hybrids</code> solver does
not make any significant difference to its behavior, since it able to
approximate the Jacobian numerically with sufficient accuracy. To
illustrate the behavior of a different derivative solver we switch to
<code>gnewton</code>. This is a traditional Newton solver with the constraint
that it scales back its step if the full step would lead “uphill”. Here
is the output for the <code>gnewton</code> algorithm,
</p>
<div class="smallexample">
<pre class="smallexample">iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03
iter = 1 x = -4.231 -65.317 f(x) = 5.231e+00 -8.321e+02
iter = 2 x = 1.000 -26.358 f(x) = -8.882e-16 -2.736e+02
iter = 3 x = 1.000 1.000 f(x) = -2.220e-16 -4.441e-15
status = success
</pre></div>
<p>The convergence is much more rapid, but takes a wide excursion out to
the point <em>(-4.23,-65.3)</em>. This could cause the algorithm to go
astray in a realistic application. The hybrid algorithm follows the
downhill path to the solution more reliably.
</p>
<hr>
<div class="header">
<p>
Next: <a href="References-and-Further-Reading-for-Multidimensional-Root-Finding.html#References-and-Further-Reading-for-Multidimensional-Root-Finding" accesskey="n" rel="next">References and Further Reading for Multidimensional Root Finding</a>, Previous: <a href="Algorithms-without-Derivatives.html#Algorithms-without-Derivatives" accesskey="p" rel="previous">Algorithms without Derivatives</a>, Up: <a href="Multidimensional-Root_002dFinding.html#Multidimensional-Root_002dFinding" accesskey="u" rel="up">Multidimensional Root-Finding</a> [<a href="Function-Index.html#Function-Index" title="Index" rel="index">Index</a>]</p>
</div>
</body>
</html>
|
