summaryrefslogtreecommitdiff
path: root/support/dktools/dk4iter.h
blob: 9ad4f544c37aedae62bf6b2dd98fc165bbcd3130 (plain)
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
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
/*
Copyright (C) 2018-2020, Dirk Krause
SPDX-License-Identifier: BSD-3-Clause
*/

/*
	WARNING: This file was generated by the dkct program (see
	http://dktools.sourceforge.net/ for details).
	Changes you make here will be lost if dkct is run again!
	You should modify the original source and run dkct on it.
	Original source: dk4iter.ctr
*/

#ifndef DK4ITER_H_INCLUDED
/** Avoid multiple inclusions. */
#define DK4ITER_H_INCLUDED 1


#line 8 "dk4iter.ctr"

/**	@file	dk4iter.h	Iteration algorithms for root finding.

This module implements the following root finding iteration algorithms:
- Bisection
- Regula falsi (primitive form, Illinois, Pegasus, Anderson-Bjoerck)
- Newton
- fix point

The function to iterate must be implemented as dk4_iter_fct_t.
This function type returns an integer value (non-zero to indicate a
successful calculation, 0 to indicate an error).
The function expects the following arguments:
- Result address<br>
  Address of a variable or an array to store the calculation result.
  Functions for Newton iteration algorithm store 2 values: Function value
  and value of the first derivative.
- X position<br>
  The x value you want a function value for.
- Address of parameter set.
  Additional parameters probably required by the function. I.e. for a
  polynomial calculation function you may specify the coefficients here.
  This parameter is optional.

Details for an iteration may be specified in an iteration context.
The dk4iter_ctx_open() creates such a context and returns a pointer.
Use dk4iter_ctx_close() to release the context when done with it.
Alternatively use a automatic/static variable of the dk4_iter_ctx_t type
and initialize it using the dk4iter_ctx_init() function.

The dk4iter_ctx_set_algorithm() function chooses the iteration method.

For Newton and fixpoint the context may specify an x interval the
iteration must not leave. Leaving the interval results in abort.
The dk4iter_ctx_set_min(), dk4iter_ctx_set_exclusive_min(),
dk4iter_ctx_set_max() and dk4iter_ctx_set_exclusive_max() functions can
be used to set closed and open interval borders.

Tolerance values (epsilon) may be specified for
- y direction unless fixpoint is used and/or
- x direction.
Use dk4iter_ctx_set_eps_y() and dk4iter_ctx_set_eps_x() to set the
tolerances.
Use 0.0 or negative tolerances to skip one check (not both!).

Once the iteration reached allowed y difference you can decide
to stop when the x difference is in allowed tolerance range too or
to continue until full machine precision is reached (new x is exactly
the same as previous x).
I do not recommend attempts to iterate to full machine precision as this
(a) may fail for some functions resulting in an oscillation and
(b) uses a larger number of iteration passes and.
Use dk4iter_ctx_set_exact() to control this.

The dk4iter_ctx_set_maxpass() function sets the maximum number of
iteration passes (iteration steps). The iteration is aborted if there
is no success within this maximum number of passes.
Although you can use 0 to set an unlimited number of passes I do not recommend
to do so.

The dk4iter_interval() function runs an iteration for a specified interval.

The dk4iter_start_point() function runs an iteration if just one
starting x value is specified.

If no context is specified the following defaults are used:

Option | Default
:----: | :------
Algorithm | DK4_ITER_ALG_RF_ANDERSON_BJOERCK for dk4iter_interval(), DK4_ITER_ALG_NEWTON for dk4iter_start_point().
Max passes | 256.
Y tolerance | 1.0e-8
X tolerance | 1.0e-8
Full machine precision | No.
Restricted x interval | None.

*/



/**	Border values for maximum number of iteration passes.
*/
enum {
								/**	Number of passes for average cases.
								*/
	DK4_ITER_PASSES_REGULAR	=	256 ,

								/**	Do not set a limit on the number of passes.
									Warning: Might result in an endless loop!
								*/
	DK4_ITER_PASSES_UNLIMITED	=	0
};



/**	Iteration algorithms.
*/
enum {
									/**	Interval bisection.
									*/
	DK4_ITER_ALG_BISECTION			=	0 ,

									/**	Primitive from of regula falsi.
									*/
	DK4_ITER_ALG_RF_PRIMITIVE ,

									/**	Regula falsi, Illinois variant.
									*/
	DK4_ITER_ALG_RF_ILLINOIS ,

									/**	Regula falsi, Pegasus variant.
									*/
	DK4_ITER_ALG_RF_PEGASUS ,

									/**	Regula falsi, Anderson-Bjoerck variant.
									*/
	DK4_ITER_ALG_RF_ANDERSON_BJOERCK ,

									/**	Newton algorithm.
										The function must place 2 elements
										in the array specified by pointer:
										Function value and derivative value.
									*/
	DK4_ITER_ALG_NEWTON ,

									/**	Fix point algorithm.
										The function must be the phi(x)
										function from the
										phi(x) = x equation.
									*/
	DK4_ITER_ALG_FIX_POINT
};



/**	Iteration result.
*/
enum {
								/**	Iteration succeeded.
								*/
	DK4_ITER_RESULT_SUCCESS		=	1 ,

								/**	Error: Too many passes without success.
								*/
	DK4_ITER_RESULT_E_PASSES	=	0 ,

								/**	Error: Infinite value or NaN in calculation.
								*/
	DK4_ITER_RESULT_E_INFINITE	=	-1 ,

								/**	Error: X left initial interval.
									Can happen with Newton and fix point
									algorithm.
								*/
	DK4_ITER_RESULT_E_OOR		=	-2 ,

								/**	Error: Not converging.
									Fix point algorithm: Interval was enlarging.
								*/
	DK4_ITER_RESULT_E_CONV		=	-3 ,

								/**	Function calculation failed.
								*/
	DK4_ITER_RESULT_E_FCT		=	-4 ,

								/**	Invalid arguments passed to function call.
								*/
	DK4_ITER_RESULT_E_ARGS		=	-5 
};



/**	Function to iterate.
	@param	d	Destination address. One element for all algorithms
	except DK4_ITER_ALG_NEWTON which saves function value and derivative
	value in a 2 elements array.
	@param	x	X position to calculate the function value for.
	@param	ps	Parameter set, may be NULL.
	@return	Non-zero value on success, 0 on error.
*/
typedef int dk4_iter_fct_t(double *d, double x, const void *ps);


/**	Iteration context.
*/
typedef struct {
	double			xmin;		/**< Minimum of x interval. */
	double			xmax;		/**< Maximum of x interval. */
	double			eps_x;		/**< Epsilon for x change. */
	double			eps_y;		/**< Epsilon for absolute y value. */
	unsigned long	maxpass;	/**< Maximum number of passes. */
	int				exact;		/**< Flag: End only if there is no x change. */
	int				algo;		/**< Iteration algorithm to use. */
	int				minmax;		/**< Flags: min(1), max(2) set. */
} dk4_iter_ctx_t;



#ifdef	__cplusplus
extern "C" {
#endif

/**	Open new iteration context.
	Allocate memory for new iteration context, set up default value.
	A context allocated by this function must be released after
	use by the dk4iter_ctx_close() function.
	@return	Valid pointer to new context on success, NULL on error.
*/

dk4_iter_ctx_t *
dk4iter_ctx_open(void);


/**	Initialize a context to default values.
	@param	ctx	Context to initialize.
*/

void
dk4iter_ctx_init(dk4_iter_ctx_t *ctx);


/**	Close iteration context.
	Release memory assigned to the context.
	@param	ctx	Context created by dk4iter_ctx_open().
*/

void
dk4iter_ctx_close(dk4_iter_ctx_t *ctx);


/**	Set epsilon value for x change.
	@param	ctx	Context to modify.
	@param	eps	Epsilon value for x change.
*/

void
dk4iter_ctx_set_eps_x(dk4_iter_ctx_t *ctx, double eps);


/**	Set maximum absolute y value.
	@param	ctx	Context to modify.
	@param	eps	Maximum absolute y value.
*/

void
dk4iter_ctx_set_eps_y(dk4_iter_ctx_t *ctx, double eps);


/**	Set maximum number of passes (iteration steps).
	@param	ctx	Context to modify.
	@param	passes	Maximum number of passes, 0 or negative values indicate
	an unlimited number of steps.
*/

void
dk4iter_ctx_set_maxpass(dk4_iter_ctx_t *ctx, unsigned long passes);


/**	Set flag for exact iteration.
	If the flag is activated, iteration is continued until there is
	no longer any change in the x value. For some functions you can
	retrieve a result in machine precision using this flag, for other
	functions the iteration may fail.
	@param	ctx	Context to modify.
	@param	flag	New flag value, 0=inactive, other=active.
	The recommended value is 0.
*/

void
dk4iter_ctx_set_exact(dk4_iter_ctx_t *ctx, int flag);


/**	Set up iteration algorithm.
	@param	ctx			Context to modify.
	@param	algorithm	Algorithm to use.
*/

void
dk4iter_ctx_set_algorithm(dk4_iter_ctx_t *ctx, int algorithm);


/**	Set interval minimum.
	The allowed interval is only used for Newton and fixed point
	algorithm.
	@param	ctx		Context to modify.
	@param	xmin	Minimum x value allowed.
*/

void
dk4iter_ctx_set_min(dk4_iter_ctx_t *ctx, double xmin);


/**	Set exlusive interval minimum (open interval border).
	The allowed interval is only used for Newton and fixed point
	algorithm.
	@param	ctx		Context to modify.
	@param	xmin	Minimum x value allowed.
*/

void
dk4iter_ctx_set_exclusive_min(dk4_iter_ctx_t *ctx, double xmin);


/**	Set interval maximum.
	The allowed interval is only used for Newton and fixed point
	algorithm.
	@param	ctx		Context to modify.
	@param	xmax	Maximum x value allowed.
*/

void
dk4iter_ctx_set_max(dk4_iter_ctx_t *ctx, double xmax);


/**	Set exclusive interval maximum (open interval border).
	The allowed interval is only used for Newton and fixed point
	algorithm.
	@param	ctx		Context to modify.
	@param	xmax	Maximum x value allowed.
*/

void
dk4iter_ctx_set_exclusive_max(dk4_iter_ctx_t *ctx, double xmax);


/**	Run iteration.
	The algorithm in the context must be one from:
	DK4_ITER_ALG_BISECTION, DK4_ITER_ALG_RF_PRIMITIVE, DK4_ITER_ALG_RF_ILLINOIS,
	DK4_ITER_ALG_RF_PEGASUS	DK4_ITER_ALG_RF_ANDERSON_BJOERCK.
	Without a context DK4_ITER_ALG_RF_ANDERSON_BJOERCK is used as default.
	@param	d	Address of result variable.
	@param	pp	Address of variable to store number of passes, may be NULL.
	@param	fct	Function to find root for.
	@param	ps	Parameter set for function, may be NULL.
	@param	a	One interval border.
	@param	b	Other interval border.
	@param	ctx	Iteration context.
	@return	DK4_ITER_RESULT_SUCCESS on success, one from
	DK4_ITER_RESULT_E_PASSES, DK4_ITER_RESULT_E_INFINITE,
	DK4_ITER_RESULT_E_OOR, DK4_ITER_RESULT_E_CONV or
	DK4_ITER_RESULT_E_FCT on error.
*/

int
dk4iter_interval(
	double					*d,
	unsigned long			*pp,
	dk4_iter_fct_t			*fct,
	void			const	*ps,
	double					 a,
	double					 b,
	dk4_iter_ctx_t	const	*ctx
);


/**	Run iteration.
	The algorithm in the context must be one from:
	DK4_ITER_ALG_NEWTON, DK4_ITER_ALG_FIX_POINT.
	Without a context, DK4_ITER_ALG_NEWTON is used as default.
	@param	d	Address of result variable.
	@param	pp	Address of variable to store number of passes, may be NULL.
	@param	fct	Function to find root for.
				For Newton iteration this function must set two values
				in the destination array: function value and derivative value.
				For fixpoint iteration the function calculates the
				phi(x) part of phi(x)=x.
	@param	ps	Parameter set for function, may be NULL.
	@param	x0	Start point.
	@param	ctx	Iteration context.
	@return	DK4_ITER_RESULT_SUCCESS on success, one from
	DK4_ITER_RESULT_E_PASSES, DK4_ITER_RESULT_E_INFINITE,
	DK4_ITER_RESULT_E_OOR, DK4_ITER_RESULT_E_CONV or
	DK4_ITER_RESULT_E_FCT on error.
*/

int
dk4iter_start_point(
	double					*d,
	unsigned long			*pp,
	dk4_iter_fct_t			*fct,
	void			const	*ps,
	double					 x0,
	dk4_iter_ctx_t	const	*ctx
);


#ifdef	__cplusplus
}
#endif

/* vim: set ai sw=4 ts=4 : */


#endif