1. No category

Universidad Politécnica de Madrid–Escuela Técnica Superior de Ingenieros Industriales
Grado en Ingeniería en Tecnologías Industriales. Curso 2015-2016-3º
Matemáticas de Especialidad–Ingeniería Eléctrica
Resolución de ecuaciones de una variable
José Luis de la Fuente O’Connor
[email protected]
[email protected]
Clase_reso_ecuaciones_2016.pdf
1/75
2/75
Índice
El problema
Método de la bisección
Iteración de punto fijo
Velocidad de convergencia
Método de Newton
Variantes del método de Newton
Métodos sin derivadas
3/75
El problema
En términos matemáticos,
dada f W R ! R; hallar un xN tal que f .x/
N D 0.
La función f se supone continua y derivable en algún conjunto abierto de R.
Los algoritmos que estudiamos son procedimientos iterativos en los que se
avanza progresivamente hacía una solución.
Se utilizará uno u otro dependiendo de la complejidad de evaluar la función y de
la disponibilidad de derivada, o de la dificultad de calcularla.
4/75
Método de la Bisección
Atribuido a Bernard Bolzano, Praga, Bohemia (hoy Rep. Checa), 1781-1848.
Idea: Si f es continua en un intervalo Œa; b en cuyos extremos cambia de
signo, f .a/f .b/ < 0, existirá un c, a < c < b, en el que f .c/ D 0.
Reduciendo Œa; b se llegará a acotar el valor de c tanto como se desee.
Origen:
5/75
Teorema Del valor intermedio Si f W R ! R es una función continua en Œa; b y
f .a/ x f .b/ o f .b/ x f .a/, existe un c, a c b, donde f .c/ D x.
Mecánica del procedimiento:
I. Comenzar con unos u D f .a/ y v D f .b/, tales que uv < 0. Determinar
c D 21 .a C b/ y w D f .c/.
Si f .c/ D 0 se ha llegado a la solución;
si no, se cumplirá que wu < 0 o wv < 0.
II. Si wu < 0, la solución estará en Œa; c; si wv < 0, en Œc; b.
III. Se estudia el nuevo intervalo donde esté la raíz, reduciéndolo secuencialmente en
dos mitades hasta que se estreche lo que se desee el intervalo de confinamiento
que contenga la solución.
f (a)
f (a)
f (x )
f (c)
f (x )
Casos que se pueden presentar al comenzar
[
]
c
a
f (c)
]
[
a
b
f (b)
[
c
]
b
f (b)
6/75
Si el intervalo con que se empieza el proceso iterativo, Œa0; b0, contiene una
solución r, usando como estimación de ésta c0 D .a0 C b0/=2, se tendrá que
e0 D jr
En cualquier iteración, ei D jr
c0 j ci j b0
bi ai
;
2
a0
2
:
i D 0; 1; 2; : : :
Teorema Al aplicar el método de la bisección a una función f W R ! R, continua en
un intervalo Œa; b, en el que f .a/f .b/ < 0, después de n iteraciones, en las que se habrán
evaluado la función n C 2 veces, se habrá obtenido un valor de la solución cn tal que su
error
b a
jr cn j nC1 ;
2
donde r es el valor real de la solución.
Definición Una solución es correcta en p posiciones decimales si el error es menor que
0;5 10 p .
7/75
Esta function lleva a cabo el método de la bisección.
function sol=Bisec_0(fun,a,b,tol)
% Método de la Bisección para resolver f(x)=0
if nargin<4, tol=sqrt(eps); end
fa=fun(a); fb=fun(b);
if sign(fa)*sign(fb)>=0, error(’ El intervalo (a,b) no contiene la solución\n’), end
while abs(b-a)/2>tol
c=(a+b)/2;
fc = fun(c);
if fc==0, break, end
if sign(fc)*sign(fa)<0
b=c; fb=fc;
else
a=c; fa=fc;
% No es necesario hacer fa=fc;
end
fprintf(’ %17.15f %17.15f\n’,a,b);
end
sol = (a+b)/2;
end
8/75
Ejemplo Resolvamos x sen.x/ 1 D 0 en 1 x 2 (radianes). Para ello,
suministramos los extremos del intervalo Œa; b y la propia función que resolver:
>> f=@(x) x*sin(x)-1;
>> Bisec_0(f,1,2}
Los resultados son los de esta tabla.
k
a
b
k
a
b
1
2
3
4
5
6
7
8
9
10
11
12
1,0000000000000
1,0000000000000
1,0000000000000
1,0625000000000
1,0937500000000
1,1093750000000
1,1093750000000
1,1132812500000
1,1132812500000
1,1132812500000
1,1137695312500
1,1140136718750
1,5000000000000
1,2500000000000
1,1250000000000
1,1250000000000
1,1250000000000
1,1250000000000
1,1171875000000
1,1171875000000
1,1152343750000
1,1142578125000
1,1142578125000
1,1142578125000
13
14
15
16
17
18
19
20
21
22
23
24
25
1,1141357421875
1,1141357421875
1,1141357421875
1,1141510009765
1,1141510009765
1,1141548156738
1,1141567230224
1,1141567230224
1.1141567230224
1.1141569614410
1.1141570806503
1.1141571402549
1.1141571402549
1,1142578125000
1,1141967773437
1,1141662597656
1,1141662597656
1,1141586303710
1,1141586303710
1,1141586303710
1,1141576766967
1.1141571998596
1.1141571998596
1.1141571998596
1.1141571998596
1.1141571700572
9/75
En la figura se representa cómo procede el método para llegar a la solución.
Intervalo inicial
Error hacia delante y error hacia atrás de un algoritmo
10/75
El proceso que se lleva a cabo para obtener la solución de un problema
mediante un algoritmo de bucle abierto se esquematiza así:
Datos que
definen el !
problema
Proceso de
solución o
algoritmo
! Solución
El problema lo constituyen los datos que se proporcionan a esa función de
transferencia que es el algoritmo.
La solución es la información que proporciona la salida.
En lo que analizamos aquí, el problema es una
ecuación de una variable;
el algoritmo, cualquiera de los que estudiamos.
11/75
Definición Si f .x/ es una función y r una raíz de ella, f .r/ D 0, supongamos que
xa es una aproximación de r obtenida por un algoritmo. El error hacia atrás del algoritmo
es jf .xa /j y el error hacia delante jr xa j.
El error hacia atrás está en el lado izquierdo del esquema anterior, o a la
entrada. Es la cantidad que tendría que cambiar el problema (la función) para
hacer que la ecuación se equilibre con la aproximación xa de la salida. Es
jf .xa /j.
El error hacia delante está en el lado derecho, o a la salida (solución del
problema). Es la cantidad que tendría que cambiar la solución aproximada para
que sea correcta, es decir, jr xa j.
8
D 0. Podemos comprobar que
Ejemplo Estudiemos x 3 2x 2 C 34 x 27
f .0/f .1/ D . 8=27/.1=27/ < 0. Apliquemos el método de la bisección.
>> f1=@(x) x^3-2*x^2+x*4/3-8/27;
>> xc=Bisec_0(f1,0,1)
Se obtiene lo siguiente:
>> xc=Bisec_0(f1,0,1)
0.500000000000000 1.000000000000000
0.500000000000000 0.750000000000000
0.625000000000000 0.750000000000000
0.625000000000000 0.687500000000000
0.656250000000000 0.687500000000000
0.656250000000000 0.671875000000000
0.664062500000000 0.671875000000000
0.664062500000000 0.667968750000000
0.666015625000000 0.667968750000000
0.666015625000000 0.666992187500000
0.666503906250000 0.666992187500000
0.666503906250000 0.666748046875000
0.666625976562500 0.666748046875000
0.666625976562500 0.666687011718750
0.666656494140625 0.666687011718750
0.666656494140625 0.666671752929688
xc =
0.666664123535156
12/75
>> xc=Bisec_0(f1,0,1,eps)
0.500000000000000 1.000000000000000
0.500000000000000 0.750000000000000
0.625000000000000 0.750000000000000
0.625000000000000 0.687500000000000
0.656250000000000 0.687500000000000
0.656250000000000 0.671875000000000
0.664062500000000 0.671875000000000
0.664062500000000 0.667968750000000
0.666015625000000 0.667968750000000
0.666015625000000 0.666992187500000
0.666503906250000 0.666992187500000
0.666503906250000 0.666748046875000
0.666625976562500 0.666748046875000
0.666625976562500 0.666687011718750
0.666656494140625 0.666687011718750
0.666656494140625 0.666671752929688
xc =
0.666664123535156
La raíz es 0;6666666666666 : : :. Aunque se aumente la precisión, el proceso es
el mismo y sólo consigue 6 dígitos significativos de la misma, ¿por qué?
Lo que ocurre no es culpa del método de la bisección sino de la incapacidad de
la aritmética de la máquina para calcular la función f con precisión suficiente
cerca de la raíz.
Cualquier otro método que se base en esta misma aritmética de maquina estará
destinado al mismo “fracaso”.
13/75
1.3 Limits of Accuracy | 45
La razón se ve en la figura.
(a)
(b)
Si la aritmética de la máquina muestra que la función es igual a cero en un valor
Figure 1.7 The shape of a function near a multiple root. (a) Plot of f (x) =
que no es una
raíz, no hay manera de que el método pueda recuperar esto.
x 3 − 2x 2 + 4/3x − 8/27. (b) Magnification of (a), near the root r = 2/3. There are
14/75
El error hacia atrás es cercano a máq 2;2 10 16, mientras que el error hacia
delante es aproximadamente 10 5. Como el error hacia atrás no puede
disminuirse por debajo de un error relativo por debajo del épsilon de la máquina,
tampoco es posible disminuir el error hacia delante.
Hay que destacar que este ejemplo es bastante especial pues la función tiene
una raíz triple en r D 2=3.
3
8
2
4
f .x/ D x 3 2x 2 C
D x
:
3 27
3
Definición Si una función continua y derivable m veces tiene en r una raíz y 0 D
f .r/ D f 0 .r/ D D f .m 1/ , pero f .m/ ¤ 0, se dice que f tiene una raíz de multiplicidad
m en r. Una raíz es simple si la multiplicidad es igual a uno y múltiple si es mayor que uno.
15/75
Iteración de punto fijo
Observemos:
>> cos(1)
ans = 0.540302305868140
>> cos(ans)
ans = 0.857553215846393
>> cos(ans)
ans = 0.654289790497779
>> cos(ans)
ans = 0.793480358742566
>> cos(ans)
ans = 0.701368773622757
>> cos(ans)
ans = 0.763959682900654
>> cos(ans)
ans = 0.722102425026708
>> cos(ans)
ans = 0.750417761763761
>> cos(ans)
ans = 0.731404042422510
>> cos(ans)
ans = 0.744237354900557
>> cos(ans)
ans = 0.735604740436347
>> cos(ans)
ans = 0.741425086610109
>> cos(ans)
ans = 0.737506890513243
>> cos(ans)
ans = 0.740147335567876
>> cos(ans)
ans = 0.738369204122323
>> cos(ans)
ans = 0.739567202212256
>> cos(ans)
ans = 0.738760319874211
>> cos(ans)
ans = 0.739303892396906
>> cos(ans)
ans = 0.738937756715344
>> cos(ans)
ans = 0.739184399771494
>> cos(ans)
ans = 0.739018262427412
>> cos(ans)
ans = 0.739130176529671
>> cos(ans)
ans = 0.739054790746917
>> cos(ans)
ans = 0.739105571926536
>> cos(ans)
ans = 0.739071365298945
>> cos(ans)
ans = 0.739094407379091
>> cos(ans)
ans = 0.739078885994992
>>
16/75
La sucesión de puntos de iterar la función coseno parece que converge a un
punto r.
Definición Un número real r es un punto fijo de f .x/ si f .r/ D r.
La iteración de punto fijo sirve para resolver problemas de punto fijo f .x/ D x.
Se puede usar para resolver f .x/ D 0 si éste se puede expresar como g.x/ D x.
Esquema general de la Iteración de Punto Fijo
Dados Un x WD x0 y una t ol . Hacer f ound WD f al se
while (not f ound ) and (k < kmax )
Hacer xk WD f .xk 1 )
if (xk xk 1 < t ol ), f ound WD true , end
k WD k C 1
end
17/75
Este código de Matlab hace el trabajo.
function xc=fpi(g,x0,k)
% fpi: Calcula la solución de f(x)=x; it. punto fijo
% Input: inline function f, punto de partida x0,
%
número de iteraciones k
% Output: Solución x
x(1)=x0;
for i=1:k
x(i+1)=g(x(i));
end
xc=x(k+1);
end
Para el caso de f .x/ D cos.x/, haciendo hasta treinta iteraciones, sería:
>> f=@(x)cos(x);
>> xc = fpi(f,1,15)
xc =
0.738369204122323
>> xc = fpi(f,1,30)
xc =
0.739087042695332
18/75
Ejemplo Resolvamos la ecuación x 3 C x
formas con el formato g.x/ D x:
x 3;
xD1
en este caso g.x/ D 1
También
1 D 0. Expresémosla de varias
x 3.
xD
p
3
1
x;
p
siendo aquí g.x/ D 3 1 x.
Si añadimos 2x 3 a los lados de la ecuación, también,
3x 3 C x 1 D 2x 3
3x 2 C 1 x D 1 C 2x 3
1 C 2x 3
x D
1 C 3x 2
siendo g.x/ D .1 C 2x 3/=.1 C 3x 2/.
ving Equations
19/75
to the diagonal
y = x. geometrías
This geometric
of a Fixed-Point
i ) across
Enheight
esta g(x
figura
se pueden
ver las line
distintas
deillustration
las y D g.x/
y los
Iteration is called a cobweb diagram.
primeros pasos de la iteración de punto fijo con ellas, partiendo de un mismo x0.
y
y
y
1
1
1
x2 x0 r
x1 1
x
x0 r x1 1
x2
x
x0 r
1
x
Figure 1.3 Geometric view of FPI. The fixed point is the intersection of g(x) and the
diagonal line. Three examples of g(x) are shown together with the first few steps of
La primera
diverge, la segunda converge y la tercera también converge, pero
FPI. (a) g(x) = 1 – x3 (b) g(x) = (1 – x)1/3 (c) g(x) = (1 + 2x3 )/(1 + 3x2 )
mucho más rápidamente. ¿Por qué?
In Figure 1.3(a), the path starts at x0 = 0.5, and moves up to the function and horizontal
to the point (0.875, 0.875) on the diagonal, which is (x1 , x1 ). Next, x1 should be substituted
Aparentemente,
algo the
tiene
que
condone
ellofor
laxpendiente
de la función, g 0.x/,
into g(x). This is done
same
wayver
it was
0 , by moving vertically to the function.
0.3300, and after moving horizontally to move the y-value to an x-value,
This del
yields
x2 ≈ fijo.
cerca
punto
we continue the same way to get x3 , x4 , . . . . As we saw earlier, the result of FPI for this g(x)
20/75
Consideremos ahora la ecuación
x3
sen.x/ D 0:
Al tratarse de una función senoidal, habrá varios puntos para los cuales
f .x/ D 0. Calcularemos el más próximo a x D 1.
Si seguimos una estrategia de punto fijo, buscaremos una función g W R ! R y
un xN tal que xN D g.x/,
N y aplicaremos un procedimiento iterativo a partir de la
relación de recurrencia xkC1 D g.xk /.
¿Qué formas de x D g.x/, y por tanto relación de recurrencia, se nos pueden
ocurrir para utilizar?
21/75
La primera podría ser x D
sería
p
3
sen.x/: La relación de recurrencia correspondiente
xkC1 D
p
3
sen.xk /:
Si comenzamos el proceso iterativo desde x0 D 1 (radianes), se tendrá que:
p
x1 D 3 sen.x0/ D 0;944
p
x2 D 3 sen.0;944/ D 0;932
p
x3 D 3 sen.0;932/ D 0;929
:::
La solución converge a xN D 0;92862.
22/75
También podemos hacer x D
sen.x/
,
x2
por lo que la relación de recurrencia sería
xkC1 D
sen.xk /
:
xk2
Partiendo de x0 D 1, se obtienen los puntos
k
0
1
2
3
4
::
:
El proceso diverge.
xk
1; 000
0; 841
1; 053
0; 783
1; 149
::
:
Veamos geométricamente qué ocurre en estos dos procesos:
El itinerario de la parte izquierda es el de xkC1 D
El que define
p
3
sen.xk /.
sen.xk /
xk2
genera lo que se denomina una tela de araña entre la recta y D x y la
función y D g.x/.
xkC1 D
23/75
Si se analiza el comportamiento de diversas relaciones de recurrencia se puede
constatar una relación directa entre el comportamiento del proceso iterativo y
las pendientes de g en el entorno de x:
N
Si jg 0.x/j
N D v < 1 y el punto de partida está cerca de x,
N el proceso converge
linealmente a velocidad v.
Si jg 0.x/j
N > 1, diverge.
En el caso de xkC1 D
p
3
sen.xk /,
.sen.x//
g .x/ D
3
En xN 0;929, g 0.0;929/ 0;23.
0
2=3
cos.x/:
Por el contrario, en xkC1 D sen.xk /=xk2 ,
cos.x/
x2
1;23.
g 0.x/ D
En xN 0;929, g 0.0;929/ 2
sen.x/
:
x3
24/75
25/75
Definición Sea ei el error en el paso i de un método iterativo. Si
ei C1
lim
D S < 1;
i !1 ei
se dice que el método converge linealmente con razón S.
Teorema
Si f W R ! R es una función continua y derivable, f .r/ D r y
S D jf .r/j < 1, la iteración de punto fijo converge linealmente a r, con razón S, para
estimaciones iniciales lo suficientemente próximas a r.
0
Velocidad de convergencia
26/75
Definición Sea una sucesión fxi g1
i D0 , xi 2 R, convergente a x . El orden de conver-
gencia de fxi g es el máximo de los números no negativos r que satisface
jxi C1 x j
0 lim
< 1:
i!1 jxi
x jr
El valor del límite, ˇ, se conoce como razón, o constante de error asintótico.
Si r D 1, la sucesión se dice que converge linealmente; si r D 2, se dice que lo
hace cuadráticamente; si r D 3, cúbicamente, etc.
Definición Un método iterativo xi C1 D f .xi /; i D 1; 2; : : :, que parte de un punto
x0 , se dice que tiene una velocidad de convergencia de orden r cuando la sucesión fxi g1
i D0
converge con orden r hacia la solución x D f .x/.
jxi C1 x j
limi!1 jxi x j
Si la convergencia es lineal y el valor de
convergencia se dice superlineal.
Ejemplo consideremos la sucesión escalar definida por
27/75
D ˇ es cero, la
k
xk D c 2 ;
donde c cumple 0 c < 1. La sucesión converge a cero. Calculemos su orden
de convergencia:
kC1
c2
jxkC1 0j
lim
D lim kC1 D 1:
2
k!1 jxk
k!1 c 2
0j
Es decir, converge cuadráticamente a 0.
La convergencia cuadrática quiere decir, grosso modo, que en las proximidades
del límite o solución el número de dígitos significativos que aporta cada paso del
proceso al valor de la solución es el doble que el anterior.
28/75
En la columna 2 de la tabla que sigue se pueden ver los distintos puntos de la
sucesión del primer ejemplo analizado para c D 0;99.
k
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
k
c2
.c D 0; 99/
0,9900000000000000
0,9801000000000000
0,9605960099999999
0,9227446944279201
0,8514577710948755
0,7249803359578534
0,5255964875255620
0,2762516676992083
0,0763149839065938
0,0058239767686636
0,0000339187054019
0,1150478576143195E-08
0,1323600954164474E-17
0,1751919485865107E-35
0,3069221884953861E-71
0,9420122979079730E-143
0,8873871694098596E-286
c2
k
.c D 2; 2/
2,200000000000000
1,483239697419133
1,217883285630907
1,103577494166543
1,050512967157732
1,024945348376065
1,012395845692812
1,006178833852518
1,003084659364561
1,001541142122759
1,000770274400054
1,000385063063246
1,000192513000995
1,000096251868287
1,000048124776146
1,000024062098581
1,000012030976918
1=k k
1,000000000000000
0,250000000000000
0,037037037037037
0,003906250000000
0,000320000000000
0,000021433470507
0,000001214265678
0,000000059604644
0,000000002581174
0,100000000000000E-10
0,350493899481392E-12
0,112156654784615E-13
0,330169095523011E-15
0,899927452978128E-17
0,228365826052116E-18
0,542101086242752E-20
0,120883864830239E-21
29/75
Consideremos la sucesión convergente a 1 que define,
k
xk D c 2 ;
con c 0.
Analicemos su orden de convergencia:
.kC1/
c2
jxkC1 1j
D lim
lim
k!1 c 2
k!1 jxk
1j
D lim
k!1
k
c 2 .kC1/
1
1
c2
.kC1/
1
1
c 2 .kC1/
C1
1
D
:
.kC1/
k!1 c 2
2
C1
Converge linealmente: Esta convergencia significa que en cada iteración se
añaden un número de dígitos constante a la solución final; concretamente,
log10 ˇ dígitos por iteración.
D lim
1
Analicemos por último la sucesión que define
1
xk D k :
k
Converge a cero. En la columna 4 de la tabla se pueden ver los primeros
puntos de esta sucesión.
Estudiemos su orden de convergencia:
1
1
jxkC1j
.k C 1/kC1
lim
D lim
D lim kC1 D 0:
1
k!1 jxk j
k!1
k!1
1
k
1
C
k
k
k
Es decir, converge superlinealmente a cero.
La convergencia superlineal significa que cada iteración añade un número
creciente de dígitos a la solución final; concretamente r dígitos más que la
iteración precedente.
30/75
31/75
La siguiente tabla resume las diferencias entre varias velocidades de
convergencia.
k
1
2
3
4
5
Error
10
10
10
10
10
2
; 10
2
; 10
2
; 10
2
; 10
2
; 10
3
;10
4
;10
3
;10
4
;10
8
;10
4
; 10 5; : : :
6
; 10 8; : : :
5
; 10 8; : : :
8
; 10 16; : : :
24
;:::
Convergencia
Lineal, con ˇ D 10 1
Lineal, con ˇ D 10 2
Superlineal: entre lineal y cuadrática
Cuadrática
Cúbica
La convergencia de una sucesión de vectores necesita de una función que los
convierta en un número.
Lo usual es usar algún tipo de norma.
Método de Newton-Raphson
De acuerdo con todo lo expuesto hasta ahora:
Sería deseable disponer de una vía sistemática y fiable de construir un modelo
x D g.x/ para hallar la solución xN de la ecuación f .x/ D 0 comenzando
desde cualquier x0, próximo a la solución, obviando que jg 0.x/j
N < 1.
Isaac Newton, Inglaterra, 1642-1727, la persona con una de las mentes más
portentosas que ha dado la humanidad, fue el primero que ideó esa vía y la
forma de llevarla a la práctica sencillamente.
32/75
33/75
Su idea consiste en reemplazar la función f .x/ en cada punto del proceso
iterativo por el modelo de ella que define su recta tangente en ese punto
(linealizar la función en un punto).
En x D x1 la ecuación de la recta tangente a f .x/ es
y D f .x1/ C f 0.x1/.x
x1/.
En x D x1, y D f .x1/ por lo que la ordenada de esta ecuación es la misma
que la de f .
La pendiente de f en x1 es la misma que la de y: f 0.x1/.
El siguiente punto del proceso iterativo lo determina la solución de
y.x/ D 0, es decir, dónde esa recta tangente corta al eje x:
0 D f .x1/ C f 0.x1/.x
x1/.
La solución de esta última ecuación es
x D x1
34/75
f .x1/
:
f 0.x1/
En la figura se describe este paso del proceso iterativo de Newton.
x2
x1
x
f (x)
La relación general de recurrencia del método de Newton, o Newton-Raphson
Joseph Raphson, Inglaterra, 1648-1715
para la solución de la ecuación f .x/ D 0 es
xkC1 D xk
f .xk /
.
f 0.xk /
35/75
Si aplicásemos el método de Newton-Raphson al problema x 3
relación de recurrencia sería:
xkC1 D xk
sen.x/ D 0, la
xk3 sen.xk /
:
3xk2 cos.xk /
Este es un código de Matlab para resolver x 3
Newton.
function raiz=Newton(fun,x,tol)
%Newton para una variable
x0=0; if nargin<3, tol=eps^0.5; end
while abs(x-x0)>tol
x0=x;
[f df] = fun(x0);
x=x0-f/df;
fprintf(’ %18.15f\n’, x);
end
raiz=x;
end
sen.x/ D 0 por el método de
function [f df]=Newt_1(x)
f = x^3-sin(x);
if nargout<2, return, end
df = 3*x*x-cos(x);
end
Si se parte de x0 D 1;4, con
36/75
>> Newton(@Newt_1,1.4)
los puntos que se obtienen con ese código son estos
xk
1,092024491974
0,958975041400
0,929997813651
0,928629313033
0,928626308746
0,928626308732
El gráfico del proceso que lleva a la solución es este:
Método de Newton
2.5
2
1.5
f(x)
k
1
2
3
4
5
6
1
Solución
0.5
0
32
0.5
0.6
0.7
0.8
0.9
1
1
x
1.1
0
1.2
1.3
1.4
1.5
37/75
Convergencia del método de Newton
Teorema Sea la función f W R ! R continúa y derivable, al menos dos veces, y un
r en el que f .r/ D 0. Si f 0 .r/ ¤ 0, el método de Newton es local y cuadráticamente
convergente a r. El error ei en el paso i de su proceso, ei D jxi rj, satisface
ei C1
D M;
i !1 e 2
i
lim
donde M D
f 00 .r/
.
2f 0 .r/
El teorema garantiza la convergencia del método de Newton sólo si se inicia
desde un punto x0 aceptable.
38/75
El método puede no funcionar si jx0
xj
N es grande:
Por ejemplo, considérese el problema clásico de hallar la solución de
arctan.x/ D 0.
Partiendo de cualquier punto del intervalo [1,39, 1,40], el método cicla
obteniéndose x1 D x0, x2 D x0, x3 D x0, : : :
Si x0 < 1;39, el procedimiento converge; si x0 > 1;40, diverge. En la figura
se representan el por qué de esto.
f (x) = arctan(x)
−x0
x0
39/75
El teorema anterior impone que f 0.r/ debe ser distinta de cero para que el
método de Newton converja cuadráticamente a r.
Teorema Si una función f , continua y derivable m C 1 veces en Œa; b, tiene en r una
raíz de multiplicidad m, el método de Newton converge linealmente a r y el error en el paso
i , ei D jxi rj satisface
ei C1
m 1
lim
D
:
i !1 ei
m
40/75
Ejemplo Raíz múltiple. Apliquemos Newton-Raphson a f1.x/ D x 2
y f2.x/ D x 2 2x C 1 D 0, partiendo de x0 D 2:
f1.x/ D x 2
2
1,25
1,025
1,0003048780488
1,0000000464611
1,0
El error es
m D 2.
e5
e4
D
1;03125 1
1;0625 1
f2.x/ D x 2
1D0
D 0;5 D
x0
x1
x2
x3
x4
x5
m 1
,
m
1D0
2x C 1 D 0
2
1,5
1,25
1,125
1,0625
1,03125
lo que verifica que la multiplicidad es
Teorema Si una función continua y derivable m C 1 veces en Œa; b tiene en r una raíz
de multiplicidad m > 1, el método de Newton modificado
xi C1 D xi
mf .xi /
f 0 .xi /
converge cuadráticamente a r.
Ejemplo Estudiar f .x/ D sen.x/ C x 2 cos.x/ x 2 x D 0. Partiendo de
x0 D 1, modificar el código de antes para admitir raíces múltiples y resolverlo.
f .x/
f 0 .x/
f 00 .x/
f 000 .x/
D sen.x/ C x 2 cos.x/ x 2 x
D cos.x/ C 2x cos.x/ x 2 sen.x/ 2x 1
D
sen.x/ C 2 cos.x/ 4x sen.x/ x 2 cos.x/
D
cos.x/ 6 sen.x/ 6x cos.x/ C x 2 sen.x/
cada derivada se evalúa en 0. Como f 000.0/ D
2
1, ésta es una raíz triple.
Comparar
function raiz_Newton(fun,x,tol)
function raiz=Newton_multiple_modificado(fun,x,tol,m)
41/75
42/75
Variantes del método de Newton-Raphson
La primeras resultan de incorporar algún mecanismo que impida que ocurran
alguno de los problemas apuntados.
Recordemos que la resolución de la ecuación de Newton no sólo define un nuevo
punto del proceso iterativo, xkC1, sino una dirección, f 0.xk /, a lo largo de la
cual se da un paso igual a xkC1 xk :
Puede que ese paso sea bueno y la función f .x/ en el nuevo punto adquiera
un valor menor que el que tenía en xk ; pero también puede que mayor y el
proceso diverja, siendo en cualquier caso buena la dirección calculada.
43/75
Intuitivamente, se puede incorporar un mecanismo de salvaguarda que permita
que, a lo largo de la dirección calculada, siempre disminuya el valor de la
función; en concreto,
Si el paso completo xkC1
jf .xkC1/j < jf .xk /j.
xk produce un aumento, disminuirlo hasta que
Ese posible mecanismo lo plasma el algoritmo que sigue.
f .xk /
f 0 .xk /
while .jf .xkC1 /j jf .xk /j/ do
xkC1 C xk
xkC1
2
xkC1 D xk
end
44/75
xk+1
(xk+1 + xk )/2
xk
xk+1
En la figura se ilustra un caso de cómo el mecanismo apuntado salva las
dificultades que surgirían de aplicar el procedimiento de Newton sin él. El punto
0
xkC1
, que sería el que determinaría el paso de Newton, no valdría. Tampoco
0
.xkC1
C xk /=2. Sí, por fin,
xkC1 D
xk C
0
xkC1
Cxk
2
2
:
Método de Newton por diferencias finitas
45/75
Hasta ahora hemos supuesto que se conoce la expresión de la derivada de la
función f .x/. No siempre es así.
Bien porque su determinación analítica es muy complicada –la función f .x/
surge de un procedimiento experimental, por ejemplo–, o,
Porque el usuario del método no desea obtenerla.
En el método de Newton por diferencias finitas se sustituye la derivada de la
función f .x/ por su definición,
f .xk C h/
h!0
h
f 0 .xk / D lim
f .xk /
;
utilizando la fórmula de recurrencia
xkC1 D xk
f .xk /
.f .xk C h/ f .xk //= h
46/75
Teorema Sea la función f W D ! R con dominio de definición en un intervalo abierto
D y derivada continua en él. Supóngase que, para todo x 2 D, jf 0 .x/j para algún
> 0. Si f .x/ D 0 tiene solución xN 2 D, existen unas constantes positivas y 0 tales
N < , la
que si fhk g es una sucesión de números reales tales que 0 < jhk j 0 y si jx0 xj
sucesión fxk g que define
xkC1 D xk
f .xk /
;
ak
con ak D
f .xk C hk /
hk
f .xk /
; k D 0; 1; : : :
converge linealmente a x.
N Si limk!1 hk D 0, la convergencia es superlineal. Si existe alguna
constante c1 tal que jhk j c1 jxk xj,
N o, de forma equivalente, una constante c2 tal que
jhk j c2 jf .xk /j, la convergencia es cuadrática. Si existe alguna constante c3 tal que
jhk j c3 jxk xk 1 j, la convergencia es al menos cuadrática cada dos pasos.
47/75
La elección de h es crítica para el buen funcionamiento del procedimiento: no
debe ser muy pequeño, de tal manera que
f l.xk C h/ ¤ f l.xk /
ni que, dado que f es continua y su derivada también, al evaluar la función en
dos puntos muy próximos, ocurra que
f l.f .xk C h// ¤ f l.f .xk //.
Una regla habitual es la de elegir
p
jhj D máx.ftip x; jxk jg;
donde tip x indica la magnitud típica de x y es la precisión de la máquina en
la que se utiliza el correspondiente código.
48/75
Para un problema bien escalado bastaría hacer
jhj D
p
maq
K
Cuando existen problemas de precisión, también es posible recurrir a la
aproximación de f 0.xk / dada por
ak D
f .xk C h/
f .xk
2h
h/
con h D
p
3
maq:
K :
¡OJO! El número de veces que se evalúa la función se duplica.
49/75
Ejemplo Calculemos la solución de f .x/ D x 2 1, partiendo de x D 2,
mediante Newton y Newton por diferencias finitas.
function newt_x2_1
% Newton de x^2-1=0
x=2.0; x0=0.0;
fx=@(x)x^2-1;
derfx=@(x)2*x;
while abs(x-x0)>eps
x0=x;
x=x0-fx(x0)/derfx(x0);
fprintf(’ %18.15f\n’,x)
end
end
function newt_x2_2
% Newton dif. finitas de x^2-1=0
x=2.0; x0=0.0; h=sqrt(eps);
fx=@(x)x^2-1;
while abs(x-x0)>eps
x0=x;
x=x0-fx(x0)/((fx(x0+h)-fx(x0))/h);
fprintf(’ %18.15f\n’,x)
end
end
50/75
Los resultados obtenidos con uno y otro código son los de las siguiente tabla.
Newton
1,250000000000000
1,025000000000000
1,000304878048780
1,000000046461147
1,000000000000001
1,000000000000000
Newton Dif. Fin.
x0
x1
x2
x3
x4
x5
1,250000000000000
1,025000001341105
1,000304878371890
1,000000046463329
1,000000000000001
1,000000000000000
Como se puede observar, son prácticamente los mismos.
En la práctica, si se acondicionan los parámetros convenientemente, el método
de Newton y el de Newton por diferencias finitas funcionan casi igual.
51/75
Método de Newton-Raphson relajado
Esta variante utiliza como dirección de búsqueda la misma en cada iteración:
f 0.x0/.
f (x)
x3 x2 x1
x0
x
Si la pendiente de f en x0 difiere significativamente de la de f en la solución,
la convergencia puede ser muy lenta o no existir.
La derivada de la función se puede evaluar cada n iteraciones.
52/75
Si se utiliza esta variante para resolver x 3 sen.x/ D 0, partiendo de x0 D 1;1,
los puntos del proceso que se obtienen son los de la tabla.
k
xk
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
1,100000000000000
1,000000000000000
0,950091669305211
0,936190919770347
0,931398691322275
0,929655744326761
0,929010358902363
0,928769834571150
0,928679981323918
0,928646384848608
0,928633818836375
0,928629118217882
0,928627359757290
0,928626701921047
0,928626455824360
0,928626363759310
La convergencia es bastante peor que la del método de Newton.
53/75
El código de Matlab que se ha utilizado es este.
x1=1.1; tol=sqrt(eps);
dx = 3*x1*x1-cos(x1); x2=1;
%
fx = @(x)x^3-sin(x);
% función
%
while abs(x2-x1)>tol
fprintf(’ %18.15f\n’, x1);
x1=x2;
x2=x1-fx(x1)/dx;
end
54/75
Una variante interesante es la conocida como composite Newton, que evalúa la
derivada cada dos iteraciones.
function raiz=Newton_composite(fun,x,tol)
% Newton composite: misma derivada cada dos iteraciones
in=1; x1=0; x2=2;
if nargin<3, tol=eps^0.5; end
if nargin==2, x1=0; x2=x; end
while abs(x2-x1)>tol
x1=x2;
if in
[f df] = fun(x1); in = 0;
else
[f] = fun(x1); in = 1;
end
x2=x1-f/df;
fprintf(’%18.15f\n’, x2);
end
raiz=x2;
end
55/75
Si se resuelve x 3 sen.x/ D 0, partiendo de x0 D 1;4, los puntos del proceso
que se obtienen con
>> raiz=Newton_composite(@Newt_1,1.4)
son los que siguen.
>> raiz=Newton_composite(@Newt_1,1.4)
1.092024491973809
1.019398480362132
0.939380250610315
0.930994424066045
0.928635243248317
0.928626375994235
0.928626308731742
0.928626308731734
raiz =
0.928626308731734
56/75
Método de Halley
Método debido a Edmund Halley, Inglaterra, 1656-1742.
Su idea consiste en reemplazar la función f .x/ en cada punto del proceso
iterativo por su aproximación por Taylor hasta segunda derivada en ese punto.
En x D x1 la aproximación hasta segunda derivada de f .x/ es
0
y D f .x1 / C f .x1 /.x
00
x1 / C f .x1 /
.x
x1 /2
.
2
El siguiente punto del proceso iterativo lo determina la solución de
y.x/ D 0, es decir, dónde esa función corta al eje x:
0
0 D f .x1 / C f .x1 /.x
00
x1 / C f .x1 /
.x
x1 /2
.
2
57/75
La solución de esta última ecuación es
x D x1 C
f 0 .x1 / ˙
Œf 0 .x1 /2 2f .x1 /f 00 .x1 /
:
f 00 .x1 /
La relación general de recurrencia del método es esta:
xkC1 D xk
p
2f .xk /f 0.xk /
2Œf 0.xk /2
f .xk /f 00.xk /
Una variante muy similar es la de Chebyshev,
xkC1 D xk
f .xk /
f .xk /f 00 .xk /
1C
f 0 .xk /
2Œf 0 .xk /2
Pafnuty Lvovich Chebyshev,
Rusia, 1821-1894.
Una codificación de este método para resolver x C ln.x/ D 0 es la que sigue.
function HalleysMethod_log
%Halleys Method; caso particular f(x) = x + ln(x)
i = 1; p0 = 1; N = 100;
% initial conditions
error = 0.000000001;
% precision required
syms ’x’
% Symbolic variables, functions, and expresions in Matlab
f(x) = x + log(x);
% The function we are root finding.
dx = diff(f);
% first and second derivatives of f(x)
ddx = diff(dx);
while i <= N
p = p0 - (f(p0)/dx(p0))*(1 - (f(p0)*ddx(p0)/dx(p0)^2))^(-1);
%Implementation of Halleys Method.
if (abs(p - p0)/abs(p)) < error
%stopping criterion when difference between iterations is below tolerance
fprintf(’Solution is %15.13f; %3.0f iterations.\n’, double(p),i)
return
end
i = i + 1;
p0 = p;
% update p0
fprintf(’%15.13f; iteration %3.0f\n’, double(p0),i)
end
fprintf(’Solution did not coverage within %d iterations; required precision of %d \n’, N, error)
end
>> tic, HalleysMethod_log, toc
0.6000000000000; iteration
2
0.5676852524497; iteration
3
0.5671434553266; iteration
4
0.5671432904098; iteration
5
Solution is 0.5671432904098;
5 iterations.
Elapsed time is 0.664269 seconds.
La ejecución:
La convergencia del método de Halley es cúbica.
58/75
59/75
Índice
El problema
Método de la bisección
Iteración de punto fijo
Velocidad de convergencia
Método de Newton
Variantes del método de Newton
Métodos sin derivadas
60/75
Métodos sin derivadas
Método de la secante
Este método utiliza como dirección de búsqueda, en vez de la que marca la
tangente en un punto, la que determina una recta secante a la función en dos
puntos sucesivos del proceso iterativo.
Si en una iteración k del proceso la ecuación de Newton es
xkC1 D xk
f .xk /
;
f 0.xk /
la idea es emplear, en vez de f 0.xk /,
f .xk /
xk
f .xk 1/
:
xk 1
La relación de recurrencia del proceso iterativo queda
xkC1 D xk
xk 1
f .xk /.
f .xk 1/
La figura ilustra esta aproximación.
f (x)
xk
f .xk /
61/75
xk+1
xk
xk−1
x
El método p
de la secante converge superlinealmente a la solución con orden
˛ D .1 C 5/=2ˇ D 1; ˇ618 –la denominada razón áurea–. De hecho, si r es la
ˇ f 00.r/ ˇ˛ 1 ˛
solución, ei C1 ˇ 2f 0.r/ ˇ ei .
Ejemplo
62/75
Resolvamos con este método y Matlab la ecuación x 3
function raiz=Newton_sec(fun,x,tol)
%Newton secante para una variable
x1=x; x2=x-0.1; if nargin<3, tol=eps^0.5; end
while abs(x2-x1)>tol
x0=x1; x1=x2;
[f df] = fun(x0,x1);
x2=x1-f/df;
fprintf(’ %18.15f\n’, x2);
end
raiz=x2;
end
sen.x/ D 0.
function [f df]=Newt_sec_1(x0,x1)
fx = @(x)x^3-sin(x);
f = fx(x1);
if nargout<2, return, end
df = (fx(x1)-fx(x0))/(x1-x0);
end
La convergencia, con x0 D 1;4 y x1 D 1;4, y la instrucción de Matlab
>> Newton_sec(@Newt_sec_1,1.4) es esta:
k
xk
1
2
3
4
5
6
7
1,065107226548895
0,978070204938512
0,937387385189276
0,929273092828631
0,928635284046901
0,928626318027714
0,928626308731868
63/75
Método de la falsa posición
Conocido como Regula Falsi. Utiliza como dirección de búsqueda una recta
secante a la función en dos puntos sucesivos del proceso iterativo donde la
función toma valores de signo opuesto.
x4
x3
x2
f (x)
x
x1
La velocidad de convergencia es superlineal de orden 1,618: la razón áurea de
nuevo.
64/75
En determinadas circunstancias desfavorables, tanto el método de la secante
como el de la falsa posición pueden presentar problemas de convergencia.
En la figura se representa un ejemplo de convergencia “lenta”.
x2
f (x)
x
x1 x3
Método de Muller
65/75
Presentado en 1956.
David Eugene Muller, EE.UU., 1924-2008.
Utiliza una interpolación cuadrática de tres puntos anteriores del proceso
iterativo para extrapolar el nuevo punto del proceso.
f (x)
p(x)
|
x0
|
x1
|
x2
|
x3
El procedimiento aproxima a x0, x1 y x2 el polinomio
x2/2 C b.x
p.x/ D a.x
66/75
x2 / C c
que pasa por .x0; f .x0//, .x1; f .x1// y .x2; f .x2//.
Los parámetros a, b y c se determinan mediante estas condiciones:
x2 /2 C b.x0
x2 /2 C b.x1
f .x0 / D a.x0
f .x1 / D a.x1
f .x2 / D c:
x2 / C c;
x2 / C c
y
Resolviendo el sistema de 2 ecuaciones con 2 incógnitas
h.x x /2 .x x /i h i hf .x / f .x /i
0
2
0
2
0
2
a
D
b
.x
x /2
.x
x /
f .x / f .x /
1
2
1
2
1
2
se obtiene la expresión de los parámetros a y b:
a
D
b
1
.x0 x2 /.x1 x2 /.x0 x1 /
.x1
.x1
x2 /
x2 /2
.x0 x2 / f .x0 /
.x0 x2 /2 f .x1 /
f .x2 /
:
f .x2 /
67/75
Para determinar el nuevo punto del proceso se aplica la fórmula
zD
2c
;
p
2
b˙ b
4ac
escogiéndose el signo, de los dos posibles, que garantice el denominador más
grande en valor absoluto:
Si b > 0, el signo positivo;
si b < 0, el negativo.
El nuevo punto x3 será entonces
x3 D x2 C z:
68/75
Este código de Matlab implementa el método.
function Muller_2(fx,x0,x1,x2) % Método de Muller
tol=eps(’double’);
fx0=fx(x0); fx1=fx(x1); fx2=fx(x2); iter=1;
while abs(fx2)>tol
c=fx2; d0=x0-x2; d1=x1-x2; det=d0*d1*(x0-x1);
b=(d0*d0*(fx1-fx2)-d1*d1*(fx0-fx2))/det;
a=(d1*(fx0-fx2)-d0*(fx1-fx2))/det;
di=sqrt(b*b-4*a*c);
isig=1; if b<0 isig=-1; end
z=(-2)*c/(b+isig*di);
x3=x2+z;
if abs(x3-x1)<abs(x3-x0)
u=x1; x1=x0; x0=u;
u=fx1; fx1=fx0; fx0=u;
end
if abs(x3-x2)<abs(x3-x1)
u=x2; x1=u;
u=fx2; fx1=u;
end
x2=x3; fx2=fx(x2);
if ~isreal(x2)
fprintf(’%17.14f+%17.14fi %23.15f %4.0f\n’,real(x2),imag(x2),fx2,iter);
else
fprintf(’%17.15f %23.15e %4.0f\n’,x2,fx2,iter);
end
iter=iter+1;
end
end
69/75
El proceso de convergencia de la resolución de la ecuación x 3
partiendo de x0 D 1;5, x1 D 1;2 y x2 D 1;0 con
sen.x/ D 0,
>> Muller_2(@Newt_1,1,1.2,1.5)
es el que expresa esta tabla.
k
xk
f .xk /
1
2
3
4
5
0,921801501077277
0,928699331903730
0,928626328365127
0,928626308731740
0,928626308731734
-1,342037938307839e-02
1,451947876007775e-04
3,903326339926849e-08
1,076916333886402e-14
-1,110223024625157e-16
Obsérvese que el número de iteraciones, para la precisión que se obtiene,
decrece apreciablemente comparándolo con el de otros métodos.
70/75
Si se trata de calcular una de las raíces complejas de x 3 C x 2 C x C 1 con este
método resultaría lo que sigue.
>> Muller_2(@Muller_2_2,0,0.1,1)
-0.21428571428571+0.65595130066457i 0.668124392614188
-0.16304914264115+1.15401500391091i 0.178872980763637
0.03069686090353+0.97963558152672i-0.016395820078345
0.00106074539714+1.00049046815007i-0.003104663615609
-0.00000339553880+1.00000042898932i 0.000005933119030
-0.00000000004554+0.99999999997470i 0.000000000141669
-0.00000000000000+1.00000000000000i 0.000000000000000
>>c=[1 1 1 1];
>> roots(c)
ans =
-1.0000
-0.0000 + 1.0000i
-0.0000 - 1.0000i
>>
1
2
3
4
5
6
7
El método de Brent
71/75
En su versión más reciente fue formulado en 1973 por Richard Peirce Brent,
Australia, 1946-. Anteriormente, en la década de 1960, fue propuesto por
Dekker y Van Wijngaarden.
Richard Peirce Brent, Australia, 1946-.
Utiliza simultáneamente el métodos de la bisección y el de Muller (o secante). Se
aplica a un intervalo Œa; b en cuyos extremos la función adopta signos distintos.
Sigue la pista a un punto xi que es el mejor en el sentido del error hacia atrás y
a un intervalo Œai ; bi  para la raíz.
72/75
Se aplica una interpolación cuadrática de Lagrange inversa, es decir, no
y D p.x/ sino x D p.y/, a los tres puntos .f .xi /; xi /, .f .ai /; ai / y
.f .bi /; bi / con el fin de reemplazar uno de ellos con aquel –único– donde
x D p.y D 0/.
1.5 Root-Finding without Derivatives | 65
Quadratic Interpolation method is attempted, and the result is used to replace one of xi , ai , bi
if
(1) these
backward
error improves
and (2) the del
bracketing
intervalde
is cut
at least incon
half. la
If not,
En esta figura
compara
la geometría
método
Muller
de
the Secant Method is attempted with the same goal. If it fails as well, a Bisection Method
interpolación
inversa.
step iscuadrática
taken, guaranteeing
that the uncertainty is cut at least in half.
y
xIQI
x0
x2 x
M
x1
x
Figure 1.13 Comparison of Muller’s Method step with Inverse Quadratic Iteration
step. The former is determined by an interpolating parabola y = p(x); the latter, by an
interpolating parabola x = p(y).
la
73/75
Si el error hacia atrás mejora y el intervalo de confinamiento de la solución se
reduce al menos a la mitad, el punto obtenido reemplaza a uno de los tres
vigentes.
Si no se cumplen los requisitos anteriores, se intenta el método de la secante
con el mismo objetivo. Si también falla, se realiza un paso del método de la
bisección.
74/75
Este método lo usa fzero de Matlab. Una implementación interesante es esta.
function value = Brent (a,b,machep,t,f)
% Se calcula un c en [a,b] tal que f(c)=0.
%
%
En [a,b] debe darse un cambio de signo en f(x).
%
La precisión de la raíz es 6*machep*abs(c)+2*t.
%
sa = a; sb = b; fa = f(sa); fb = f(sb);
c = sa; fc = fa; e = sb-sa; d = e;
while 1
if abs(fc)<abs(fb), sa=sb; sb=c; c=sa; fa=fb; fb=fc; fc=fa; end
tol = 2.0*machep*abs(sb) + t;
m = 0.5*(c-sb);
if abs(m)<=tol || fb==0.0, break, end
if abs(e)<tol || abs(fa)<=abs(fb)
e = m; d = e;
% Bisección, decrece poco
else
s = fb/fa;
if sa==c, p=2.0*m*s; q=1.0-s;
% Interpolación lineal
else
q = fa/fc; r = fb/fc;
% Interpolación cuadrática
p = s*(2.0*m*q*(q-r)-(sb-sa)*(r-1.0)); %
inversa
q = (q-1.0)*(r-1.0)*(s-1.0);
end
if 0.0<p, q=-q;
else
p=-p;
end
s = e; e = d;
if 2.0*p<3.0*m*q-abs(tol*q) && p<abs(0.5*s*q)
d = p/q;
% Interpolación aceptada
else
e = m; d = e;
% Bisección
end
end
sa = sb; fa = fb;
% Completa paso
if tol<abs(d), sb=sb+d;
% Nuevo punto
elseif 0.0<m
sb = sb + tol;
else
sb = sb - tol;
end
fb = f(sb);
if (0.0<fb&&0.0<fc) || (fb<=0.0&&fc<=0.0)
c = sa; fc = fa; e = sb-sa; d = e;
end
end
value = sb;
end
75/75
Probemos Brent y fzero con el problema anterior x 3
>> f1=@(x) x^3-2*x^2+x*4/3-8/27;
>> tic, x=Brent(0,1,eps,eps,f1), toc
x =
0.666670143350873
Elapsed time is 0.000297 seconds.
>> tic, fzero(f1,[0 1],optimset(’Display’,’iter’)), toc
Func-count
x
f(x)
Procedure
2
1
0.037037
initial
3
0.888889
0.0109739
interpolation
4
0.843776
0.00555553
interpolation
5
0.798619
0.00229748
interpolation
6
0.798619
0.00229748
bisection
7
0.755766
0.00070733
interpolation
8
0.737384
0.000353658
interpolation
9
0.71944
0.000146977
interpolation
10
0.71944
0.000146977
bisection
11
0.702418
4.56958e-05
interpolation
12
0.695012
2.27741e-05
interpolation
13
0.687828
9.47672e-06
interpolation
14
0.687828
9.47672e-06
bisection
15
0.681016
2.95457e-06
interpolation
16
0.67804
1.47117e-06
interpolation
17
0.675159
6.124e-07
interpolation
18
0.675159
6.124e-07
bisection
19
0.672426
1.91084e-07
interpolation
20
0.671232
9.51209e-08
interpolation
21
0.670075
3.95999e-08
interpolation
22
0.670075
3.95999e-08
bisection
23
0.668979
1.23591e-08
interpolation
24
0.668499
6.15183e-09
interpolation
25
0.668035
2.56116e-09
interpolation
26
0.668035
2.56116e-09
bisection
27
0.667595
7.99392e-10
interpolation
28
0.667402
3.97894e-10
interpolation
29
0.667216
1.65654e-10
interpolation
30
0.667216
1.65654e-10
bisection
2x 2 C 34 x
31
0.667039
5.17053e-11
32
0.666962
2.5736e-11
33
0.666887
1.07147e-11
34
0.666887
1.07147e-11
35
0.666816
3.34444e-12
36
0.666785
1.66456e-12
37
0.666755
6.93112e-13
38
0.666755
6.93112e-13
39
0.666727
2.16271e-13
40
0.666714
1.07692e-13
41
0.666702
4.4964e-14
42
0.666702
4.4964e-14
43
0.666691
1.39888e-14
44
0.666686
6.99441e-15
45
0.666681
2.88658e-15
46
0.666681
2.88658e-15
47
0.666676
9.99201e-16
48
0.666674
3.33067e-16
49
0.666673
3.33067e-16
50
0.666673
3.33067e-16
51
0.66667 -1.11022e-16
52
0.666671
1.11022e-16
53
0.66667
0
Zero found in the interval [0, 1]
ans =
0.666670143350873
Elapsed time is 0.193985 seconds.
8
27
D 0.
interpolation
interpolation
interpolation
bisection
interpolation
interpolation
interpolation
bisection
interpolation
interpolation
interpolation
bisection
interpolation
interpolation
interpolation
bisection
interpolation
interpolation
interpolation
bisection
interpolation
interpolation
bisection