Post Reply 
(50g) Simpson's rule for f(x,y)
11-17-2018, 09:27 PM (This post was last modified: 11-18-2018 11:14 PM by peacecalc.)
Post: #1
(50g) Simpson's rule for f(x,y)
Hello friends,

like Eddie Shore showed us HERE an algorithm for integration a function with two variables with the simpson rule and a matrix. He implemented this for the HP 71B.

I do the same thing for the HP 50G but not so elegant, it is brute force:

I implementated this formular:

\[ F = \int_a^b\int_c^d f(t,s)dtds \\ \sim \frac{ha}{3} \left( \frac{hi}{3}\left(
f(a,c) + f(b,c) + \sum_{j=1}^{k-1}\left( 2f(t_j,c) + 4f(t_j,c)\right) +
f(a,d) + f(b,d) + \sum_{j=1}^{k-1}\left( 2f(t_j,d) + 4f(t_j,d)\right) + \\
\sum_{i=1}^{m-1}\left(2\left(f(a,s_i)+f(b,s_i) +
\sum_{j=1}^{k-1}\left( 2f(t_j,s_i) + 4f(t_j,s_i)\right)\right) + \\ 4\left(f(a,s_i)+f(b,s_i) +
\sum_{j=1}^{k-1}\left( 2f(t_j,s_i) + 4f(t_j,s_i)\right)\right)\right)\right)\right) \]

That looks horrible, but I used the stack to sum up all function values and multiplied them afterwards with 2 or 4. And the indices in formular above has to be disdinguish between even or odd numbers (only values with odd indices has to be multiplied be 4 and even with 2). I used in the FOR loops no integer values but the values for the variables (the hp 50g is very happy to use a real variable in the FOR loop.

For instance I used my little program for estimate antiderivatives with harmonic sphere function multiplied with a light function to geht the coeffecients.The one angle goes from 0 to pi, the other one from 0 to 2pi. With N = 15 the hp 50g has to calculate 30*60 = 180 function values and it takes 2 minutes at average. That seems to be very long, but it is faster as you take the built-in function \[ \int \]. I have the impression that the built in function works then (when you have more variables) with recursion.

Code:

%%HP: T(2)A(R)F(.);
« \-> UA OA UI OI N    @@lower bound outer antiderivative
                @@higher bound outer antiderivative
                @@lower bound inner antiderivative
                @@higher bound inner antiderivative
                @@counts of partitial areas (N >= 2)
                @@otherwise the formular here is wrong 
                @@in between the bounds (same for both)
  
   «
    « DUP 2. * \-> TU TO SI H H2 @@ inner loop defined as a subprogram
      « TO SI FCSIM              @@load the function with two variables
                             @@with two input (numeric) values, each for                          @@every variable
    TU SI FCSIM              @@same again 

    0. DUP                  @@reserve two stackplaces for summation
                             @@the function values
    TU H + TO H2 - FOR K          @@init the for-loop
    K SI FCSIM + SWAP          @@save the value in the one stack place
                             @@and sum it up with the former ones
    K H + SI FCSIM + SWAP     @@save the value in the other stack place
                             @@and sum it up with the former ones
    H2 STEP                      @@increment the loop variable K with the 
                             @@double value 
    TO H - SI FCSIM +          @@calulate the last value which is multi-
                             @@plied later on with 4
    4. * SWAP 2. * + + +         @@multiply the sum of values in their 
                             @@stackplaces and sum it up
      »
    » 

    OI UI - N 2. * / DUP 2. *   @@calculate the delta H for the inner loop
                            @@and the double value of that
    OA UA - N 2. * / DUP 2. * @@calculate the delta H for the outer loop
                            @@and the double value of that
    
    \-> SFU HF HF2 HT HT2    @@local variables for the subprogram and the
                                @@delta values    
    
    « TICKS -2. SF -3. SF -105. SF @@ticks for time mesure (can be left out)
                               @@switch to numerical necessary for speed  
    UI OI UA HF SFU EVAL           @@load the subprogram with the boundaries
                               @@of the inner antiderivative, 
                                @@the value of outer antiderivative and the 
                                      @@delta h of the inner antiderivative
    UI OI OA HF SFU EVAL           @@same one line before 

    0. DUP                    @@reserve two stackplaces for summation
                               @@the function values
    UA HT + OA HT2 - FOR I      @@init the for-loop
    UI OI I HF SFU EVAL + SWAP    @@save the value in the one stack place
                             @@and sum it up with the former ones 
    UI OI I HT + HF SFU EVAL + SWAP @@save the value in the other stackplace                                                     
                                                   @@and sum it up with 
        HT2 STEP                            @@increment the loop variable K with the 
                                    @@double value 
    UI OI OA HT - HF SFU EVAL +       @@calulate the last value which is                                        
                                                       @@multiplied later on with 4
    4. * SWAP 2. * + + + HF * HT * 9. /   @@multiply the sum of values in 
                                                            @@their stackplaces and
                                        @@multiply it the delta h 
                                        @@for both antiderivatives
          TICKS ROT - 8192. /    @@can be left out only for time mesure 
    -2. CF -3. CF -105. CF    @@switch back to symbolic mode 
    »
  »
»
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
(50g) Simpson's rule for f(x,y) - peacecalc - 11-17-2018 09:27 PM



User(s) browsing this thread: 1 Guest(s)