1/** @file  test-iter.c  Testing iteration functions.
2
3By default - if SERIALCIRCUIT is defined to 1 - the program calculates
4the current through a serial connection of a resistor and a diode if the
5summary voltage, the resistor value and the diodes saturation current
6are specified.
7
8As command line arguments specify
9- summary voltage in Volts,
10- resistor value in Ohms and
11- Saturation current in Ampere (usually in range 1.0e-18 to 1.0e-12).
12
13If the values are not specified on the command line, the program prompts
14for values.
15
16When defining SERIALCIRCUIT to 0 you can do tests with mathematical
17equations   e^x=x^2-2 (if MATH_FCT1 defined to 1) or 0=sqrt(x)-x
18(if MATH_FCT1 defined to 0).
19
20You can play with the eps variable which is used for both x and y
21tolerance here.
22
23Change full to 1 to calculate to the full machine precision
24(does not necessarily work for all equations and all algorithms).
25*/
26
27#include <stdio.h>
28#include <stdlib.h>
29#include <math.h>
30
31#include <libdk4c/dk4iter.h>
32
33
34/** Switch between an electronic circuit (resistor and diode in serial)
35    or mathematical function.
36*/
37#define SERIALCIRCUIT 0
38
39/** For mathematical function switch between two versions.
40*/
41#define MATH_FCT1       1
42
43/** Attempt to iterate to full machine precision.
44*/
45#define FULL_MACHINE_PRECISION  1
46
47
48
49/** Epsilon.
50*/
51static const double eps     = 1.0e-8;
52
53
54
55/** Flag: Attempt to iterate to full machine precision.
56*/
57static const int    full    =   FULL_MACHINE_PRECISION;
58
59
60
61/** Variables for the serial circuit.
62*/
63static	double	i0	=	0.0;
64static	double	uges	=	0.0;
65static	double	r	=	0.0;
66
67
68
69/** This is a constant. 26mV.
70*/
71static const double UT      = 0.026;
72
73
74/** Exit status code.
75*/
76static int exval = EXIT_FAILURE;
77
78
79
80static
81int
82f(double *dptr, double x, const void *pset)
83{
84#if SERIALCIRCUIT
85    *dptr = r * x + UT * log(1.0 + x / i0) - uges;
86    return 1;
87#else
88#if MATH_FCT1
89    *dptr = exp(x) + x * x - 2.0;
90    return 1;
91#else
92    *dptr = sqrt(x) - x;
93    return 1;
94#endif
95#endif
96    return 1;
97}
98
99
100
101static
102int
103f_dfdx(double *dptr, double x, const void *pset)
104{
105#if SERIALCIRCUIT
106    *(dptr++)   = r * x + UT * log(1.0 + x / i0) - uges;
107    *dptr       = r + UT / (1.0 + x / i0);
108    return 1;
109#else
110#if MATH_FCT1
111    *(dptr++)   = exp(x) + x * x - 2.0;
112    *dptr       = exp(x) + 2.0 * x;
113    return 1;
114#else
115    *(dptr++)   = sqrt(x) - x;
116    *dptr       = 1.0/(2.0*sqrt(x)) - 1.0;
117    return 1;
118#endif
119#endif
120}
121
122
123
124static
125int
126f_fix(double *dptr, double x, const void *pset)
127{
128#if SERIALCIRCUIT
129    double  logarg;
130    int     back        = 0;
131
132    logarg = (i0 + x) / i0;
133    if (0.0 < logarg) {
134        *dptr = -1.0 * ((UT * log(logarg) - uges) / r);
135        back = 1;
136    }
137    return back;
138#else
139#if MATH_FCT1
140    double  logarg;
141    int     back        = 0;
142
143    logarg = 2.0 - x * x;
144    if (0.0 < logarg) {
145        *dptr = log(logarg);
146        back = 1;
147    }
148    return back;
149#else
150    *dptr = sqrt(x);
151    return 1;
152#endif
153#endif
154}
155
156
157
158static
159void
160show_results(const char *verfname, int res, double x, unsigned long pa)
161{
162    fputs(verfname, stdout);
163    fputc(' ', stdout);
164    switch (res) {
165        case DK4_ITER_RESULT_SUCCESS : {
166            printf("Success   x = %lg   (%lu steps)\n", x, pa);
167        } break;
168        case DK4_ITER_RESULT_E_PASSES : {
169            fputs("ERROR: Too many iteration steps!\n", stdout);
170        } break;
171        case DK4_ITER_RESULT_E_INFINITE : {
172            fputs("ERROR: Infinite/undefined calculation result!\n", stdout);
173        } break;
174        case DK4_ITER_RESULT_E_OOR : {
175            fputs("ERROR: Leaving allowed interval!\n", stdout);
176        } break;
177        case DK4_ITER_RESULT_E_CONV : {
178            fputs("ERROR: No convergence!\n", stdout);
179        } break;
180        case DK4_ITER_RESULT_E_FCT : {
181            fputs("ERROR: Failed to evaluate function!\n", stdout);
182        } break;
183        case DK4_ITER_RESULT_E_ARGS : {
184            fputs("ERROR: Unusable arguments!\n", stdout);
185        } break;
186    }
187}
188
189
190
191static
192void
193calculation(void)
194{
195    dk4_iter_ctx_t  ctx;
196    double          a   =   0.0;
197    double          b;
198    double          x   =   0.0;
199    unsigned long   pa  =   0UL;        /* Number of passes */
200    int             res =   0;          /* Operation result */
201
202#if SERIALCIRCUIT
203    b = uges / r;
204#else
205#if MATH_FCT1
206    b = 0.7;
207#else
208    a = 0.8;
209    b = 1.2;
210#endif
211#endif
212
213    /*  Halbierung
214    */
215    dk4iter_ctx_init(&ctx);
216    dk4iter_ctx_set_eps_x(&ctx, eps);
217    dk4iter_ctx_set_eps_y(&ctx, eps);
218    dk4iter_ctx_set_maxpass(&ctx, DK4_ITER_PASSES_REGULAR);
219    dk4iter_ctx_set_exact(&ctx, full);
220    dk4iter_ctx_set_algorithm(&ctx, DK4_ITER_ALG_BISECTION);
221    res = dk4iter_interval(&x, &pa, f, NULL, a, b, &ctx);
222    show_results("Bisection       ", res, x, pa);
223
224    dk4iter_ctx_init(&ctx);
225    dk4iter_ctx_set_eps_x(&ctx, eps);
226    dk4iter_ctx_set_eps_y(&ctx, eps);
227    dk4iter_ctx_set_maxpass(&ctx, DK4_ITER_PASSES_REGULAR);
228    dk4iter_ctx_set_exact(&ctx, full);
229    dk4iter_ctx_set_algorithm(&ctx, DK4_ITER_ALG_RF_PRIMITIVE);
230    res = dk4iter_interval(&x, &pa, f, NULL, a, b, &ctx);
231    show_results("Regula falsi    ", res, x, pa);
232
233    dk4iter_ctx_init(&ctx);
234    dk4iter_ctx_set_eps_x(&ctx, eps);
235    dk4iter_ctx_set_eps_y(&ctx, eps);
236    dk4iter_ctx_set_maxpass(&ctx, DK4_ITER_PASSES_REGULAR);
237    dk4iter_ctx_set_exact(&ctx, full);
238    dk4iter_ctx_set_algorithm(&ctx, DK4_ITER_ALG_RF_ILLINOIS);
239    res = dk4iter_interval(&x, &pa, f, NULL, a, b, &ctx);
240    show_results("Illinois        ", res, x, pa);
241
242    dk4iter_ctx_init(&ctx);
243    dk4iter_ctx_set_eps_x(&ctx, eps);
244    dk4iter_ctx_set_eps_y(&ctx, eps);
245    dk4iter_ctx_set_maxpass(&ctx, DK4_ITER_PASSES_REGULAR);
246    dk4iter_ctx_set_exact(&ctx, full);
247    dk4iter_ctx_set_algorithm(&ctx, DK4_ITER_ALG_RF_PEGASUS);
248    res = dk4iter_interval(&x, &pa, f, NULL, a, b, &ctx);
249    show_results("Pegasus         ", res, x, pa);
250
251    dk4iter_ctx_init(&ctx);
252    dk4iter_ctx_set_eps_x(&ctx, eps);
253    dk4iter_ctx_set_eps_y(&ctx, eps);
254    dk4iter_ctx_set_maxpass(&ctx, DK4_ITER_PASSES_REGULAR);
255    dk4iter_ctx_set_exact(&ctx, full);
256    dk4iter_ctx_set_algorithm(&ctx, DK4_ITER_ALG_RF_ANDERSON_BJOERCK);
257    res = dk4iter_interval(&x, &pa, f, NULL, a, b, &ctx);
258    show_results("Anderson-Bjoerck", res, x, pa);
259
260    dk4iter_ctx_init(&ctx);
261    dk4iter_ctx_set_eps_x(&ctx, eps);
262    dk4iter_ctx_set_eps_y(&ctx, eps);
263    dk4iter_ctx_set_maxpass(&ctx, DK4_ITER_PASSES_REGULAR);
264    dk4iter_ctx_set_exact(&ctx, full);
265    dk4iter_ctx_set_algorithm(&ctx, DK4_ITER_ALG_NEWTON);
266    res = dk4iter_start_point(&x, &pa, f_dfdx, NULL, b, &ctx);
267    show_results("Newton          ", res, x, pa);
268
269    dk4iter_ctx_init(&ctx);
270    dk4iter_ctx_set_eps_x(&ctx, eps);
271    dk4iter_ctx_set_eps_y(&ctx, eps);
272    dk4iter_ctx_set_maxpass(&ctx, DK4_ITER_PASSES_REGULAR);
273    dk4iter_ctx_set_exact(&ctx, full);
274    dk4iter_ctx_set_algorithm(&ctx, DK4_ITER_ALG_FIX_POINT);
275    res = dk4iter_start_point(&x, &pa, f_fix, NULL, b, &ctx);
276    show_results("Fixpoint        ", res, x, pa);
277
278}
279
280
281
282static
283int
284read_a_double(const char *prompt, double *da)
285{
286    char    buf[256];
287    int     back    = 0;
288    fputs(prompt, stdout); fputc('\n', stdout);
289    if (NULL != fgets(buf, sizeof(buf), stdin)) {
290        if (0 != sscanf(buf, "%lg", da)) {
291            back = 1;
292        }
293        else {
294            fputs("ERROR: Not a floating point number: ", stderr);
295            fputs(buf, stderr);
296            fputc('\n', stderr);
297        }
298    }
299    else {
300        fputs("ERROR: Failed to read input!\n", stderr);
301    }
302    return back;
303}
304
305
306
307int main(int argc, char *argv[])
308{
309    if (4 == argc) {
310        if (0 != sscanf(argv[1], "%lg", &uges)) {
311            if (0 != sscanf(argv[2], "%lg", &r)) {
312                if (0 != sscanf(argv[3], "%lg", &i0)) {
313                    exval = EXIT_SUCCESS;
314                    calculation();
315                }
316            }
317        }
318    }
319    else {
320#if SERIALCIRCUIT
321	    if (0 != read_a_double("Uges in V: ", &uges)) {
322		    if (0 != read_a_double("R in Ohm: ", &r)) {
323			    if (0 != read_a_double("Is in A: ", &i0)) {
324				    if ((0.0 < uges) && (0.0 < i0) && (0.0 < r)) {
325					    exval = EXIT_SUCCESS;
326					    calculation();
327				    }
328			    }
329		    }
330	    }
331#else
332    exval = EXIT_SUCCESS;
333    calculation();
334#endif
335    }
336	return exval;
337}
338
339/* vim: set ai sw=4 ts=4 expandtab : */
340