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