( sim_equns)
( 15th July 2015)

CR ." sim_equns: "

( Solving simultaneous equations using Pivotal Elimination)


( 1 - move equation with greatest 1st coefficient to the 'top' - zeroth posn)
( 2 - divide the zeroth posn equation by the first coefficient - to get '1')
( 3 - a{n,m}=a{n,m} - a{0,m} * a{n,0} for n>0 [row no.] )
( 4 - divide the n=1 row by first coefficient to get '1')
( 5 - a{n,m}=a{n,m} - a{0,m} x a{n,0} for n>1 [row no.] )
( 6 - repeat until get '1' in each diagonal position)
( 7 - solve for xn working from last equation using substitution)


( divide row n by nth element of row - first non-zero value)
: norm_row     ( n...)
  DUP DUP matrix @              ( n,divisor...)
  #columns @ 0 DO 2DUP SWAP I matrix @   ( n,divisor,divisor,val...)
         SWAP F/                ( n,divisor,norm value...)
         ROT                    ( divisor,norm value,n...)
         DUP >R I matrix ! R>   ( divisor,n...)
         SWAP                   ( n,divisor...)
      LOOP 2DROP ;
 
: pivot
  #rows @ 1- 0 DO  I DUP n ! norm_row 
        #rows @ I 1+ DO 
                 I J matrix @ a{n,m} !   ( ...)         ( multiplier for normalised row)
                 #columns @ 0 DO  n @  I          ( n,m...)      ( start m from value of n)
                      matrix @           ( v1...)        ( element in normalised row)
                      a{n,m} @           ( v1,v2...)     ( multiplier)
                      F*                 ( product...)   ( now subtract this from element on row n)
                      FMNF               ( negate - to use F+!)
                      J I matrix         ( -product,address...)
                      F+!                ( ...)
                     LOOP
                 LOOP
        LOOP
        #rows @ 1- norm_row ;

: test pivot dsp_matrix ;

: dsp_xn #rows @ 0 DO  CR ." x(" I . ." ) = " I x(i) @ F. LOOP ;

: zero_xn #rows @ 0 DO 0 I x(i) ! LOOP ;

: calc_xn
   zero_xn  
   #zrows  #zcolumns matrix @  #zrows x(i) !      ( bottom equation - involving no extra terms)
   0  #zrows 1- DO  I #zcolumns matrix @          ( rhs column value for row I)
           I 1+   #zcolumns 1-  DO J I matrix @   I x(i) @   F*   F-
                            -1 +LOOP 
                             I x(i) !  
         -1 +LOOP ;

: check ( ...)       ( put values into original equations)
   CR ." Check by substitution into original equations:" CR
   3 FLD !
   filename read_matrix    ( re-read original matrix)
   #rows @ 0 DO
     0  #zcolumns 0 DO  J I matrix @ DUP F. ."  * "
                            I x(i) @  DUP F. 
                            F*   F+  
                            I #columns @ 2- < IF ."  + " ENDIF
                    LOOP
                    ."  = " F.
             I #zcolumns matrix @  ."  ( cf: " F. ." )" CR
             LOOP ;

: solve_equns ( ...)   pivot  CR ." Modified Matrix:" CR dsp_matrix
                      calc_xn CR  dsp_xn  CR check ;

( Select stored file from here ------------)
: sel_file ( ...addr)
  filename prefile$  STRCPY        ( 'file' in filename)
  100 dirnames 1+                  ( max no. objects)
  CR show_names
  CR ." Input file number: "  QUERY INTERPRET
  DUP  10 < IF filename zero$ STRCAT ENDIF
  obj$ STR$ DROP                   ( numeric string suffix - nos 01 to max)
  filename obj$ STRCAT
  CR  filename .STR                ( ...)
  filename ;       ( ...addr)

: read_file   sel_file  read_matrix 
  CR ." stored matrix:"
  dsp_matrix 
  CR ." Solve? [Y;N]: "  Y/N
  IF solve_equns ENDIF ;  
