+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (50g) Simpson's rule for f(x,y) (/thread-11818.html)
(50g) Simpson's rule for f(x,y) - peacecalc - 11-17-2018 09:27 PM
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
»
»
»RE: (50g) Simpson's rule for f(x,y) - Valentin Albillo - 11-17-2018 10:31 PM
.
Hi, peacecalc:
(11-17-2018 09:27 PM)peacecalc Wrote: 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 didn't see Eddie's post at the time but he's wrong on one count, namely (my bolding):
Eddie W. Shore Wrote:On the HP 71B, matrices cannot be typed directly, elements have to be stored and recalled on element at a time. The program presented does not use modules.
That's not correct. HP-71B's BASIC language allows for filling in all elements of an arbitrary size matrix at once by including the values in one or more DATA statements and then reading them all into the matrix using a single READ statement, no extra ROM modules needed. Thus, this lengthy initialization part in Eddie's code:
14 DIM I(5,5)
20 I(1,1) = 1
21 I(1,2) = 4
22 I(1,3) = 2
23 I(1,4) = 4
24 I(1,5) = 1
25 I(2,1) = 4
26 I(2,2) = 16
27 I(2,3) = 8
28 I(2,4) = 16
29 I(2,5) = 4
30 I(3,1) = 2
31 I(3,2) = 8
32 I(3,3) = 4
33 I(3,4) = 8
34 I(3,5) = 2
35 I(4,1) = 4
36 I(4,2) = 16
37 I(4,3) = 8
38 I(4,4) = 16
39 I(4,5) = 4
40 I(5,1) = 1
41 I(5,2) = 4
42 I(5,3) = 2
43 I(5,4) = 4
44 I(5,5) = 1
can be replaced by this much shorter, much faster version (OPTION BASE 1 assumed):
14 DATA 1,4,2,4,1,4,16,8,16,4,2,8,4,8,2,4,16,8,16,4,1,4,2,4,1
20 DIM I(5,5) @ READ I
20 DIM I(5,5) @ READ I
where the READ I fills in all data into the matrix with a single statement, no individual assignments or loops needed, thus it's much faster and uses less program memory.
Notice that this also works for arbitrary numerical expressions in the DATA, i.e.: the following hipothetical code would work Ok:
10 DATA 5,-3,2.28007e20,X,2*Z+SIN(Y),FNF(X+Y,X-Y),FNZ(2+FNF(C+D),3/FNF(6,8)),43
20 DIM M(2,4) @ READ M
20 DIM M(2,4) @ READ M
Anyway, Simpson's rule is suboptimal for integration purposes, either one-dimensional or multi-dimensional. There are much better methods providing either significantly increased accuracy for the same number of function evaluations or the same accuracy with fewer evaluations.
V.
RE: (50g) Simpson's rule for f(x,y) - peacecalc - 11-19-2018 08:33 PM
Hello Valentin,
I second your statement:
Quote:Anyway, Simpson's rule is suboptimal for integration purposes, either one-dimensional or multi-dimensional. There are much better methods providing either significantly increased accuracy for the same number of function evaluations or the same accuracy with fewer evaluations.
But as I wrote with double integrals the build-in solution for the hp 50g is slower then the brute force simpson rule.