[Back to MATH SWAG index]  [Back to Main SWAG index]  [Original]

{  Calculate definite integral using Romberg method
   ------------------------------------------------
   Features:
     You can calculate *ANY* functions with any range (except infinites)
     with a *PRECISE* calculation
   Drawbacks:
     May have errors if the real result is 0. But the error is not
     significant. It may display something like 6.1200032e-18 or any
     small insignificant numbers.

   Program made by: Roby Johanes
   http://www.geocities.com/SiliconValley/Park/3230
   Finished in September 1997

   Latest update : December 1997 for submission to SWAG

   Note from author:
*  Your computer should have a numeric coprocessor intact otherwise this
   program would execute very slowly since it consumes lots of CPU power
*  If you want to have more precise calculation, make the mj (integration
   improvement) and mk (integration order) bigger.

   A little intro...
   -----------------------------
   What is Romberg integration?

   It is a method of massive recursive calculations to calculate the value
   of an integral. As many methods are being developed, this method seems
   to generalize them all. I don't want to explain it in-depth mathematic,
   but this is the general idea:

   Integral is actually the area below the function. The first idea to
   calculate definite integral without even integrating the function
   itself is using trapezium formulae. It is:

       (b - a)
   A = -------  *  [ f(b) + f(a) ]
          2

   This method is then developed into Simpson integral, which is a slight
   improvement of calculation. The first Simpson is the well-known 1/3
   Simpson's rule. Later on that rule is improved into 3/8 Simpson's rule
   which offers more precise calculations.

   The last improvement made was Boyle's rule. It is a very good
   approximation in calculating integrals. If you are interested in this
   area, I suggest you to refer to any Numerical Methods books.

   Other approach is, beside improving methods, by doing iteration using
   existing method. So, the integrating area is then divided (usually
   vertically, and I did here in my program) into n equal parts so that
   each part is small enough in order to minimize the calculation error.
   The result of each individual parts are accumulated into the final
   result. However, the drawback is that the error is also accumulated.
   But, iteration is widely used throughout the world because this
   method is the easiest way to accomplish the calculation. Also, if each
   part is small enough, the error is negligible.

   Then, there come recursive method. This method emerges after computer
   become ready to assist scientists in calculations (silly). Iterations
   can be done in the computer but costs a lot of computing resources.
   The main goal of recursive method is to reduce the calculating time.
   (oh really? By the way, this is the theorem...)
   -----------------
   Recursive trapezium approximation formulae:
                          2^(J-1)
     T(J) = T(J-1)/2 + h *   ä   f(X    )       for j > 0
                            k=1     2k-1

     Where h = (b-a) / 2^J  and Xi = a + ih

     If J = 0, then T(0) = h/2 * ( f(a) + f(b) )
   -----------------
   As we may see that the trapezium method is calculated recursively as
   above. The J here is called the order. The higher the order, the more
   precise the calculation.

   1/3 Simpson's rule, 3/8 Simpson's Rule, and Boyle's Rule are also
   implemented into recursive. However, Romberg sees more. He then gene-
   ralize them all into one formulae:
   -----------------
   Romberg formulae:

               4^k * R(J, K-1) - R(J-1, K-1)
   R (J, K) = -------------------------------   with 1 <= K <= J
                         4^k - 1

   When K = 0, R(J, 0) = Trapezium(J)
   -----------------
   Actually, R(J,0) is Trapezium Jth order
             R(J,1) is 1/3 Simpson's Rule Jth order
             R(J,2) is 3/8 Simpson's Rule Jth order
             R(J,3) is Boyle's Rule Jth order
   So, you can derive any better improvement methods. If Boyle's Rule is
   the third improvement method, then R(J,4) is the fourth, R(J,5) is the
   fifth and so on.

   I won't give the mathematical prove here, but read Numerical Methods
   books instead. I also don't remember the mathematical prove of the
   order of error, but I can tell you this: The order of error is
   O(n ^ 2k). So, the higher the K (improvement number), the error is
   much less (squared).

   Therefore, Romberg integration is a very powerful method in calculating
   integrals. You can give any number to J and K (in program mj and mk)
   but remember that your computer has a limited calculating power. Beside,
   you need a tremendous stack since the calculations involve massive
   recursive calls.

------------------------
 * USING THE PROGRAM *
------------------------

You can modify the constants mj, mk, a, and b and try it yourself!
Also, you can modify the evaluated function. The default is 1/sqrt(x).
You can say any functions as you wish as long as your Pascal provides
a mean to express it.

The program itself is explanatory. All implementations refer to the
theorem above.

}
{$N+,E-,X+}
{$M 65520,0,655360}
uses crt;
const
  mj : integer = 10;       { Integration order }
  mk : integer = 4;        { Integration improvement }
  a  : real    = 0.25;     { Lower bound }
  b  : real    = 4.0;      { Upper bound }

function f(x : real):real;
begin
  f:= 1/sqrt(x);
end;

{
  Recursive trapezium approximation formulae:
                          2^(J-1)
     T(J) = T(J-1)/2 + h *   ä   f(X    )       for j > 0
                            k=1     2k-1

     Where h = (b-a) / 2^J  and Xi = a + ih

     If J = 0, then T(0) = h/2 * ( f(a) + f(b) )
}
function trapeze(j: integer): real;
var k, n   : integer;
    h, sum : real;
begin
  n := 1 shl j;     { n = 2^J }
  h := (b-a)/n;
  if j>0 then
  begin
    sum:=0.0;
    for k:=1 to (n shr 1) do sum:=sum+f(a + (2*k-1)*h);
    trapeze:=trapeze(j-1)/2.0 + h*sum;
  end
  else trapeze:= (h/2.0) * (f(a)+f(b));
end;

{
  Romberg formulae:

              4^k * R(J, K-1) - R(J-1, K-1)
  R (J, K) = -------------------------------   with 1 <= K <= J
                        4^k - 1

  When K = 0, R(J, 0) = Trapezium(J)
}
function romberg(j, k: integer): real;
var n : real;
begin
  n := exp(k*ln(4.0));  { n = 4^k }
  if k>0 then
     romberg:=(n*romberg(j,k-1)-romberg(j-1,k-1)) / (n-1.0)
  else
     romberg:=trapeze(j);
end;

begin
  clrscr;
  writeln(romberg(mj, mk):4:8);
  readkey;
end.

[Back to MATH SWAG index]  [Back to Main SWAG index]  [Original]