//
// Εργασία Αριθμητικής Ανάλυσης
// Χριστοπούλου Μαρία
// Α.Μ : 200600237
//
// Άσκηση 1α - Επίλυση Εξίσωσης με τη μέθοδο της διχοτόμησης
//
function bisection(a,b)
    // Επίλυση Εξίσωσης με τη μέθοδο διχοτόμησης
    //
    // Δεδομένα
    //    Η εξίσωση που δίνεται στην εργασία είναι:
    // f(x) = x^5-5*x^3+9*x-6
    //    με παράγωγο:
    // f'(x) = 5*x^4-15*x^2+9
    // 
    // Παράμετροι Συνάρτησης
    //    a : πίνακας 1-επί-1 doubles, το κάτω όριο του διαστήματος
    //    b : πίνακας 1-επί-1 doubles, το πάνω όριο του διαστήματος
    //    f : η δοσμένη συνάρτηση
    //    c : η ρίζα της εξίσωσης
    //
    // Περιγραφή συνάρτησης
    //    Η συνάρτηση υπολογίζει την ρίζα της εξίσωσης f(x)=0 με τη μέθοδο της διχοτόμησης με ακρίβεια 4 δεκαδικών ψηφίων.
    //    Από τη γραφική παράσταση της f(x) επιλέγουμε το κατάλληλο διάστημα, οπου περιέχεται μια ρίζα της εξίσωσης, και εισάγουμε το κάτω και άνω όριο στη συνάρτηση.
    //    Σε κάθε επανάληψη υπολογίζουμε το μέσο του διαστήματος [a,b] και ελέγχουμε αν αυτό είναι ρίζα της εξίσωσης.
    //    Αν δεν είναι, τότε ορίζεται νέο διάστημα ελέγχου, όπου το ένα άκρο είναι το μέσο του προηγούμενου διαστήματος. Οι τιμές των άκρων πρέπει να είναι ετερόσημες (ώστε να εξασφαλίζεται η ύπαρξη ρίζας).
    //    Το πρόγραμμα τυπώνει τον αριθμό επαναλήψεων, τις τιμές των άκρων του διαστήματος, το μέσο του διαστήματος, καθώς και την τιμή της συνάρτησης στο μέσο.
    //
    // Αναφορές
    // http://mathworld.wolfram.com/Bisection.html
    //
    
    // Ορισμός ακρίβειας
    tol = 1.e-4;
    // Αριθμός μέγιστων επαναλήψεων
    nmax = 100;
    // Counter
    aux = 1;
    //
    [lhs,rhs] = argn();
    // Έλεγχος εισόδου
    if (rhs ~= 2) then
        error("Wrong number of inputs. Input parameters should be 2.")
    end
    // Έλεγχος τύπου μεταβλητής
    if (length(a)~=1) then
        error("a is a scalar, not a vector.")
    end
    if (length(b)~=1) then
        error("b is a scalar, not a vector.")
    end
    // Έλεγχος a<b
    if ((myfun(a)*myfun(b)>0) | (a>b)) then
        error("It must be a<b, f(a)*f(b)<0 or f(a)*f(b)>0")
     end
     //
      mprintf("\n\nIter\t  A\t    C\t    B\t  F(a)\t   F(c)\t    F(b)\n");
     // Bisection Method
     while(abs(a-b)>= tol)
         // Υπολογισμός μέσου του διαστήματος [a,b] 
         c = (a+b)/2;
         //
         mprintf("\n%d\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t %.4f",aux,a,c,b,myfun(a),myfun(c),myfun(b));
         // Αριθμός επαναλήψεων
         aux = aux + 1;
         // Εύρεση διαστήματος, όπου εμπεριέχεται η ρίζα
         if (myfun(a)*myfun(c)> 0) then
             a = c;
         elseif (myfun(b)*myfun(c)>0)
             b = c;
         end
     end
     mprintf("\n\nRoot is %.4f\n",c);
     
     
endfunction
