/**************************************************** * Program to demonstrate the Aitken Steffenson * * iteration subroutine * * ------------------------------------------------- * * Reference: BASIC Scientific Subroutines, Vol. II * * By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * * * * C++ version by J-P Moreau, Paris. * * (www.jpmoreau.fr) * * ------------------------------------------------- * * SAMPLE RUN: * * (Example: find a real root of f(x) = x-2*SIN(x)) * * * * Input the initial guess: * * * * X0 = -2 * * * * Convergence criterion: 1e-6 * * * * Convergence factor: -1 * * * * Maximum number of iterations: 25 * * * * * * The calculated zero is X = -1.895494 * * * * The associated Y value is Y = -0.000000 * * * * The number of iterations was: 3 * * * * ------------------------------------------------- * * Note: this algorithm fails with function (x+1)^5. * ****************************************************/ #include #include double c,e,x,x0; int m, n; //************************************************* // Function subroutine double Y(double x) { return (x-2.0*sin(x)); c = 1.0-2.0*cos(x); c = -1.0/c; } //************************************************* /********************************************** * Aitken Steffenson iteration subroutine * * ------------------------------------------- * * This subroutine calculates the zeroes of a * * function Y(x) by iterations, and employs * * Aitken acceleration to speed up convergence.* * An initial guess is required, x0, and two * * convergence factors, c and e. e relates to * * the accuracy of the estimate, and c is used * * to aid convergence. Also required is the * * maximum number of iterations, m. c=-1 is a * * normal value, if divergence occurs, smaller * * and/or positive values should be tried. * * The root is returned in x, the number of * * iterations in n. * **********************************************/ void Steffenson() { //Labels: e50,e100,e200 double x1,x2,xk,yy; int m1; n=0; e50: m1=0; x=x0; // Get yy e100: yy = Y(x); yy=x+c*yy; // Enough points for acceleration ? if (m1>0) goto e200; x1=yy; x=x1; n=n+1; m1=m1+1; goto e100; e200: x2=yy; // Perform acceleration // Guard against a zero denominator xk=x2-2.0*x1+x0; if (xk==0) xk += 0.001; xk=(x1-x0)*(x1-x0)/xk; x0=x0-xk; // Test for convergence if (n>=m) return; if (fabs(x-x0)