forked from mnhrdt/ccmath
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathC05-roots
207 lines (155 loc) · 8.2 KB
/
C05-roots
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
Chapter 5
ROOTS and OPTIMA
Summary
The functions in this library segment are used to
find the roots of nonlinear equations and systems
of such equations and to locate local extrema of
functions. The areas covered are:
o Roots of Equations
o Optimization
-------------------------------------------------------------------------------
Notes on Contents
Functions in this section of the library are used to extract the roots of
nonlinear equations and to find local optima.
o Roots of Equations:
solnl -------- solve a system of nonlinear equations.
solnlx ------- solve a nonlinear system using an initial
Jacobian.
secrt -------- solve a nonlinear equation using the secant
method.
plrt --------- find the roots (real and complex) of a polynomial.
polyc -------- evaluate a polynomial for complex arguments.
o Optimization:
optmiz ------- perform an unconstrained search for the
optimal parameter vector.
optsh -------- perform a golden section search for optima in
one dimension.
-------------------------------------------------------------------------------
General Technical Comments
The general root finding functions employ a matrix update algorithm
developed by Broyden. The second form offers an enlarged domain of
convergence at the expense of requiring the Jacobian matrix of the system at
the initial point of the search. The secant method is a simple and standard
root extraction technique, provided for convenience. Roots of a polynomial
with real coefficients are computed by plrt.
The nonlinear optimal search algorithm uses the BFGS matrix update
method. The robustness of this procedure has been confirmed in numerous tests
of optimization techniques. The golden section method is useful and efficient
for bounded search in one-dimension.
Functions for finding the roots of a nonlinear system of equations
and for nonlinear optimization are based on efficient quasi-Newton matrix
update algorithms. This approach has been shown to combine high execution
efficiency (ie. a small number of function evaluations) with a large
convergence domain. Simplicity of use was also an important criterion in the
selection of the algorithms. The functions implemented do not require user
supplied procedures for computation of derivatives. Fortunately, this does
not entail a performance sacrifice, since the quasi-Newton procedures exhibit
excellent convergence properties when operating with numerical derivatives.
The algorithms employed for finding roots and optima are iterative.
They require an initial estimate of the solution as input. While the
implemented algorithms have robust convergence, refinement of these initial
estimates is often the most effective way to improve algorithm performance.
A few hours spent estimating appropriate scales for variables and refining
initial point estimates may save days of agonizing search for a "golden
algorithm". Nonlinear analysis is a field so rich that vast bodies of code
are no substitute for a thoughtful understanding of the problem at hand.
-------------------------------------------------------------------------------
FUNCTION SYNOPSES
-------------------------------------------------------------------------------
Roots of Equations:
-------------------------------------------------------------------------------
solnl
Solve a system of nonlinear equations.
int solnl(double *x,double *f,double (*fvec[])(),int n,double test)
x = pointer to array containing initial solution estimate,
converted to solution vector at exit
f = pointer to array containing final values of the system functions
{ f[k]=(*fvec[k])(x) }
fvec = array of pointers to the system functions
{ (*fvec[k])(x)=0 for k=0,1,2,--,n-1 }
n = dimension of system
test = convergence threshold (f~*f <test)
return: 1 -> normal exit
0 -> convergence failure
-------------------------------------------------------------------
solnx
Variant of solnl, which permits input of an initial system
Jacobian. (Has an enlarged convergence domain.)
int solnx(double *x,double *f,double (*fvec)(),double *jm,int n, \
double *test)
x = pointer to array containing initial solution estimate,
converted to solution vector at exit
f = pointer to array containing final values of the system
functions { f[k]=(*fvec[k])(x) }
fvec = array of pointers to the system functions
{ (*fvec[k])(x)=0 for k=0,1,2,--,n-1 }
jm = pointer to array of dimension n*n containing an initial system
Jacobian ( row order: JM[i,j] = df[i]/dx[j] )
n = dimension of system
test = convergence threshold (f~*f <test)
return: 1 -> normal exit
0 -> convergence failure
--------------------------------------------------------------------
plrt
Find the real and complex roots of a polynomial with real
coefficients.
typedef struct complex {double re,im;} Cpx;
int plrt(double *cof,int n,Cpx *root,double ra,double rb)
cof = pointer to array containing polynomial coefficients, with
the leading coefficient cof[n] != 0
root = pointer to array of complex structures, containing the
polynomial roots at exit (dimension n)
n = degree of polynomial
ra,rb = initial root estimates, with
( rb>0 -> real = ra, imag = rb and
rb<0 -> real roots ra+rb, ra-rb )
return: 0 -> normal exit
m>0 -> convergence failure for roots k<m
The polynomial is given by
Pn(x) = Sum(i=0 t0 n) cof[i]*x^i .
--------------------------------------------------------------------
polyc
Evaluate a polynomial with real coefficients at complex values
of the argument z.
typedef struct complex {double re,im;} Cpx;
Cpx polyc(Cpx z,double *cof,int n)
z = complex structure containing the value of the argument
cof = pointer to array containing the polynomial coefficients
(dim = n) [see plrt above]
n = degree of polynomial
return: p = complex structure containing the value of Pn(z)
-------------------------------------------------------------------
secrt
Compute a function root using the secant method.
double secrt(double (*func)(),double x,double dx,double test)
func = pointer to function
x = initial estimate of root position
dx = displacement to 2nd root estimate x'=x-dx
test = convergence threshold ( |x'-x|<test )
return: y = root of equation (*func)(y)=0
-------------------------------------------------------------------------------
Optimization:
-------------------------------------------------------------------------------
optmiz
Unconstrained optimal search algorithm of BFGS type.
int optmiz(double *x,int n,double (*func)(),double de, \
double test,int max)
func = pointer to function to be minimized
( min. f = (*func)(x) )
x = pointer to array of function parameters containing
an initial estimate at entry and the solution at exit
n = dimension of parameter vector x
de = interval used for derivative computation
test = convergence threshold ( |f'-f|<test )
max = iteration limit for line searches
return: m>0 -> success (iteration count)
m=0 -> convergence failure
---------------------------------------------------------------------
optsh
Conduct a golden section line search for a function minima.
double optsch(double (*func)(),double a,double b,double test)
func = pointer to function to be minimized
a,b = bounds on function argument (a<=x<=b)
test = convergence threshold (length of section
containing minima <test)
return: y = location of minimum of f = (*func)(y)