(Comp.sys.hp48) Option: Item: 1224 by colbach@nessie.cs.id.ethz.ch [Philippe Colbach] Subj: NUMERICAL LIBRARY (v1.2) Date: 29 Jun 1993 NUMERICAL LIBRARY ================= Title: NUM ROMID: 962 Bytes: 20593 CHKSM: # 4245h ERRATUM (documentation follows) ======= - I didn't mentioned in the documentation below that must have 2 arguments and (or <1>) ^^^^^^^^^^^ Combining [IPOL] with [SOLV], [SYST] and [INTG]: | [SOLV] | [SYST] | [INTG] ---------------------------------------------------------------------- Interpolation | predefined | predefined | predefined | | | Derivative | predefined | UNDEFINED | predefined | | ^^^^^^^^^ | Argument of | F(X) | F(X,Y) | F(X) | | ^^^^^^ | Arguments of | F(X) & Fx(X) | UNDEFINED | F(X) & Fx(X) | F(X) & 1 | ^^^^^^^^^ | F(X) & 1 Examples -------- !!!! Fcn has 1 argument: Fcn(F(X)), Fcn(F(Y)) or Fcn(F(X,Y)) !!!! derFcn has 2 arguments: derFcn(F(X),Fx(X)) or derFcn(F(X),1)) ^^^^^^^^^^^ ^^^^^ ^^^ 1) derFcn = d Fcn / d X = Fcn' --------------------------- Fcn = Fcn ( F(X) ) => d Fcn / d X = derFcn ( F(X) , Fx(X) ) -> Fx(X) = dF/dX ^^^^^^^ In fact, derFcn(F(X),Fx(X)) = Fx(X)*derFcn(F(X)) (command SYNTAX) (calculation FORMULA) ^^^^^^ ^^^^^^^ f.e.: Fcn(F(X))=Fcn(SIN(X)) => derFcn(F(X),Fx(X)) = derFcn(SIN(X),COS(X)) = COS(X)*derFcn(SIN(X)) Fcn = Fcn ( X ) => d Fcn / d X = derFcn ( X , 1 ) -> F(X) = X => Fx = 1 Fcn = Fcn ( X^2+2*X ) => d Fcn / d X = derFcn ( X^2+2*X , 2*X+2 ) -> F(X) = X^2+2*X => Fx = 2*X+2 try this: 'Fcn(F(X))' 'X' [derive] -> 'derFcn(F(X),derF(X,1))' -> derF(X,1)=Fx(X) 2) derFcn = function as is ----------------------- derFcn = derFcn(F(X),1) = 1 * derFcn(F(X)) ^^^ ^^^ ATTENTION: derFcn(F(X),0) = 0 * derFcn(F(X)) = 0 ^^^ ^^^ ^^^ 3) [SOLV] & [INTG] --------------- Both commands accepts arguments like: * 'Fcn(X^2+2*X)*SIN(X+pi)' * 'derFcn(X,1)' * 'derFcn(SIN(X^2),2*COS(X^2))' = 'derFcn(F(X),Fx(X))' = d Fcn / d X -> F(X)=SIN(X^2) => Fx(X)=2*X*COS(X^2) ^^^^^^^^^^^ * 'derFcn(SIN(X^2),1)' = 'derFcn(F(X),1)' = function as is -> = 1*derFcn(SIN(X^2)) ^^^^^^^^^^^^^^ try this: 'Fcn(SIN(X^2))' 'X' [derive] -> 'derFcn(SIN(X^2),2*COS(X^2))' ^^^^^^^^ Refer to [SOLV] and [INTG] for further information (-> NUMERICAL LIBRARY PART 01/02 published on 16 Jun 1993) 4) [SYST] ------ Systems of equations solved by [SYST] must be DERIVABLE: d Fcn(F(X,Y)) / d X = derFcn(F(X,Y),Fx(X,Y)) d Fcn(F(X,Y)) / d Y = derFcn(F(X,Y),Fy(X,Y)) The derivative of is defined by , but you can't solve systems including because its derivative is undefined! Refer to [SYST] for further information -------------------------------- DOCUMENTATION ============= FIRST Menu: [SOLV] Solving Equation [SYST] Solving System of Equations [INTG] Integrating Function [DIFF] Computing Differential Equation [ASIM] Continual Analogical Simulation SECOND Menu: [IPOL] Interpolation [CFIT] Curve Fitting [TRANS] Transforming Coordinates THIRD Menu: [REVä] Reversing äDAT Matrix [ORDä] Ordering äDAT Matrix [DRWä] Plotting äDAT Matrix (autoscaled) [SCLä] Plotting äDAT Matrix (scale plot) [LINä] Connecting scattered Points by Lines DOCUMENTATION ============= Special characters (used in this documentation) ----------------------------------------------- SQRT Square Root 131 S Integral 132 ä Sigma 133 pi Pi 135 a Alpha 140 [right-shift] [A] d Delta 146 [right-shift] [D] l Lambda 150 [right-shift] [L] r Rho 151 [right-shift] [R] Omega Omega 157 [right-shift] [O] Flag State: - 2 CF symbolic representation (IMPORTANT: [INTG]) ----------- (-2 SF: numerical representation) - 4 CF unused system flag -55 CF LASTARG disabled (probably WRONG LAST ARGUMENTS left on the stack in case of errors otherwise) (-55 SF: LASTARG enabled) -60 SF ALPHA in case of alpha mode (-60 will be set by [Y=SF(xy)dx], [Y=SS[>dt], [CFIT] and [TRANS]) (-60 CF: 1ALPHA mode) ALPHA: Alpha lock by pressing [a] twice 1ALPHA: Alpha lock by pressing [a] once The following programs have been tested for these flag states. [SOLV] SOLVING EQUATION ======================== INPUT: 2: 'F(X)' or 'EQ(X)' 1: { GuessX1 | GuessX2 GuessX3 | } or GuessX OUTPUT: 1: X: Root DISPLAY: running: current X (after pressing any key) root found: root interpretation ERROR: no special error message Both (internal or user-defined) and can be used as argument. You can enter either one guess or a list of maximum 3 guesses. [SOLV] uses the internal procedure of the [SOLVR] menu. No variable ('EQ' or 'X') is left in the current directory. N.B. - 'X' has to be the unknown variable (only 2 inputs!) ^^^ ^^^^^^^^^^^^^^^^ - Changing variable: 'F(Z)' { Z X } | -> 'F(X)' |, WHERE command The interpolation command (-> [IPOL]) is represented by the 'Fcn(X)' and 'derFcn(X)' functions. Refer to [IPOL] for further information. EXAMPLE: 2: 'SIN(X)' 1: { 3 4 } RAD Mode [SOLV] -> 1: X: 3.14159265359 DISPLAY: Sign Reversal [SYST] SOLVING SYSTEM OF EQUATIONS =================================== INPUT: 4: 'F1(X,Y)' or 'EQ1(X,Y)' 3: 'F2(X,Y)' or 'EQ2(X,Y)' 2: GuessX 1: GuessY OUTPUT: 1: X: RootX 2: Y: RootY DISPLAY: 1) message displayed during derivation 2) current X and Y displayed during computation ERROR: no special error message Both (internal or user-defined) and can be used as argument. But you can't enter lists of guesses. [SYST] uses the Newton Algorithm: both curve functions are derived by the program (message displayed during derivation). No variable is left in the current directory. N.B. 'F1(X,Y)' (or 'EQ1(X,Y)') and 'F2(X,Y)' (or 'EQ2(X,Y)') must be DERIVABLE! 'X' and 'Y' have to be the unknown variables. ^^^^^^^^^ ^^^^^^^^^^^ ^^^^^^^^^^^^^^^^^ The interpolation command (-> [IPOL]) is represented by the 'Fcn(X)' and 'derFcn(X)' functions. The computation time of systems including interpolation 'Fcn(X)' and its derivative is very long because 1) the computation of 'Fcn(X)' and 'derFcn(X)' takes a long time and 2) the derivation curve 'derFcn(X)' has a wave form. Refer to [IPOL] command for further explanations. EXAMPLE: 4: 'SQ(X)+SQ(Y)=1' circle 3: 'Y=SQRT(X)' 2: .5 1: .5 As no root interpretation procedure is included, the program doesn't interrupt the computation at a sign reversal or in case of no root! The program can be interrupted any time by pressing [ON]. [SYST] -> 2: X: .618033988751 1: Y: .786151377759 [INTG] INTEGRATING FUNCTION ============================ INPUT: 2: 'F(X)' 1: { LowerLimit UpperLimit } OUTPUT: 1: S: Integral DISPLAY: message displayed in case of numerical computation (computation time can be long!) ERROR: no special error message [INTG] first tries an algebraic integration, then a numerical one (message displayed). No variable is left in the current directory. In case of numerical computation, accuracy is defined by current display format (STD, FIX, ENG, SCI; STD allows the best accuracy, but increases the length of calculation time). LowerLimit and Upperlimit can be algebraic expressions including pi, e or user-defined constants (f.e. '2*pi' or 'a+e',a=1). N.B. - Flag -2 must be cleared! ^^^^^^^ ^^^^^^^ - 'X' has to be the unknown variable (only 2 inputs!). ^^^ ^^^^^^^^^^^^^^^^ Interpolation: 'Fcn(X)' and 'derFcn(X)' (increases the lenght of calculation time!). Combining [INTG] and [SOLV]: << -> A << 'LN(X)' 1 A 2 ->LIST INTG 1 - >> >> 'FCN' STO 5 FIX 'FCN(X)' { 2 3 } [SOLV] EXAMPLE: 2: 'SIN(X)' 2: 'SIN(ABS(X))' 1: { 0 '2*pi' } 1: { 0 '2*pi' } RAD and STD mode selected (algebraic computation) (numerical computation) [INTG] -> 1: S: 0 1: S: -3.67374857952E-12 EXAMPLE: 2: 'LN(X)' 1: { 0 e } (algebraic computation) [INTG] -> 1: S: 0 This calculation would be impossible by using numerical computation! [DIFF] COMPUTING DIFFERENTIAL EQUATION ======================================= [DIFF] activates a special menu: [Y=SF(xy)dx] Computing Y'=f(X,Y,l), Y'=dY/dX (label 1 & 2 & 3) [Y=SS[>dxdx] Computing Y"=f(X,Y,d), Y"=d^2Y/dX^2 (label 4 & 5 & 6) | [Y=SF(xy)dx] | [Y=SF(xy)dx] | [Y=SS[>dt] | (first algorithm) | (second algorithm) | ------------------------------------------------------------------------------- computes | Y' | Y' | Y" | | |^^^^ variables | X, Y | X, Y, l | X, Y, d | | ^^^ | integr.var. | dx | dx | dx | | | algorithm | Runge-Kutta | (Heun) | Analogical | | | order | 4th | 2nd | 1st |^^^^^ | | accuracy: | flat: very good | flat: good | flat: suffic. | ^^^^^^^^^ | | | +/-90: bad | +/-90: good | +/-90: very bad | | ^^^^ | init.val.= | X0, Y0 | X0, Y0 | Y0, d0 | | | comp.step dH= | +/-dX | +/-SQRT(dX^2+dY^2) | +/-dX | | | stor.step dS= | +dS | +dS | +dS | | | dX= | con=dH | var=dH*COS(a) | con=dH | | ^^^^^^^^^^^^^ | sim.interval= | [ Xb Xe ] | No | [ X0 Xe ] | | | RAD mode | no | yes | no | | | äDAT | appends | appends | overwrites | ^^^^^^^ | ^^^^^^^ | [Y=SF(xy)dx] Computing Y'=f(X,Y,l), Y'=dY/dX ============================================== First EXAMPLE vvvvvvvvvvvvv [Y=SF(xy)dx] -> |Enter Differential EQ | Y'=f(X,Y.l): Algebraic enter mode is activated. 1ALPHA in case of alpha mode. l (lambda; [right-shift] [L]) represents the length of the curve arc. EDIT menu is activated (-> [->STK] command). ^^^^^^^^^^^^^^^^^^^^^^^ |'TAN(SQ(l)/2) l = [right-shift] [L] |Enter Initial Values | X0 and Y0: l0 (lambda0)=0 by default! EDIT menu is now desactivated ([left-shift] [+/-] ^^^^^^^^^^^^^^ to reactivate). | 0 0 |Enter Computation step | +/-dH: | .02 |Enter Storage Step | +dS: |----:---------------------/ /--------------------------> curve arc [l] l0 ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ computation step dH=0.2 ^ ^ ^ ^ ^ ^ ^ ^ storage step dS=3*dH=0.6 0 1 2 n-2 n-1 n n=number of steps Choosing a computation step |dH| < than dS allows a greater accuracy of stored points. Data points are stored into 'äDAT'. If 'äDAT' already exists, points are appended! ^^^^^^^^ | .06 Every third ( ABS(IP(dS/dH))=ABS(IP(.06/.02))=3 ) point is stored into 'äDAT'. |>ae]-90 90[ >ae[0 360] | |dx=h dx=h*cosa | Algorithm formulae |dy=h*tana dy=h*sina | |l=undefined l=Eh > l=current length of curve arc * |>4th Order >2nd Order > Order of computation ** |[ ] /\ [ ] [ ] /\ [ ] Press [/\] to choose algorithm * First algorithm is unable to compute the current length of the curve arc. The computation time would have been too long! ** The first algorithm comptes at a better accuracy in case of FLAT curves, but fails in case of GREAT gradients (around +/-90 degrees) because dX=dH is constant (choose the second algorithm: dX=dH*COS(a)!). A 4th order computation in case of a variable dX would have been too difficult to program and the calculation too long, too! IMPORTANT: During computation of the second algorithm, RAD mode is automatically selected for programming reasons. ^^^^^^^^^^^^^^^^ ^^^^^^^^ | second algorithm (able to compute l) |Enter Number of Steps | No: | | 80 80 points ( [ Xn Yn ] ) are stored into 'äDAT' (240 points are computed). The result curve is then drawn as a scale curve. If an error occures during computation, the program interrupts without drawing a curve. Second EXAMPLE vvvvvvvvvvvvvv 'äDAT' PURGE Data points would be appended to the former computation [Y=SF(xy)dx] -> |'-SIN(X) Y'=-SIN(X) Differential equation | 0 1 X0=0 and Y0=1 Initial values | .02 dH=.02 Computation step | .1 dS=.1 Storage step |------------------------------------------------------> [X] ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ computation step dH=.02 ^ ^ ^ ^ ^ storage step dS=5*dH=1 | First algorithm (dH=dX) IMPORTANT: During computation of the first algorithm, RAD mode is *NOT* selected automatically as during the second one: ^^^^^^^^^^^^^^^^^^^^^^^^^^ [left-shift] [1] to select RAD. |Enter Storage Interval | [Xb Xe]: | 0 6.3 X0 <= Xb < Xe vvvvvvvvvvvvv |-----:---------------:----------:------------> [X] X0 Xb Xe |--------------------------> computation interval |----------> storage interval X0 >= Xb > Xe vvvvvvvvvvvvv [-X] <------:--------------:-----------:-------| Xe Xb X0 <--------------------------| computation interval <--------------| storage interval These data points are used by the second example of [IPOL]. [Y=SS[>dt] Computes Y"=f(X,Y,d), Y"=d^2Y/dX^2 ================================================= The [>-block represents the analogical algorithm. [Y=SS[>] uses the [ASIM] programs to compute the differential equation of the 2nd degree. [ASIM] does continual simulation => time T = independent variable! But for the computation of differential equations with [Y=SS[>dt], T = X ! [Y=SS[>dt] defines 3 blocks and combines them to an analogical system: Y" Y' Y ---[INTEGRATOR 50>-------[INTEGRATOR 50>--- | | | | | | X --dt] -> |Enter Differential EQ | Y"=f(X,Y,d): | 2 Y"=2 => Y'=2*X => Y=X^2 |Enter Initial Values | Y0 and d0: d0=Y'0 | 0 0 | Compiling System [ASIM] compiler |Simulation Interval | [Tb Te]: Simulation interval [Tb Te] = computation interval [Xb Xe] (N.B. [ASIM] does continual simulation -> time T = independent variable!)! Tb=Xb=X0! Tb=Xb can be > than Te=Xe -> negative computation step dT=dX. ^^^^^^^^ ^^^^^ ^^^^^ ^^^^^ Xb=X0 < Xe vvvvvvvvvv |--------:-----------------------:--------------> [X] X0 Xe |-----------------------> computation interval = storage interval Xb=X0 > Xe vvvvvvvvvv [-X] <--------:----------------------:-------| Xe X0 <----------------------| computation interval = storage interval | 0 1 |Enter Computation Step | +dT: Again: dT=dX! (remember that!) | .01 |Enter Storage Step | +dS: |---------------------------------------------------> [X] ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ computation step dX=.01 ^ ^ ^ storage step dS=10*dH=.1 | .1 100 points are computed, 10 are stored in 'äDAT'. Former 'äDAT' points are cleared (data points are NOT APPENDED! -> [ASIM]). Result ('äDAT' variable): X Ycalc Ycorr [[ 0.0 0.000 ] 0.000 [ 0.1 0.011 ] 0.010 [ 0.2 0.042 ] 0.040 [ 0.3 0.093 ] 0.090 [ 0.4 0.164 ] 0.160 [ 0.5 0.255 ] 0.250 [ 0.6 0.366 ] 0.360 [ 0.7 0.497 ] 0.490 [ 0.8 0.648 ] 0.640 [ 0.9 0.819 ] 0.810 [ 1.0 1.010 ] 1.000 [ASIM] CONTINUAL ANALOGICAL SIMULATION ======================================= [ASIM] activates special menu: [()-[> ] DEFINE, RECALL or CLEAR Blocks [()-()-] CHANGE Parameter of Blocks [<< >> ] COMPILE analogical Structure [ E [> ] STORE Blocks during Computation [ o-[] ] SIMULATE analogical Structure [|__/\_] DRAW Block Curve These commands allow to define blocks and to combine them to an analogical structure. DEFINITION (Schaum's outline series: Theory of feedback and control systems): A system is an arrangement of physical components connected or related in such a manner as to form and/or act as an entire unit. Control systems are classified into two general categories: OPEN-LOOP and CLOSED-LOOP systems. The distinction is determined by the CONTROL ACTION, which is that quantity responsible for activating the system to produce the output: - an OPEN-LOOP control system is one in which the control action is independent of the output. - an CLOSED-LOOP control system is one in which the control action is somehow dependent on the output. Both OPEN-LOOP and CLOSED-LOOP control systems can be combined. ATTENTION: An CLOSED-LOOP control system MUST ALWAYS contain at least ONE INTEGRATOR (block type 50,51 or 52)! Control systems simulated by the [ASIM]-Program MUST ALWAYS include at least ONE INTEGRATOR, even in case of an OPEN-LOOP! Otherwise the [<< >> ]-compiler errors ("Integrator missing" or "Short Circuit in **"; refer to [<< >> ]). EXAMPLES vvvvvvvv INPUT OUTPUT-< ]-INPUT -[I>- INTEGRATOR [C>- CONSTANT | | INPUT-[ >-OUTPUT INPUT (t)- CURRENT TIME -[Q> QUIT (t)---[ >---[ >---[I>---[ >---[ > (t)---[ >---[ >---[I>---[ > | | | | -------[ >---- ----< ]---- -------< ]---- -[Q> (t)---[ >---[ >---[ >---[I> | | | | | (t)---[ >---[ >---[I>---[ >- --[ >--- | | | --< ]--- | [C>---[I>---[ >---[ >---[ > | | | | ---------- [C>---------- N.B. - In fact, simulating an analogical structure means to compute a (system of) differential equation(s). ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - [()-[> ], [<< >> ] and [ E [> ] create several variables in the current menu (named 'Omega***'). Before defining an analogical system, create a new directory! (Character 'Omega' = [right-shift] [O]) ^^^^^^^^ ^^^^^^^^^^^^^ [()-[> ] DEFINE, RECALL or CLEAR Blocks ======================================= DEFINE BLOCK CLEAR BLOCK RECALL BLOCK vvvvvvvvvvvv vvvvvvvvvvv vvvvvvvvvvvv INPUT: 4: Block Number 3: Block Type 2: { Entries } 2: Block Number 1: { Parameters } 1: 0 1: Block Number OUTPUT: none none 4: Block Number 3: Block Type 2: { Entries } 1: { Parameters } DEFINE BLOCK ------------ Every block has several inputs (from other blocks or parameters), but only ONE output. A block is defined by its block number (its "name"), its block type (defining its action), its entries (from other blocks) and its parameters. There are special block types without any entry (56 and 60). Block number 100 represents the current time during computation. [()-[> ] always owerwrites former block definitions. 50 different block types are predefined (refer to block list at the end of this documentation), some types are undefined (f.e. 23). The creation of user-defined block types is possible: A user-defined block type consists of a program stored under global name 'OmegaNNN' (NNN > 100). This program must compute the exact number of block entries (E) and parameters (P) you indicated in the { Entires } and { Parameters } lists during definition. The count of E's and P's is not tested as for the predefined block types! The program must only have ONE output E0! RECALL BLOCK ------------ To recall a block definition, just enter its number. Block number, block type, entry list and parameter list are put on stack. CLEAR BLOCK ----------- To clear an already defined block, enter its number and a zero. In fact, blocks are not cleared, but replaced by the block type 60 (Null). This to avoid a block entry from an undefined (cleared) block (refer to [<< >> ] errors). Note: - user-defined block type program form: << -> E1 E2 P1 P2 P2 << ... >> >> Block entries (E) BEFORE parameters (P)! EXAMPLE: block type: 101 action: E0=E1^P1+E2^P2 E: 2 P: 2 'Omega101': << -> E1 E2 P1 P2 << E1 P1 ^ E2 P2 ^ + >> >> << -> E1 E2 P1 P2 'E1^P1+E2^P2' >> >> << ROT SWAP ^ 3 ROLLD ^ + >> is of course possible! !!!! - A negative block entry means that the sign of the input from this !!!! block is automatically changed (refer to the example). - { } means NO entry or parameter (E=0 or/and P=0 in block list). - Character Omega = [right-shift] [O] !!!! - AFTER every change, the system has to be compiled BEFORE it can be !!!! computed as new system! VARIABLES: OmegaBLC block entries OmegaPAR block parameters OmegaTYP block types OmegaINT Integrators OmegaDER d/dt blocks OmegaINIT initial values of integrators ERRORS: "Wrong Entry Count" Refer to block list under column E "Wrong Parameter Count" Refer to block list under column P "Type ** undefined" Refer to block list under column NBR & ACTION "Block ** undefined" Trying to recall an undefined block [()-()-] CHANGE Parameter of Blocks =================================== INPUT: 2: Block Number 1: { Parameters } or Parameter (Initial Value of Integrator) OUTPUT: none This command allows to change the parameters of an ALREADY defined block. To simply change the initial value of an integrator (block number 50,51,52), just enter its block number and the initial value (not included in a list). N.B. AFTER every change, the system has to be compiled BEFORE it can be simulated! ERRORS: "Wrong Parameter Count" Refer to block list under column P "Block ** undefined" Trying to change the parameters of an undefined block [<< >> ] COMPILE analogical Structure ===================================== INPUT: none OUTPUT: none Before a system can be simulated, it has to be compiled. The compiled system (program) is stored into global name 'OmegaPRG'. IMPORTANT: AFTER every change (block definition or just a parameter), the system has to be compiled BEFORE it can be computed! VARIABLES: OmegaORD order of block output computation OmegaPRG compiled system ERRORS: "Integrator missing" No integrator at all was defined "Short Circuit in **" Closed-loop (containing block **) without integrator "Block ** undefined" Block entry from an undefined block [ E [> ] STORE Blocks during Computation ======================================== INPUT: 1: { Block Numbers } or { 0 } or BlockNumber OUTPUT: none Not every block output is important. You can pick out all the outputs you are interested in. { 0 } means that ALL blocks are stored during computation. Remember that block number 100 represents the current time. [ E [> ] creates a program stored into global name 'OmegaSTO'. [ E [> ] changes don't have to be compiled, but the storage list has to be defined before the system can be simulated! To create a 'äDAT' matrix compatible with [IPOL], [CFIT] and [TRANS] etc (size = { n 2 }), enter { 100 BlockNumber } as storage list (X: Time, Y: Bock) or just BlockNumber (means { 100 BlockNumber }). VARIABLES: OmegaSTP list of blocks to be stored during simulation OmegaSTO program to store blocks during simulation { 100 5 3 } [ E [> ] -> OmegaSTP = { 100 5 3 } OmegaSTO = << OmegaTIM OmegaVAL 5 GET OmegaVAL 3 GET 3 ->ARRY ä+ OmegaTIM 7 DISP >> You can of course change that program or even replace it by an user-defined program (f.e. which prints the current block values) EXAMPLE: OmegaSTO = << " " 1 OmegaTIM 0 FIX ->STR REPL 5 OmegaVAL 5 GET 3 FIX ->STR REPL 10 OmegaVAL 3 GET 1 FIX ->STR REPL PR1 DROP OmegaTIM 7 DISP ^^^^ >> prints without storing results! N.B. OmegaSTO must not have an output! ^^^^^^^^^^^^^^^^^^^^^^^ ERRORS: no special error message [ o-[] ] SIMULATE analogical Structure ====================================== INPUT: none OUTPUT: none [ o-[] ] -> |Simulation Interval: | [Tb Te]: Enter the time interval (f.e. from 0 [?] to 100 [?], unit [?] depends always on the current system): Tb < Te vvvvvvv |------:--------------------:-----------> [T] Tb Te |--------------------> computation interval = storage interval Tb > Te vvvvvvv [-T] <-------:--------------------:----------| Te Tb <--------------------| computation interval = storage interval |0 100 |Enter Computation Step: | +dT: Enter the time step of computation (f.e. 0.1 [?]): |.1 |Enter Storage Step: | +dS: Enter the time step to store the blocks indicated by the command [ E [> ]: |--------------------------------------------------------> [T] ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ computation step dT=.1 ^ ^ ^ storage step dS=10*dT=1 |1 | SIMULATING Simulation begins ... Blocks outputs are stored into 'äDAT' matrix. During simulation, the program uses special variables (which can be used by user-defined block types): - OmegaTIM: current simulation time - OmegaCNT: step count 1st step=initial time Tb 2nd step=Tb+dT, dT=time step number of steps = (Te-Tb)/dT - OmegaDT: time step dT ERRORS: "Structure not compiled" Trying to simulate a non-compiled system. Refer to [<< >> ] command. "Storage List missing" Trying to simulate without indicating the blocks to be stored during computation. Refer to [ E [> ] command. [|__/\_] DRAW Block Curve ========================= INPUT: 1: Block Number or 0 OUTPUT: none Just enter the number of the block of which you want to draw the result curve (Plot Type: BARPLOT). If you entered { 100 BlockNumber } or BlockNumber as storage "list" (means 2 "blocks" only!) (refer to [ E [> ], 0 as BlockNumber produces a SCATTER drawing. ERRORS: "Block ** undefined" Trying to draw the result curve of an undefined block "Block ** not stored" Trying to draw the result curve of a block which was not included in the storage list (refer to [ E [> ]) BLOCK TYPE DEFINITIONS at the end of documentation ---------------------- FIRST EXAMPLE ------------- __ \___ affluent R [m^3/h] \___________ weir W [m^3/h] \ area A [m^2] ___ \~~~~~~~~~~~~~| | \ ^ \ depth H | | | \ [m] | | | initial depth of water H0=20m \ | | | \_________| | v | |dam: heigth=21m, width=4m Area(H) [m2] -------- ------|Omega101>--------> -- [m3/h] -- [m3/s] -- [m/s] ---- [m] | -------- (R) (+) |15>------->|20>------->|18>--->(1)| 50 >-----| *1* ---> -- -- -- (dH) ---- (H) | |(-) *2* *3* *4* | *5* | | (20)=(H0) | *BlockNumber* | | --------- | Wire(H) *7* | |BlockType> | [m3/h] -------- | --------- ------------------------------------ outlet ). Block Definition: NBR TYPE ENTRIES PARAMETERS 1 101 { 100 } { } [()-[> ] 2 15 { 1 -7 } { } [()-[> ] Note input signs! 3 20 { 2 } { 3.6 } [()-[> ] 4 18 { 3 6 } { } [()-[> ] 5 50 { 4 } { 20 1 } [()-[> ] 6 102 { 5 } { } [()-[> ] 7 103 { 5 } { 21 4 } [()-[> ] User-defined block types: Block type 101: (Affluent) << -> E1 << ... >> >> 'Omega101' STO Block type 102: (Area: function of depth H) << -> E1 << ... >> >> 'Omega102' STO Block type 103: (Wire: function of depth H, dam height and dam width) << -> E1 P1 P2 << ... >> >> 'Omega103' STO Compiling: [<< >> ] Storage List: { 1 5 7 } [ E [> ] Simulating: [ o-[] ] -> |0 24 Time Interval |.1 Comp. Step |.4 Storage Step Draw Curve: 5 [|__/\_] Depth Curve SECOND EXAMPLE -------------- Differential Equation Y'''=f(Y,Y',Y'',t), Y'''=d^3Y/dt^3 Y''(0)=A3 Y'(0)=A2 Y(0)=A1 Analogical System: (A3,1) (A2,1) (A1,1) | Y'' | Y' | Y ---[INTEGRATOR 50>-------[INTEGRATOR 50>-------[INTEGRATOR 50>--- | *1* | *2* | *3* | | | | | T --] Integrator 1: E0=Y'' 2 50 { 1 } { A2 1 } [()-[>] Integrator 2: E0=Y' 3 50 { 2 } { A1 1 } [()-[>] Integrator 3: E0=Y 4 101 { 1 2 3 100 } { ? } [()-[>] Differential: E0=Y''' !!!! You should ENTER (not: number; not important) the INTEGRATORS in THAT !!!! particular order (E0=Y'' -> E0=Y' -> E0=Y). If you "connect" several !!!! integrators in series, you should always ENTER the INTEGRATORS in "flow !!!! direction" order. This is due to programming and algorithm reasons. Omega101: << -> Y2 Y1 Y0 T | Pn | << ... >> >> Compile and simulate (from T0=0 to T=TE) the system ... NOTE: this method corresponds to the [Y=SS[>dt] algorithm (dT=const; 1st order computation). This analogical system fails in case of great gradients (around +/-90 degrees). Refer to [Y=SS[>dt] for more information about differential equations. ACCURACY -------- Computing differential equations either with the Runge-Kutta (4th order) or Heun (2nd order) algorithm is very accurate, but integrating these equations step by step with an analogical system is not! Calculation blocks as [COSINE>, [LOGARITHM>, [MULTIPLICATION> or [LIMIT> (blocks 1-22, 25-33 and 36-47) are as accurate as the HP48 itself, but accuracy of [INTEGRATION>-blocks depends on their algorithm formula => on dt: (I0) | E----[I>----O, E=F(t) (open-loop) or E=F(t,O) (closed-loop) - at t=0, I=I0 m=n - In+1 = In + E * dt (first order) = I0 + ä dI, dI=E*dt m=0 m=n n - dt ~ 0 => ä dI ~ S dI, ä Sigma, S Integral M=0 0 => choosing dt<<1 -> greater accuracy but: increases length of time to calculate! [IPOL] [CFIT] [TRANS] [REVä] [ORDä] [DRWä] [SCLä] [LINä] ============================================================= ****************************** * I M P O R T A N T * ****************************** All these commands need 'äDAT' matrix under special form: - 'äDAT' size = { n 2 } - First column of 'äDAT' = Xn - Second column of 'äDAT' = Yn - [IPOL] ONLY: X-column in GROWING order, Xn > Xn-1! ^^^^^^^^^^^ ^^^^^^^ [DIFF] commands create 'äDAT' in growing order if: -dX=dH > 0 (first algorithm of [Y=SF(xy)dx]) -dX=dH*COS(a) > 0 (second algorithm of [Y=SF(xy)dx]) -dT=dX > 0 ([Y=SS[>dt]) Otherwise: [REVä] to reverse 'äDAT' [ORDä] to order 'äDAT' [IPOL] INTERPOLATION ===================== INPUT: 1: RealX OUTPUT: 2: d: Derivative 1: Y: RealY DISPLAY: no special message displayed ERROR: "Bad äDAT Data" äDAT either doesn't exist or is not a matrix "Wrong äDAT Size" äDAT size is different from { n 2 } [IPOL] lets you interpolate data points stored in GROWING order in 'äDAT': [Y] ^ (m) (1) | x x (2) (m-1) x| x :| x x | :| x : :| : --------:-----/ /--------:-------> [X] :| : : : -------><-----------------><-------- 2nd 5th 2nd polynomial degree (of interpolation) 1st 4th 1st derived polynomial degree Interpolation between (Xn,Yn) and (Xn+1,Yn+1): (1) F(Xn)=Yn (2) F(Xn+1)=Yn+1 (3) F'(Xn)=(Yn+1 - Yn-1)/(Xn+1 - Xn-1) (4) F'(Xn+1)=(Yn+2 - Yn)/(Xn+2 - Xn) (5) F"(Xn)=0 (6) F"(Xn+1)=0 6 conditions => IPOL(X), polynomial of 5th degree => derivation DER(X), polynomial of 4th degree and: (5) & (6) => DER'(Xn)=DER'(Xn+1)=0 => wave form DER(X) curve: __._---_.__ __.__/ | | \__.__-- __.__/ | | | | __/ | | | | | | | | | | -----:-----:-----:-----:-----:-----> [X] First EXAMPLE vvvvvvvvvvvvv Let's create an apropriate 'äDAT' matrix: 'äDAT' PURGE ([Y=SF(xy)dx] appends data points to former computation) [DIFF] [Y=SF(xy)dx] -> | l Y'=l (lambda) | 0 1 X0=0, Y0=1 | .1 dH=.1 | .1 dS=.1=dH (every point is stored) | second algorithm (able to compute l) | 10 10 points Differential equation Y'=l [Gradient = Function(CurveArc)] => Y=COSH(X) !! That's why I proposed the initial values X0=0 and Y0=1! [IPOL] interpolates the data points from the 'äDAT' matrix: .5 [IPOL] -> d(.5)/dX= 0.510410896178 | SINH(0.5)= 0.521095305494 Y= 1.127241126120 | COSH(0.5)= 1.127625965210 These data points are also used for the example of [CFIT]. Second EXAMPLE vvvvvvvvvvvvvv Use the data points from the second example of [Y=SF(xy)dx]. As Y'=-SIN(X), Y=COS(X). Let's try to find the root X=pi/2 by using [IPOL] and [SOLV]. As [IPOL] has TWO tagged outputs, we would have to create a program (function form) which has ONLY the (real) output Y: << -> X << X IPOL SWAP DROP 0 + >> >> 'Fcn' STO But during execution of [SOLV] ([SYST] and [INTG]), the interpolation functions 'Fcn(X)' and 'derFcn(X)' are predefined! Refer to [SOLV] and [SYST]. 2: 'Fcn(X)' 1: { 1 2 } [SOLV] -> 1: X: 1.57077308547 and: pi/2= 1.57079632680 Note: -Program with Y output (as real): << -> X << | function form X IPOL | calculates 2: d: and 1: Y: SWAP DROP | drops d: 0 + | Y: (tagged) -> Y (real) >> >> | -Program with dY/dX output (as real): << -> X << | function form X IPOL | calculates 2: d: and 1: Y: DROP | drops Y: 0 + | d: (tagged) -> d (real) >> >> | !!!! Internal ROOT and DRAW commands don't accept tagged arguments; neither do !!!! [SOLV] and [INTG]! [CFIT] CURVE FITTING ===================== INPUT: none OUTPUT: 1: Correlation ERROR: "Bad äDAT Data" äDAT either doesn't exist or is not a matrix "Wrong äDAT Size" äDAT size is different from { n 2 } [CFIT] fits the data points from the 'äDAT' matrix. HP48 Curve Fitting program offers only 4 base modles, but [CFIT] allows user modles: For this example, we use the data points from the first example of [IPOL]. [CFIT] -> | Curve Fitting Modles |LIN: B + f(X)*M = g(Y) |PWR: B * f(X)^M = g(Y) | |[LIN] [] [] [] [] [PWR] [CFIT] proposes 2 base modles: LIN and PWR - f(X) can be ANY function of X like f(X)=COS(X)^2+SIN(X)^2 - in definition of g(Y), the variable Y can only appear ONE time because Y is isolated by the internal [ISOL] command: f.e. g(Y)=SQRT(1+SQ(Y)) ^ |[LIN] |Modle: LINFIT |B + f(X) *M= g(Y) |Enter Function f(X): ^ Algebraic enter mode is activated. 1ALPHA in case of alpha mode. EDIT menu is activated (-> [->STK] command). |'COSH(X) Just pressing would mean that g(X)=X! |Modle: LINFIT |B + f(X) *M= g(Y) |Enter Function g(Y): ^ |' Pressing without entering any function element means that g(Y)=Y. |Computing LR-Elements |Enter äLINE Name: ALPHA mode is activated. You can enter a function name for the curve. |APRX |Defining äLINE FCN 1: Correlation: .999999996453 Test: 0.5 [APRX] -> 1.12739143401 0.5 [IPOL] -> 1.12724112612 0.5 [COSH] -> 1.12762596521 [TRANS] TRANSFORMING COORDINATES ================================= INPUT: none OUTPUT: none ERROR: "Bad äDAT Data" äDAT either doesn't exist or is not a matrix "Wrong äDAT Size" äDAt size is different from { n 2 } [TRANS] transforms data points from the 'äDAT' matrix: [[ .. .. ] [[ .. .. ] [ Xn Yn ] -> [ f(Xn,Yn,an,rn) g(Xn,Yn,an,rn) ] [ .. .. ]] [ .. .. ]] 'äDAT' PURGE [DIFF] [Y=SF(xy)dx] -> |'TAN(l) (lambda: [right-shift] [NXT]) | 0 -1 X0=0 and Y0=-1 | .1 dH=.1 | .1 dS=.1=dH (every point is stored) | algorithm able to compute l (lambda) | 63 63 points [TRANS] -> |Enter X-Transformation |X1=f(X,Y,a,r) !!!! a = ATAN(Y/X) [right-shift] [A] !!!! r = SQRT(SQ(X)+SQ(Y)) [right-shift] [R] !!!! [X Y] = cartesian coordinates !!!! [a r] = polar coordinates Algebraic enter mode is activated. 1ALPHA in case of alpha mode. EDIT menu is activated (-> [->STK] command). |'X*SIN(2*a) a, Alpha Just typing would mean that X1=X! |Enter Y-Transformation |Y1=g(X,Y,a,r) |'Y*SIN(2*a) a, Alpha Just typing would mean that Y1=Y! |Transforming (X,Y) Data points are transformed and drawn as a scale curve. [REVä] REVERSING äDAT MATRIX =============================== INPUT: none OUTPUT: none ERROR: "Bad äDAT Data" äDAT either doesn't exist or is not a matrix "Wrong äDAT Size" äDAt size is different from { n 2 } ACTION: reverses 'äDAT' matrix EXAMPLE: [[ 1 1 ] [REVä] [[ 4 8 ] [ 2 4 ] [ 3 6 ] [ 3 6 ] [ 2 4 ] [ 4 8 ]] [ 1 1 ]] [ORDä] ORDERING äDAT MATRIX ============================== INPUT: none OUTPUT: none ERROR: "Bad äDAT Data" äDAT either doesn't exist or is not a matrix "Wrong äDAT Size" äDAt size is different from { n 2 } ACTION: orders 'äDAT' matrix (X-column in GROWING order; Xn > Xn-1) ^^^^^^^^^^^^^^^^^^^^^^^^^ EXAMPLE: [[ 1 -2 ] [ORDä] [[ -9 2 ] [ 6 7 ] [ 1 -2 ] [ 3 4 ] [ 3 4 ] [ -9 2 ]] [ 6 7 ]] [DRWä] DRAWING äDAT MATRIX (autoscaled) ========================================== INPUT: none OUTPUT: none ERROR: "Bad äDAT Data" äDAT either doesn't exist or is not a matrix "Wrong äDAT Size" äDAt size is different from { n 2 } ACTION: erases PICT and draws a scatter plot of 'äDAT' (AUTOSCALED) autoscaled => äMIN=PMIN & äMAX=PMAX ^^^^^^^^^^ ^^^^^^^^^^ [SCLä] DRAWING äDAT MATRIX (scale plot) ========================================== INPUT: none OUTPUT: none ERROR: "Bad äDAT Data" äDAT either doesn't exist or is not a matrix "Wrong äDAT Size" äDAt size is different from { n 2 } ACTION: erases PICT and draws a SCALE 'äDAT' plot ^^^^^ [LINä] CONNECTING SCATTERED POINTS BY LINES ============================================= INPUT: none OUTPUT: none ERROR: "Bad äDAT Data" äDAT either doesn't exist or is not a matrix "Wrong äDAT Size" äDAt size is different from { n 2 } ACTION: uses current plot parameters to connect 'äDAT' data points by lines EXAMPLE: [DRWä] [LINä] or [SCLä] [LINä] BLOC DEFINITIONS ================ NBR TYPE ACTION E P 0* NULL ENTRY E0=0 - - 1 ABSOLUTE VALUE E0=ABS(E1) 1 0 2 SIGN E0=SIGN(E1) 1 0 3 NEGATIVE VALUE E0=-E1 1 0 4 SIGN TRANFER E0=E1*SIGN(E2) 2 0 5 CLIP NEGATIVE E0=IFTE(E1>0,E1,0) 1 0 6 CLIP POSITIVE E0=IFTE(E1>0,0,E1) 1 0 7 FREE SPACE E0= * [see below] 1 2 8 LIMIT E0= * [see below] 1 2 9 SCALE E0= * [see below] 1 3 10 PHASE E0= * [see below] 1 3 11 RELAY 1 E0=IFTE(E1>=0,E2,E3) 3 0 12 RELAY 2 E0=IFTE(E1>0,E2,0) 2 0 13 MINIMUM E0=MIN(E1,E2) 2 0 14 MAXIMUM E0=MAX(E1,E2) 2 0 15 ADDITION 2 E0=E1+E2 2 0 16 ADDITION 3 E0=E1+E2+E3 3 0 17 MULTIPLICATION E0=E1*E2 2 0 18 DIVISION E0=E1/E2 2 0 19 ADDITION CONSTANT E0=E1+P1 1 1 20 MULTIPLICATION CONSTANT E0=E1*P1 1 1 21 ADD+MUL 2 E0=E1*P1+E2*P2 2 2 22 ADD+MUL 3 E0=E1*P1+E2*P2+E3*P3 3 3 23* 24* 25 INVERSE E0=INV(E1) 1 0 26 SQUARE E0=SQ(E1) 1 0 27 SQUARE ROOT E0=SQRT(E1) 1 0 28 POWER E0=E1^E2 2 0 29 POLYNOM E0=P1+P2*E1+P3*E1^2 1 3 30 COMMON LOGARITHM E0=P1*LOG(E1+P2)+P3 1 3 31 COMMON EXPONENTIAL E0=P1*ALOG(E1*P2)+P3 1 3 32 NATURAL LOGARITHM E0=P1*LN(E1+P2)+P3 1 3 33 NATURAL EXPONENTIAL E0=P1*EXP(E1*P2)+P3 1 3 34* 35* 36 SINE E0=P1*SIN(E1*P2+P3) 1 3 37 COSINE E0=P1*COS(E1*P2+P3) 1 3 38 TANGENT E0=P1*TAN(E1*P2+P3) 1 3 39 ARC SINE E0=P1*ASIN(E1*P2+P3) 1 3 40 ARC COSINE E0=P1*ACOS(E1*P2+P3) 1 3 41 ARC TANGENT E0=P1*ATAN(E1*P2+P3) 1 3 42 SINE HYPERBOLIC E0=P1*SINH(E1*P2+P3) 1 3 43 COSINE HYPERBOLIC E0=P1*COSH(E1*P2+P3) 1 3 44 TANGENT HYPERBOLIC E0=P1*TANH(E1*P2+P3) 1 3 45 ARC SINE HYPERBOLIC E0=P1*ASINH(E1*P2+P3) 1 3 46 ARC COSINE HYPERBOLIC E0=P1*ACOSH(E1*P2+P3) 1 3 47 ARC TANGENT HYPERBOLIC E0=P1*ATANH(E1*P2+P3) 1 3 48* 49* 50 INTEGRATION 1 E0=P1+Integral[E1*P2]dt 1 2 51 INTEGRATION 2 E0=P1+Integral[E1*P2+E2*P3]dt 2 3 52 LIMIT INTEGRATION E0=P1+Integral[E1]dt,P2<=E0<=P3 1 3 53 DERIVATIVE ** [see below] E0=dE1/dt 1 0 54* 55* 56 RANDOM E0=RAND 0 0 57 CONSTANT E0=P1 0 1 58 QUIT quit if E1>E2 2 0 59* 60 NULL E0=0 0 0 61* 62* 63* 64* 65* 66* 67* 68* 69* 70* 71* 72* 73* 74* 75* 76* 77* 78* 79* 80* 81* 82* 83* 84* 85* 86* 87* 88* 89* 90* 91* 92* 93* 94* 95* 96* 97* 98* 99* 100* CURRENT TIME E0=current time 0 0 nn* undefined block type 0* can only be used to determine a null entry (En=0) in entry list 100* can only be used as time input ** Block 53: derivative through the last 3 points Note: - Input E---[D>---O Output O=dE/dt - Input E---[I>---[D>---O Output O is NOT equal to E for algorithm reasons: TIME t INPUT E INTEGRAL I OUTPUT O tn En In On-1 (using In-3, In-2 and In-1) tn+1 En+1 In+1 On (using In-2, In-1 and In) On(tn+1)=En(tn)! *** > 100 USER DEFINED BLOCKS Block 7 ------- ^ ^ / ^ | | / \ | | / | 1:1 / \ 1:1 | | 1:1 / | / \ | P1 | / | P1 / \ |P1 ----:=====|=====:---> --|--:=====:--------> ----:=====|=:-------> / | P2 | / P2 P2 | \ / | |/ | \ / 1:1 | / 1:1 |1:1 \ / | /| | \ | | | \ P1 < P2 P1 < P2 P1 > P2 Block 8 ------- P2- ------ ^ ^ | / P2- ------------ ------- - -P2 | / 1:1 | / \ | | / | / 1:1 1:1 \ | P1 |/ | / \| P1 ------:---|----:----> ---P1- ------:---|----:----> /| P2 | P2 |\ 1:1 / | | P1 | \ / | --|--:---:----------> | \ 1:1 ------- -P1 | P2 | \ | | -P1 - ------ P1 < P2 P1 < P2 P1 > P2 Block 9 ------- P2>0 ^ <-------> ^ | ------ P1 | | - P3 / ====:-----|-:-------> --------- - P3 | / \ | P1+P2 \| | / P3:P2 \ | \ | / \ | |\ P3:P2 | / P3:P2 \ | | \ | / \| P1+P2 | \ ==|=====:------:----> \ --------:-|---:=====> | P1 P1+P2 |\ | P1 | <------> P3- --------- <-----> P2>0 P2<0 Block 10 -------- ^ | ---|P3--------- | | | | | | | | | ===:---|-----------:===> P1 | P2 | P1 must be < than P2!