NUMERICAL ANALYSIS: INTERPOLATION AND EXTRAPOLATION
CONVERGENCE ACCELERATION DURING THE 20TH CENTURY 1 INTRODUCTION IN NUMERICAL ANALYSIS MANY METHODS PRODUCE SEQUENCES, FOR INSTANCE ITERATIVE METHODS FOR SOLVING SYSTEMS OF EQUATIONS, METHODS INVOLVING SERIES EXPANSIONS, DISCRETIZATION METHODS (THAT IS METHODS DEPENDING ON A PARAMETER SUCH THAT THE APPROXIMATE SOLUTION TENDS TO THE EXACT ONE WHEN THE PA- RAMETER TENDS TO ZERO), PERTURBATION METHODS, ETC SOMETIMES, THE CONVERGENCE OF THESE SEQUENCES IS SLOW AND THEIR EECTIVE USE IS QUITE LIMITED CONVERGENCE ACCELERATION METHODS CONSIST OF TRANS- FORMING A SLOWLY CONVERGING SEQUENCE (SN) INTO A NEW SEQUENCE (TN) CONVERGING TO THE SAME LIMIT FASTER THAN THE INITIAL ONE AMONG SUCH SEQUENCE TRANSFORMATIONS, THE MOST WELL KNOWN ARE CERTAINLY RICHARDSON’S EXTRAPOLATION ALGORITHM AND AITKEN’S 2 PROCESS ALL KNOWN METHODS ARE CONSTRUCTED BY EXTRAPOLATION AND THEY ARE OFTEN CALLED EXTRAPOLATION METHODS THE IDEA CONSISTS OF INTERPOLATING THE TERMS SN; SN+1; : : : ; SN+K OF THE SEQUENCE TO BE TRANSFORMED BY A SEQUENCE SATISFYING A CERTAIN RELATIONSHIP DEPENDING ON PARAMETERS THIS SET OF SEQUENCES IS CALLED KERNEL OF THE TRANSFORMATION AND EVERY SEQUENCE OF THIS SET IS TRANSFORMED INTO A CONSTANT SEQUENCE BY THE TRANSFORMATION INTO CONSIDERATION FOR EXAMPLE, AS WE WILL SEE BELOW, THE KERNEL OF AITKEN’S 2 PROCESS IS THE SET OF SEQUENCES SATISFYING ?N; A0(SN ?S) + A1(SN+1 ?S) = 0, WHERE A0 AND A1 ARE PARAMETERS SUCH THAT A0 + A1 ?= 0 IF AITKEN’S PROCESS IS APPLIED TO SUCH A SEQUENCE, THEN THE CONSTANT SEQUENCE (TN = S) IS OBTAINED THE PARAMETERS INVOLVED IN THE DENITION OF THE KERNEL ARE UNIQUELY DETERMINED BY THE INTERPOLATION CONDITIONS AND THEN THE LIMIT OF THE INTERPOLATING SEQUENCE OF THE KERNEL IS TAKEN AS AN APPROXIMATION OF THE LIMIT OF THE SEQUENCE TO BE ACCELERATED SINCE THIS LIMIT DEPENDS ON THE INDEX N, IT WILL BE DENOTED BY TN EECTIVELY, THE SEQUENCE (SN) HAS BEEN TRANSFORMED INTO A NEW SEQUENCE (TN) THIS PAPER, WHICH IS BASED ON [31], BUT INCLUDES NEW DEVELOPMENTS OBTAINED SINCE 1995, PRESENTS MY PERSONAL VIEWS ON THE HISTORICAL DEVELOPMENT OF THIS SUBJECT DURING THE 20TH CENTURY I DO NOT PRETEND TO BE EXHAUSTIVE NOR EVEN TO QUOTE EVERY IMPORTANT CONTRIBUTION (IF A REFERENCE DOES NOT ANISI; N = 0; 1; : : : ; WHERE THE NUMBERS ANI ARE CONSTANTS INDEPENDENT OF THE TERMS OF THE SEQUENCE (SN) SUCH A LINEAR TRANSFORMATION IS USUALLY CALLED A SUMMATION PROCESS AND ITS PROPERTIES ARE COMPLETELY DETERMINED BY THE MATRIX A=(ANI) FOR PRACTICAL REASONS, ONLY A NITE NUMBER OF THE COECIENTS ANI ARE DIERENT FROM ZERO FOR EACH N AMONG SUCH PROCESSES ARE THOSE NAMED AFTER EULER, CESARO AND HOLDER IN THE CASE OF LINEAR METHODS, THE CONVERGENCE OF THE SEQUENCE (TN) TO S FOR ANY CONVERGING SEQUENCE (SN) IS GOVERNED BY THE TOEPLITZ SUMMABILITY THEOREM; SEE [115] FOR A REVIEW EXAMPLES OF SUCH PROCESSES ARE TN = N + SI OR TN = K + N+K X I=N SI: IN THE SECOND CASE, THE SEQUENCE (TN) ALSO DEPENDS ON A SECOND INDEX, K, AND THE CONVERGENCE HAS TO BE STUDIED EITHER WHEN K IS XED AND N TENDS TO INNITY, OR WHEN N IS XED AND K TENDS TO INNITY N : IN A 1927 PAPER [120] HE USED THE SAME TECHNIQUE TO SOLVE A 6TH ORDER DIERENTIAL EIGENVALUE PROBLEM HIS PROCESS WAS CALLED (H2; H4)-EXTRAPOLATION RICHARDSON EXTRAPOLATION CONSISTS OF COMPUTING THE VALUE AT 0, DENOTED BY T (N) K , OF THE INTERPOLATION POLYNOMIAL OF THE DEGREE AT MOST K, WHICH PASSES THROUGH THE POINTS (XN; SN); : : : ; (XN+K; SN+K) USING THE NEVILLE–AITKEN SCHEME FOR THESE INTERPOLATION POLYNOMIALS, WE IMMEDIATELY OBTAIN T (N) K+1 = XN+K+1T (N) K ?XNT (N+1) K XN+K+1 ?XN WITH T (N) = SN LET US MENTION THAT RICHARDSON REFERRED TO A 1926 PAPER BY NIKOLAI NIKOLAEVICH BOGOLYUBOV (BORN IN 1909) AND NIKOLAI MITROFANOVICH KRYLOV (1879–1955) WHERE THE PROCEDURE (OFTEN CALLED THE DEFERRED APPROACH TO THE LIMIT) CAN ALREADY BE FOUND [11] IN 1955, WERNER ROMBERG (BORN IN 1909) WAS THE RST TO USE REPEATEDLY AN ELIMINATION APPROACH FOR IMPROVING THE ACCURACY OF THE TRAPEZOIDAL RULE [121] HE HIMSELF REFERS TO THE BOOK OF LOTHAR COLLATZ (1910–1990) OF 1951 [50] THE PROCEDURE BECAME WIDELY KNOWN AFTER THE RIGOROUS ERROR ANALYSIS GIVEN IN 1961 BY FRIEDRICH L BAUER [3] AND THE WORK OF EDUARD L STIEFEL (1909–1978) [138] ROMBERG’S DERIVATION OF HIS PROCESS WAS HEURISTIC IT WAS PROVED BY PIERRE-JEAN LAURENT IN 1963 [81] THAT THE PROCESS COMES OUT FROM THE RICHARDSON PROCESS BY CHOOSING XN LAURENT ALSO GAVE CONDITIONS ON THE CHOICE OF THE SEQUENCE (XN) IN ORDER THAT THE SEQUENCES ( TEND TO S EITHER WHEN K OR N TENDS TO INNITY WEAKER CONDITIONS WERE GIVEN BY MICHEL CROUZEIX AND ALAIN L MIGNOT IN [52, PP 52–55] AS WE SHALL SEE BELOW, EXTENSIONS OF ROMBERG’S METHOD TO NONSMOOTH INTEGRANDS LEADS TO A METHOD CALLED THE E-ALGORITHM APPLICATIONS OF EXTRAPOLATION TO THE NUMERICAL SOLUTION OF ORDINARY DIERENTIAL EQUATIONS WERE STUD- IED BY HC BOLTON AND HI SCOINS IN 1956 [12], ROLAND BULIRSCH AND JOSEF STOER IN 1964–1966 [47] AND WILLIAM B GRAGG [65] IN 1965 THE CASE OF DIERENCE METHODS FOR PARTIAL DIERENTIAL EQUATIONS WAS TREATED BY GURI IVANOVICH MARCHUK AND VV SHAIDUROV [91] STURM–LIOUVILLE PROBLEMS ARE DISCUSSED IN [117] FINALLY, WE MENTION THAT HEINZ RUTISHAUSER (1918–1970) POINTED OUT IN SN+2 ?2SN+1 + SN ; N = 0; 1; : : : THE METHOD WAS STATED BY ALEXANDER CRAIG AITKEN (1895–1967) IN 1926 [1], WHO USED IT TO AC- CELERATE THE CONVERGENCE OF BERNOULLI’S METHOD FOR COMPUTING THE DOMINANT ZERO OF A POLYNOMIAL AITKEN POINTED OUT THAT THE SAME METHOD WAS OBTAINED BY HANS VON NAEGELSBACH (1838) IN IN HIS STUDY OF FURSTENAU’S METHOD FOR SOLVING NONLINEAR EQUATIONS [104] THE PROCESS WAS ALSO GIVEN BY JAMES CLERK MAXWELL (1831–1879) IN HIS TREATISE ON ELECTRICITY AND MAGNETISM OF [95] HOWEVER, NEITHER NAEGELSBACH NOR MAXWELL USED IT FOR THE PURPOSE OF ACCELERATION MAXWELL WANTED TO ND THE EQUILIBRIUM POSITION OF A POINTER OSCILLATING WITH AN EXPONENTIALLY DAMPED SIMPLE HARMONIC MOTION FROM THREE EXPERIMENTAL MEASUREMENTS IT IS SURPRISING THAT AITKEN’S PROCESS WAS KNOWN TO TAKAKAZU SEKI (1642–1708), OFTEN CONSIDERED THE GREATEST JAPANESE MATHEMATICIAN IN HIS BOOK KATSUYO SANPO, VOL IV, HE USED THIS PROCESS TO COMPUTE THE VALUE OF , THE LENGTH OF A CHORD AND THE VOLUME OF A SPHERE THIS BOOK WAS WRITTEN AROUND 1680 BUT ONLY PUBLISHED IN 1712 BY HIS DISCIPLE MURAHIDE ARAKI PARTS OF IT CAN BE FOUND IN [73] LET US MENTION THAT THE JAPANESE CHARACTERS CORRESPONDING TO TAKAKAZU HAVE ANOTHER PRONOUNCIATION WHICH IS KOWA THIS IS THE REASON WHY THIS MATHEMATICIAN IS OFTEN CALLED, ERRONEOUSLY AS IN [29,31] SEKI KOWA WHAT MAKES AITKEN’S PROCESS SO POPULAR IS THAT IT ACCELERATES THE CONVERGENCE OF ALL LINEARLY CONVERGING SEQUENCES, THAT IS SEQUENCES SUCH THAT ?A ?= LIM N?? SN+1 ?S SN ?S = A: IT CAN EVEN ACCELERATE SOME LOGARITHMIC SEQUENCES (THAT IS CORRESPONDING TO A = 1) WHICH ARE THOSE WITH THE SLOWEST CONVERGENCE AND THE MOST DICULT TO ACCELERATE AITKEN’S 2 PROCESS IS EXACT (WHICH MEANS THAT ?N; TN = S) FOR SEQUENCES SATISFYING, A1(SN+1 ?S) = 0; ?N; A0A1 ?= 0; A0 + A1 ?= 0 SUCH SEQUENCES FORM THE KERNEL OF AITKEN’S PROCESS THE IDEA NATURALLY AROSE OF NDING A TRANSFORMATION WITH THE KERNEL A0(SN ?S) + · · · + AK(SN+K ?S) = 0; ?N; A0AK ?= 0; A0 + · · · + AK ?= 0 A PARTICULAR CASE OF K = 2 WAS ALREADY TREATED BY MAXWELL IN HIS BOOK OF 1873 AND A PARTICULAR CASE OF AN ARBITRARY VALUE OF K WAS STUDIED BY TH O’BEIRNE IN : WHEN K = 1, SHANKS TRANSFORMATION REDUCES TO THE AITKEN’S 2 PROCESS IT CAN BE PROVED THAT EK(SN)=S; ?N IF AND ONLY IF (SN) BELONGS TO THE KERNEL OF THE TRANSFORMATION GIVEN ABOVE THE SAME RATIOS OF DETERMINANTS WERE OBTAINED BY RJ SCHMIDT IN 1941 [127] IN HIS STUDY OF A METHOD FOR SOLVING SYSTEMS OF LINEAR EQUATIONS THE DETERMINANTS INVOLVED IN THE DENITION OF EK(SN) HAVE A VERY SPECIAL STRUCTURE THEY ARE CALLED HANKEL DETERMINANTS AND WERE STUDIED BY HERMANN HANKEL (1839–1873) IN HIS THESIS IN 1861 [72] SUCH DETERMINANTS SATISFY A VE-TERM RECURRENCE RELATIONSHIP THIS RELATION WAS USED BY O’BEIRNE AND SHANKS TO IMPLEMENT THE TRANSFORMATION BY COMPUTING SEPARATELY THE NUMERATORS AND THE DENOMINATORS OF THE EK(SN)’S HOWEVER, NUMERICAL ANALYSTS KNOW IT IS DICULT TO COMPUTE DETERMINANTS (TOO MANY ARITHMETICAL OPERATIONS ARE NEEDED AND ROUNDING ERRORS DUE TO THE COMPUTER OFTEN LEAD TO A COMPLETELY WRONG RESULT) A RECURSIVE PROCEDURE FOR COMPUTING THE EK(SN)’S WITHOUT COMPUTING THE DETERMINANTS INVOLVED IN THEIR DENITION WAS NEEDED THIS ALGORITHM WAS OBTAINED IN 1956 BY PETER WYNN IT IS CALLED THE -ALGORITHM [147] IT IS AS FOLLOWS ONE STARTS WITH (N) ?1 = 0; (N) 0 = SN AND THEN (N) K+1 = (N+1) K?1 + (N+1) K ?(N) K : NOTE THAT THE NUMBERS (N) K ’S LL OUT A TWO-DIMENSIONAL ARRAY THE -ALGORITHM IS RELATED TO SHANKS TRANSFORMATION BY (N) 2K = EK(SN) AND (N) 2K+1 = 1=EK(SN): THUS, THE ’S WITH AN ODD LOWER INDEX ARE ONLY AUXILIARY QUANTITIES THEY CAN BE ELIMINATED FROM THE ALGORITHM, THUS LEADING TO THE SO-CALLED CROSS RULE DUE TO WYNN [153] WHEN IMPLEMENTING THE -ALGORITHM OR USING WYNN’S CROSS RULE, DIVISION BY ZERO CAN OCCUR AND THE ALGORITHM MUST BE STOPPED HOWEVER, IF THE SINGULARITY IS CONNED, A TERM THAT WILL AGAIN BE USED IN SECTION 16, THAT IS IF IT OCCURS ONLY FOR SOME ADJACENT VALUES OF THE INDEXES JUMP OVER IT BY USING SINGULAR RULES AND CONTINUE THE COMPUTATION IF A DIVISION BY A NUMBER CLOSE TO ZERO ARISES, THE ALGORITHM BECOMES NUMERICALLY UNSTABLE DUE TO THE CANCELLATION ERRORS A SIMILAR SITUATION HOLDS FOR THE OTHER CONVERGENCE ACCELERATION ALGORITHMS THE STUDY OF SUCH PROBLEMS WAS INITIATED BY WYNN IN 1963 [151], WHO PROPOSED PARTICULAR RULES FOR THE -ALGORITHM WHICH ARE MORE STABLE THAN THE USUAL RULES THEY WERE EXTENDED BY FLORENT CORDELLIER IN 1979 [51,151] PARTICULAR RULES FOR THE -ALGORITHM WERE OBTAINED BY REDIVO ZAGLIA [155] THE CONVERGENCE AND ACCELERATION PROPERTIES OF THE -ALGORITHM HAVE ONLY BEEN COMPLETELY DE- SCRIBED ONLY FOR TWO CLASSES OF SEQUENCES, NAMELY TOTALLY MONOTONIC AND TOTALLY OSCILLATING SEQUENCES [154,15,16] SHANKS’ TRANSFORMATION AND THE -ALGORITHM HAVE CLOSE CONNECTIONS TO PADE APPROXIMANTS, CONTIN- UED FRACTIONS AND FORMAL ORTHOGONAL POLYNOMIALS; SEE, FOR EXAMPLE [18] 23 SUBSEQUENT DEVELOPMENTS THE SHANKS TRANSFORMATION AND THE -ALGORITHM SPARKED THE REBIRTH OF THE STUDY OF NONLINEAR AC- CELERATION PROCESSES THEY NOW FORM AN INDEPENDENT CHAPTER IN NUMERICAL ANALYSIS WITH CONNECTIONS TO OTHER IMPORTANT TOPICS SUCH AS ORTHOGONAL AND BIORTHOGONAL POLYNOMIALS, CONTINUED FRACTIONS, AND PADE APPROXIMANTS THEY ALSO HAVE APPLICATIONS TO THE SOLUTION OF SYSTEMS OF LINEAR AND NONLINEAR EQUATIONS, THE COMPUTATION OF THE EIGENVALUES OF A MATRIX, THE SOLUTION OF SYSTEMS OF LINEAR AND NONLINEAR EQUATIONS, AND MANY OTHER TOPICS, SEE [40] AMONG OTHER ACCELERATION METHODS WHICH WERE OBTAINED AND STUDIED, ARE THE W-PROCESS OF SAMUEL LUBKIN [89], THE METHOD OF KJELL J OVERHOLT [110], THE -ALGORITHM OF WYNN [148], THE G-TRANSFORMATION OF HL GRAY, TA ATCHISON AND GV MCWILLIAMS [70], THE -ALGORITHM OF CLAUDE BREZINSKI [14], THE TRANSFORMATIONS OF BERNARD GERMAIN– BONNE [63] AND THE VARIOUS TRANSFORMATIONS DUE TO DAVID LEVIN [85] TO MY KNOWLEDGE, THE ONLY KNOWN ACCELERATION THEOREM FOR THE -ALGORITHM WAS OBTAINED BY NAOKI OSADA [108] SIMULTANE- OUSLY, SEVERAL APPLICATIONS BEGAN TO APPEAR FOR EXAMPLE, THE -ALGORITHM PROVIDES A QUADRATICALLY CONVERGENT METHOD FOR SOLVING SYSTEMS OF NONLINEAR EQUATIONS AND ITS DOES NOT REQUIRE THE KNOWL- EDGE OF ANY DERIVATIVE THIS PROCEDURE WAS PROPOSED SIMULTANEOUSLY BY BREZINSKI [13] AND ECKHART GEKELER [61] IT HAS IMPORTANT APPLICATIONS TO THE SOLUTION OF BOUNDARY VALUE PROBLEMS FOR ORDINARY DIERENTIAL EQUATIONS [44] MANY OTHER ALGORITHMS ARE GIVEN IN THE WORK OF ERNST JOACHIM WENIGER [145], WHICH ALSO CONTAINS APPLICATIONS TO PHYSICS, OR IN THE BOOK OF BREZINSKI AND MICHELA RE- DIVO ZAGLIA [40] WHERE APPLICATIONS TO VARIOUS DOMAINS OF NUMERICAL ANALYSIS CAN BE FOUND THE AUTHORS OF THIS BOOK PROVIDE FORTRAN SUBROUTINES THE BOOK OF ANNIE CUYT AND LUC WUYTACK MUST ALSO BE MENTIONED [53] THE -ALGORITHM HAS BEEN APPLIED TO STATISTICS, SEE THE WORK OF ALAIN BERLINET [9], AND TO THE ACCELERATION OF THE CONVERGENCE OF SEQUENCES OF RANDOM VARIABLES, CONSID- ERED BY HELENE LAVASTRE [82] APPLICATIONS TO OPTIMIZATION WERE PROPOSED BY LE FERRAND [84] AND BOUCHTA RHANIZAR [118] INSTEAD OF USING A QUITE COMPLICATED ALGORITHM, SUCH AS THE -ALGORITHM, IT CAN BE INTERESTING TO USE A SIMPLER ONE (FOR INSTANCE, AITKEN’S 2 PROCESS) ITERATIVELY SUCH A USE CONSISTS OF AP- PLYING THE ALGORITHM TO (SN) TO PRODUCE A NEW SEQUENCE (TN), THEN TO APPLY THE SAME ALGORITHM TO (TN), AND SO ON FOR EXAMPLE, APPLYING THE ITERATED 2 PROCESS TO THE SUCCESSIVE CONVERGENTS