Numeric Methods of Integration

I will keep this light and understandable…​ I promise!

Background

Some of the first applications of mechanical and later electronic computers were for doing tedious mathematic computations. This was long before user input and databases were a thing.

Modern engineering and science applications depend heavily on calculus, especially integrals, to compute answers to challenging problems.

Many of us slugged through calculus classes learning integration techniques. But sometimes that is either too difficult or actually impossible, and it is necessary to approximate the integrals using something called numerical methods which uses the sum (addition) of many little approximations to build up a number which is close to the actual result.

This blog article shows how to approximate integrals over a range [low to high], by dividing the range into a large number of little steps.

In general, the more steps you take, the longer the program will take to execute. But the complexity of the function also can impact the speed.

The Problem

I wanted to solve an Engineering problem which involved integrals and wondered if I could do it with EWB. Would it

  • be accurate enough

  • be fast enough

  • would it be so slow as to time out and give user warnings

In many of these calculation problems there is the speed versus accuracy. You can get quick results, or you can get accurate results, but not always both.

I would need to determine what is the approximate number of calculations per second and how many calculations are needed for a good result.

Example

My first exposure to this challenge was during the time of Commodore Pets when I was in grade nine.

Those were 1 MHz 8 bit computers with incredibly slow math processing in its BASIC interpretter.

I was doing a project which needed to calculate the ability of solar panels to heat a house based on the lattitude and other factors.

The calculations took several minutes to come up with an estimate. Those equations looked foreign and magical.

A more recent example would be estimating the size of the population who would contract Covid-19 during a given month given all the input conditions. The formulas for this sort of thing were well known, but at the time we only had estimates for the parameters so they could give approximate boundaries of high and low estimates.

Mere JavaScript running on even cell phones is thousands of times faster now and can do so much!

Deliverables

You can look at this program to see how to:

  • convert an equation in a text string to a JavaScript function

  • perform simple integrals (using Simpson method)

  • time a function

  • calculate the percentage difference between two numbers

Results

I include four functions trials.

f(x) is the function in word processing terms.

source code is the pure javascript of the function.

range is the numbers we are integrating over [low, high]

Steps is the number of tiny steps. Larger n will give more accurate results but will be slower.

Exact answer is the perfect answer.

Time is how long the strategy took (Chrome on M1 Mac processor)

f(x) source code range steps exact answer time (chrome)

x^3

x*x*x

[0, 1]

100

0.25

0 ms

1/x

1/x

[1, 100]

1000

4.605170..

1 ms

x

x

[0, 5000]

5,000,000

12,500,000

1517 ms

x

x

[0, 6000]

6,000,000

18,000,000

1823 ms

I gave the timing results using Chrome, with its fast V8 JavaScript interpretter engine. EWB’s IDE uses the pitiful Internet Explorer which is not only not HTML5 compliant, but is much slower. It is so slow that it timed out on these tests.

Chrome seems to execute functions in about 0.0003 ms, giving up to about 3000 steps per millisecond.

Functions

The function to evaluate is written in pure JavaScript. So if you want to use log or power, sin or cos, you would need to write them as Math.power(x), Math.log(x), etc.

In our examples here, the math is easily written with just involving numbers and x.

A Bit of Math

The integrate() function looks a little strange, what is that funky math?

As I said, I use the Simpson’s method here. It doesn’t just break up the graph into n small samples, it uses rounding of the curve to try to improve accuracy of the approximation.

You can read about that in calculus web sites if you are interested. But in general, it usually produces more accurate results.

Converting Text of f(x) into a JavaScript Function

I could have coded each of the functions as a pure Pascal function which I could compile at design time.

Or I could have included an interpretter which parsed and computing the equations. That would be helpful for portability to other Pascal environments. But it would be slow because the interpretters are usually not as fast as native code and speed is a concern.

Instead, I decided I wanted to pass a string of the JavaScript code and use CreateObject() to execute the bare JavaScript, so we could change the function at run-time.

I would caution against letting users enter their own functions this way. They could call any JavaScript function and could impact the security of your system.

However, since I control all the strings, I’m not worried.

To accomplish this, I pass a string with the user-defined function to TestOutput where I save it to a global funciton which can later be evaluated.

  // replacec x with fn_val which will be global x
  fn_exec :=  strreplace( fntext, 'x','fn_val', True);

Then I define a function which looks at that global string and evaluates it with CreateObject().

CreateObject returns a class, but we need the result as a double. We use casting, first to a variant to remove the class, then to double.

function fn ( x : double ) : double;
begin
  fn_val := x;

  result := double( variant(  CreateObject( fn_exec )));
end;

Actual output

Here is the actual output of my program.

I did the 1/x function twice, once with 1,000 steps, once with 10,000 steps. The difference was improving the accuracy by 0.019%, but taking several times longer to compute.

x^3 : [ 0, 1] (100) = 0.2500000008333334 should be  0.25 off by 0% , took 0 ms
1/x : [ 1, 100] (1000) = 4.701057474047187 should be  4.60517 off by 0.021% , took 0 ms
1/x : [ 1, 100] (10000) = 4.615037676928028 should be  4.60517 off by 0.002% , took 4 ms
x : [ 0, 5000] (5000000) = 12500000.000000333 should be  12500000 off by 0% , took 1490 ms
x : [ 0, 6000] (6000000) = 18000000.00000033 should be  18000000 off by 0% , took 1786 ms

Trial and error is usually sufficient. You can try a certain number of steps, and then try 1/10 as many, and see if the accuracy benefits of a large number of steps outweighs the slower calculation time.

unit integratefn;

interface

uses WebCore, Webdom, WebUI, WebForms, WebCtrls, WebEdits;

type

   TForm1 = class(TDialog)
      MultiLineEdit1: TMultiLineEdit;
      procedure Form1Show(Sender: TObject);
   private
      { Private declarations }

   public
      { Public declarations }
      procedure TestOutput( fnname, fntext : string; left, right : double ; n : integer ; answer : double );

   end;

var
   Form1: TForm1;

implementation

type
   TRealFunction = function( x : double ) : double;

var
   external fn_exec : string;
   external fn_val : double;


function integrate(  f : TRealFunction ; a, b : Double ; n : integer ):double;
var
   integral, h, cur: double;
   k: integer;
   i : integer;
   frac : double;
   sum1, sum2 : double;
begin
   integral := 0;
   h := (b-a)/n;

   sum1 := f(a + h/2);
   sum2 := 0;

   for i := 0 to n-1 do begin

      sum1 := sum1 + f(a + h * i + h/2);
      sum2 := sum2 + f(a + h * i);
   end;
   result := (h / 6) * (f(a) + f(b) + 4*sum1 + 2*sum2);
end;

function fn ( x : double ) : double;
begin
  fn_val := x;

  result := double( variant(  CreateObject( fn_exec )));
end;

// fnname is textual description
// fntext is JavaScript of function
// function over [a,b]
// n is number of steps
// answer is the expected true answer
procedure TForm1.TestOutput( fnname, fntext : string; left, right : double ; n : integer ; answer : double );
var
  s : string;
  res : double;
  start, stop : datetime;
  difftime : integer;
  pcentdiff : double;

begin
  if left > right then
      exit;

  // replacec x with fn_val which will be global x
  fn_exec :=  strreplace( fntext, 'x','fn_val', True);

  // record begin time
  start := now;
  res := integrate( fn, left, right, n );
  stop := now;
  difftime := Integer( stop )  - Integer( start );

  pcentdiff :=  round( 1000 *( res - answer) / answer) / 1000;

  s := fnname + ' : [ '+ FloatToStr( left ) + ', '+ FloatToStr( right ) + '] ('+IntToStr( n)+') = ' + FloatToStr( res ) +
    ' should be  ' + FloatToStr(answer) +
    ' off by ' + FloatToStr( pcentdiff ) + '% , took ' + IntToStr( difftime ) + ' ms';

  multilineedit1.Lines.add( s );
end;




procedure TForm1.Form1Show(Sender: TObject);
begin
  TestOutput( 'x^3', 'x*x*x', 0,    1,    100, 0.25 );
  TestOutput( '1/x', '1/x'  , 1,  100,   1000, 4.605170 );
  if isie then exit;
  TestOutput( '1/x', '1/x'  , 1,  100,   10000, 4.605170 );
  TestOutput( 'x'  , 'x'    , 0, 5000, 5000000, 12500000 );
  TestOutput( 'x'  , 'x'    , 0, 6000, 6000000, 18000000 );

{
  ƒ(x) = x3,  where   x   is     [0,1],       with           100 approximations.   The exact result is     0.25               (or 1/4)
  ƒ(x) = 1/x, where   x   is   [1,100],     with        1,000 approximations.   The exact result is     4.605170+     (natural log of 100)
  ƒ(x) = x,   where   x   is   [0,5000],   with 5,000,000 approximations.   The exact result is   12,500,000
  ƒ(x) = x,   where   x   is   [0,6000],   with 6,000,000 approximations.   The exact result is   18,000,000
}
end;

end.