(#391). BÚSQUEDA DE SOLUCIONES (V): MÉTODO DE LA POSICIÓN FALSA

[MONOTEMA] Una variante del método de la secante es el método de la posición falsa (Regula Falsi), que básicamente ejecuta el procedimiento de la secante pero cambiando levemente el algoritmo para que entre dos puntos de dos iteraciones sucesivas siempre se encuentra la solución, de forma análoga a como ocurre con el método de la bisección. Lo haremos, como siempre, desde el punto de vista práctico, y siguiendo a Burden, Faires & Burden (2017).

Función de partida

Emplearemos la misma función que en los ejemplos anteriores, donde se relaciona la inversión en publicidad con los beneficios.

f_ceronegativo:ev(funcion,x=-0.03297097167558977), numer;
f_ceropositivo:ev(funcion,x=3.03297097167559), numer;
plot2d([[discrete, [[-0.0329709716755897,f_ceronegativo],
[3.03297097167559, f_ceropositivo]]], funcion],[x,-1,4],[y,-3,3],
[xlabel, “Inversión en publicidad (miles de euros)”],
[ylabel, “Beneficios (miles de euros)”],
[style, points, lines], [color, red, green], [legend, false]);

Método de la posición falsa

(1) Objetivo: Aproximarse numéricamente a la solución p de una función , tal que

(2) Condiciones: La función debe ser continua en el intervalo considerado.

(3) Descripción rápida: Partimos del método de la secante:

Para calcular p2 empleamos la fórmula anterior (tras proveer las semillas p0 y p1). Sin embargo, para calcular p3, usaremos los valores (p2,p1,f(p2),f(p1)) sólo si f(p2) y f(p1) tienen signo opuesto. Si ambos tiene el mismo signo, entonces en lugar del par (p1,f(p1)), escogemos (p0,f(po)). Y se sigue ese mismo criterio para el resto de iteraciones.

(4) Convergencia: Un criterio suficiente (pero no necesario) para que lo haga es que:

Si se cumple la fórmula anterior cuando está en un intervalo alrededor de la raíz entonces la solución converge en cualquier punto de ese intervalo. Para emplear este criterio, hemos de añadir el requerimiento de que la función sea dos veces derivable.

Normalmente, el método de la posición falsa converge de forma más lenta que el de Newton y el de la secante.

(5) Estimación del error: Hay diversas opciones, aunque una de las más recomendables es el error relativo:

Implementación en Maxima

Vamos a implementar el algoritmo partiendo de los dos valores con los que inicíabamos el método de la bisección: [2.5,3.5], que actúan como las primeras semillas en el método de la posición falsa, y vamos a adaptar en Maxima la propuesta de Ramírez (2008):

falsapos(arg):=block([a:arg[1],b:arg[2],f:arg[3],xm], xm:(a*ev(f,x=b)-b*ev(f,x=a))/(ev(f,x=b)-ev(f,x=a)),
if (ev(f,x=xm)*ev(f,x=a)<0) then [a,xm,f] else [xm,b,f] );
ini:[2.5,3.5,-(x^2)+(3*x)+0.1]; for i:1 thru 10 do (print(ini[1]),ini:falsapos(ini));

Y nos genera estos resultados:

2.5
2.95
3.021739130434783
3.031481481481481
3.03277399056109
3.03294493097818
3.032967529289182
3.032970516620634
3.032970911521146
3.032970963723678

Para trabajar un poco estos procedimientos, los estudiantes pueden intentar programar una rutina con Maxima donde se codifiquen los procedimientos de Newton, secante y posición falsa, y comparar las aproximaciones.

Conclusión

En este post hemos explicado brevemente en qué consiste el método de la posición falsa. Burden, Faires & Burden (2017) no recomiendan por lo general este método, pero no deja de ser interesante para trabajar las rutinas de programación y para entender los diferentes procedimientos numéricos de búsquedas de raíces.

Todos los posts relacionados




(#390). BÚSQUEDA DE SOLUCIONES (IV): MÉTODO DE LA SECANTE

[MONOTEMA] Una modificación del método de Newton permite una aproximación a este sin computar la derivada. Se trata del método de la secante, que explicaremos a continuación. Lo haremos, como siempre, desde el punto de vista práctico, y siguiendo a Burden, Faires & Burden (2017).

Función de partida

Emplearemos la misma función que en los ejemplos anteriores, donde se relaciona la inversión en publicidad con los beneficios.

f_ceronegativo:ev(funcion,x=-0.03297097167558977), numer;
f_ceropositivo:ev(funcion,x=3.03297097167559), numer;
plot2d([[discrete, [[-0.0329709716755897,f_ceronegativo],
[3.03297097167559, f_ceropositivo]]], funcion],[x,-1,4],[y,-3,3],
[xlabel, “Inversión en publicidad (miles de euros)”],
[ylabel, “Beneficios (miles de euros)”],
[style, points, lines], [color, red, green], [legend, false]);

Método de la secante

(1) Objetivo: Aproximarse numéricamente a la solución p de una función , tal que

(2) Condiciones: La función debe ser continua en el intervalo considerado.

(3) Descripción rápida: Partimos del método de Newton, que recordemos permite construir la sucesión donde:

Pero esa derivada puede ser compleja de calcular, por lo que podemos aproximarla de la siguiente forma:

 De este modo, obtenemos:

(4) Convergencia: Un criterio suficiente (pero no necesario) para que lo haga es que:

Si se cumple la fórmula anterior cuando está en un intervalo alrededor de la raíz entonces la solución converge en cualquier punto de ese intervalo. Para emplear este criterio, hemos de añadir el requerimiento de que la función sea dos veces derivable.

Normalmente, el método de la secante converge de forma más lenta que el de Newton, pero más rápido que otros, como el de la bisección.

(5) Estimación del error: Hay diversas opciones, aunque una de las más recomendables es el error relativo:

Implementación en Maxima

Vamos a implementar el algoritmo partiendo de los dos valores con los que inicíabamos el método de la bisección: [2.5,3.5], que actúan como las primeras semillas en el método de la secante:

f(i):=-x[i]^2+3*x[i]+0.1;
x[0]:2.5;
x[1]:3.5;
for i:2 thru 8 do
    (
    x[i]:x[i-1]-((x[i-1]-x[i-2])*f(i-1)/(f(i-1)-f(i-2))),
    print(i-2,x[i-2]));
datos: makelist([[i-2], x[i-2]], i, 2, 8);
matriz_resultados:apply(matrix,datos);

Le hemos dicho al programa que nos retorne en el output siguiente:

b390_2

Como se puede observar, el método converge rápido, y también podemos establecer criterios de parada con programaciones similares a las realizadas con los otros métodos.

Otra forma de hacerlo

El método de la secante, como su propio nombre indica, traza rectas secantes a la función y que pasan por el eje X. La idea es ir aproximando secantes entre puntos de la función hasta converger en la solución. Para ello, es útil considerar la ecuación de la recta que une los puntos (a, f(a)) y (b, f(b)):

Podemos verlo gráficamente si estipulamos que a y b sean 2.5 y 3.5, respectivamente, que son los valores empleados como intervalos en el método de la bisección:

funcion:-(x^2)+(3*x)+0.1;
a:2.5;
fa:(ev(funcion,x=a));
b:3.5;
fb:(ev(funcion,x=b));
secante:((fb-fa)*(x-a)/(b-a))+fa, ratsimp;
f_ceronegativo:ev(funcion,x=-0.03297097167558977), numer;
f_ceropositivo:ev(funcion,x=3.03297097167559), numer;
plot2d([[discrete, [[-0.0329709716755897,f_ceronegativo],
[3.03297097167559, f_ceropositivo]]], funcion, secante],[x,-1,4],[y,-3,3],
[xlabel, “Inversión en publicidad (miles de euros)”],
[ylabel, “Beneficios (miles de euros)”],
[style, points, lines, lines], [color, red, green, orange], [legend, false]);

b390_3Esta primera recta secante, corta a la función en f(a) y f(b), y ya nos quedamos bastante cerca de la solución. Simplemente hay que seguir trazando secantes iterando entre puntos consecutivos, es decir, el siguiente paso sería realizar el mismo cálculo pero tomando f(b) en lugar de f(a), y f(m) en lugar de f(b), siendo f(m) el valor de la función para el punto m (el punto de corte de la primera secante con el eje X). En Maxima, una forma de computarlo se puede encontrar en Ramírez (2008). Nosotros vamos a adaptar ese código de la siguiente manera:

secante(arg):=block([a:arg[1],b:arg[2],f:arg[3],xsec], xsec:(a*ev(f,x=b)-b*ev(f,x=a))/(ev(f,x=b)-ev(f,x=a)),
[b,xsec,f] );

ini:[2.5,3.5,-(x^2)+(3*x)+0.1];

for i:1 thru 10 do (print(float(ini[1])),ini:secante(ini));

El resultado es:

2.5
3.5
2.95
3.021739130434783
3.033284564740307
3.0329698187459
3.032970971557676
3.032970971675589

Conclusión

En este post hemos explicado brevemente en qué consiste el método de la secante, una de las variantes del método de Newton, que provee una convergencia más lenta, pero que tiene la ventaja de no tener que calcular derivadas en cada iteración.

Todos los posts relacionados




(#389). BÚSQUEDA DE SOLUCIONES (III): METODO DE NEWTON

[MONOTEMA] Seguimos en la búsqueda de los ceros de una función con aproximaciones numéricas. Tras ver el método de la bisección y el método del punto fijo,  vamos a continuar con el método de Newton, también conocido como Newton-Raphson. Lo haremos, como siempre, desde el punto de vista práctico, y siguiendo a Burden, Faires & Burden (2017).

Función de partida

Emplearemos la misma función que en los ejemplos anteriores, donde se relaciona la inversión en publicidad con los beneficios.

f_ceronegativo:ev(funcion,x=-0.03297097167558977), numer;
f_ceropositivo:ev(funcion,x=3.03297097167559), numer;
plot2d([[discrete, [[-0.0329709716755897,f_ceronegativo],
[3.03297097167559, f_ceropositivo]]], funcion],[x,-1,4],[y,-3,3],
[xlabel, “Inversión en publicidad (miles de euros)”],
[ylabel, “Beneficios (miles de euros)”],
[style, points, lines], [color, red, green], [legend, false]);

Método de Newton

(1) Objetivo: Aproximarse numéricamente a la solución p de una función , tal que

(2) Condiciones: La función debe ser continua y dos veces diferenciable en el intervalo considerado.

(3) Descripción rápida: Partimos del primer polinomio de Taylor para  expandido alrededor de  y evaluado en  . Suponemos que es una buena aproximación a , es decir, que la semilla inicial está relativamente cerca de la solución. Cuando , entonces, obviando el término de error de la serie de Taylor:

A partir de esa fórmula podemos construir la sucesión donde:

(4) Convergencia: Un criterio suficiente (pero no necesario) para que lo haga es que:

Si se cumple la fórmula anterior cuando está en un intervalo alrededor de la raíz entonces la solución converge en cualquier punto de ese intervalo. La convergencia es no lineal (cuadrática), y esta es una de las razones por las que es uno de los procedimientos más empleados.

(5) Estimación del error: Hay diversas opciones, aunque una de las más recomendables es el error relativo:

Implementación en Maxima

Vamos a implementar el algoritmo partiendo de 3 semillas diferentes, para ver que en función de la elección de ellas nos puede llevar a una de las dos soluciones posibles.  Para ello, empleamos un bucle con 10 iteraciones, tomando como semillas: 3, 0.1 y 8. Esto hace que evaluemos 3 funciones diferentes (f1, f2, f3) a la hora de analizar su convergencia.

f1(i):= -x1[i]^2+3*x1[i]+0.1;
derivada_f1(i):=-(2*x1[i])+3;
x1[0]:3;
for i:1 thru 10 do(
x1[i]:x1[i-1]-(f1(i-1)/derivada_f1(i-1)));
f2(i):= -x2[i]^2+3*x2[i]+0.1;
derivada_f2(i):=-(2*x2[i])+3;
x2[0]:0.2;
for i:1 thru 10 do(
x2[i]:x2[i-1]-(f2(i-1)/derivada_f2(i-1)));
f3(i):= -x3[i]^2+3*x3[i]+0.1;
derivada_f3(i):=-(2*x3[i])+3;
x3[0]:8;
for i:1 thru 10 do(
x3[i]:x3[i-1]-(f3(i-1)/derivada_f3(i-1)));
datos: makelist([i,”Iteración”,x1[i],f1(i-1), x2[i],f2(i-1),
    x3[i], f3(i-1)], i, 1,10);
matriz_resultados:apply(matrix,datos);

Le hemos dicho al programa que nos retorne en el output siguiente:

 b389_2

El método converge bastante rápido, incluso cuando la semilla se aleja de la raíz, y como puede observarse en función de la elección de la semilla converge hacia una de las dos soluciones.

Criterio de parada

Al igual que hicimos en el método de la bisección, podemos estipular que nos baste con una tolerancia de error, es decir, pedirle al programa que pare las iteraciones cuando se cumpla una condición determinada, por ejemplo, que el error relativo sea menor que una cota que fijemos. De este modo, si fijamos que el error relativo sea como máximo 0.0001:

f1(i):= -p[i]^2+3*p[i]+0.1;
derivada_f1(i):=-(2*p[i])+3;
p[0]:8;
for i:1 while Error_rel[i]>0.0001 do(
p[i]:p[i-1]-(f1(i-1)/derivada_f1(i-1)),
Error_abs[i]:abs(p[i]-p[i-1]),
Error_rel[i]: Error_abs[i]/abs(p[i]),
    print(i,p[i],Error_rel[i]));

Le hemos dicho al programa que fije un error relativo mayor que 0.0001, por lo que parará en esa iteración, y de esta manera sabremos que en la siguiente iteración el error relativo estará por debajo de esa cota. Si ejecutamos esas instrucciones en Maxima vemos que para en la iteración número 5, por lo que incluso con semillas alejadas de la raíz la convergencia puede ser muy rápida.

Conclusión

En este post hemos explicado brevemente en qué consiste el método de Newton, un algoritmo muy empleado en diversas aplicaciones estadísticas, como para obtener los estimadores que minimizan una función de error (hacen que su derivada sea cero, por lo que es realmente obtener la aproximación numérica a la raíz de una función). El método de Newton tiene, además, diversas variantes que comentaremos en futuras entradas.

Todos los posts relacionados




(#388). BÚSQUEDA DE SOLUCIONES (II): METODO DEL PUNTO FIJO

[MONOTEMA] Continuamos buscando los ceros de una función con aproximaciones numéricas. Tras ver el método de la bisección vamos a continuar con el método del punto fijo. Lo haremos, como siempre, desde el punto de vista práctico, y siguiendo a Burden, Faires & Burden (2017).

Función de partida

Emplearemos la misma función que en el método de la bisección, donde se relaciona la inversión en publicidad con los beneficios.

Lo primero que vamos a hacer es representar esa función en Maxima, tal y como hicimos en el post anterior:

funcion: -x^2+3*x+0.1;
plot2d(funcion,[x,0,4],[y,-3,3],
[xlabel, “Inversión en publicidad (miles de euros)”], [ylabel, “Beneficios (miles de euros)”], [color, green]);

Pero ahora vamos a expandir el eje de las X, porque nos interesa por motivos ilustrativos representar el corte en el tramo negativo:

funcion: -x^2+3*x+0.1;
plot2d(funcion,[x,-1,4],[y,-3,3],
[xlabel, “Inversión en publicidad (miles de euros)”], [ylabel, “Beneficios (miles de euros)”],
    [color, green]);

b388_2Y vemos que tenemos dos ceros (las dos veces que corta al eje X). Dijimos que la parte negativa no estaba dentro del dominio de la función, pero nos interesa representarla para ver lo que ocurre con la iteración del método del punto fijo.

De este modo tenemos dos ceros de la función  y .

Podemos representar esos puntos en el gráfico anterior de la siguiente manera:

f_ceronegativo:ev(funcion,x=-0.03297097167558977), numer;
f_ceropositivo:ev(funcion,x=3.03297097167559), numer;
plot2d([[discrete, [[-0.0329709716755897,f_ceronegativo],
                [3.03297097167559, f_ceropositivo]]], funcion],[x,-1,4],[y,-3,3],
     [xlabel, “Inversión en publicidad (miles de euros)”],
     [ylabel, “Beneficios (miles de euros)”],
[style, points, lines], [color, red, green], [legend, false]);

Método del punto fijo

(1) Objetivo: Aproximarse numéricamente a la solución p de una función , tal que

(2) Condiciones: La función debe ser continua en el intervalo considerado.

(3) Descripción rápida: Si la función , entonces podemos construir la función . De este modo, hemos de escribir nuestra función original como . Esto significa para que haya un cero , es decir, que haya una x despejada en uno de los lados de la igualdad.

(4) Convergencia: No siempre converge. Un criterio suficiente (pero no necesario) para que lo haga es que  , donde es una constante positiva. Si el método converge cuadráticamente.

(5) Estimación del error: Hay diversas opciones, aunque una de las más recomendables es el error realtivo:

Implementación en Maxima

El primer paso, por tanto es crear la expresión a partir de la función original. Suele haber varias formas de hacerlo, en este caso:

Por tanto, tenemos que:

Dibujemos con Maxima todas las funciones:

f(x):= -x^2+3*x+0.1;
g1(x):=((x^2)/3)-(0.1/3);
g2(x):=(-0.1/(-x+3));
plot2d([f(x),g1(x),g2(x), x],[x,-1,4],[y,-3,5]
, [xlabel, “Inversión en publicidad (miles de euros)”],
     [ylabel, “Beneficios (miles de euros)”],
[style, lines, lines, lines], [color, green, red, blue, black],
     [legend, “f(x)”, “g1(x)”, “g2(x)”, “y=x”]);

b388_4Como se puede apreciar, la recta corta a  y en el valor de x que corresponde a las soluciones de . No vamos a escoger porque tenemos una singularidad en x=3, y por tanto la función no es continua en un entorno cercano a la solución positiva.

Vamos a programar con Maxima el algoritmo, que es bastante sencillo. Para ello, empleamos un bucle con 15 iteraciones, comenzando en 3 como semilla inicial.

g1(i):=(((x[i])^2)/3)-(0.1/3);
x[0]:3;
for i:1 thru 15 do(
x[i]:g1(i-1),print(i,”_Iteración”, x[i-1],g1(i-1)));
datos: makelist([i,”Iteración”,x[i-1],g1(i-1)], i, 1,15);
matriz_resultados:apply(matrix,datos);

Le hemos dicho al programa que nos retorne en el output siguiente:

Como puede observarse, cada imagen de x es el valor de prueba de la siguiente iteración, hasta aproximarnos a la solución. Pero, atención, la solución es negativa, la que no nos interesaba por no tener sentido en nuestro ejemplo.

¿Por qué no converge a la solución positiva pese a que hemos escogido una semilla muy cercana a la solución real?

Si cogemos un intervalo en las proximidades de la solución y lo comparamos con , veremos que no se cumple esa condición. Por ejemplo, en el intervalo [2.5, 3.5] que es el que usamos en el método de la bisección, podemos coger x=3 y comprobar lo que sucede:

En cualquier otro punto de ese intervalo no se cumple el criterio suficiente para la convergencia. No obstante, podría haber convergencia (ya que el criterio no es necesario), pero hemos comprobado que no es así.

El método del punto fijo se va hacia la solución en el eje negativo de las X, donde sí que en un intervalo cercano se cumple esa condición. Por ejemplo, en el intervalo [-1,0].

Aceleración de Aitken

La convergencia se puede acerlerar imlpementando el método de Aitken, donde se construye una nueva sucesión a partir de los 2 primeros términos de la original (la calculada por el método numérico, en este caso del punto fijo):

Podemos programar en Maxima el método de Aitken en conjunción con el de punto fijo, de la siguiente forma:

g1(i):=(((p[i])^2)/3)-(0.1/3);
p[0]:3;
for i:1 thru 15 do(
p[i]:g1(i-1),
p_aitken[i+2]:p[i-1]-(((p[i]-p[i-1])^2)/(p[i+1]-2*p[i]+p[i-1])));
datos: makelist([i,”Iteración”,p[i-1],g1(i-1), p_aitken[i+2]], i, 1,15);
matriz_resultados:apply(matrix,datos);

Al implementar ese código en Maxima, se puede ver que la convergencia es más rápida.

Conclusión

En este post hemos dado unas pinceladas sobre la implementación del método de punto fijo, encontrando que pueden existir problemas de convergencia para determinadas funciones. Asimismo, la aceleraciónde Aitken puede emplearse para intentar aproximarse a la solución con menor coste computacional.

Todos los posts relacionados




(#387). BÚSQUEDA DE SOLUCIONES (I): MÉTODO DE LA BISECCIÓN

[MONOTEMA] Vamos a comenzar una serie de posts sobre la búsqueda de soluciones  numéricas (también llamadas raíces o ceros)  para una función. El objetivo es que los alumnos de marketing se familiaricen con técnicas matemáticas y de computación que les ayuden a resolver problemas más complejos en el futuro. Para ello vamos a emplear el software Maxima, de libre acceso.

Función de partida

Imaginemos que conocemos la forma de una función que relaciona la inversión en publicidad con los beneficios. En entradas anteriores hemos explicado que esa relación es probablemente no lineal, y que puede adoptar diversas formas. Vamos a partir de la siguiente función:

 Lo primero que vamos a hacer es representar esa función en Maxima:

funcion: -x^2+3*x+0.1;
plot2d(funcion,[x,0,4],[y,-3,3],
     [xlabel, “Inversión en publicidad (miles de euros)”], [ylabel, “Beneficios (miles de euros)”]);

En esta hipotética situación, el analista de marketing puede computar el máximo y el cero de la función, de la siguiente forma:

derivada: diff(funcion, x);
maximo:solve(derivada,x), numer;
cero: solve(funcion,x),numer;
f_maximo: ev(funcion,x=1.5), numer;
f_cero:ev(funcion,x=3.03297097167559), numer;
plot2d([[discrete, [[1.5,f_maximo], [3.03297097167559, f_cero]]], funcion],[x,0,4],[y,-3,3],
     [xlabel, “Inversión en publicidad (miles de euros)”], [ylabel, “Beneficios (miles de euros)”],
[style, points, lines], [color, red, green], [legend, false]);

De este modo tenemos que el máximo es , y el cero es cuando (hay otra solución al ser una ecuación de segundo grado pero no nos interesa porque se refiere a una inversión en publicidad negativa, lo cual no es posible).

Podemos representar esos puntos en el gráfico anterior de la siguiente manera:

b387_3Es evidente que nos interesa conocer el valor de inversión que maximiza los beneficios, pero también nos puede resultar atractivo conocer el punto a partir del cual incrementar más la inversión produce pérdidas. Y este es el objetivo en el que nos vamos a centrar.

Para ello implementar el método de la bisección, que puede consultarse en Burden, Faires & Burden (2017).

Método de la bisección

(1) Objetivo: Aproximarse numéricamente a la solución p de una función , tal que 

(2) Condiciones: La función  debe ser continua en el intervalo  considerado, donde  y son de distinto signo, es decir, .

(3) Descripción rápida: Escoger un intervalo  a priori donde e ir evaluando soluciones   en la mitad de ese intervalo y de los siguientes   hasta que se llegue a un cierto número de iteraciones  o a un valor prefijado de error . De este modo, se trata de ir evaluando la función en los valores de medios de cada intervalo sucesivo (donde la anchura de cada intervalo es la mitad del anterior) hasta aproximarse a la solución. Por el teorema de Bolzano, sabemos que vamos a encontrar esa raíz.

(4) Convergencia: Siempre converge a una solución, aunque su velocidad puede ser lenta.

(5) Estimación del error: Hay diversas opciones, aunque una de las más recomendables es el error realtivo:

Implementación en Maxima

Vamos a programar con Maxima el algoritmo de la bisección. Para ello, empleamos un bucle con 15 iteraciones. La inspección gráfica nos permite estipular un intervalo inicial [2.5,3.5]:

a[0]:2.5;
b[0]:3.5;
for i:1 thru 15 do(
p[i]:a[i-1]+((b[i-1]-a[i-1])/2),
f[i]: -(p[i])^2+3*(p[i])+0.1,
g[i]: -(a[i-1])^2+3*(a[i-1])+0.1,
if g[i]*f[i]<0
then (b[i]:p[i],
        a[i]:a[i-1])
else (a[i]:p[i],
                b[i]:b[i-1]),
Error_abs[i]:abs(p[i+1]-p[i]),
Error_rel[i]: Error_abs[i]/abs(p[i+1]));

datos: makelist([i,”Iteración”,a[i],b[i],p[i],f[i], Error_abs[i-1], Error_rel[i-1]], i, 1,15);
matriz_resultados:apply(matrix,datos);

Le hemos dicho al programa que nos retorne en el output siguiente:

b387_4_1Como puede apreciarse en las dos primeras columnas de resultados tenemos los extremos de cada intervalo del algoritmo. En la tercera está el valor de p, que siempre es la mitad del valor de esos extremos. la cuarta columna evalúa la función en cada p considerado, mientras que las columnas quinta y sexta corresponden, respectivamente, al error absoluto y al relativo.

Podemos ver aquí como las iteraciones nos acercan a la solución:

matriz_t: transpose(matriz_resultados);

p_list: matriz_t[5];

fp_list: matriz_t[6];

funcion: -x^2+3*x+0.1;
plot2d([[discrete, p_list,fp_list],funcion],[x,0,4],[y,-3,3],
     [xlabel, “Inversión en publicidad (miles de euros)”],
     [ylabel, “Beneficios (miles de euros)”],
[style, points, lines], [color, red, green], [legend, false]);

b387_5 Y también es interesante comprobar cómo disminuye el error relativo con el número de iteraciones:

Error_rel_list: matriz_t[8];
Iteraciones_list:matriz_t[1];

plot2d([discrete, Iteraciones_list,Error_rel_list],[x,1,15],[y,0,0.09],
     [xlabel, “Iteraciones”],
     [ylabel, “Error relativo”],
[style, lines], [color, red], [legend, false]);

b387_6 Criterio de parada

Podemos estipular que nos baste con una tolerancia de error, es decir, pedirle al programa que pare las iteraciones cuando se cumpla una condición determinada, por ejemplo, que el error relativo sea menor que una cota que fijemos. De este modo, si fijamos que el error relativo sea como máximo 0.001:

a[0]:2.5;
b[0]:3.5;
for i:1 while Error_rel[i]>0.001
    do(
p[i]:a[i-1]+((b[i-1]-a[i-1])/2),
f[i]: -(p[i])^2+3*(p[i])+0.1,
g[i]: -(a[i-1])^2+3*(a[i-1])+0.1,
if g[i]*f[i]<0
then (b[i]:p[i],
        a[i]:a[i-1])
else (a[i]:p[i],
                b[i]:b[i-1]),
Error_abs[i]:abs(p[i+1]-p[i]),
Error_rel[i]: Error_abs[i]/abs(p[i+1]),
        print(i,p[i], f[i],Error_rel[i]));

Le hemos dicho al programa que fije un error relativo mayor que 0.001, por lo que parará en esa iteración, y de esta manera sabremos que en la siguiente iteración el error relativo estará por debajo de esa cota.

Alternativamente podríamos conocer de forma conservadora las iteraciones necesarias para una precisión determinada empleado la siguiente fórmula:

Conclusión

En este post hemos explicado con carácter eminentemente práctico el método de la bisección para buscar una aproximación numérica a la solución de una ecuación. Es importante señalar que los errores absoluto y relativo que hemos definido en las programaciones con Maxima se refieren siempre a las estimaciones de la solución en sus respectivas iteraciones, pero no al error absoluto y relativo tomado cuando se conoce valor real de la solución.

En próximos posts, comentaremos otros métodos diferentes.

Todos los posts relacionados