doc/share/romberg.usg

   1
   2 ROMBERG is now MACSYM;ROMBRG FASL and is autoloading.  Just call
   3 ROMBERG and it will be there. - JPG
   4 \x1f
   5 GJC@MIT-MC 01/12/81 13:50:27
   6 To: (FILE [SHARE;ROMBRG USAGE]) at MIT-MC
   7 The old conventions of what are efficient ways to call ROMBERG,
   8 or what are ways which work in at all, are now repaired.
   9
  10 Given F(N):=ROMBERG(expression,VAR,0,N)$
  11
  12 If F is translated/compiled then the expression
  13 will also be translated/compiled! Also, the
  14 old "efficient" way, which could ONLY BE USED
  15 ON TRANSLATED FUNCTIONS, (what a pain in the neck
  16 right?), now works on all functions.
  17
  18 This change is made possible due to the increased
  19 power of the macsyma->lisp translator, and some
  20 reworking of the numerical code.
  21
  22 [Anonymous functions, i.e. LAMBDA expressions, are
  23 now translated into lisp functions.]
  24
  25 The entry point to ROMBERG is now ROMBERG_SUBR,
  26 which only takes 3 arguments, the first being a function.
  27 [This is never seen by the user.]
  28
  29 ROMBERG(F,A,B) => ROMBERG_SUBR(F,A,B)
  30 ROMBERG(X^2,X,A,B) => ROMBERG_SUBR(LAMBDA([X],X^2),A,B)
  31
  32 The arrow "=>" signifies a transformation which takes
  33 place through the same mechanisms as the *new* macsyma
  34 "macro" feature. To see how this feature can be
  35 useful in numerical applications see SHARE;SIMPSN MACRO.
  36
  37 Proper use of MODE_DECLARE is still important of course.
  38
  39 p.s. These same extensions apply to the syntax of the
  40 INTERPOLATE command. However, the restrictions on the
  41 use of the PLOT2 command unfortunately still hold.
  42 This takes lots of work on PLOT2 to fix.
  43 \x1f
  44 GJC@MIT-MC 03/20/80 19:09:30
  45 To: (FILE [SHARE;ROMBRG USAGE]) at MIT-MC
  46 The source for ROMBERG is now MAXSRC;ROMBRG > and NUMER;ROMBRG MMACRO.
  47 The three argument version of ROMBERG now works for macsyma functions,
  48 not just translated functions. Problems with arguments not getting
  49 properly floated, e.g. ROMBERG(X,X,1/10,2/10) have been fixed, problems
  50 with the arguments to ROMBERG not getting evaluated in compiled code
  51 have been fixed.
  52 \x1f
  53 CFFK@MIT-MC 12 AUG 1979 1743-EDT
  54 To: (FILE [SHARE;ROMBRG USAGE]) at MIT-MC
  55 At SK's suggestion, ROMBERG (BROMBERG) will now exit if an absolute
  56 error condition is satisfied.  This is governed by the new variable
  57 ROMBERGABS (BROMBERABS) - default value 0.0 (0.0B0).
  58
  59 Assuming that successive estimates produced by ROMBERG are Y[0], Y[1],
  60 Y[2] etc., then ROMBERG will return after N iterations if
  61 (roughly speaking)
  62 (ABS(Y[N]-Y[N-1]) <= ROMBERGABS OR
  63  ABS(Y[N]-Y[N-1])/(IF Y[N]=0.0 THEN 1.0 ELSE Y[N]) <= ROMBERGTOL)
  64 is TRUE.  (The condition on the number of iterations given by ROMBERGMIN
  65 must also be satisfied.)  (The old exit condition, was on the relative
  66 error if ABS(Y[N]) > 1.0, and on the absolute error otherwise.  Pretty
  67 random!)
  68
  69 Thus if ROMBERGABS is 0.0 (the default) you just get the relative error
  70 test.  The usefulness of the additional variable comes when
  71 you want to perform an integral, where the dominant contribution
  72 comes from a small region.  Then you can do the integral over
  73 the small dominant region first, using the relative accuracy
  74 check, followed by the integral over the rest of the region
  75 using the absolute accuracy check.
  76
  77 Example:  Suppose you want to compute
  78    Integral(exp(-x),x,0,50)
  79 (numerically) with a relative accuracy of  1 part in 10000000.
  80
  81   /* Define the function.  N is a counter, so we can see how many
  82      function evaluations were needed. */
  83 F(X):=(MODEDECLARE(N,INTEGER,X,FLOAT),N:N+1,EXP(-X))$
  84 TRANSLATE(F)$
  85   /* First of all try doing the whole integral at once */
  86 BLOCK([ROMBERGTOL:1.E-6,ROMBERABS:0.],N:0,ROMBERG(F,0,50)); ==> 1.00000003
  87 N; ==> 257  /* Number of function evaluations*/
  88   /* Now do the integral intelligently, by first doing
  89      Integral(exp(-x),x,0,10) and then setting ROMBERGABS to 1.E-6*(this
  90      partial integral).  */
  91 BLOCK([ROMBERGTOL:1.E-6,ROMBERGABS:0.,SUM:0.],
  92       N:0,SUM:ROMBERG(F,0,10),ROMBERGABS:SUM*ROMBERGTOL,ROMBERGTOL:0.,
  93       SUM+ROMBERG(F,10,50));  ==> 1.00000001  /* Same as before */
  94 N;  ==> 130
  95
  96 So if F(X) were a function that took a long time to compute, the
  97 second method would be about 2 times quicker.
  98 \x1f
  99 CFFK@MIT-MC 05/04/78 16:52:43
 100 To: (FILE [SHARE;ROMBRG USAGE]) at MIT-MC
 101 There is a new option ROMBERGMIN whose default value is 0, that
 102 govern the minimum number of function evaluations that ROMBERG will
 103 make.  ROMBERG will evaluate its first arg. at least 2^(ROMBERGMIN+2)+1
 104 times.  This is useful for integrating oscillatory functions, when
 105 the normal converge test might sometimes wrongly pass.
 106
 107 There are 2 ways of using this function:
 108
 109 1) An inefficient way that looks like the definite integral version of
 110 INTEGRATE:
 111         ROMBERG(<integrand>,<variable of integration>,<lower limit>,
 112                 <upper limit>);
 113 Examples:
 114         ROMBERG(SIN(Y),Y,1,%PI);
 115                 TIME= 39 MSEC.          1.5403023
 116         F(X):=1/(X^5+X+1);
 117         ROMBERG(F(X),X,1.5,0);
 118                 TIME= 162 MSEC.         - 0.75293843
 119
 120 2) An efficient way that is more like RJF's ROMBERG function:
 121         ROMBERG(<function name>,<lower limit>,<upper limit>);
 122 Example:
 123         F(X):=(MODEDECLARE([FUNCTION(F),X],FLOAT),1/(X^5+X+1));
 124         TRANSLATE(F);
 125         ROMBERG(F,1.5,0);
 126                 TIME= 13 MSEC.          - 0.75293843
 127 The first argument must be a TRANSLATEd or compiled function, which
 128 returns a floating point number.  (If it is
 129 compiled it must be declared to return a FLONUM.)  If the first argument
 130 is not already TRANSLATEd, ROMBERG will not attempt to TRANSLATE it but
 131 will give an error.
 132
 133 The accuracy of the integration is governed by the global variables
 134 ROMBERGTOL (default value 1.E-4) and ROMBERGIT (default value 11).
 135 ROMBERG will return a result if the relative difference in successive
 136 approximations is less than ROMBERGTOL.  It will try halving the
 137 stepsize ROMBERGIT times before it gives up.
 138
 139 ROMBERG may be called recursively and thus can do double and triple
 140 integrals.
 141 Example:
 142         INTEGRATE(INTEGRATE(X*Y/(X+Y),Y,0,X/2),X,1,3);
 143                                 13/3 (2 LOG(2/3) + 1)
 144         %,NUMER;
 145                                 0.81930233
 146
 147         DEFINE_VARIABLE(X,0.0,FLOAT,"Global variable in function F")$
 148         F(Y):=(MODE_DECLARE(Y,FLOAT), X*Y/(X+Y) )$
 149         G(X):=ROMBERG('F,0,X/2)$
 150         ROMBERG(G,1,3);
 151                                 0.8193023
 152
 153 The advantage with this way is that the function F can be used for other
 154 purposes, like plotting. The disadvantage is that you have to think up
 155 a name for both the function F and its free variable X.
 156
 157 Or, without the global:
 158         G1(X):=(MODE_DECLARE(X,FLOAT), ROMBERG(X*Y/(X+Y),Y,0,X/2))$
 159         ROMBERG(G1,1,3);
 160                                 0.8193023
 161 The advantage here is shortness.
 162
 163         Q(A,B):=ROMBERG(ROMBERG(X*Y/(X+Y),Y,0,X/2),X,A,B)$
 164         Q(1,3);
 165                                 0.8193023
 166 It is even shorter this way, and the variables do not need to be declared
 167 because they are in the context of ROMBERG.
 168
 169 Use of ROMBERG for multiple integrals can have great disadvantages,
 170 though.  The amount of extra calculation needed because of the geometric
 171 information thrown away by expressing multiple integrals this way can
 172 be incredible.  The user should understand and use the ROMBERGTOL and
 173 ROMBERGIT switches.
 174