NumericalIntegration

Using BASIC for Numerical Integration
By Tom Nally This article first appeared in the Liberty BASIC Newsletter Issue #105, February, 2003. It is reprinted here with the author's permission. In this article, all references to Liberty BASIC (LB) can be applied to Just BASIC (JB) as well.

Introduction: Numerical Integration
--- (WikiWriter, USA) One of the fun things that we learn during the first course in calculus is to find the area under a curve. By "area under a curve", we mean the enclosed area bounded by the following:
 * a function, f(x), on the top;
 * the x-axis on the bottom;
 * a real number, "A", creating a boundary on the left; and
 * another real number, "B", creating a boundary on the right.



Figure 01 above shows the area under the curve between A and B shaded in light blue.

In this article, I'm going to discuss a procedure for finding the area under a curve using a Liberty BASIC computer program that features "numerical integration". This procedure is called "numerical integration" because it peforms the integration (that is, calculates the area under the curve) without strictly relying on calculus. Rather, the procedure performs this duty by slicing the area up into hundreds of thin "elements", finding the area of each element, then adding them all together.

Before we get to that, however, we are going to find the area under a curve using traditional calculus. We are doing that for two reasons:

> 1. We want to demonstrate what traditional calculus looks like, so you can see how it differs from numerical integration; and > 2. We want to develop an answer to a calculus problem using traditional calculus as a way of checking the answer that we get using numerical integration.

If you haven't had a calculus course, the traditional calculus might not make much sense to you. Don't worry about that. Understanding traditional calculus is not a prerequisite for understanding numerical integration. So, just humor us for a few minutes, and we'll get to the LB program shortly.

Finding an Area using Traditional Calculus
--- Before we can find the area under a curve, we have to start by identifying a function. What is a "function", anyway? If I recall correctly, a function is a relationship between a set of x values and a set of y values such that each x value has one and only one y value. Any single y value can be associated with more than one x value, but the reverse cannot be true. That's how it was defined for me, anyway.

In this traditional calculus example, we will use the following function:



Here is what that same function looks like when graphed:



Let's find the area under this curve between code format="vb" x = 2 and x =8 code The first step in finding the area under a function in calculus is to find the "indefinite integral" of the function. The indefinite integral of a function is another function which is found by applying the rules of integration to the first function. For instance, the indefinite integral of (a*x^b) with respect to the variable "x" is (a/(b+1))*(x^(b+1)).

For the function in which we are interested, the indefinite integral is shown in Figure 04, below, using calculus notation:



To find the area, or "definite integral" of the function, we have to evaluate the indefinite integral between our two limits, code format="vb" x = 2 and x =8 code Figure 05 shows what that process usually looks like using calculus notation:



The result of our integration using traditional calculus is code format="vb" Area = 27.9 code With that out of the way, we will find the area using numerical integration, and come up with the same answer. OK, it might not be exactly the same, but it will be so close that the difference will probably be insignicant.

Finding the Area by Numerical Integration
--- First, let's establish the concept around which we are going to build our Liberty Basic program.



If you look at Figure 06 above, you will see the same function that was shown in Figure 01. In this one, however, a particular "element" of the area under the function is highlighted. Look at that highlighted element closely for a few moments. It resembles a shape with which we are quite familiar: the trapezoid. Of course, it is probably not a true trapezoid, because the function, f(x), is curved at the top of the element. However, as xa and xb get closer and closer together, the function between f(xa) and f(xb) becomes nearly a straight line! This collection of observations is the key to finding the area of this trapezoidal element, as well as the key to finding the area under the function.

So, what is the area of this trapezoidal element? That question was treated extensively in December, 2002, in the Yahoo LibertyBasic news group during the discussion of the area inside of a polygon. On several occasions, I offered that the area of the trapezoid was as follows: code format="vb" TrapArea = (xb - xa)*f(xa) + (1/2)*(xb - xa)*(f(xb) - f(xa)) code While that is correct, this formula has a simpler form that we will use here forward: code format="vb" TrapArea = (xb - xa)*(f(xa) + f(xb))/2 code If you look at this second formula, it's easy to see why that's the area of the trapezoid. (xb - xa) is the base of the trapezoid. Then, (f(xa)+f(xb))/2 is the average of the heights of the two sides of the trapezoid. Essentially, then, we are multiplying the base of the trapezoid by the average height between the left and right sides. Pretty straight-forward if you axe me (as we say in New Orleans).

Finding the area under our curve, then, involves slicing the shaded area into many small trapezoidal elements, and using the formula above to find the area of each one. When all the areas of the trapezoids are added together, then the result will be the area under the curve.

Friends, arithmetic such as this is child's play for our computers running Liberty BASIC.

The Numerical Integration Program
--- Before we get into all of the slicing and dicing of the shaded area, the first thing that we will do is define a Liberty BASIC function to hold our mathematical function. This will help make the program a little more modular: if we want to change mathematical functions, all we have to do is alter the contents of the Liberty BASIC function.

Let's make our Liberty BASIC function look like the code shown below. (Even though I'm discussing this function first, it will actually appear at the end of the code.) code format="vb" function f(x) f = 0.15*x^2 - 0.6*x + 3.45 end function code Let's get right into the display of the code, with comments included. code format="vb" '' 'Liberty BASIC program to find 'the area under a curve. ' 'by Tom Nally 'Released as Open Source ' ''' ' 'Note that the function, f(x), is defined 'within a Liberty BASIC function statement 'at the end of the program. ' 'First, dimension an area to hold the 'individual areas of each trapezoid

Dim TrapArea(10000)

'Establish the left and right boundaries 'for the area under the curve. ' A = 2 B = 8

'We will cut our area into 1000 'extremely thin elements, each of 'which will have the same width. 'we will call this width "dx".

NumElements = 1000

dx = (B - A) / NumElements

'Setup is complete. Time to find 'the area of each trapezoid. We will' 'use a for...next loop for this 'operation.

for i = 0 to (NumElements - 1) xa = A + i*dx xb = xa + dx   fxa = f(xa) fxb = f(xb) TrapArea(i) = (xb - xa)*(fxa + fxb)/2 next i

'The areas of all the trapezoidal elements 'have now been found. Now we will sum 'them all together in a second loop.

TotalArea = 0

for i = 0 to (NumElements - 1) TotalArea = TotalArea + TrapArea(i) next i

print "TotalArea = "; TotalArea

End

 function f(x) f = 0.15*x^2 - 0.6*x + 3.45

'If you want to find the area 'under a different function 'just change the definition 'above!

end function code If you copy this code into Liberty BASIC and run it, a single line of output will be produced: code format="vb" TotalArea = 27.9000054 code While this is not exactly the same result as we saw when we found the area under the curve using traditional calculus, you have to admit that it is pretty close to the exact result of 27.9.

Three Final Notes About Numerical Integration and the LB Program
--- >(1) You may slice the area under the curve into as many or as few elements as you wish by assigning different values to the variable NumElements. The more elements used, the more accurate the integration will be. Be aware, however, that high levels of accuracy may not be more useful than lesser levels of accuracy. Even dividing the area up into 10 elements (rather than 1000) produces an Area of 27.954, compared to the true area of 27.9. >(2) Finding areas by integration, whether by traditional calculus or by numerical integration, can sometimes produce values less than zero. Negative areas might occur when (a) all or part of the function lies below the X-axis, or (b) we integrate from the right-hand boundary to the left-hand boundary, instead of vice versa. This is not a flaw in the methods. This is the way it is supposed to work. >(3) Look again at the differential element highlighted in Fig 06. Yes, it's shaped like a trapezoid. But, as xb gets closer and closer to xa, the difference between f(xa) and f(xb) gets very small compared to either f(xb) or f(xa). Because f(xa) and f(xb) become nearly identical, the element might also be characterized as a tall, slender rectangle rather than a trapezoid. If that's the case, then the area of the element in the BASIC program can be written like this code format="vb" TrapArea(i) = (xb - xa)*(fxa) code instead of the way shown above: code format="vb" TrapArea(i) = (xb - xa)*(fxa + fxb)/2 code If you make that change in the BASIC program you will get these results... code format="vb" TotalArea = 27.8838054 code Compare this to the true value of 27.9. For some problems, this may be more than enough accuracy. The TrapArea formula which you select will depend on the tradeoffs you wish to make between accuracy and computing time. --- Tom Nally Steelweaver52@aol.com

--- Remember that all references to Liberty BASIC (LB) pertain equally as well to Just BASIC (JB).

--- Back to Nally's Code