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