Date post: | 12-Jul-2015 |
Category: |
Education |
Upload: | ariel-belino |
View: | 1,624 times |
Download: | 1 times |
Solución de ecuaciones no lineales
Curso de Programación numérica
TemarioMétodos cerrados:
Métodos gráficos
Método de bisección
Método de la posición falsa
Métodos abiertos
Iteración simple de punto fijo
Método de Newton-Raphson
Método de la secante
Raíces de polinomios
Método de Müller
Método de Bairstow
Métodos gráficos
Los métodos gráficos consisten en graficar la función f(x) y observar donde la función cruza el eje x.
Ejemplo 1( ) ( ) 0401
38.667 146843.0 =−−= − xex
xf
x f(x)
4 34.114889388 17.65345264
12 6.06694996316 -2.26875420820 -8.400624408
Encontrar la raíz de:
-15
-10
-5
0
5
10
15
20
25
30
35
40
0 5 10 15 20 25
Ejemplo 2
x f(x)
0.00 1.000.25 1.330.50 -0.890.75 0.311.00 -1.531.25 -0.891.50 0.441.75 -0.462.00 1.872.25 0.412.50 0.212.75 0.313.00 -1.903.25 -0.063.50 -0.903.75 0.054.00 1.594.25 -0.014.50 1.454.75 -0.485.00 -1.02
-2.50
-2.00
-1.50
-1.00
-0.50
0.00
0.50
1.00
1.50
2.00
2.50
0.00 1.00 2.00 3.00 4.00 5.00 6.00
Grafica de: f(x) = sen 10x + cos 3x
Ejemplo 2 (cont.)
-0.02
0.00
0.02
0.04
0.06
0.08
0.10
0.12
4.18 4.20 4.22 4.24 4.26 4.28 4.30 4.32
Grafica de: f(x) = sen 10x + cos 3x
x f(x)
4.20 0.084.21 0.054.22 0.024.23 0.004.24 -0.014.25 -0.014.26 -0.014.27 0.014.28 0.044.29 0.074.30 0.11
TareaUtilice Excel para los siguientes problemas.
Determine las raíces reales de: f(x) = –0.5x2 + 2.5x + 4.5
Gráficamente. Confirme utilizando la fórmula cuadrática.
Determine las raíces reales de: f(x) = 5x3 – 5x2 + 6x – 2
Gráficamente.
Método de la bisecciónSe trata de encontrar los ceros de
f(x) = 0
Donde f es una función continua en [a,b] con f(a) y f(b) con signos diferentes.
y = f(x)
x
y
a
b
f(b)
f(a)
Método de la bisección
De acuerdo con el teorema del valor medio, existe p ∈ [a,b] tal que f(p) = 0.
El método consiste en dividir a la mitad el intervalo y localizar la mitad que contiene a p.
El procesos se repite hasta la lograr la precisión deseada.
Método de la bisección
y = f(x)
x
y
a
b
f(b)
f(a)
p1=(a+b)/2
f(p1)
p
Mitad del intervalo que contiene a p
Primera iteración del algoritmo
Método de la bisección
y = f(x)
x
y
a =p1
b
f(b)
f(a)
p2=(a+b)/2
f(p2)p
Mitad del intervalo que contiene a p
Segunda iteración del algoritmo
Método de la bisección
Algoritmo bisección
Entradas: extremos a,b; número de iteraciones ni; tolerancia tol
1. p=a; i=1; eps=1;2. mientras f(p)≠0 y i≥ ni eps>tol 2.1. pa = p; 2.2. p = (a+b)/2 2.3. si f(p)*f(a)>0 entonces a=p; 2.4. sino 2.5. si f(p)*f(b)>0 entonces b=p; 2.6. i = i + 1; eps = |p-pa|/p;
Bisección en Cdouble biseccion(double a, double b, double error, int ni){ double p,pa,eps; int i; p = a; i = 1; eps = 1; while(f(p) != 0 && i<ni && eps > error){ pa = p; p = (a+b)/2; if(f(p)*f(a)>0) a = p; else if(f(p)*f(b)>0) b = p; i++; eps = fabs(p-pa)/p; } return p;}
Ejemplo
Función de ejemplo
Función en C:
double f(x){
return sqrt(x*x + 1) - tan(x);
}
)tan(12 xx −+
Tarea
Haga funciones en C para encontrar la solución de las siguientes ecuaciones utilizando la función biseccion():
1. ex – x2 + 3x – 2 = 0 para 0 <= x <= 1
2. ( ) ( ) 040138.667 146843.0 =−−= − xe
xxf
Error en el método de bisecciónPara el método de bisección se sabe que la raíz esta dentro del intervalo, la raíz debe situarse dentro de ±∆x / 2, donde ∆x = xb – xa.
La solución en este caso es igual al punto medio del intervalo
xr = (xb + xa) / 2
Deberá expresarse por
xr = (xb + xa) / 2 ± ∆x / 2
Error aproximado
2abanterior
rnuevor
xxxx
−=− 2abnuevo
r
xxx
+=
%100ab
aba xx
xx
+−
=ε
%100nuevor
anteriorr
nuevor
ax
xx −=ε
sustituyedo
Número de iteracionesEl error absoluto en la primera iteración es:
0000 xxxE aba ∆=−=
El error absoluto en la iteración n-ésima es:
nna
xE
2
0∆=
Si el error deseado es Ead,El número de iteraciones será:
( )
∆=∆=ad
ad
E
xExn
0
2
0
log2log
/log
Volumen del abrevadero
h
r
L
rh
αβ
=
rh
senβ
α2sectorarea r=
−=−= −
rh
sen 1
22πβπα
( )
−== − rhsenrr /
2sectorarea 122 πα
( ) 2212 /2
triangularareasectorareaA hrhrhsenr −−
−=−= −π
22
2alturabase
2triangulararea hrh −=×=
( )
−−
−== − 2212 /
2hrhrhsenrLLAV
π
Tarea17. Un abrevadero de longitud L tiene una sección transversal en forma de semicírculo con radio r (véase la figura) Cuando se llena de agua hasta una distancia h de la parte superior, el volumen V de agua es
V = L [ 0.5 πr2 – r2 arcsen(h/r) – h(r2 – h2)1/2 ]
Escriba un programa en C amigable para el usuario que lea los datos de este problema y encuentre la profundidad h del abrevadero. Utilice el método de bisección para encontrar la solución.
h
r
L
Resumen
Requiere que se conozca el intervalo en donde está la raíz.
Los valores de la función en los extremos deben tener signos diferentes.
Converge lentamente, a cada paso el intervalo se divide en 2.
Método de falsa posición
xl
xu
f(xl)
f(xu)
xr
Este método considera cual límite del intervalo está más próximo a la raíz.
De la figura
( ) ( )ur
u
lr
l
xx
xf
xx
xf
−=
−
( ) ( )( ) ( )ul
uluur xfxf
xxxfxx
−−−=
Despejando
f(xr)
Ejemplo en Excel
xl xu xr f(xl) f(xu) f(xr)12.0000000 16.0000000 14.9113077 6.0669500 -2.2687542 -0.254277512.0000000 14.9113077 14.7941976 6.0669500 -0.2542775 -0.027257212.0000000 14.7941976 14.7817001 6.0669500 -0.0272572 -0.002907612.0000000 14.7817001 14.7803676 6.0669500 -0.0029076 -0.000310012.0000000 14.7803676 14.7802255 6.0669500 -0.0003100 -0.0000330
( ) ( ) 040138.667 146843.0 =−−= − xe
xxfEncontrar la raíz de:
TareaEncuentre la raíz real de f(x) = (0.8 – 0.3x)/x, por el método de falsa posición. Utilice valores iniciales de 1 y 3, calcule el error porcentual verdadero en cada iteración. Encuentre la raíz analiticamente.
Falsa posición en C
double falsaPosicion(double xl, double xu, double ee, int imax){ double error,fl,fu,fr,xr,xrOld; int iter=0,il=0,iu=0; fl = f(xl); fu = f(xu); do{ xrOld = xr; xr = xu - fu*(xl-xu)/(fl-fu); fr = f(xr); iter++; if(xr!= 0) error=fabs((xr-xrOld)/xr*100);
if(fl*fr<0){ xu = xr; fu = f(xu); iu = 0; il++; if(il>=2) fl/=2; } else{ xl = xr; fl = f(xl); il = 0; iu++; if(iu>=2) fu/=2; else; error = 0; } }while(error>ee && iter<=imax); return xr;}
Iteración de punto fijo
Un punto fijo de una función g(x) es un número p tal que g(p) = p.
Dado un problema f(x) = 0, se puede definir una función g(x) con un punto fijo en p de diferentes maneras. Por ejemplo g(x) = x – f(x).
TeoremaSi g ∈ C [a, b] y g(x) ∈ C [a, b] para toda x ∈ C [a, b], entonces g tiene un punto fijo en [a, b].
Si además g’(x) existe en (a, b) y una constante positiva k<1 existe con
|g’(x)| <= k, pata toda x ∈ (a, b),
Entonces el punto fijo en [a, b] es único.
y = g(x)
y
a b
a
b
p
p=g(p)
x
y = x
y = g(x)
y
a b
a
b
p
p=g(p)
x
y = x|g’(x)>1|g’(x)<=1
Algoritmo de punto fijo
Obtener una solución a p = g(p) dada un aproxiamción inicial p0.
ENTRADA aproximación inicial p0; tolerancia TOL; número máximo de iteraciones N0.
1. Tome i = 1.2. Mientras i <= N0 hacer 3. p = g(p0) 4. Si |p – p0| < TOL entonces5. Regresar p6. i = i +17. p0 = p8. Fin mientras9. Imprime ‘El procedimiento fracasó después de N0 iteraciones’
Gráfica del algoritmo de punto fijo
y = g(x)
y
x
y = x
p0
p1= g(p0)
p3 p2p1
p2= g(p1)
p3= g(p2)
y = g(x)
y
x
y = x
p0
p1= g(p0)
p2p1
p2= g(p1)
p3= g(p2)
Casos de no convergencia
y = g(x)
y
x
y = x
y = g(x)
y
x
y = x
EjemploSea la función: x3 + 4x2 –10 = 0 tiene una raíz en [1, 2]
Puede despejarse en:
a. x = g1(x) = x – x3 – 4x2 +10
b. x = g2(x) = ½(10 – x3)½
c. x = g3(x) = (10/(4 + x))½
d. x = g4(x) = x – (x3 + 4x2 – 10)/(3x2 + 8x)
Iteraciones de punto fijo
(b)
1.51.2869537671.4025408031.3454583741.3751702521.3600941921.3678469671.3638870031.3659167331.3648782171.3654100611.3651378201.3652772081.3652058501.3652423831.3652295781.3652300281.365230012
(c)
1.51.3483997241.3673763711.3649570151.3652647481.3652255941.3652305751.3652299411.3652300221.3652300121.3652300131.365230013
(a)
1 1.52 -0.8753 6.7324218754 -469.720012005 1.02754555E86 -1.084933870E247 1.277055591E728 -2.082712908E2169 NaN101112131415202530
(d)
1.51.3733333331.365262014 1.3652300131.365230013
Funciones graficadas en MathLab
a) b)
c) d)
Teorema de punto fijo
Si g ∈ C [a, b] y g(x) ∈ C [a, b] para toda x ∈ C [a, b], además supongamos que existe g’(x) en (a, b) y una constante positiva k<1 cuando
|g’(x)| <= k, pata toda x ∈ (a, b),
Entonces, para cualquier punto p0 en [a, b] la sucesión definida por
pn = g(pn–1), n >=1
Converge en el único punto fijo p en [a, b].
Corolario
Si g satisface las hipótesis de teorema del punto fijo, las cotas de error que supone utilizar pn para aproximar a p están dadas por
| pn – p| <= kn max(p0 – a, b – p0)
Y por
| pn – p| <= kn | p1 – p0|/ (1 – k), para toda n>=1
Análisis del ejemploCaso (a)
g1(x) = x – x3 – 4x2 +10
g1’(x) = 1 – 3x2 – 8x
g1’(1) = – 11, g1’(2) = – 28
No se cumple |g1’(x)| <1
Caso (b)
g2(x) = ½(10 – x3)½
g2’(x) = – 3/4x2(10 – x3)–½
g2’(1) = – 0.25, g1’(2) = – 2.1213
No se cumple |g1’(x)| <1
Caso (c)g3(x) = (10/(4 + x))½ g3’(x) = (– 5/3.16)(4 + x)–1.5
<= (– 5/3.16)(5)–1.5 <= 0.15
Para toda x en [1, 2]
Caso (d)
g4(x) = x – (x3 + 4x2 – 10)/(3x2 + 8x)Se cumple |g4’(x)| es aún menor que en el caso (c) para toda x en [1, 2]
Programa en Matlabfunction y = PuntoFijo(f_name, p0, tol, ni)%f_name - nombre de la funcion%p0 - valor inicial de la raiz%tol – tolerancia%ni – número de iteracionesi = 1;while i<=ni p = feval(f_name,p0); if(abs(p0-p)<tol) y = p; break; end i = i + 1; p0 = p;endfprintf('No se encontro solucion.');
Función en Cdouble PuntoFijo(double p0, double tol, int ni){ int i = 1; double p; while(i<=ni){ p = f(p0); if(fabs((p0-p)/p)<tol) return p; i++; p0 = p; }
std::cout << "NO solucion en :" << ni << “ iteraciones.\n"; return p;}
Método de Newton-Raphson
f (xn)Pendiente = f ’ (xn)
xnxn+1
La ecuación de la recta tangente es:
y – f(xn) = f ’ (xn)(x – xn)
Cuando y = 0, x = xn+1 o sea
0 – f(xn) = f ’ (xn)(xn+1– xn)
o
x xf x
f xn nn
n+ = −1
( )
'( )
f(x)
Algoritmo NewtonPara obtener una solución a f(x) = 0 dada una aproximación p0.
ENTRADA aproximación inicial p0; tolerancia tol; número máximo de iteraciones N0.
1. i = 1
2. Mientras i<=N0 hacer
2.1. p = p0 – f(p0)/f’(p0)
2.2. Si |p – p0|< tol entonces regrese p
2.3. i = i + 1
2.4. p0 = p
3. fracaso en encontrar la raíz en N0 iteraciones
Ejemplo
f(x) = x – cos(x) f’(x) = 1 + sen(x)
Tomando p0 = 0, se obtiene
pn f(pn) f’(pn) pn+1
0 -1 1 11 0.459698 1.8414 0.75036390.7503639 0.0189 1.6819 0.73911280.7391128 0.00005 1.6736 0.7390851 0.7390851 3E-10 1.6736 0.7390851
pn+1 = pn – (pn – cos(pn))/(1 + sen(pn))
Ejercicio
Encontrar la solución de
x3 + 4x2 – 10 = 0
En el intervalo [1, 2] con el método de Newton
Código en C
double Newton(double x0, double ee, int ni){ int i = 0; double x,fx,dfx; while(i<ni){ fx = f(x0); dfx = df(x0); x = x0-fx/dfx; if(fabs((x-x0)/x)<ee) return x; i++; x0 = x; } std::cout << "No solución en "<< i << " pasos\n"; return x;}
Ejemplo: cuenta de ahorrosEl valor acumulado de una cuenta de ahorros puede calcularse con la ecuación de anualidad vencida
A = P[(1 + i )n - 1 ] / i
En esta ecuación A es el monto de la cuenta, P es la cantidad que se deposita periódicamente e i es la tasa de interés por periodo para los n periodos de depósito. A un ingeniero le gustaría tener una cuenta de ahorros con un monto de $ 750,000 dólares al momento de retirarse dentro de 20 años, y puede depositar $ 1,500 dólares mensuales para lograr dicho objetivo. ¿Cuál es la mínima tasa de interés a que puede invertirse ese dinero, suponiendo que es un interés compuesto mensual?
Escriba un programa en C para este problema, el programa deberá pedir todos los datos necesarios y utilizar el método de Newton para calcular el interés a que debe invertirse el dinero.
SoluciónPara estimar el valor inicial de i podemos desarrollar el binomio (1 + i)n para aproximarlo a la segunda potencia. El resultado es
( )( ) Pnn
nPAi
1
20 −
−=
Se sugiere validar los datos de entrada. El capital a obtener debe ser mayor que el depósito por el número de abonos, es decir
A > nP
Ejemplos resuelto en ExcelA = $750,000.00 i f(i) f'(i) i n+1P = $1,500.00 0.009065551 4784.893234 2361961.89 0.00703974n = 240 0.007039738 1297.701361 1175049.3 0.00593536
0.005935357 255.8695592 730982.881 0.00558532A(calculado)= $750,000.00 0.005585323 20.97312565 612780.041 0.0055511
0.005551096 0.18919948 601739.117 0.005550780.005550782 1.58807E-05 601638.103 0.005550780.005550782 -1.99179E-10 601638.094 0.00555078
i = 0.56%
A = $350,000.00 i f(i) f'(i) i n+1P = $20,000.00 0.166666667 15099.14998 450849.857 0.13317625n = 10 0.133176249 3212.297411 266179.386 0.12110808
0.121108082 346.4384394 209573.765 0.11945502A(calculado)= $350,000.00 0.11945502 6.113559001 202191.641 0.11942478
0.119424783 0.002029206 202057.423 0.119424770.119424773 2.47383E-10 202057.379 0.119424770.119424773 6.54836E-11 202057.379 0.11942477
i = 11.94%
Método alternativo para evaluar la(derivada (método de la secante
Es posible calcular la derivada en xn usando:
( ) ( ) ( )h
xfhxfxf nn
n
−+='
O utilizando
( ) ( ) ( )h
hxfxfxf nn
n
−−='
Algoritmo Newton2Para obtener una solución a f(x) = 0 dada una aproximación p0.
ENTRADA aproximación inicial p0; tolerancia tol; número máximo de iteraciones N0.
1. i = 12. h = 0.0013. Mientras i<=N0 hacer
2.1. y = f(p0) 2.2. y_deriv =(f(p0+h)-y)/h 2.3. p = p0 – y/y_deriv
2.4. Si |p – p0|< tol entonces regrese p 2.5. i = i + 1 2.6. p0 = p3. fracaso en encontrar la raíz en N0 iteraciones
Código en C
double Newton(double x0, double ee, int ni){ int i = 0; double x,fx,dfx,h; h = 0.0001; while(i<ni){ fx = f(x0); dfx = (f(x0+h)-fx)/h; x = x0-fx/dfx; if(fabs((x-x0)/x)<ee) return x; i++; x0 = x; } std::cout << "No solución en "<< i << " pasos\n"; return x;}
Programa en Matlab
function x = Newt_n(f_name, xO)% Iteración de Newton sin gráficosx = xO; xb = x-999;n=0; del_x = 0.01;while abs(x-xb)>0.000001 n=n+1; xb=x ; if n>300 break; end y=feval(f_name, x) ; y_driv=(feval(f_name, x+del_x) - y)/del_x; x = xb - y/y_driv ; fprintf(' n=%3.0f, x=%12.5e, y=%12.5e, ', n,x,y) fprintf(' yd = %12.5e \n', y_driv)endfprintf('\n Respuesta final = %12.6e\n', x) ;
Calcula derivada con incrementos
Raíz cuadrada con NewtonPara extraer la raíz cuadrada de un número se puede resolver la ecuación
f(x) = x2 – c = 0
La derivada es
f’(x) = 2x
La fórmula de recurrencia de Newton es
xn+1 = xn – (xn2 – c)/(2xn)
= xn/2 + c/(2xn)
= (xn + c/xn)/2
Ejemplo: raíz cuadrada de 5 con x0 = 1.
xn xn+11 3.000000
3.000000 2.3333332.333333 2.2380952.238095 2.2360692.236069 2.236068
DesventajasEn algunos casos la convergencia es muy lenta, considere
f(x) = xn – 1
Se obtiene la siguiente secuencia empezando en x = 0.5
iteración x0 0.51 51.652 46.4853 41.83654 37,65285
..
-- 1.000000
Desventajas (cont.)
x0x1
x2 x
f(x)
x0 x1x2 x
f(x)
x0
x1
x
f(x)
raíz cerca de punto de inflexiónmínimo local
varias raíces
x
f(x)
la iteración en un mínimo
x0 x1
Ejemplo
Resolver utilizando Excel
sen x - e-x = 0 para 0<= x <= 1 y 3<= x <= 4 y 6<= x <= 7
Resultadosh= 0.1
xn f(xn) f'(xn) xn+1
0.00000000 -1.00000000 1.94995999 0.51283104
0.51283104 -0.10815190 1.41522716 0.58925121
0.58925121 0.00099615 1.33011566 0.58850229
0.58850229 -0.00004224 1.33095756 0.58853402
0.58853402 0.00000178 1.33092188 0.58853269
h= 0.01
xn f(xn) f'(xn) xn+10.00000000 -1.00000000 1.99499996 0.501253140.50125314 -0.12524617 1.47731614 0.586032670.58603267 -0.00347081 1.38411969 0.588540270.58854027 0.00001043 1.38133294 0.588532710.58853271 -0.00000004 1.38134134 0.58853274
Tarea #14
La carga en un circuito RLC serie esta dada por
( )
−= − t
LR
LCeqtq LRt
2)2/(
0 21
cos
suponga q0/q = 0.01, t = 0.05 s, L = 5H y C = 10-6 F.
Encuentre el valor de la Resistencia R usando el método de Newton. Haga un programa en C para este problema.
Convergencia en el punto fijo
El algoritmo de punto fijo es de tipo lineal. Se puede demostrar que el error verdadero en la iteración i+1 es:
Et,i+1= g’(ξ)Et,i
donde
Et,i = xr - xi
Convergencia en Newton Raphson
El algoritmo de Newton es de tipo cuadrático. Se puede demostrar que el error verdadero en la iteración i+1 es:
Et,i+1= (- f ’’(xr)/2f ’(xr))E2t,i
Esto significa que el número de decimales exactos se duplica con cada iteración.
Raíces múltiples
En el caso de que un polinomio tenga raíces múltiples, la función tendrá pendiente igual a cero cuando cruce el eje x.
Tales casos no pueden detectarse en el método de bisección si la multiplicidad es par.
En el método de Newton la derivada en la raíz es cero.
Generalmente el valor de la función tiende a cero más rápido que la derivada y puede utilizarse el método de Newton
Ejemplo
n xn f(xn) f'(xn) xn+10 0.50000000 -0.62500000 2.75000000 0.727272731 0.72727273 -0.16904583 1.31404959 0.855917672 0.85591767 -0.04451055 0.63860849 0.925616943 0.92561694 -0.01147723 0.31413077 0.962153414 0.96215341 -0.00291894 0.15568346 0.980902605 0.98090260 -0.00073639 0.07748373 0.990406366 0.99040636 -0.00018496 0.03865069 0.995191767 0.99519176 -0.00004635 0.01930234 0.997593008 0.99759300 -0.00001160 0.00964539 0.998795789 0.99879578 -0.00000290 0.00482125 0.99939771
10 0.99939771 -0.00000073 0.00241026 0.99969881
Polinomio: f(x) = (x – 3) (x – 1) (x – 1)