home *** CD-ROM | disk | FTP | other *** search
/ HAM Radio 1 / HamRadio.cdr / tech / design2 / romberg.pas < prev    next >
Pascal/Delphi Source File  |  1986-06-16  |  3KB  |  83 lines
  1. Program Rom;          { Romberg Integration   07/31/84  J. E. Crawford  }
  2.  
  3.     label 10;
  4.  
  5.     var
  6.          R: array[1..10,1..10] of real;
  7.          a,b,total: real;
  8.          i,j,k,m,order,nosteps: integer;
  9.  
  10. function YtoX(y,x: real):real;
  11.     begin
  12.     YtoX:= exp(x*ln(y));
  13.     end;
  14.  
  15. function fnc(x:real):real;
  16.     begin                             { Evaluate Error Function Integral }
  17.     fnc:= 1/sqrt(2*Pi)*exp(-x*x/2);   { The Function to be Integrated is }
  18.     end;                              { Placed Here                      }
  19.  
  20. function integral(a,b: real; nosteps: integer): real;
  21.  
  22.  
  23. { Reference    Applied Numerical Methods for Digital Computation  }
  24. {                     James, Smith, Wolford                       }
  25. {                     1977   Harper & Row Publishers              }
  26.  
  27. {  This function does simple trapezoidal integration of nosteps strips }
  28.     var
  29.          h,f1,f2,sum: real;
  30.          ii: integer;
  31.  
  32.     {  a, b are limits of integration }
  33.     {  f1 = function evaluated at first endpoint }
  34.     {  f2 = function evaluated at last endpoint  }
  35.     {  sum = 2.0*intermediated functional values }
  36.     {  h = strip thickness                       }
  37.  
  38. begin
  39. h:= (b-a)/nosteps;
  40. sum:= 0.0;
  41. f1:= fnc(a);
  42. for ii:= 1 to nosteps - 1  do
  43.     begin
  44.     sum:= sum + fnc(a+ii*h);
  45.     end;
  46. f2:= fnc(b);
  47. integral:= h/2.0*(f1 + 2.0*sum + f2);   { Rule for trapezoidal integration }
  48. end;
  49.  
  50. begin
  51. for i:= 1 to 10 do       { Initialize array to zero }
  52.     begin
  53.     for j:= 1 to 10 do
  54.          begin
  55.          R[i,j]:= 0.0;
  56.          end;
  57.     end;
  58. write(' Enter lower and upper limits  ');
  59. readln(a,b);
  60. nosteps:= 2;                 { # Total Steps Depends on Desired Accuracy }
  61. R[1,1]:= Integral(a,b,nosteps);     { First integration using nosteps }
  62. nosteps:= nosteps*2;
  63. R[1,2]:= Integral(a,b,nosteps);     { 2nd integration using 2*nosteps }
  64. order:= 1;
  65. while order < 10 do                 { Places limit on how many iterrations }
  66.  
  67.     begin
  68.     m:= order + 1;                 { used as 2nd index of first matrix }
  69.     For i:= 1 to order do          { i is order of extrapolation       }
  70.          begin
  71.          R[i+1,m-1]:= (YtoX(4.0,i)*R[i,m]-R[i,m-1])/(YtoX(4.0,i) - 1.0);
  72.          m:= m - 1;
  73.          end;
  74.                             { The Next Line Controls Accuracy Desired }
  75.  
  76.    if (abs(R[order+1,1] - R[order,2])/R[order+1,1] < 1.e-7 ) then  goto 10;
  77.    nosteps:= nosteps*2;
  78.    R[1,order + 2]:= Integral(a,b,nosteps);  { Does integration for more strips }
  79.   order:= order + 1;          { Increase available order of extrapolation used }
  80. end;
  81. 10: writeln('Integral = ',R[order + 1,1]);
  82. end.
  83.