Approximate Integration - execution of the program

This program will produce the definite integral of the function f(x)=x² between the limits x=A and x=B (B>A). On running the program, the display shows A? (the lower limit of integration). To begin with, enter 0 and then press [EXE]. The display then shows B? (the upper limit of integration). Enter 3 and press [EXE]. The display then shows D?. This is the number of pairs of strips that the interval between x=A and x=B is to be divided into. Enter 2 and press [EXE]. The display then shows M and 9.

This value is correct because the integral of f(x)=x² between the limits x=0 and x=a is (1/3)a³ and with a=3, this is (1/3)3³=9. This is a useful way to check that the code has been implemented correctly.

Let's try something harder. It is very difficult to obtain an expression for the integral of the function f(x)=exp(-x³) with respect to x but an estimate of the definite integral can be obtained numerically. Changing the function in the program from x² to exp(-x³) can be done by following the notes in this previous post. For the new function key in [SHIFT] [e^x] [(-)] [ALPHA] [X] [x³] [)].

Integrating this new function between the limits of x=0 and x=1 we obtain the following results as the number of pairs of strips is increased from 1 to 9:-

D=1 M=0.816311175
D=2 M=0.807843586
D=3 M=0.807572494
D=4 M=0.807530216
D=5 M=0.807518914
D=6 M=0.807514894
D=7 M=0.807513181
D=8 M=0.807512351
D=9 M=0.807511911

As can be seen, the results show that an estimate of this integral is 0.808 to 3 significant figures (as there is no change in the fifth significant figure for D>5) and this value agrees with other estimates.

Approximate Integration - further comments on the program

One of the other things to note about this program is that I have had to reuse several of the memory variables. The fx-50F comes with six memory variables (A, B, C, D, X and Y) and one independent memory (M) but this was not sufficient for the purposes of the program. For example, I calculate the value of the x-ordinate and store it in memory variable X, but then immediately calculate the y-ordinate, x², and store that in memory X too. This can be done because the x-ordinate is not needed after the y-ordinate has been calculated.

Approximate Integration - comments on the program

An important part of the program I presented previously is defining the function that is being integrated. I chose the simple function f(x)=x² and the statement in the code :X²→X: achieves this. In order to integrate a different funtion the program must be edited. For example, if I wanted to have f(x)=cos(x), then I would replace :X²→X: with :cos(X)→X:.

To edit a program it is best to ensure that your calculator is in insert mode by keying [SHIFT] [INS]. This will insert text before the vertical cursor and not overwrite existing text. To edit the program key [MODE] [6] [1] and then your program number. Use the replay button to move the blinking vertical cursor | to after X², so that it looks like :X²|→X:. Then press [DEL] twice and replace the function by keying [COS] [AlPHA] [X] [)], say.

It is a bit tedious, but I haven't found a way round this on the fx-50F.

Approximate Integration - the program

The logic for the approximate integration program was presented previously. Here is the logic converted into a program for the fx-50F:-

?→A:?→B:?→D:0→M:(B-A)÷(2D)→Y:1→C:While C≤D:1→B:Lbl 0:
B=1=>A+2(C-1)Y→X:B=2=>A+2CY→X:B=3=>A+(2C-1)Y→X:X²→X:
YX÷3→X:B=3=>4X→X:XM+:B+1→B:B≤3=>Goto 0:C+1→C:WhileEnd:
M▲

This is pretty much how the program looks when keyed into the calculator in edit mode except it all appears on one continuous line. Special Program Commands (? → => : Lbl While WhileEnd = ≤ ▲ Goto) are input using the [SHIFT] [P-CMD] keys. The cursor key (marked Replay) is used to switch between various screens of commands. Memory variables A, B, C, D, M, X and Y are input using [ALPHA] [A], [B] etc keys. The ":" Separator Code command is used frequently and can be alternatively input using the [EXE] key.

This program is 136 steps in length.

Approximate Integration - the logic

I previously gave a formula (Simpson's Rule) for finding the approximate integral of a function f(x) between x=a and x=b (b>a). This consisted of a sum of n bracketed terms corresponding to the approximate area under the curve of n pairs of strips, with each strip being of width h.

This sum provides the logic for calculating the area under the curve and hence a value for the definite integral as follows:-

Input A (the lower limit of the integration)
Input B (the upper limit of the integration)
Input D (the number of pairs of strips)
Store 0 in M (this will hold the sum of the area)
Store (B-A)/2D in Y (the value of h, the width of a strip)
Store 1 in C (C indicates which pair of strips is being processed)
While C is less than or equal to D
   Store 1 in B (B now indicates one of three ordinates)
   Label 0
   If B=1 then
      Store A+2(C-1)Y in X (left ordinate)
   EndIf
   If B=2 then
      Store A+2CY in X (right ordinate)
   EndIf
   If B=3 then
      Store A+(2C-1)Y in X (middle ordinate)
   EndIf
   Store f(X) in X
   Store YX/3 in X
   If B=3 then
      Store 4X in X
   EndIf
   Store M+X in M (the summation)
   Store B+1 in B
   If B is less than or equal to 3 then
      Goto Label 0
   EndIf
   Store C+1 in C
WhileEnd
Display M

Approximate Integration

Sometimes it is necessary to obtain a numerical estimate of the definite integral of a function f(x) between x=a and x=b (where b>a). This integral is represented by the area bounded by the curve y=f(x), the x-axis and the ordinates at x=a and x=b. One reasonably accurate method of estimating this area, and hence this integral, is to use Simpson's Rule. Here the area is approximated by

(h/3)(y₀ + 4y₁ + y₂) + (h/3)(y₂ + 4y₃ + y₄) + ... + (h/3)(y_2n-2 + 4y_2n-1 + y_2n)

The area is divided into 2n strips of width h=(b-a)/2n and the ordinates are given by y₀=f(a), y₁=f(a+h), y₂=f(a+2h) etc. The method uses a quadratic of the form AX²+BX+C to approximate the function between three adjacent points on the curve (i.e. between two adjacent strips).

Prime Factors - increasing efficiency

Previously I mentioned that for simplicity I have just incremented the prime factor variable X in my program by one each time, so when X is 3, the next prime factor to be tried is 4. As I said then, clearly 4 is not a prime number and so it is ineffecient to be using this number in the search for prime factors (the same is also true of 6, 8, 9, 10...etc.).

So, how could this be improved? Well, the answer is to add some logic like this after incrementing X:-

If X=4 then
   X=X+1
If-End

So this will mean that 4 is skipped and the next prime number that is tried is 5. A similar approach can be used when X=6 (the next prime being 7), but at X=8 we need to add 3 as the next prime number is 11.

I will leave it to the reader to make these adjustments to their code.