Τετάρτη 15 Αυγούστου 2018

Επίλυση Γραμμικού συστήματος

Ο κώδικας είναι γραμμένος στα αγγλικά και αναρτήθηκε πρώτη φορά στο RosettaCode.

Για τη λύση έχει επιλεχθεί η απαλοιφή Gauss όπως ήταν και το απαιτούμενο για το έργο στο RosettaCode.

Μας ενδιαφέρει η ακρίβεια. Για το σκοπό αυτό θα χρησιμοποιήσουμε αριθμούς τύπου Decimal. Πληροφορίες για τον τύπο αριθμών Decimal υπάρχουν εδώ. Ο αριθμός 10@ είναι τύπου decimal λόγω του χαρακτήρα @ που ακολουθεί τα ψηφία. Έτσι το α=1212.121@ θα κάνει την α τύπου Decimal.

Από την γραμμή εντολών μπορούμε να ορίσουμε γενικές μεταβλητές τύπου Decimal.
>α=1212.121@
>? Τύπος$(α)
Decimal

Ακολουθούν δυο παραδείγματα και στα δύο αυτά υπάρχει ο πίνακας a() με δυο σετ τιμών. Διαλέγουμε πριν εκτελέσουμε κάποιο από τα δύο (σκιάζουμε την γραμμή με τα στοιχεία με έναν χαρακτήρα"\" - συνήθως βάζουμε δύο- \\.

Σε καθένα θα γίνει η εύρεση των λύσεων (για απλοποίηση δεχόμαστε ότι τα δεδομένα που βάζουμε έχουν λύση).

Στο τέλος θα γίνει επαλήθευση ως προς ένα βαθμό στρογγυλοποίησης. Η Round() κάνει στρογγυλοποίηση στο μισό, δηλαδή αν ζητάμε στο 3 δεκαδικό τότε αν το τέταρτο είναι από 5 και πάνω, κερδίζουμε ένα χιλιοστό.

Οι Decimal αριθμοί έχουν μέγιστο 28 ψηφία, και μπορούν να μετακινούν το σημείο δεκαδικών, έτσι αν έχουμε έναν αριθμό με ακέραια ψηφία 24 δεν μπορούμε να έχουμε 5 δεκαδικά, αλλά το πολύ τέσσερα. Ομοίως αν έχουμε έναν αριθμό με ακέραια ψηφία 4 δεν μπορούμε να έχουμε 25 δεκαδικά, αλλά το πολύ 24.

Αν παρελπίδα δώσουμε σε μεταβλητή που δέχεται decimal αριθμό που δεν χωράει τότε βγαίνει λάθος. Αν ο Decimal είναι σε πίνακα, τότε δεν βγαίνει λάθος αλλά αλλάζει μορφή γίνεται Double.

Μας ενδιαφέρει το πρόγραμμα να εξαγει αποτελέσματα στην οθόνη και στο πρόχειρο (για να τα επικολλήσουμε κάπου, όπως στο Rosettacode). Για το σκοπό αυτό υπάρχει ένα αντικείμενο Document, το οποίο είναι αλφαριθμητικό με παραγράφους. Σε αυτό τροφοδοτούμε στοιχεία με το "=" τα οποία προστίθενται στο τέλος ( λειτουργία append). Επειδή θέλουμε να βάλουμε και τους πίνακες μέσα έχουμε μια βοηθητική συνάρτηση η οποία έχει εκτός από τον πίνακα που μας ενδιαφέρει (θέλουμε να υποστηρίζει πίνακες μιας και δυο διαστάσεων) και άλλες παραμέτρους. Δέχεται προαιρετικά  το μέγεθος περιθωρίου αριστερά (left margin),  επίσης προαιρετικά το μέγεθος πεδίου αριθμού, και τέλος επίσης προαιρετικά ένα αλφαριθμητικό που μπορεί να είναι κενό ή να έχει έναν αριθμό. Αν είναι κενό τότε αφήνει την υποδιαστολή ως έχει, αν δεν είναι κενό τότε αν είναι μηδέν δεν δείχνει την υποδιαστολή και στρογγυλεύει το νούμερο στο ακέραιο (με το μισό και πάνω υπέρ της μονάδας). Αν δώσουμε το "2" τότε θα στρογγυλέψει στα δυο δεκαδικά και ακόμα αν δεν έχουμε δεκαδικά θα βάλει δυο μηδενικά!


Στο πρώτο πρόγραμμα χρησιμοποιούμε αποθήκευση σκέτου αριθμού, ενώ στο δεύτερο αποθηκεύουμε κλάσματα.

Μια εντολή form 80,50 μπορεί να μεγαλώσει την ανάλυση της κονσόλας, για να βλέπετε όλα τα αποτελέσματα. (ή Φόρμα 80,80  με ή χωρίς τόνο στην εντολή).



module checkit {
      Dim Base 1, a(6, 6), b(6)
      a(1,1)= 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.63, 0.39, 0.25, 0.16, 0.10, 1.00, 1.26, 1.58, 1.98, 2.49, 3.13, 1.00, 1.88, 3.55, 6.70, 12.62, 23.80, 1.00, 2.51, 6.32, 15.88, 39.90, 100.28, 1.00, 3.14, 9.87, 31.01, 97.41, 306.02
      \\ remove \\ to feed next array
      \\ a(1,1)=1.1,0.12,0.13,0.12,0.14,-0.12,1.21,0.63,0.39,0.25,0.16,0.1,1.03,1.26,1.58,1.98,2.49,3.13, 1.06,1.88,3.55,6.7,12.62,23.8, 1.12,2.51,6.32,15.88,39.9,100.28,1.16,3.14,9.87,31.01,97.41,306.02
      for i=1 to 6 : for j=1 to 6 : a(i,j)=val(a(i,j)->Decimal) :Next j:Next i
      b(1)=-0.01, 0.61, 0.91, 0.99, 0.60, 0.02
      for i=1 to 6 : b(i)=val(b(i)->Decimal) :Next i
      function GaussJordan(a(), b()) {
            cols=dimension(a(),1)
            rows=dimension(a(),2)
            \\ make augmented matrix
            Dim Base 1, a(cols, rows)
            \\ feed array with rationals
            Dim Base 1, b(Len(b()))
            for diag=1 to rows {
                        max_row=diag
                        max_val=abs(a(diag, diag))
                        if diag<rows Then {
                              for ro=diag+1 to rows {
                                    d=abs(a(ro, diag))
                                    if d>max_val then max_row=ro : max_val=d
                              }
                        }
            \\         SwapRows diag, max_row
                        if diag<>max_row then {
                              for i=1 to cols {
                                    swap a(diag, i), a(max_row, i)
                              }
                              swap b(diag), b(max_row)
                        }
                        invd= a(diag, diag)
                        if diag<=cols then {
                              for col=diag to cols {
                                    a(diag, col)/=invd
                              }
                        }
                        b(diag)/=invd
                        for ro=1 to rows {
                              d1=a(ro,diag)
                              d2=d1*b(diag)
                              if ro<>diag Then {
                                         for col=diag to cols {a(ro, col)-=d1*a(diag, col)}
                                          b(ro)-=d2
                              }
                        }        
                  }
            =b()
      }
      Function ArrayLines$(a(), leftmargin=6, maxwidth=8,decimals$="") {
            \\ defualt  no set  decimals, can show any number
            ex$={
            }
           const way$=", {0:"+decimals$+":-"+str$(maxwidth,"")+"}"
            if dimension(a())=1 then {
                  m=each(a())
                  while m {ex$+=format$(way$,array(m))}
                  Insert 3, 2 ex$=string$(" ", leftmargin)
                  =ex$ : Break
            }
            for i=1 to dimension(a(),1) {
                  ex1$=""
                  for j=1 to dimension(a(),2 ) {
                              ex1$+=format$(way$,a(i,j))
                  }
                  Insert 1,2 ex1$=string$(" ", leftmargin)
                  ex$+=ex1$+{
                  }
            }
            =ex$
      }
      mm=GaussJordan(a(), b())
            c=each(mm)
            while c {
                  print array(c)
            }
      \\ check accuracy
      link mm to r()
      \\ prepare output document
      Document out$={Algorithm using decimals
                  }+"Matrix A:"+ArrayLines$(a(),,,"2")+{
                  }+"Vector B:"+ArrayLines$(b(),,,"2")+{
                  }+"Solution: "+{
                  }
      acc=25
      for i=1 to dimension(a(),1)
            sum=a(1,1)-a(1,1)
            For j=1 to dimension(a(),2)
                  sum+=r(j)*a(i,j)
            next j
            p$=format$("Coef. {0::-2},  rounding to {1} decimal, compare {2:-5}, solution: {3}", i, acc, round(sum-b(i),acc)=0@, r(i) )
            Print p$
            Out$=p$+{
            }
      next i
      Report out$
      clipboard out$
}
checkit




Εδώ στο δεύτερο πρόγραμμα αποθηκεύουμε σε πίνακα δυο στοιχείων κάθε στοιχείο του πίνακα Α, όπου το στοιχείο στη θέση 0 είναι ο αριθμητής και στη θέση 1 ο παρανομαστής. Με αυτό τον τρόπο έχουμε πιο αργό αλγόριθμο αλλά έχει μεγαλύτερη ακρίβεια! (ένα δεκαδικό)
 Ο αλγόριθμος για τα κλάσματα παίρνει βελτιώσεις, αλλά για την ώρα μένουμε εδώ (για να μην δυσκολέψει το πράγμα στη κατανόηση).



Module Checkit2 {
      Dim Base 1, a(6, 6), b(6)
      \\ a(1,1)= 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.63, 0.39, 0.25, 0.16, 0.10, 1.00, 1.26, 1.58, 1.98, 2.49, 3.13, 1.00, 1.88, 3.55, 6.70, 12.62, 23.80, 1.00, 2.51, 6.32, 15.88, 39.90, 100.28, 1.00, 3.14, 9.87, 31.01, 97.41, 306.02
      a(1,1)=1.1,0.12,0.13,0.12,0.14,-0.12,1.21,0.63,0.39,0.25,0.16,0.1,1.03,1.26,1.58,1.98,2.49,3.13, 1.06,1.88,3.55,6.7,12.62,23.8, 1.12,2.51,6.32,15.88,39.9,100.28,1.16,3.14,9.87,31.01,97.41,306.02
      for i=1 to 6 : for j=1 to 6 : a(i,j)=val(a(i,j)->Decimal) :Next j:Next i
      b(1)=-0.01, 0.61, 0.91, 0.99, 0.60, 0.02
      for i=1 to 6 : b(i)=val(b(i)->Decimal) :Next i
      \\ modules/function to use rational nymbers
      Module Global subd(m as array, n as array) { ' change m
            link m to m()
            link n to n()
            if m(0)=0 then return m, 0:=-n(0), 1:=n(1) : exit
            if n(0)=0 then exit
             return m, 0:=m(0)*(n(1)/m(1))-n(0), 1:=n(1)
      }
      Function Global Inv(m as array){
            link m to m()
            if m(0)=0@ then =m : exit
            =(m(1), m(0))
      }
      Function Global mul(m as array, n as array){' nothing change
             link m to m()
            link n to n()
            if n(0)=0 or n(1)=0 then =(0@,0@) : exit
           =((m(0)/n(1))*n(0),m(1))
      }
      Module Global mul(m as array, n as array) { ' change m
             link m to m()
            link n to n()
            if n(0)=0 or n(1)=0 then m=(0@,0@) : exit
             return m, 0:=(m(0)/n(1))*n(0)
      }
      Function Global Res(m as array) {
            link m to m()
            if m(0)=0@ then =0@: exit
            =m(0)/m(1)
      }
      \\  GaussJordan  get arrays byvalue
      function GaussJordan(a(), b()) {
            Function copypointer(m) { Dim a() : a()=m:=a()}
            \\ we can use : def copypointer(a())=a(0),a(1)
            cols=dimension(a(),1)
            rows=dimension(a(),2)
            Dim Base 1, a(cols, rows)
            for i=1 to cols : for j=1 to rows : a(i, j)=(a(i, j), 1@) : next j : next i
            def d as decimal
            for j=1 to rows : b(j)=(b(j), 1@) : next j
            for diag=1 to rows {
                        max_row=diag
                        max_val=abs(Res(a(diag, diag)))
                        if diag<rows Then {
                              for ro=diag+1 to rows {
                                    d=abs(Res(a(ro, diag)))
                                    if d>max_val then max_row=ro : max_val=d
                              }
                        }
            \\         SwapRows diag, max_row
                        if diag<>max_row then {
                              for i=1 to cols {
                                    swap a(diag, i), a(max_row, i)
                              }
                              swap b(diag), b(max_row)
                        }
                        invd= Inv(a(diag, diag))
                        if diag<=cols then {
                              for col=diag to cols {
                                    mul a(diag, col), invd
                              }
                        }
                        mul b(diag), invd
                         for ro=1 to rows {
                              \\ work also d1=(a(ro,diag)(0), a(ro,diag)(1))
                              d1=copypointer(a(ro, diag))
                              if ro<>diag Then {
                                         for col=diag to cols {subd a(ro, col), mul(d1, a(diag, col))}
                                          subd b(ro), mul(d1, b(diag))
                              }
                        }        

                  }
                  dim base 1, ans(len(b()))
                  for i=1 to cols {
                        ans(i)=res(b(i)) \\ : Print b(i)  ' print pairs
                  }
                  =ans()
      }
      Function ArrayLines$(a(), leftmargin=6, maxwidth=8,decimals$="") {
            \\ defualt  no set  decimals, can show any number
            ex$={
            }
           const way$=", {0:"+decimals$+":-"+str$(maxwidth,"")+"}"
            if dimension(a())=1 then {
                  m=each(a())
                  while m {ex$+=format$(way$,array(m))}
                  Insert 3, 2 ex$=string$(" ", leftmargin)
                  =ex$ : Break
            }
            for i=1 to dimension(a(),1) {
                  ex1$=""
                  for j=1 to dimension(a(),2 ) {
                              ex1$+=format$(way$,a(i,j))
                  }
                  Insert 1,2 ex1$=string$(" ", leftmargin)
                  ex$+=ex1$+{
                  }
            }
            =ex$
      }
      mm=GaussJordan(a(), b())
            c=each(mm)
            while c {
                  print array(c)
            }
      \\ check accuracy
      link mm to r()
      for i=1 to dimension(a(),1)
            sum=a(1,1)-a(1,1)
            For j=1 to dimension(a(),2)
                  sum+=r(j)*a(i,j)
            next j
            Print round(sum-b(i),26), b(i)
      next i
      \\ check accuracy
      Document out$={Algorithm using pair of decimals as rational numbers
                  }+"Matrix A:"+ArrayLines$(a(),,,"2")+{
                  }+"Vector B:"+ArrayLines$(b(),,,"2")+{
                  }+"Solution: "+{
                  }
      acc=26
      for i=1 to dimension(a(),1)
            sum=a(1,1)-a(1,1)
            For j=1 to dimension(a(),2)
                  sum+=r(j)*a(i,j)
            next j
            p$=format$("Coef. {0::-2},  rounding to {1} decimal, compare {2:-5}, solution: {3}", i, acc, round(sum-b(i),acc)=0@, r(i) )
            Print p$
            Out$=p$+{
            }
      next i
      Report out$
      clipboard out$
}
Checkit2




Παρακάτω ακολουθούν αποτελέσματα με αλλαγές στο acc, και στην επιλογή τιμών στον πίνακα a().





Algorithm using decimals
Matrix A:
          1,10,     0,12,     0,13,     0,12,     0,14,    -0,12
          1,21,     0,63,     0,39,     0,25,     0,16,     0,10
          1,03,     1,26,     1,58,     1,98,     2,49,     3,13
          1,06,     1,88,     3,55,     6,70,    12,62,    23,80
          1,12,     2,51,     6,32,    15,88,    39,90,   100,28
          1,16,     3,14,     9,87,    31,01,    97,41,   306,02

Vector B:
         -0,01,     0,61,     0,91,     0,99,     0,60,     0,02
Solution: 
Coef.  1,  rounding to 26 decimal, compare  True, solution: -0,0597391027501962649904316335
Coef.  2,  rounding to 26 decimal, compare  True, solution: 1,8501896672627829700670299288
Coef.  3,  rounding to 26 decimal, compare  True, solution: -1,9727833018116428175300387318
Coef.  4,  rounding to 26 decimal, compare  True, solution: 1,4697587750651240151384675034
Coef.  5,  rounding to 26 decimal, compare  True, solution: -0,5538741847821888403564152897
Coef.  6,  rounding to 26 decimal, compare  True, solution: 0,0723048745759411900531809852

Algorithm using pair of decimals as rational numbers
Matrix A:
          1,10,     0,12,     0,13,     0,12,     0,14,    -0,12
          1,21,     0,63,     0,39,     0,25,     0,16,     0,10
          1,03,     1,26,     1,58,     1,98,     2,49,     3,13
          1,06,     1,88,     3,55,     6,70,    12,62,    23,80
          1,12,     2,51,     6,32,    15,88,    39,90,   100,28
          1,16,     3,14,     9,87,    31,01,    97,41,   306,02

Vector B:
         -0,01,     0,61,     0,91,     0,99,     0,60,     0,02
Solution: 
Coef.  1,  rounding to 26 decimal, compare  True, solution: -0,0597391027501962649904316335
Coef.  2,  rounding to 26 decimal, compare  True, solution: 1,8501896672627829700670299288
Coef.  3,  rounding to 26 decimal, compare  True, solution: -1,9727833018116428175300387317
Coef.  4,  rounding to 26 decimal, compare  True, solution: 1,4697587750651240151384675034
Coef.  5,  rounding to 26 decimal, compare  True, solution: -0,5538741847821888403564152897
Coef.  6,  rounding to 26 decimal, compare  True, solution: 0,0723048745759411900531809852



Algorithm using decimals
Matrix A:
          1,00,     0,00,     0,00,     0,00,     0,00,     0,00
          1,00,     0,63,     0,39,     0,25,     0,16,     0,10
          1,00,     1,26,     1,58,     1,98,     2,49,     3,13
          1,00,     1,88,     3,55,     6,70,    12,62,    23,80
          1,00,     2,51,     6,32,    15,88,    39,90,   100,28
          1,00,     3,14,     9,87,    31,01,    97,41,   306,02

Vector B:
         -0,01,     0,61,     0,91,     0,99,     0,60,     0,02
Solution: 
Coef.  1,  rounding to 25 decimal, compare  True, solution: -0,01
Coef.  2,  rounding to 25 decimal, compare  True, solution: 1,6027903945021139442641548525
Coef.  3,  rounding to 25 decimal, compare  True, solution: -1,6132030599055614189052834829
Coef.  4,  rounding to 25 decimal, compare  True, solution: 1,2454941213714367443882298102
Coef.  5,  rounding to 25 decimal, compare  True, solution: -0,4909897195846576129526569211
Coef.  6,  rounding to 25 decimal, compare  True, solution: 0,0657606961752320046201065486


Algorithm using pair of decimals as rational numbers
Matrix A:
          1,00,     0,00,     0,00,     0,00,     0,00,     0,00
          1,00,     0,63,     0,39,     0,25,     0,16,     0,10
          1,00,     1,26,     1,58,     1,98,     2,49,     3,13
          1,00,     1,88,     3,55,     6,70,    12,62,    23,80
          1,00,     2,51,     6,32,    15,88,    39,90,   100,28
          1,00,     3,14,     9,87,    31,01,    97,41,   306,02

Vector B:
         -0,01,     0,61,     0,91,     0,99,     0,60,     0,02
Solution: 
Coef.  1,  rounding to 26 decimal, compare  True, solution: -0,01
Coef.  2,  rounding to 26 decimal, compare  True, solution: 1,6027903945021139442641548522
Coef.  3,  rounding to 26 decimal, compare  True, solution: -1,6132030599055614189052834817
Coef.  4,  rounding to 26 decimal, compare  True, solution: 1,2454941213714367443882298085
Coef.  5,  rounding to 26 decimal, compare  True, solution: -0,4909897195846576129526569203
Coef.  6,  rounding to 26 decimal, compare  True, solution: 0,0657606961752320046201065485

Δεν υπάρχουν σχόλια:

Δημοσίευση σχολίου

You can feel free to write any suggestion, or idea on the subject.