martes, 9 de mayo de 2017

Regresión Polinomial

post original: https://metodosnumericosgera.wordpress.com/2017/04/06/regresion-polinomial/

Antecedentes

La aproximación polinomial es describir el comportamiento de una nube de datos a manera de una curva, en lugar de una recta. Se busca que esta recta difiera lo menos posible a la nube de datos.
Podríamos, por ejemplo, intentar aproximar la nube de datos a un polinomio de grado 2:
null
De tal forma que se minimice la suma de los cuadrados de los residuos, Sr:
null
Que, de igual manera, podría tratarse de un polinomio de m-ésimo grado:
null
El error estándar de estimación se calcula igual que en una regresión lineal:
null
El coeficiente de determinación y correlación también se calculan de la misma manera:
null
null

Justificación y propósito

Hemos hablado antes sobre la regresión lineal. Pero, ¿qué sucede cuando la nube de datos no está cerca a comportarse de manera lineal? Podríamos darnos cuenta si esto ocurre si al hacer la aproximación lineal tenemos un error estándar demasiado grande. Es aquí donde sería mejor usar regresión polinomial.

Diagrama de flujo

RegresionPolinomial.png

Codigo en C++


  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
#include <iostream>
#include <fstream>
#include <math.h>
#include <cmath>

using namespace std;
void regresionPolinomial(double x[], double y[], int n){

 const int grado = 2;            //Cambiar el grado para cambiar el grado del polinomio resultante <------- ESTE SI SE CAMBIA DEPENDIENDO DEL PROBLEMA

 double matrix[grado+1][grado+1];
 double solutions[grado+1];

 double sumX = 0;
 double sumXY = 0;

 cout << "------------SISTEMA DE ECUACIONES-------------" << endl;
 //Inicializamos el sistema de ecuaciones
 for (int i = 0; i < grado+1; ++i)
 {

  sumXY = 0;
  for (int k = 0; k < n; ++k)
  {
   sumXY += pow(x[k],i)*y[k];
  }
  solutions[i] = sumXY; //Get the solutiona of the ecuation system

  for (int j = 0; j < grado+1; ++j)
  {
   if (i == 0 && j == 0)
   {
    matrix[i][j] = n;   //Start the matrix[0,0] as n 
   }else{
    sumX = 0;
    for (int k = 0; k < n; ++k)
    {
     sumX += pow(x[k],(j+i));
    }

    matrix[i][j] = sumX;        //Get the other values of the matrix
   }

   cout << matrix[i][j] << ", ";
  }
  cout << "= " << solutions[i] << endl;
 }

 //Ahora resolvemos el sistema de ecuaciones con algun método conocido
 cout << "-------------------GAUSS SEIDEL---------------------" << endl;
 int size = grado +1;
 double EPSILON = 0.0000001; // error iterativo como tolerancia

 double S[size]; //Array with the solutions
 double OS[size];  //Array with old solutions
 double E[size]; //Array with th error of each variable
 for (int i = 0; i < size; ++i)
 {
  S[i] = 0;  //We initialize each solution at 0;
  OS[i] = 0; //Pseudo minus infinity at each old solution;

  E[i] = 0;
 }

 bool errorCheck = false;

 int counter=0;
 
 while(!errorCheck){
  for (int i = 0; i < size; ++i)
  {
   double sum = 0;
   for (int j = 0; j < size; ++j)
   {
    if (j != i)
    {
     sum += matrix[i][j]*S[j];
    }
   }
   OS[i] = S[i];  //The old value is the actual value
   S[i] = (solutions[i] - sum)/matrix[i][i];  //The actual value changes to the substitution of the variable

   E[i] = (S[i]-OS[i]);  //Error iterativo absoluto
   //E[i] = (S[i]-OS[i])/S[i];  // Error iterativo relativo
   //E[i] = ((S[i]-OS[i])/S[i])*100;  // Error iterativo porcentual
  }

  for (int i = 0; i < size; ++i)
  {
   if(abs(E[i]) > EPSILON){
    break;
   }
   errorCheck = true;
  }
  counter++;
 }

 for (int i = 0; i < size; ++i)
 {
  cout  << " a" << i << ": " << S[i] << endl;
 }
 cout << "it took " << counter << " times to terminate" << endl;
 //-------------------Aqui termina GaussSeidel------------------------

 cout << "------------------POLINOMIO---------------" << endl;
 for (int i = size-1; i >= 1; --i)
 {
  cout << S[i] << "X^"<< i << " + ";
 }
 cout << S[0] << endl;

 cout << "------------------ERRORES------------------" << endl;
 double e2[n]; //Calcula el error para cada X
 for (int i = 0; i < n; ++i)
 {
  double Ycalculada = 0;
  for (int j = size-1; j >= 1; --j)
  {
   Ycalculada += S[j]*(pow(x[i], j));
  }
  Ycalculada += S[0];
  e2[i] = pow(y[i]-Ycalculada,2);
 }

 double sr = 0;
 double sumY = solutions[0];
 double st = 0;
 for (int i = 0; i < n; ++i)
 {
  sr += e2[i];
  st += pow(y[i] - (sumY/n),2);
 }
 cout << "Sr: " << sr << endl;

 double sxy = sqrt(sr/(n-size));
 cout << "St: " << st << endl;

 double r2 = (st-sr)/st;
 double r = sqrt(r2);
 cout << "Error Estandar (Sx/y): " << sxy << endl;
 cout << "Coeficiente de determinacion (r2): " << r2 << endl;
 cout << "Coeficiente de correlacion (r): " << r << endl;
 cout << endl;
}

int main(){
 //CON ESTE SI DA TODO BIEN
 double x[] = {0,1,2,3,4,5};
 double y[] = {2.1,7.7,13.6,27.2,40.9,61.1};
 int dataNum = 6;
 
 regresionPolinomial(x,y,dataNum);
}

Ejemplos

Valores:
xy
02.1
17.7
213.6
327.2
440.9
561.1
Aproximación con polinomio de grado 2:
y = 1.86071x^2 + 2.35929x + 2.47857
Gráfica:
null
___________________________________________________________________
Valores:
xy
10.269.81
11.482.75
11.581.75
12.580.38
13.185.89
13.475.32
13.669.81
15.078.54
15.281.29
15.399.20
15.686.35
16.4110.23
16.5106.55
17.085.50
17.190.02
Aproximación con polinomio de grado 4:
y = -0.0025011X^4 + 0.0836578X^3 + -0.21331X^2 + -11.8438X^1 + 157.18
Gráfica:
null

Gauss Seidel y Jacobi

post original: https://metodosnumericosgera.wordpress.com/2017/03/15/gauss-seidel-y-jacobi/

Tanto el método Gauss-Seidel y Jacobi, son algoritmos que resuelven sistemas de ecuaciones mayores o iguales a segundo grado.
Al igual que los métodos anteriores de Gauss y Gauss-Jordan (ya que resuelven sistemas ecuaciones) necesitan un sistema de ecuaciones para su funcionamiento. Es decir, tener la misma cantidad de funciones que de incógnitas. No obstante, para que los métodos lleguen a una solución, dichas funciones no deben ser múltiplos entre sí, en otras palabras, no deben ser paralelas, de tal manera que las funciones se intersectan para que exista una solución para las incógnitas.
También es importante considerar que para ambos métodos, las ecuaciones del sistema de ecuaciones inicial necesitan estar acomodadas de manera que los coeficientes más grandes estén en la diagonal.
Esto significa que si comenzaramos con el siguiente sistema de ecuaciones:
image03
Lo primero que tendríamos que hacer es reacomodarlo de la siguiente manera:
image00
Los dos métodos son prácticamente iguales, en código fuente solo es cuestion de cambiar una cuantas líneas ya que utilizan el mismo principio y los mismos pasos para encontrar la solución. La unica diferencias es, como se muestra en la siguiente imagen, es la manera en que converge los valores de la incógnita. En Gauss-Seidel, dado los valores semillas para las incógnitas, el método evalúa x1 y con ella evalúa x2 los cuales son utilizadas para la siguiente incógnita, y así sucesivamente para las demás en la misma iteración. En cambio, Jacobi evalúa las incógnitas con y el valor nuevo lo guarda para la iteración siguiente.
null

Diagramas de flujo

Gauss-Seidel

Gauss-Seidell

Jacobi

JAcobi
Es un método iterativo porque utiliza la medición del error para detener las iteraciones en vez de realizar un proceso con un número de pasos definidos por el grado del sistema de ecuaciones.

Codigo Fuente en C++

Gauss-Seidel


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include 
#include 
using namespace std;

double EPSILON = 0.00001; // error iterativo como tolerancia

int x=0;
    int a=x; // copy of x
int y=0;
    int b=y; // copy of y
int z=0;
    int c=z; // copy of z

// augmented matrix
int rowA[]={2,-1,1,7}; // 2x -y +1 = 7
int rowB[]={1,2,-1,6}; // x +2y -z = 6
int rowC[]={1,1,1,12}; // x + y + z = 12

double anteriorX=-999999; // pseudo minus infinity
double anteriorY=-999999; // pseudo minus infinity
double anteriorZ=-999999; // pseudo minus infinity

int times=0;

void gaussSeidel(){
    while(true){
        times++;
        x=(rowA[3]-(rowA[1]*b)-(rowA[2]*c))/rowA[0];
        y=(rowB[3]-(rowB[0]*x)-(rowB[2]*c))/rowB[1]; // we use previously obtained value of x
        z=(rowC[3]-(rowC[0]*x)-(rowC[1]*y))/rowC[2]; // we use previously obtained values of x and y
        a=x;
        b=y;
        c=z;
        
        cout << "x: " << x << ", anteriorX: " << anteriorX << endl;
        cout << "y: " << y << ", anteriorY: " << anteriorY << endl;
        cout << "z: " << z << ", anteriorZ: " << anteriorZ << endl << endl;

        if(abs((x-anteriorX)/x)<=EPSILON && abs((y-anteriorY)/y)<=EPSILON && abs((z-anteriorZ)/z)<=EPSILON){
            break;
        }
        anteriorX=x;
        anteriorY=y;
        anteriorZ=z;
    }
}

int main(){
    gaussSeidel();
    cout << "FINAL ANSWER..." << endl;
    cout << "x: "<< x << endl;
    cout << "y: " << y <<endl;
    cout << "z: " << z << endl;
    cout << "it took " << times << " times to terminate";
    return 0;
}

Jacobi


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
#include <iostream>
#include <cmath>
using namespace std;

double EPSILON = 0.00001; // error iterativo como tolerancia

double x=0;
    double a=x; // copy of x
double y=0;
    double b=y; // copy of y
double z=0;
    double c=z; // copy of z

// augmented matrix
int rowA[]={2,-1,1,7}; // 2x -y +1 = 7
int rowB[]={1,2,-1,6}; // x +2y -z = 6
int rowC[]={1,1,1,12}; // x + y + z = 12

double anteriorX=-999999; // pseudo minus infinity
double anteriorY=-999999; // pseudo minus infinity
double anteriorZ=-999999; // pseudo minus infinity

int times=0;

void gaussSeidel(){
    while(true){
        times++;
        if(times>2){
            x=(rowA[3]-(rowA[1]*anteriorY)-(rowA[2]*anteriorZ))/rowA[0];
            y=(rowB[3]-(rowB[0]*anteriorX)-(rowB[2]*anteriorZ))/rowB[1];
            z=(rowC[3]-(rowC[0]*anteriorX)-(rowC[1]*anteriorY))/rowC[2];
        }else{
            x=(rowA[3]-(rowA[1]*b)-(rowA[2]*c))/rowA[0];
            y=(rowB[3]-(rowB[0]*a)-(rowB[2]*c))/rowB[1];
            z=(rowC[3]-(rowC[0]*a)-(rowC[1]*b))/rowC[2];
        }
        a=x;
        b=y;
        c=z;
        
        cout << "x: " << x << ", anteriorX: " << anteriorX << endl;
        cout << "y: " << y << ", anteriorY: " << anteriorY << endl;
        cout << "z: " << z << ", anteriorZ: " << anteriorZ << endl << endl;

        if(abs((x-anteriorX)/x)<=EPSILON && abs((y-anteriorY)/y)<=EPSILON && abs((z-anteriorZ)/z)<=EPSILON){
            break;
        }
        anteriorX=x;
        anteriorY=y;
        anteriorZ=z;
    }
}

int main(){
    gaussSeidel();
    cout << "FINAL ANSWER..." << endl;
    cout << "x: "<< x << endl;
    cout << "y: " << y <<endl;
    cout << "z: " << z << endl;
    cout << "it took " << times << " times to terminate";
    return 0;
}


Criterio de convergencia

La convergencia del método depende de algo llamado diagonal dominante, de acuerdo a esta pagina. Los valores de la diagonal deben ser mayores o iguales a la suma de los demás elementos de la fila. Es una condición que asegura la convergencia más no indica que obligatoriamente debe cumplir esa condición para que converja

Estos métodos pueden utilizar algo conocido como relajación, el cual, es una modificación al método de Gauss-Seidel donde para cada incógnita evaluada, se modifica con un promedio con peso obtenido con los resultados anteriores y con los resultados de la iteración presente. Esto, con el motivo de mejorar la convergencia del método. La fórmula para relajar es la siguiente, donde lambda, es un valor entre 0 y 2 previamente asignada. Para acelerar aún más la convergencia, se recomienda usar valores para lamba en un intervalo (1,2], esto se le llama sobrerrelajación.
image

Pruebas y Resultados

Consideremos los siguiente sistemas de ecuaciones de 3 incógnitas:
(A)
2x + z = 1

2x + 3y + z = 4

-3x + y - 2z = -3
(B)
-2x + 3y + 2z = 2

3x + z = 16

x - 4y = 11
(C)
3x + y - 2z = -10

2x - 3y + z = -4

x + 2y - 3z = -16
Las respuestas para cada uno de estos casos son:
CasoMétodoRespuestaIteraciones
AGauss-Seidelx = -2,
y = 1,
z = 5
38
AJacobix = -2,
y = 1,
z = 5
76
BGauss-Seidelx = 3,
y = -2,
z = 7
10
BJacobix = 3,
y = -2,
z = 7
30
CGauss-Seidelx = 1,
y = 5,
z = 9
17
CJacobix = 1,
y = 5,
z = 9
41
Como podemos observar, la velocidad de convergencia de Gauss-Seidel es aproximadamente 2 veces más rápida que la de Jacobi.

Casos de falla

Primeramente, si el sistema no tiene solución, el método nunca va a converger; además, si alguna variable es cero, se generaría un error de división entre cero. El programa no tiene ninguna manera de saber si fallo, al menos que se presente una división entre cero, ya que el método seguirá iterando indefinidamente dado que los valores se siguen alejando a la solución. Además, si no cumple con la propiedad de diagonal dominante (descrita anteriormente) el método difícilmente, si es que lo hace, va a converger.

Conclusiones

Estos métodos son muy deficientes en comparación a los otros métodos de los que se han visto para resolver sistemas de ecuaciones. Primeramente, este método tiene dos problemas: o simplemente no converge o lo hace muy lento. Esto se debe a que es un método iterativo, que, a comparación de los otros que no lo son y trabajan con matrices, realizan un número de procesos ya definidos por su cantidad de incógnitas donde es seguro que converge, si es que tiene solución, y es más sencillo monitorear su proceso y observar cual es el error en caso de que falle.

Método de Cramer

Post original: https://metodosnumericosgera.wordpress.com/2017/03/08/el-metodo-de-cramer/

El método de Cramer o regla de Cramer, es un método que sirve para resolver sistemas de ecuaciones lineales mediante el cálculo de determinante de matrices. Un determinante es un número real que corresponde a una matriz de n x n.
Para funcionar, primeramente, se necesita un sistema de ecuaciones, después, se toman los coeficientes que están antes de la igualación para formar una matriz, y los que están después de la igualación para forma una matriz de una columna (o vector), por ejemplo: dada la ecuación
formula1
Se crean las siguientes matrices:
formula2
Después, se calcula el determinante de la primera matriz; luego, la primera columna se reemplaza con la columna del vector b y se calcula su determinante detX. Se toma de nuevo la matriz original A y se reemplaza la segunda columna o columna de y por el vector para sacar su determinante; se repite el proceso para todas las columnas. Finalmente, para sacar los valores de las incógnitas, se toma el determinante detX y se divide por el determinante de la matriz A, esto calcula el valor de x; y para las demás se repite el proceso, pero con su respectiva determinante.
Este método falla cuando el determinante de la matriz original es cero, ya que, al evaluar las incógnitas mediante la división de determinantes se indeterminará. Esto puede significar que el sistema no tiene solución o muchas soluciones, y por lo tanto se debe utilizar otro método.

Diagrama de flujo

Cramer
Para un código en scilab, a lo más, el número de iteraciones es igual al tamaño de la matriz, como usa una matriz cuadrada se puede utilizar ya sea el tamaño de columnas o renglones para iterar. Incluso, las determinantes y las incógnitas se pueden calcular en una línea mediante la expresión A \ b, donde A es la matriz y b es el vector.

Código Fuente en SciLab


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
//Declara la matriz cuadrada de coeficientes como A
//Declara el vector columna de constantes como b

// Direct Solution A \ b

//Step to step Solution
Ap = A
[rX cX] = size(A)
//rX toma el size en renglones de la matriz
//cX toma el size en columnas de la matriz
D = det(A)

//Dado que rX y cX son iguales, por ser de una matriz cuadrada,
//se puede usar cualquiera de los dos para el ciclo
for i=1:1:rX
    Ap(:,i) = b         //Cambiar columna x1 al vector b
    DX = det(Ap)
    disp (DX/D)
    Ap(:,i) = A(:,i)    //Revertir cambio del primer paso
end

Pruebas con 2 incógnitas

formula3

Pruebas con 3 incógnitas

formula5

Conclusiones

El método de Cramer es muy sencillo de implementar, en especial en scilab, y funciona en la gran mayoría de los casos. Sin embargo, el proceso para evaluar los determinantes incrementa exponencialmente dependiendo del tamaño de la matriz, lo que lo hace ineficiente, en términos de procesamiento. Por lo tanto, solamente se recomienda para sistemas de ecuaciones de tres incognitas o menos.

Bibliografía:

Gilberto, E. (2001). Matrices and Linear Algebra with SCILAB.

Método de la Secante

El método de secante es un algoritmo que busca la raíz de una función. Podría considerarse como una variante del método Newton-Raphson por todas las similitudes que tienen, ambas son métodos abiertos, funcionan a partir de una fórmula que proviene de la pendiente de una gráfica, y tienen las mismas ventajas y desventajas. A diferencia de Newton-Raphson, el método secante necesita dos valores de entrada para funcionar y no se necesita evaluar la derivada de la función en la cual se quiere encontrar la raíz.
Al igual que el método Newton-Raphson, el método secante funciona entorno a una formula, en este caso, la fórmula de la secante; la cual, dado dos puntos como valores de entrada, el método calcula una secante para acercarse a la raíz. Gráficamente, funciona igual que Newton-Raphson, por excepción de que, en vez de calcular una tangente dado un punto, este calcula una secante dado dos puntos: xy x1. El punto en la gráfica, donde la secante de la iteración anterior cruza eje x, se convierte en x(o nuevo x0), y se calculara una nueva secante con xy x2; este proceso se repetirá hasta que uno de los puntos evaluado en la función sea cero.
Por ser método abierto, no tiene requisitos estrictos para que funcione, como es el caso del método bisección; sin embargo, dada una función donde su raíz sea un punto de inflexión, el método tendrá que hacer muchas más iteraciones, o incluso a nunca terminar, a comparación con los métodos bisección y Newton-Raphson. Además, los dos puntos dados no deben formar una línea paralela al eje x al calcular la secante, ya que esto provocaría que el método se itere indefinidamente.
secant

El método se detiene se detiene si y solamente si la función evaluada en uno de los puntos donde la secante cruza el eje x, es cercana a cero; o puede que nunca termine si la secante es paralela al eje x.

Codigo Fuente en C++


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#include <iostream>
#include <cmath>

using namespace std;
double EPSILON = 0.0000000000001;

double f (double x) {
    //return exp (-x) - x;
    //return log (x);
    //return pow(x,2);
    //return (sin(x) + 2 * x - 1);
    return pow(x-1, 3) + 0.512;
    //return sin(x);
}

void finish( double x, int count ){
    cout << "One root is: " << x << endl; 
    cout << "Iterations: " << count; 
    return;
}

int main() {
    double x0, x1, xn, rb;  //rb -> first root evaluated
    int count = 1;          //counter
    x0 = 4;               // -> Entry points 
    x1 = 1;                 // -> Entry points
    
    //Check entry points
    if ( f(x0) == 0 ) {
        finish( x0, 0 );
    }
    if ( f(x1) == 0 ) {
        finish( x1, 0 );
    }
    //end check entry points
    
    //Start loop
    do {
        //Evaluate entry points
        double fX0 = f( x0 );
        double fX1 = f( x1 );
        
        //Evaluate new x
        xn = x1 - ( fX1 ) * ( x0 - x1) / ( fX0 - fX1 );
        
        //Get Eaprox% and Print
        double Eap = abs( (xn - x1)  / xn ) * 100;
        cout << "Current iteration: " << count << "\t";
        cout << "X new Value: " << xn << "\t";
        cout << "Error Iterativo Porcentual: " << Eap << endl;
        //End Get Eaprox% and print
        
        //Check new X
        if ( abs( f(xn) ) <= EPSILON ) {
            break;
        }
        
        //Changing points
        x0 = x1;
        x1 = xn;
        
        count++;
    } while (true);
    
    finish( x1, count );
    
    return 0;
}

Pruebas y Resultados

PruebaFunciónInputOutputIteraciones
1X2[-2, -1]-3.507X10-732
2X2[-2, 2]INF
3sin(x) + 2x – 1[-2, 5]0.3354187
4sin(x) + 2x – 1[-5, 5]0.3354186
5e-x – x[5, 6]0.5671437
6Sin (x)[4, -1]3.141594
7Sin (x)[4, -0.2]-3.81498x10-75
8Sin (x)[4, -2] INF
9(x – 1)+ 0.512[4, 1]0.2183
10(x – 1)+ 0.512[10, 3]0.215
En la prueba 2, el método falla porque la secante evaluada es paralela al eje x, y por lo tanto, nunca cambia en las posiciones de x en cada iteración. Lo mismo sucede en la prueba 8, en la gráfica de seno, los puntos forman una línea casi horizontal.
La gran cantidad de iteraciones en la prueba 9, es debido a que en la gráfica, en el rango de aproximadamente de -1 a 2, la pendiente es casi 0, por lo que la secante evaluada es casi horizontal pero lo suficientemente inclinada para encontrar un valor de x en la segunda iteración, el cual era de -153, después vuelve a 0.9 y se aleja de nuevo de la raíz.

Conclusiones

El método de la secante, es una variación del método de Newton-Raphson; obtiene algunas desventajas a cambio de hacer la formula principal en algo más sencillo. Esto es debido a que no necesita evaluar la derivada de la función. Sin embargo, los resultados no solamente dependen de qué punto de la función empieza el algoritmo (como en el caso de Newton-Raphson), sino que ahora depende de dos puntos; los cuales, si están muy separados o la línea generada entre los dos puntos es muy cercana a una línea horizontal, provocaría que el algoritmo haga más iteraciones o tenga mayor error en comparación con Newton-Raphson.
A pesar de lo dicho anteriormente, no es posible concluir cuál de los dos métodos es mejor, el error y las iteraciones de cada uno varia con respecto a la función la cual se quiera encontrar su raíz. Desde mi punto de vista, Newton-Raphson es más fácil de usar, dado que solo ocupa un dato de entrada y si utilizas la formula siguiente para sacar la derivada.
f( x + ALMOST_CERO ) – f ( x ) ) / ( ALMOST_CERO ). Para entender mejor su uso, vea el codigo del blog de Newton-Raphson.