DataMuseum.dk

boolean procedure Romberg special(f,n,a,b,delta,max ord,rom);
_______ _________
  value                             n,a,b,delta             ;
  _____
  real                                a,b,delta             ;
  ____
  integer                           n,          max ord     ;
  _______
  array                                                 rom ;
  _____
  procedure                       f                         ;
  _________
    comment Romberg special integrates simultaneously n functions
    _______
    (f) depending upon a real variable from a to b.
    The parameters are:
    Input
      f:     A procedure which computes the values of the n
               _________
             functions. It must have the form: f(variable,result),
             with real variable and array result [1:n].
                  ____              _____
      n:     Number of functions integrated simultaneously.
      a:     Lower limit.
      b:     Upper limit.
      delta: A relative tolerance.
    Input/Output
      max ord: On entry the maximum number of iterations.
               On exit the number of iterations performed.
    Output
      rom:     An array containing the n integration results.
                  _____
      Romberg special: is true when the desired accuracy is
                          ____
                       obtained, otherwise it is false;
                                                 _____
begin
_____
  real step,x;
  ____
  integer i,j,k,p;
  _______
  array I1,I2,sum,error,f0[1:n],trapez[1:n,1:max ord+1];
  _____
  step:=b-a;
  f(a,I1); f(b,I2);
  for i:=1 step 1 until n do
  ___      ____   _____   __
    begin
    _____
      trapez[i,1]:=(I1[i]+I2[i])⨯step/2;
      I1[i]:=0
    end;
    ___
  for k:=1 step 1 until max ord do
  ___      ____   _____         __
    begin
    _____
      for i:=1 step 1 until n do sum[i]:=error[i]:=0;
      ___      ____   _____   __
      step:=step/2;
      p:=2∧k;
          |
      for j:=p-1 step -2 until 1 do
      ___        ____    _____   __
        begin
        _____
          x:=j/p;
          x:=x⨯a+(1-x)⨯b;
          f(x,f0);
          for i:=1 step 1 until n do
          ___      ____   _____   __
            begin
            _____
              I2[i]:=sum[i]+f0[i];
              error[i]:=error[i]+(f0[i]-(I2[i]-sum[i]));
              sum[i]:=I2[i]
            end for i
            ___
        end summation of function values and errors;
        ___
      for i:=1 step 1 until n do
      ___      ____   _____   __
        I2[i]:=trapez[i,k+1]:=trapez[i,k]/2+(I2[i]+error[i])⨯step;
          comment this was a trapezoidal integration with 2∧k
          _______                                          |
          subintervals;
      p:=1;
      for j:=k step -1 until 1 do
      ___      ____    _____   __
        begin
        _____
          p:=4⨯p;
          for i:=1 step 1 until n do
          ___      ____   _____   __
            I2[i]:=trapez[i,j]:=(I2[i]⨯p-trapez[i,j])/(p-1)
        end the extrapolation from the k first approximants;
        ___
      for i:=1 step 1 until n do
      ___      ____   _____   __
        begin
        _____
          if abs(I2[i]-I1[i])>delta⨯abs(I2[i]) then goto continue
          __                                   ____ ____
        end;
        ___
      goto finis;
      ____
continue:
      for i:=1 step 1 until n do I1[i]:=I2[i]
      ___      ____   _____   __
    end for k, i.e. a whole iteration;
    ___
finis:
  Romberg special:=k_max ord;
                    <
  max ord:=k;
  for i:=1 step 1 until n do rom[i]:=I2[i]
  ___      ____   _____   __
end Romberg special;
___
  [end] [stop]