Function: romberg (expr, x, a, b)

Function: romberg (expr, x, a, b)

Function: romberg (expr, a, b)

Romberg integration. There are two ways to use this function. The first is an inefficient way like the definite integral version of integrate: romberg (<integrand>, <variable of integration>, <lower limit>, <upper limit>).

Examples:

(%i1) showtime: true$

(%i2) romberg (sin(y), y, 0, %pi);

Evaluation took 0.00 seconds (0.01 elapsed) using 25.293 KB.

(%o2)                   2.000000016288042

(%i3) 1/((x-1)^2+1/100) + 1/((x-2)^2+1/1000) + 1/((x-3)^2+1/200)$

(%i4) f(x) := ''%$

(%i5) rombergtol: 1e-6$

(%i6) rombergit: 15$

(%i7) romberg (f(x), x, -5, 5);

Evaluation took 11.97 seconds (12.21 elapsed) using 12.423 MB.

(%o7)                   173.6730736617464

The second is an efficient way that is used as follows:

romberg (<function name>, <lower limit>, <upper limit>);

Continuing the above example, we have:

(%i8) f(x) := (mode_declare ([function(f), x], float), ''(%th(5)))$

(%i9) translate(f);

(%o9)                          [f]

(%i10) romberg (f, -5, 5);

Evaluation took 3.51 seconds (3.86 elapsed) using 6.641 MB.

(%o10)                  173.6730736617464

The first argument must be a translated or compiled function. (If it is compiled it must be declared to return a flonum.) If the first argument is not already translated, romberg will not attempt to translate it but will give an error.

The accuracy of the integration is governed by the global variables rombergtol (default value 1.E-4) and rombergit (default value 11). romberg will return a result if the relative difference in successive approximations is less than rombergtol. It will try halving the stepsize rombergit times before it gives up. The number of iterations and function evaluations which romberg will do is governed by rombergabs and rombergmin.

romberg may be called recursively and thus can do double and triple integrals.

Example:

(%i1) assume (x > 0)$

(%i2) integrate (integrate (x*y/(x+y), y, 0, x/2), x, 1, 3)$

(%i3) radcan (%);

                    26 log(3) - 26 log(2) - 13

(%o3)             - --------------------------

(%i4) %,numer;

(%o4)                   .8193023963959073

(%i5) define_variable (x, 0.0, float, "Global variable in function F")$

(%i6) f(y) := (mode_declare (y, float), x*y/(x+y))$

(%i7) g(x) := romberg ('f, 0, x/2)$

(%i8) romberg (g, 1, 3);

(%o8)                   .8193022864324522

The advantage with this way is that the function f can be used for other purposes, like plotting. The disadvantage is that you have to think up a name for both the function f and its free variable x. Or, without the global:

(%i1) g_1(x) := (mode_declare (x, float), romberg (x*y/(x+y), y, 0, x/2))$

(%i2) romberg (g_1, 1, 3);

(%o2)                   .8193022864324522

The advantage here is shortness.

(%i3) q (a, b) := romberg (romberg (x*y/(x+y), y, 0, x/2), x, a, b)$

(%i4) q (1, 3);

(%o4)                   .8193022864324522

It is even shorter this way, and the variables do not need to be declared because they are in the context of romberg. Use of romberg for multiple integrals can have great disadvantages, though. The amount of extra calculation needed because of the geometric information thrown away by expressing multiple integrals this way can be incredible. The user should be sure to understand and use the rombergtol and rombergit switches.

Option variable: rombergabs

Default value: 0.0

Assuming that successive estimates produced by romberg are y[0], y[1], y[2], etc., then romberg will return after n iterations if (roughly speaking)

 (abs(y[n]-y[n-1]) <= rombergabs or

 abs(y[n]-y[n-1])/(if y[n]=0.0 then 1.0 else y[n]) <= rombergtol)

is true. (The condition on the number of iterations given by rombergmin must also be satisfied.) Thus if rombergabs is 0.0 (the default) you just get the relative error test. The usefulness of the additional variable comes when you want to perform an integral, where the dominant contribution comes from a small region. Then you can do the integral over the small dominant region first, using the relative accuracy check, followed by the integral over the rest of the region using the absolute accuracy check.

Example: Suppose you want to compute

'integrate (exp(-x), x, 0, 50)

(numerically) with a relative accuracy of 1 part in 10000000. Define the function. n is a counter, so we can see how many function evaluations were needed. First of all try doing the whole integral at once.

(%i1) f(x) := (mode_declare (n, integer, x, float), n:n+1, exp(-x))$

(%i2) translate(f)$

Warning-> n is an undefined global variable.

(%i3) block ([rombergtol: 1.e-6, romberabs: 0.0], n:0, romberg (f, 0, 50));

(%o3)                   1.000000000488271

(%i4) n;

(%o4)                          257

That approach required 257 function evaluations. Now do the integral intelligently, by first doing 'integrate (exp(-x), x, 0, 10) and then setting rombergabs to 1.E-6 times (this partial integral). This approach takes only 130 function evaluations.

(%i5) block ([rombergtol: 1.e-6, rombergabs:0.0, sum:0.0],

  n: 0, sum: romberg (f, 0, 10), rombergabs: sum*rombergtol, rombergtol:0.0,

      sum + romberg (f, 10, 50));

(%o5)                   1.000000001234793

(%i6) n;

(%o6)                          130

So if f(x) were a function that took a long time to compute, the second method would be about 2 times quicker.

Option variable: rombergit

Default value: 11

The accuracy of the romberg integration command is governed by the global variables rombergtol and rombergit. romberg will return a result if the relative difference in successive approximations is less than rombergtol. It will try halving the stepsize rombergit times before it gives up.

Option variable: rombergmin

Default value: 0

rombergmin governs the minimum number of function evaluations that romberg will make. romberg will evaluate its first arg. at least 2^(rombergmin+2)+1 times. This is useful for integrating oscillatory functions, when the normal converge test might sometimes wrongly pass.

Option variable: rombergtol

Default value: 1e-4

Function: tldefint (expr, x, a, b)

Equivalent to ldefint with tlimswitch set to true.