Announcement

Collapse
No announcement yet.

Crout's method for solving equations

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Crout's method for solving equations

    https://forum.powerbasic.com/forum/u...ions#post24270

    Kerry,

    more of a blast than you might think.
    The original version was written for a BBC micro in the 1980s.

    Yes, it's the latest version.

    I don't understand the rest of your comments..
    "definitions of the fields"?

    The idea is that you have a number of interconnected but unknown quantities and equations that link them.
    You need to calculate the unknown quantities.
    For example, you buy 3 apples, 2 oranges and 2 bananas and the cost is $20.
    Then you buy 3 apples, 3 oranges and 4 bananas, and the cost is $30
    Next you buy 2 apples, 1 orange and 1 banana, and the cost is $11.
    How much is one banana?

    The program calculates that and gives 1 apple costs $2, 1 orange costs $4 and one banana costs $3.



    "when you say 'get the array data - can I type it in?"
    That's what the example shows, the array data being typed in.

    "How come this statement is before the dim statement?"
    It isn't.
    The first "LINE INPUT "array dimension?",d$ " gets the required array dimension and the array is dimensioned to that size in the next 2 lines.

    The next "LINE INPUT d$ " gets the required amount of data to fill the array.

  • #2
    Thanks Paul

    Big help!

    Kerry
    [I]I made a coding error once - but fortunately I fixed it before anyone noticed[/I]
    Kerry Farmer

    Comment


    • #3
      Kerry,
      if you're typing in the numbers every week then you can come up with a better user interface.
      When that was written the calculation was important and the user wasn't relevant!

      Comment


      • #4
        I looked at that old code and my mind was boggled by all the d$, d%, d##, s##,s$, o$=0 etc, etc.
        So I made it a bit easier for me to follow and an easier way to load the array:

        '
        Code:
        '************************************************************************************************************'
        '* PBWin10/ PBCC6                                                                                           *'
        '* solve simultaneous linear equations using Crout's method.                                                *'
        '*                                                                                                          *'
        '* based on code for PBCC1 by Paul Dixon, 2002:                                                             *'
        '*                                                                                                          *'
        '*  Equations in comma delimited file (crout.txt)                                                           *'
        '*  Equation 1a+2b+3c+4d+5e = 6                                                                             *'
        '*  Represented by 1,2,3,4,5,6                                                                              *'
        '*  File must have at least as many equations lines as variables                                            *'
        '*  any additional equations discarded.                                                                     *'
        '*                                                                                                          *'
        '* original Paul Dixon comments:                                                                            *'
        '*    REM program in PBCC1.0 to solve arrays using Crout's method                                           *'
        '*    REM The original array is destroyed in the process and the result is left                             *'
        '*    REM IN column zero of that array. The ARRAY is populated FROM 1 but dimmed                            *'
        '*    REM FROM 0 so this column is always free to take the results. Row0 is also unused.                    *'
        '*    REM The initial array is of the form:                                                                 *'
        '*    REM                                                                                                   *'
        '*    REM                     col0 col1 col2 col3                                                           *'
        '*    REM            (row 0) unused                                                                         *'
        '*    REM  equation 1(row 1)        1    2    3   representing 1a(1) + 2a(2) = 3                            *'
        '*    REM  equation 2(row 2)        4    5    6   representing 4a(1) + 5a(2) = 6                            *'
        '*    REM                                                                                                   *'
        '*    REM The above figures result in a(1)=-1,a(2)=2 which can easily be substituted back into the initial  *'
        '*    REM equations to verify they are correct. i.e. 1x-1 + 2x2 = 3  and 4x-1 + 5x2 = 6                     *'
        '*    REM Although only 2 dimensions are shown, the method is very efficient for large arrays.              *'
        '*   REM e.g. on a  400MHz K6-III, solving 200 simultaneous equations, each with 200 terms takes            *'
        '*   REM about 0.6seconds                                                                                   *'
        '************************************************************************************************************'
        #COMPILE EXE
        #DIM ALL
        
        FUNCTION PBMAIN() AS LONG
            ? SolveIt()
        #IF %DEF(%PB_CC32)
              WAITKEY$
        #ENDIF
        END FUNCTION
        
        FUNCTION SolveIt() AS STRING
            LOCAL strTemp AS STRING
            LOCAL lngVariables,lngLine,x,y,z,lngCol,lngRow AS LONG
            LOCAL extArr(),extSum,extTemp AS EXT
        '== Get data ==
            OPEN "crout.txt" FOR INPUT AS #1
            WHILE NOT EOF(#1)
                INCR lngLine
                LINE INPUT #1, strTemp
                IF lngLine = 1 THEN
                   lngVariables = PARSECOUNT(strTemp,",") -1
                   IF lngVariables < 3 THEN
                       FUNCTION = "Bad file format"
                       EXIT FUNCTION
                   END IF
                   DIM extArr(0 TO lngVariables+1, 0 TO lngVariables)
                END IF
                IF lngLine > lngVariables THEN DECR lngLine: EXIT LOOP ' drop any extra equations
                IF PARSECOUNT(strTemp,",") <> lngVariables + 1 THEN
                     FUNCTION = "Incorrect number of values in line "  & STR$(lngLine)
                    EXIT FUNCTION
                END IF
        
                FOR x = 1 TO lngVariables + 1
                     extArr(x,lngLine) = VAL(PARSE$(strTemp,",",x))
                NEXT
            WEND
            CLOSE #1
            IF lngLine < lngVariables THEN
                FUNCTION = "Insufficient equations: "  & STR$(lngLine)  & " equations for " & STR$(lngVariables) & " variables"
                EXIT FUNCTION
            END IF
        '== Solve ==
            FOR y = 1 TO lngVariables
                FOR lngRow = y TO lngVariables
                    extSum=0
                    FOR z=0 TO y-1
                        extSum+=extArr(z,lngRow)*extArr(y,z)
                    NEXT
                    extArr(y,lngRow)=extArr(y,lngRow)-extSum
                NEXT
        
                extTemp=extArr(y,y)
                FOR lngCol=y TO lngVariables+1
                    extSum=0
                    FOR z=0 TO y-1
                        extSum+=extArr(lngCol,z)*extArr(z,y)
                    NEXT
                    extArr(lngCol,y)=(extArr(lngCol,y)-extSum)/extTemp
                NEXT
            NEXT
        
            FOR lngRow = lngVariables TO 1 STEP -1
                extSum=0
                FOR z = lngRow TO lngVariables
                    extSum+=extArr(z,lngRow)*extArr(0,z)
                NEXT
                extArr(0,lngRow)=extArr(lngVariables+1,lngRow)-extSum
            NEXT
        '== Display result
            strTemp = ""
            FOR lngRow = 1 TO lngVariables
                strTemp +=  "a(" & FORMAT$(lngRow) & ")=" & STR$(extArr(0,lngRow)) & $CRLF
            NEXT
            FUNCTION = strTemp
        END FUNCTION
        '
        Test DATA using https://handymath.com/cgi-bin/matrix5.cgi
        Eq.1: 1V + 2W + 3X + 4Y + 5Z = 15 V = -0.02131
        Eq.2: 2V + 3W + -4X + 6Y + 7Z = 20 W = 2.82295
        Eq.3: 3V + 4W + 2X + 1Y + 5Z = 18 X = 0.22787
        Eq.4: 2V + 4W + 3X + 7Y + 5Z = 23 Y = 0.79180
        Eq.5: 5V + 6W + 3X + 6Y + 7Z = 30 Z = 1.10492

        Contents of Crout.txt:
        1,2,3,4,5,15
        2,3,-4,6,7,20
        3,4,2,1,5,18
        2,4,3,7,5,23
        5,6,3,6,7,30

        Output:
        ---------------------------
        PowerBASIC
        ---------------------------
        a(1)=-2.13114754098361E-2
        a(2)= 2.82295081967213
        a(3)= .227868852459016
        a(4)= .791803278688524
        a(5)= 1.10491803278688
        ---------------------------
        OK
        ---------------------------
        Last edited by Stuart McLachlan; 24 Aug 2020, 08:09 PM.

        Comment


        • #5
          Stuart,
          the original was written for a BBC micro when variables were typed by a suffix.
          I wouldn't write it that way now.

          Comment


          • #6
            Originally posted by Paul Dixon View Post
            Stuart,
            the original was written for a BBC micro when variables were typed by a suffix.
            I wouldn't write it that way now.
            I'm sure you wouldn't. I realise it was very old code.

            When I tried to work out what it was doing I kept getting lost so I did some find and replace as I worked through it. I then thought that posting the update would be useful for others (especially Kerry)

            Comment


            • #7
              Wow thanks Stuart - very helpful indeed.
              [I]I made a coding error once - but fortunately I fixed it before anyone noticed[/I]
              Kerry Farmer

              Comment


              • #8
                Interested in this topic. I studied Crout's method via enough YouTube videos to develop a sideways wobble to my head .

                I don't remember being taught Matrix Algebra so was lucky to find this site which helped me see how to use some of the MAT statements available in PB as an alternative to Crout's method to solve simultaneous equations as shown in test code below..

                Code:
                #Compile Exe        ' PBCC604 , PBWin1004
                #Dim All
                #If %Def(%PB_CC32)
                 #Console On
                 #Break On
                #EndIf
                #Include "win32api.inc"
                
                #If %Def(%PB_Win32)
                  Macro StdOut = MsgBox
                #EndIf
                '------------------/
                
                Function PBMain () As Long
                 Local Eqs, Terms, r, c As Long
                 Local sTemp, sErr As String
                 Local A(), Ainv(), X(), B() As Double
                
                  OPEN "mat1.csv" FOR INPUT AS #1                  ' Open data file
                    FileScan #1, Records To Eqs : If Eqs < 2 Then sErr = "File Data incompatable" : Goto Fin
                    Line Input #1, sTemp                          ' Read a line
                      Terms = ParseCount(sTemp,",") -1            ' (need at least three terms, same number of equations as variables)
                      If Terms <> Eqs Then sErr = "File Data incompatable" : Goto Fin
                      Redim A(1 To Eqs,1 To Eqs), Ainv(1 To Eqs,1 To Eqs), X(1 To Eqs), B(1 To Eqs)
                        Seek #1, 1
                        For r = 1 To Eqs                          ' load terms into arrays (matrices)
                          Line Input #1, sTemp
                          For c = 1 TO Eqs
                             A(r,c) = VAL(PARSE$(sTemp,",",c))    '
                          NEXT
                          B(r) = VAL(PARSE$(sTemp,",",Eqs+1))     ' Last item per line '=' B() constants
                        Next
                    CLOSE #1
                    sTemp = ""
                  ' Solve equations using PB MAT Statements
                  ' X = Ainv*B
                  MAT Ainv() = Inv(A())
                  MAT X() = Ainv() * B()
                  For r = 1 To Eqs
                   sTemp += $CRLF + Using$("X_(#_) _= ##.#### ", r, X(r))
                  Next
                
                  StdOut "Solved Variables" + $CRLF + sTemp
                Fin:
                   If Len(sErr) Then StdOut sErr
                
                #If %Def(%PB_CC32)
                  If Con.Handle Then
                    ? "Any key to exit"
                    WaitKey$
                  End If
                 #EndIf
                End Function
                '------------------/PBMain
                #IF 0
                Sample Data sets for "mat1.csv"
                (Five Vars)
                1,2,3,4,5,15
                2,3,-4,6,7,20
                3,4,2,1,5,18
                2,4,3,7,5,23
                5,6,3,6,7,30
                (Output
                B(1) = -0.0213
                B(2) = 2.8230
                B(3) = 0.2279
                B(4) = 0.7918
                B(5) = 1.1049)

                (Three Vars)
                3,2,2,20
                3,3,4,30
                2,1,1,11
                (Output
                B(1) = 2.0000
                B(2) = 4.0000
                B(3) = 3.0000)

                (Two Vars)
                1,2,3
                2,3,-4
                (Output
                X(1) = -17.0000
                X(2) = 10.0000)

                #ENDIF
                Last edited by Dave Biggs; 28 Aug 2020, 06:53 PM.
                Rgds, Dave

                Comment


                • #9
                  My Crout's code was written long before PB Matrix operations existed so so I didn't have them to use.
                  I think it was written in 1983, that might be before PB existed.

                  Comment


                  • #10
                    Wow! Your original code was 19 years old when posted as PBCC1 code in 2002. And still going strong!
                    Rgds, Dave

                    Comment


                    • #11
                      Originally posted by Dave Biggs View Post
                      Interested in this topic. I studied Crout's method via enough YouTube videos to develop a sideways wobble to my head .

                      I don't remember being taught Matrix Algebra so was lucky to find this site which helped me see how to use some of the MAT statements available in PB as an alternative to Crout's method to solve simultaneous equations as shown in test code below..
                      Thank you, I knew it could be simplified with MAT but I haven't touched a matrix for over 50 years so I didn't want to spend a lot of time re-learning about them.
                      That's a really good link!

                      Your method is going in my code library.

                      Comment


                      • #12
                        Originally posted by Stuart McLachlan View Post


                        Your method is going in my code library.
                        already in mine - thanks Kerry
                        [I]I made a coding error once - but fortunately I fixed it before anyone noticed[/I]
                        Kerry Farmer

                        Comment


                        • #13
                          Originally posted by Kerry Farmer View Post

                          already in mine - thanks Kerry
                          Watcha thanking yourself for?

                          Comment


                          • #14
                            Originally posted by Stuart McLachlan View Post

                            Watcha thanking yourself for?
                            oops - quoted wrong person!!! Never get old!!!
                            [I]I made a coding error once - but fortunately I fixed it before anyone noticed[/I]
                            Kerry Farmer

                            Comment


                            • #15
                              What's a coma between friends

                              Rgds Dave
                              Rgds, Dave

                              Comment


                              • #16
                                Originally posted by Dave Biggs View Post
                                What's a coma between friends

                                Rgds Dave
                                "Let's eat grandma!"



                                Comment


                                • #17
                                  Great thread! Thank you

                                  Comment


                                  • #18
                                    Originally posted by Stuart McLachlan View Post


                                    Your method is going in my code library.
                                    Stuart

                                    My code library is somewhat haphazard. I rely to a degree on the PB code backup system.

                                    Do you have a better system? And are you willing to share it?

                                    Thanks

                                    Kerry

                                    [I]I made a coding error once - but fortunately I fixed it before anyone noticed[/I]
                                    Kerry Farmer

                                    Comment


                                    • #19
                                      Originally posted by Kerry Farmer View Post

                                      Stuart

                                      My code library is somewhat haphazard. I rely to a degree on the PB code backup system.

                                      Do you have a better system? And are you willing to share it?

                                      Thanks

                                      Kerry
                                      Nope, I'm very disorganised, but I've got a good search tool

                                      I've just got a folder called PBMisc with lots of files with long descriptive names and the "Everything" tool from https://www.voidtools.com/

                                      This one went into the folder as "LinearEquationMatrix.bas"

                                      (compilable are *.bas, all other snippets are *.inc)

                                      I can find this one almost instantly by entering "*equation*.bas" (or even just "equation") in Everything

                                      Everything is a brilliant tool once you learn its capabilities. With bookmarkable search/filter capabilities its easy to locate anything

                                      For example, a bookmark in the top menu of "*\PBMisc\*.bas | *\PBMisc\*.inc !Backup" displays all the bas and inc files in PBMisc and its child folders, excluding backup files in about a tenth of a second and sorts them by name, path, date or size on request.

                                      Comment


                                      • #20
                                        Thanks Stuart - good thoughts
                                        [I]I made a coding error once - but fortunately I fixed it before anyone noticed[/I]
                                        Kerry Farmer

                                        Comment

                                        Working...
                                        X