Implementando el método de Falsa Posición en MAPLE

  El método de Falsa Posición (FP) es segundo en efectividad con respecto al método de bisección. La principal ventaja es que no parte en dos el intervalo de búsqueda de la raíz sino que esta división se basa en la siguiente fórmula:

$x_r=x_2-\frac{f(x_2)(x_1-x_2)}{f(x_1)-f(x_2)}$

en donde $x_r$ es el punto en donde una diagonal, que se crea a partir de los valores del intervalo de búsqueda, corta el eje horizontal. Créanlo o no, la idea del método de FP es mas efectiva que el de bisección. El esquema de abajo muestra gráficamente  lo que acabo de explicar:

Un comando en MAPLE

En MAPLE existe un comando muy sencillo. Solo debemos darle la función, el intervalo de búsqueda y decir qué tanta exactitud deseamos. Esto es como sigue:

> with(Student[NumericalAnalysis]):

> ftest:= x**2+8*x-3:

> FalsePosition(ftest, x = [-10, -8], tolerance = 1.5e-4);

-8.358663165

Debe notarse que podemos mejorar la exactitud haciendo más pequeño el número que se asigna a tolerance.

Una mejora al método de FP

Resulta que si usamos el método de FP para hallar una raíz como en la gráfica anterior, lo que veremos es que tardaremos en obtener la raíz y luego cierta exactitud. Esto se debe (como puede verse en la figura) a que:

• si la función tiene una parte asintótica (casi horizontal) y luego se hace vertical,
• si el intervalo de búsqueda se extiende más hacia la parte asintótica (casi horizontal) que hacia el otro extremo por donde está la raíz

el método de FP avanzará lentamente hacia la raíz, incluso más que el método de bisección. Entonces, surge el método de Falsa Posición Modificado (FPm). La idea detrás, es que si el programa del método de FPm detecta que si de una iteración a otra $f(x_1)$ cambia muy poco entonces se tiene el "estancamiento" que mencionaba antes, y la solución es que $f(x_2)$ sea dividido entre 2 para reducir el intervalo de búsqueda. Lo mismo se haría si el "estancamiento" ocurriera del lado de $x_2$.

Implementando un programa para método de FPm

Hay varias ideas acerca de como programar el método de FPm. Chapra y Canale dan un algoritmo en la séptima edición de su libro Métodos Numéricos para Ingenieros en donde al igual que en el método de bisección revisan en qué lado de la partición del intervalo está la raíz y cuentan del número de veces consecutivas que la raíz se ubica en la partición de $x_1$ y así también con $x_1$. Si, digamos, la raíz aparece dos veces (o mas) en la partición del lado de $x_2$ eso quiere decir que hay un posible "estancamiento" y por consecuencia $f(x_1)$ debe ser dividido entre 2.

Yo implementé en MAPLE el algoritmo propuesto por Chapra y Canale usando un procedure. El programa es como sigue:

> FlsPosmet:= proc(FUNC::algebraic,X1::float,X2::float,XACC::float)

     local JMAX, f1, f2, fr, x1, x2, xr, xrold, i1, i2, DX, J, ep:

  #Este es el número máximo de iteraciones

  JMAX:=40: 

  x1:= X1:

  x2:= X2:

  #Evaluar la función problema en los extremos proporcionados

  f2:= eval(FUNC,x=X2):

  f1:= eval(FUNC,x=X1):

  #Revisar si hay una raíz

  if f1*f2>=0 then 

    print("Bracketing Needed");

  end if:

  #Inicializar xr

  xr:= 0:
  #Inicializar contadores de repeticiones (modificación)
  i2:= 0: i1:= 0:

  #Cálculo de la raíz

  for J from 1 to JMAX do

  #Calcular el punto de cruce de la diagonal (falsa posición)

    xrold:= xr:

    xr:= x2 - f2*(x1 - x2)/(f1 - f2):

    fr:= eval(FUNC,x=xr):

  #Calcular el error porcentual

    ep:= abs((xr - xrold)/xr)*100:

  #Revisar de qué lado de xr está la raíz

    if f1*fr<0 then

      x2:= xr:

      f2:= fr:

      i2:= 0:

      i1:= i1 + 1:

      if i1>=2 then

        f1:= f1/2:

      end if

    else

      x1:= xr:

      f1:= fr:

      i1:= 0:

      i2:= i2 + 1:

      if i2>=2 then

        f2:= f2/2:

      end if

    end if:

    if ep<XACC or J=JMAX then break;

    end if:

  end do:

  return xr;

  end proc:

y para echar a andar el procedure usaré la misma función anterior, el mismo intervalo de búsqueda y la misma tolerancia. Es como sigue:

> FlsPosmet(ftest,-10.,-8.,1.5e-4);

-8.358898939

Los resultados son bastante parecidos. No puedo decir exactamente porque no son iguales. Es por lo que sigue:

• no sabemos si el comando de MAPLE corresponde al método de FP o al de FPm,
• no sabemos si el comando de MAPLE divido entre dos el intervalo, como sí se hace en nuestro procedure FlsPosmet,

entre otras cosas. Sin embargo, tenemos que contentarnos con que ambos van convergiendo a la misma raíz.

Otra idea de modificación

Chapra y Canale proponen que esto se haga con dos contadores: uno del lado de $x_1$ y el otro del lado de $x_2$ (los contadores son i1 y i2). Pero esa no es la única manera.

Se me ocurre que en vez de contar las veces consecutivas en que la raíz aparece en uno de los extremos del intervalo sería mejor cuantificar el cambio entre iteraciones de la función. Es decir, que si entre dos iteraciones sucesivas la raíz se encuentra del lado de $x_1$, entonces revisemos cuanto cambió $f(x_1)$ . Y que si ese cambio es menor al 1%, por ejemplo, multiplicar por 0.7 a $f(x_1)$ y si ese cambio está entre 1-2% por ejemplo multiplicar por 0.5 a $f(x_1)$. Esto podría ayudar a salir más rápido del "estancamiento" al programa.

Bueno. Solo es una idea que dejo al lector que trate de implementarla. Tal vez en un futuro la presente aquí.

Si tienes alguna duda ponla en los comentarios y trataré de ayudar.

Ildebrando.