%% options copyright owner = Dirk Krause copyright year = 2018-xxxx SPDX-License-Identifier: BSD-3-Clause %% header /** @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
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
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 : */ %% module #include "dk4conf.h" #include #if DK4_HAVE_ASSERT_H #ifndef ASSERT_H_INCLUDED #include #define ASSERT_H_INCLUDED 1 #endif #endif #if DK4_HAVE_STDLIB_H #ifndef STDLIB_H_INCLUDED #include #define STDLIB_H_INCLUDED 1 #endif #endif #if DK4_HAVE_LIMITS_H #ifndef LIMITS_H_INCLUDED #include #define LIMITS_H_INCLUDED 1 #endif #endif #if DK4_HAVE_STDINT_H #ifndef STDINT_H_INCLUDED #include #define STDINT_H_INCLUDED 1 #endif #endif #include "dk4mem.h" #include "dk4iter.h" #include "dk4math.h" $!trace-include /** Initial function value condition. */ enum { CONDITION_ILLEGAL = 0 , /** fa < 0, fb >= 0. */ CONDITION_FA_LESS_ZERO , /** fa <= 0, fb > 0. */ CONDITION_FA_LEQ_ZERO , /** fa > 0, fb <= 0. */ CONDITION_FA_GREATER_ZERO , /** fa >= 0, fb < 0. */ CONDITION_FA_GEQ_ZERO }; /** Note which interval border was changed in previous step. */ enum { /** No interval border change yet. */ CHANGED_NONE = 0, /** Previous step changed interval border a. */ CHANGED_A = -1, /** Previous step changed interval border b. */ CHANGED_B = 1 }; /** Flags which restrictions apply to x values. */ enum { /** Minimum specified. */ MINMAX_MINIMUM = 0x0001 , /** Maximum specified. */ MINMAX_MAXIMUM = 0x0002 , /** Specified minimum is exclusive. */ MINMAX_MIN_EXCL = 0x0004 , /** Specified maximum is exclusive. */ MINMAX_MAX_EXCL = 0x0008 }; void dk4iter_ctx_init(dk4_iter_ctx_t *ctx) { #if DK4_USE_ASSERT assert(NULL != ctx); #endif if (NULL != ctx) { DK4_MEMRES(ctx,sizeof(dk4_iter_ctx_t)); ctx->xmin = 0.0; ctx->xmax = 0.0; ctx->eps_x = 1.0e-8; ctx->eps_y = 1.0e-8; ctx->maxpass = (unsigned long)(DK4_ITER_PASSES_REGULAR); ctx->exact = 0; ctx->algo = DK4_ITER_ALG_RF_ANDERSON_BJOERCK; ctx->minmax = 0; } } dk4_iter_ctx_t * dk4iter_ctx_open(void) { dk4_iter_ctx_t *back = NULL; back = dk4mem_new(dk4_iter_ctx_t,1,NULL); if (NULL != back) { dk4iter_ctx_init(back); } return back; } void dk4iter_ctx_close(dk4_iter_ctx_t *ctx) { #if DK4_USE_ASSERT assert(NULL != ctx); #endif if (NULL != ctx) { dk4mem_free(ctx); } } void dk4iter_ctx_set_eps_x(dk4_iter_ctx_t *ctx, double eps) { #if DK4_USE_ASSERT assert(NULL != ctx); #endif if (NULL != ctx) { ctx->eps_x = eps; } } void dk4iter_ctx_set_eps_y(dk4_iter_ctx_t *ctx, double eps) { #if DK4_USE_ASSERT assert(NULL != ctx); #endif if (NULL != ctx) { ctx->eps_y = eps; } } void dk4iter_ctx_set_maxpass(dk4_iter_ctx_t *ctx, unsigned long passes) { #if DK4_USE_ASSERT assert(NULL != ctx); #endif if (NULL != ctx) { ctx->maxpass = passes; } } void dk4iter_ctx_set_exact(dk4_iter_ctx_t *ctx, int flag) { #if DK4_USE_ASSERT assert(NULL != ctx); #endif if (NULL != ctx) { ctx->exact = flag; } } void dk4iter_ctx_set_algorithm(dk4_iter_ctx_t *ctx, int algorithm) { #if DK4_USE_ASSERT assert(NULL != ctx); #endif if (NULL != ctx) { ctx->algo = algorithm; } } void dk4iter_ctx_set_min(dk4_iter_ctx_t *ctx, double xmin) { #if DK4_USE_ASSERT assert(NULL != ctx); #endif if (NULL != ctx) { ctx->xmin = xmin; ctx->minmax |= MINMAX_MINIMUM; } } void dk4iter_ctx_set_max(dk4_iter_ctx_t *ctx, double xmax) { #if DK4_USE_ASSERT assert(NULL != ctx); #endif if (NULL != ctx) { ctx->xmax = xmax; ctx->minmax |= MINMAX_MAXIMUM; } } void dk4iter_ctx_set_exclusive_min(dk4_iter_ctx_t *ctx, double xmin) { #if DK4_USE_ASSERT assert(NULL != ctx); #endif if (NULL != ctx) { ctx->xmin = xmin; ctx->minmax |= (MINMAX_MINIMUM | MINMAX_MIN_EXCL); } } void dk4iter_ctx_set_exclusive_max(dk4_iter_ctx_t *ctx, double xmax) { #if DK4_USE_ASSERT assert(NULL != ctx); #endif if (NULL != ctx) { ctx->xmax = xmax; ctx->minmax |= (MINMAX_MAXIMUM | MINMAX_MAX_EXCL); } } static int find_condition(double fa, double fb) { int back = CONDITION_ILLEGAL; if ((0.0 > fa) && (0.0 <= fb)) { back = CONDITION_FA_LESS_ZERO; } else { if ((0.0 < fa) && (0.0 >= fb)) { back = CONDITION_FA_GREATER_ZERO; } else { if ((0.0 >= fa) && (0.0 < fb)) { back = CONDITION_FA_LEQ_ZERO; } else { if ((0.0 <= fa) && (0.0 > fb)) { back = CONDITION_FA_GEQ_ZERO; } } } } return back; } /** Calculate gamma value for PEGASUS variant. @param fg Function value at border last set. @param fx Function value at recent x position. @return Gamma value. */ static double dk4iter_gamma_pegasus(double fg, double fx) { double back; back = fg / (fg + fx); if (0 == dk4ma_is_finite(back)) { back = 0.5; } else { if (0.0 >= back) { back = 0.5; } #if 0 /* 2018-07-06 Can not happen, fg and fx have same sign. */ else { if (1.0 < back) { back = 1.0; } } #endif } return back; } /** Calculate gamma value for ANDERSON-BJOERCK variant. @param fg Function value at border last set. @param fx Function value at recent x position. @return Gamma value. */ static double dk4iter_gamma_anderson_bjoerck(double fg, double fx) { double back; back = 1.0 - fx / fg; if (0 == dk4ma_is_finite(back)) { back = 0.5; } else { if (0.0 >= back) { back = 0.5; } #if 0 /* 2018-07-06 Can not happen as fx and fy have same sign. */ else { if (1.0 < back) { back = 1.0; } } #endif } return back; } 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 ) { dk4_iter_ctx_t mctx; /* Copy of context */ double fa = 0.0; /* Function value for border a */ double fb = 0.0; /* Function value for border b */ double x = 0.0; /* Text x value */ double fx = 0.0; /* Function value for x */ double afx = 0.0; /* Absolute function value for x */ double xo = 0.0; /* Previous step x value */ double gamma = 0.5; /* Correction factor */ unsigned long passno = 0UL; /* Number of current pass */ int res = 0; /* Operation result */ int cond = 0; /* Condition of interval borders */ int cc = 0; /* 1=continue, 0=finished, -1=abort */ int pc = CHANGED_NONE; /* Previous border change */ int nc = CHANGED_NONE; /* Next change */ int back = DK4_ITER_RESULT_E_ARGS; $? "+ dk4iter_interval" /* Check function call arguments */ #if DK4_USE_ASSERT assert(NULL != d); assert(NULL != fct); #endif if ((NULL == d) || (NULL == fct)) { $? "! d or fct" goto finished; } if (a == b) { $? "! zero length interval" goto finished; } /* Copy or initialize context */ if (NULL != ctx) { $? ". use context" DK4_MEMCPY(&mctx,ctx,sizeof(dk4_iter_ctx_t)); } else { $? ". no context, use defaults" dk4iter_ctx_init(&mctx); } /* Check algorithm specification */ if (DK4_ITER_ALG_BISECTION != mctx.algo) { if (DK4_ITER_ALG_RF_PRIMITIVE != mctx.algo) { if (DK4_ITER_ALG_RF_ILLINOIS != mctx.algo) { if (DK4_ITER_ALG_RF_PEGASUS != mctx.algo) { if (DK4_ITER_ALG_RF_ANDERSON_BJOERCK != mctx.algo) { goto finished; $? "! illegal algorithm" } } } } } /* Calculate borders at beginning, check */ res = (*fct)(&fa, a, ps); if ((0 == res) || (!(dk4ma_is_finite(fa)))) { back = DK4_ITER_RESULT_E_FCT; $? "! function calculation for a" goto finished; } res = (*fct)(&fb, b, ps); if ((0 == res) || (!(dk4ma_is_finite(fb)))) { back = DK4_ITER_RESULT_E_FCT; $? "! function calculation for b" goto finished; } cond = find_condition(fa, fb); if (CONDITION_ILLEGAL == cond) { $? "! invalid interval" goto finished; } /* Prepare iteration loop */ passno = 0UL; x = xo = 0.0; cc = 1; pc = CHANGED_NONE; /* Run iteration loop */ do { /* Keep previous x position */ xo = x; /* Increase pass number, but avoid wrapping */ if (ULONG_MAX > passno) { passno++; } $? ". pass %lu", passno /* Calculate x position */ $? ". a = %lg fa = %lg b = %lg fb = %lg", a, fa, b, fb switch (mctx.algo) { case DK4_ITER_ALG_RF_PRIMITIVE : case DK4_ITER_ALG_RF_ILLINOIS : case DK4_ITER_ALG_RF_PEGASUS : case DK4_ITER_ALG_RF_ANDERSON_BJOERCK : { /* Regula falsi */ x = (a * fb - b * fa) / (fb - fa); } break; default : { /* Bisection */ x = (a + b) / 2.0; } break; } if (0 == dk4ma_is_finite(x)) { back = DK4_ITER_RESULT_E_INFINITE; cc = -1; } else { $? ". x = %lg", x /* Calculate function for x, if x is usable */ res = (*fct)(&fx, x, ps); if (0 == res) { $? "! function calculation for x" back = DK4_ITER_RESULT_E_FCT; cc = -1; } if (0 == dk4ma_is_finite(fx)) { $? "! function calculation for x" back = DK4_ITER_RESULT_E_INFINITE; cc = -1; } afx = fabs(fx); if (0 == dk4ma_is_finite(afx)) { $? "! absolute value for fx" back = DK4_ITER_RESULT_E_INFINITE; cc = -1; } /* Check whether we are done, if we could continue */ if (1 == cc) { $? ". y = %lg", fx $!trace-code if (1UL < passno) { $? ". x-delta = %lg", fabs(x - xo) $!trace-code } if (isgreater(mctx.eps_y,0.0)) { /* Must check y */ if (afx < mctx.eps_y) { $? ". y small enough" if ((x == xo) && (1UL < passno)) { cc = 0; $? ". OK exact match" } else { if (0 == mctx.exact) { if(isgreater(mctx.eps_x,0.0)) { if (isless(fabs(x - xo),mctx.eps_x)) { if (1UL < passno) { cc = 0; $? ". OK small step" } } } else { cc = 0; $? ". OK no x step restriction" } } } } } else { $? ". no y check required" /* No y check */ if ((x == xo) && (1UL < passno)) { cc = 0; $? ". OK exact match" } else { if (0 == mctx.exact) { if(isgreater(mctx.eps_x,0.0)) { if (isless(fabs(x - xo),mctx.eps_x)) { if (1UL < passno) { cc = 0; $? ". OK x step small enough" } } } else { /* Neither eps_x nor eps_y, so we wait for x=xo. */ } } } } } /* Continue for usable y values only */ if (1 == cc) { /* Find direction for next change */ nc = CHANGED_NONE; switch (cond) { case CONDITION_FA_LESS_ZERO : { if (0.0 > fx) { nc = CHANGED_A; } else { nc = CHANGED_B; } } break; case CONDITION_FA_LEQ_ZERO : { if (0.0 >= fx) { nc = CHANGED_A; } else { nc = CHANGED_B; } } break; case CONDITION_FA_GREATER_ZERO : { if (0.0 < fx) { nc = CHANGED_A; } else { nc = CHANGED_B; } } break; case CONDITION_FA_GEQ_ZERO : { if (0.0 <= fx) { nc = CHANGED_A; } else { nc = CHANGED_B; } } break; } /* Apply interval border change */ switch (nc) { case CHANGED_A: { $? ". use x as a" if (CHANGED_A == pc) { switch (mctx.algo) { case DK4_ITER_ALG_RF_ILLINOIS : { fb = 0.5 * fb; if (0 == dk4ma_is_finite(fb)) { back = DK4_ITER_RESULT_E_INFINITE; cc = -1; $? "! fb" } } break; case DK4_ITER_ALG_RF_PEGASUS : { gamma = dk4iter_gamma_pegasus(fa, fx); fb *= gamma; if (0 == dk4ma_is_finite(fb)) { back = DK4_ITER_RESULT_E_INFINITE; cc = -1; $? "! fb" } } break; case DK4_ITER_ALG_RF_ANDERSON_BJOERCK : { gamma = dk4iter_gamma_anderson_bjoerck( fa, fx ); fb *= gamma; if (0 == dk4ma_is_finite(fb)) { back = DK4_ITER_RESULT_E_INFINITE; cc = -1; $? "! fb" } } break; } } a = x; fa = fx; pc = CHANGED_A; } break; case CHANGED_B: { $? ". use x as b" if (CHANGED_B == pc) { switch (mctx.algo) { case DK4_ITER_ALG_RF_ILLINOIS : { fa = 0.5 * fa; if (0 == dk4ma_is_finite(fa)) { back = DK4_ITER_RESULT_E_INFINITE; cc = -1; $? "! fa" } } break; case DK4_ITER_ALG_RF_PEGASUS : { gamma = dk4iter_gamma_pegasus( fb, fx ); fa *= gamma; if (0 == dk4ma_is_finite(fa)) { back = DK4_ITER_RESULT_E_INFINITE; cc = -1; $? "! fa" } } break; case DK4_ITER_ALG_RF_ANDERSON_BJOERCK : { gamma = dk4iter_gamma_anderson_bjoerck( fb, fx ); fa *= gamma; if (0 == dk4ma_is_finite(fa)) { back = DK4_ITER_RESULT_E_INFINITE; cc = -1; $? "! fa" } } break; } } b = x; fb = fx; pc = CHANGED_B; } break; default: { /* ERROR: Must not happen */ back = DK4_ITER_RESULT_E_CONV; $? "! bug" cc = -1; } break; } } } /* Check number of passes */ if ((1 == cc) && (0UL < mctx.maxpass) && (passno >= mctx.maxpass)) { back = DK4_ITER_RESULT_E_PASSES; $? "! too many passes" cc = -1; } } while (1 == cc); /* Success */ if (0 == cc) { *d = x; if (NULL != pp) { *pp = passno; } back = DK4_ITER_RESULT_SUCCESS; } finished: $? "- dk4iter_interval %d", back return back; } /** Run Newton iteration. @param d Destination (address of result variable). @param pp Address of variable to store number of passes on success. @param fct Iteration function, returns two values into the array at address d: the function value and the first derivative value. @param ps Parameter set, may be NULL if fct does not use it. @param x0 Start point. @param ctx Iteration context, may be NULL. @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, DK4_ITER_RESULT_E_FCT, or DK4_ITER_RESULT_E_ARGS on error. */ static int dk4iter_newton( double *d, unsigned long *pp, dk4_iter_fct_t *fct, void const *ps, double x0, dk4_iter_ctx_t const *ctx ) { double v[2]; double xn = 0.0; unsigned long passno = 0UL; int back = DK4_ITER_RESULT_E_ARGS; int cc = 1; int res = 0; $? "+ dk4iter_newton" #if DK4_USE_ASSERT assert(NULL != d); assert(NULL != fct); #endif res = (*fct)(v, x0, ps); if (0 == res) { cc = -1; $? "! function calculation" back = DK4_ITER_RESULT_E_FCT; } if (!((dk4ma_is_finite(v[0])) && (dk4ma_is_finite(v[1])))) { cc = -1; $? "! infinite values" back = DK4_ITER_RESULT_E_INFINITE; } while (1 == cc) { if (ULONG_MAX > passno) { passno++; } $? ". begin pass %lu", passno xn = x0 - v[0] / v[1]; if (dk4ma_is_finite(xn)) { $? ". x = %lg", xn if (0 != (MINMAX_MINIMUM & (ctx->minmax))) { if (0 != (MINMAX_MIN_EXCL & (ctx->minmax))) { if (xn <= ctx->xmin) { cc = -1; $? "! out of range, too small" back = DK4_ITER_RESULT_E_OOR; } } else { if (xn < ctx->xmin) { cc = -1; $? "! out of range, too small" back = DK4_ITER_RESULT_E_OOR; } } } if (0 != (MINMAX_MAXIMUM & (ctx->minmax))) { if (0 != (MINMAX_MAX_EXCL & (ctx->minmax))) { if (xn >= ctx->xmax) { cc = -1; $? "! out of range, too large" back = DK4_ITER_RESULT_E_OOR; } } else { if (xn > ctx->xmax) { cc = -1; $? "! out of range, too large" back = DK4_ITER_RESULT_E_OOR; } } } if (1 == cc) { res = (*fct)(v, xn, ps); if (0 == res) { cc = -1; $? "! function calculation" back = DK4_ITER_RESULT_E_FCT; } if (!((dk4ma_is_finite(v[0])) && (dk4ma_is_finite(v[1])))) { cc = -1; $? "! infinite values" back = DK4_ITER_RESULT_E_INFINITE; } if (1 == cc) { $? ". y = %lg", v[0] $? ". dy/dx = %lg", v[1] if (isgreater(ctx->eps_y,0.0)) { $? ". check y" if (fabs(v[0]) < ctx->eps_y) { $? ". y in range" if (xn == x0) { cc = 0; $? ". OK exact value" } else { if (0 == ctx->exact) { $? ". check step" if(isgreater(ctx->eps_x,0.0)) { if (isless(fabs(xn - x0),ctx->eps_x)) { cc = 0; $? ". OK small step" } } else { cc = 0; $? ". OK no x step constraint" } } } } } else { $? ". no y check" if (xn == x0) { cc = 0; $? ". OK exact match" } else { if (0 == ctx->exact) { if (isgreater(ctx->eps_x,0.0)) { if (isless(fabs(xn - x0),ctx->eps_x)) { cc = 0; $? ". OK x step small enough" } } else { /* Neither eps_x nor eps_y, so we wait for x=x0. */ } } } } } x0 = xn; } } else { cc = -1; $? "! infinite value" back = DK4_ITER_RESULT_E_INFINITE; } if ((1 == cc) && (0UL < ctx->maxpass) && (passno >= ctx->maxpass)) { cc = -1; back = DK4_ITER_RESULT_E_PASSES; } } if (0 == cc) { *d = xn; if (NULL != pp) { *pp = passno; } back = DK4_ITER_RESULT_SUCCESS; } $? "- dk4iter_newton %d", back return back; } /** Run fix point iteration. @param d Destination (address of result variable). @param pp Address of variable to store number of passes on success. @param fct Iteration function, the phi part of phi(x)=x. @param ps Parameter set, may be NULL if fct does not use it. @param x0 Start point. @param ctx Iteration context, may be NULL. @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, DK4_ITER_RESULT_E_FCT, or DK4_ITER_RESULT_E_ARGS on error. */ static int dk4iter_fix_point( double *d, unsigned long *pp, dk4_iter_fct_t *fct, void const *ps, double x0, dk4_iter_ctx_t const *ctx ) { #if 0 double a; /* Border from previous steps */ double b; /* Border from previous steps */ #endif double xn = 0.0; /* New x value */ unsigned long passno = 0UL; /* Iteration step number */ int res = 0; /* Function evaluation result */ int cc = 1; /* Flag: Can continue */ int back = DK4_ITER_RESULT_E_ARGS; $? "+ dk4iter_fix_point" #if DK4_USE_ASSERT assert(NULL != d); assert(NULL != fct); #endif while (1 == cc) { if (ULONG_MAX > passno) { passno++; } $? ". begin pass %lu", passno /* Calculate new x, check result */ res = (*fct)(&xn, x0, ps); if (0 == res) { cc = -1; $? "! function evaluation" back = DK4_ITER_RESULT_E_FCT; } if (!(dk4ma_is_finite(xn))) { cc = -1; $? "! infinite result" back = DK4_ITER_RESULT_E_INFINITE; } if (1 == cc) { $? ". x = %lg", xn /* Check whether specified interval is exceeded */ if (0 != (MINMAX_MINIMUM & (ctx->minmax))) { if (0 != (MINMAX_MIN_EXCL & (ctx->minmax))) { if (xn <= ctx->xmin) { cc = -1; $? "! out of range, too small" back = DK4_ITER_RESULT_E_OOR; } } else { if (xn < ctx->xmin) { cc = -1; $? "! out of range, too small" back = DK4_ITER_RESULT_E_OOR; } } } if (0 != (MINMAX_MAXIMUM & (ctx->minmax))) { if (0 != (MINMAX_MAX_EXCL & (ctx->minmax))) { if (xn >= ctx->xmax) { cc = -1; $? "! out of range, too large" back = DK4_ITER_RESULT_E_OOR; } } else { if (xn > ctx->xmax) { cc = -1; $? "! out of range, too large" back = DK4_ITER_RESULT_E_OOR; } } } if (1 == cc) { if (1 == cc) { /* Check whether we are finished */ if (xn == x0) { cc = 0; $? ". OK exact match" } else { if ((0 == ctx->exact) && (0.0 < ctx->eps_x)) { if (isless(fabs(xn-x0),ctx->eps_x)) { cc = 0; $? ". OK change small enough" } #if TRACE_DEBUG else { $? ". change %lg", fabs(xn-x0) } #endif } else { $? "! No condition to check %lg", (xn-x0) } } } } x0 = xn; } /* Stop if too many passes */ if (1 == cc) { if ((0UL < ctx->maxpass) && (passno >= ctx->maxpass)) { cc = -1; back = DK4_ITER_RESULT_E_PASSES; } } } if (0 == cc) { *d = xn; if (NULL != pp) { *pp = passno; } back = DK4_ITER_RESULT_SUCCESS; } $? "- dk4iter_fix_point %d", back return back; } 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 ) { dk4_iter_ctx_t mctx; int back = DK4_ITER_RESULT_E_ARGS; /* Check function call arguments */ #if DK4_USE_ASSERT assert(NULL != d); assert(NULL != fct); #endif if ((NULL == d) || (NULL == fct)) { goto finished; } /* Copy or initialize context */ if (NULL != ctx) { DK4_MEMCPY(&mctx,ctx,sizeof(dk4_iter_ctx_t)); } else { dk4iter_ctx_init(&mctx); mctx.algo = DK4_ITER_ALG_NEWTON; } /* Check x0 in interval */ if (0 != (MINMAX_MINIMUM & (mctx.minmax))) { if (0 != (MINMAX_MIN_EXCL & (mctx.minmax))) { if (x0 <= mctx.xmin) { goto finished; } } else { if (x0 < mctx.xmin) { goto finished; } } } if (0 != (MINMAX_MAXIMUM & (mctx.minmax))) { if (0 != (MINMAX_MAX_EXCL & (mctx.minmax))) { if (x0 >= mctx.xmax) { goto finished; } } else { if (x0 > mctx.xmax) { goto finished; } } } /* Call function for specified algorithm. */ switch (mctx.algo) { case DK4_ITER_ALG_NEWTON : { back = dk4iter_newton(d, pp, fct, ps, x0, &mctx); } break; case DK4_ITER_ALG_FIX_POINT : { back = dk4iter_fix_point(d, pp, fct, ps, x0, &mctx); } break; } finished: return back; } /* vim: set ai sw=4 ts=4 : */