Práctica 4. Resolución de ecuaciones (II)
Esta práctica tratará el método de la secante, el método de Newton-Raphson y los métodos de iteración funcional.
Contents
Método de la secante
El método de la secante busca aproximaciones a una raíz de una función de forma iterativa. Dadas dos iteraciones, se halla la recta que pasa por los puntos de la gráfica de la función cuyas abscisas son esas dos iteraciones. La siguiente iteración es el punto de corte de dicha recta con el eje OX.
Recordamos el pseudocódigo que permite implementar el método de la secante:
1. Entrar f , a, b, tol 2. Mientras |b - a| >= tol 3. Hacer c =b-f(b)*(b-a)/(f(b)-f(a)) 4. Si f (c) = 0 entonces c es raíz. Fin 5. Hacer a = b y b = c. 6. Ir a 2
% Implementación del método de la secante para una función f definida. a=1; b=2.618; tol=0.0001; while abs(b-a)>=tol c=b-f(b)*(b-a)/(f(b) - f(a)); % fprintf('a= %d, b=%d, c=%d, f(c)=%d\n',a,b,c,f(c)) if f(c)==0 disp('Encontrada raíz.') break end a=b; b=c; end
Podemos observar la salida del algoritmo mostrando c y f(c).
c f(c)
c =
0.4429
ans =
8.8148e-09
Ej 1. Dibuje la función f(x)=x-5*log(x) en el intervalo [0,2]. Encuentre una raíz de f mediante el método de la secante, determinando antes unos valores apropiados para a y b que definan un intervalo de longitud 0.2. (Sol. Por ejemplo con a=1.2, b=1.4 se llega a la aproximación 1.295856)
Ej 2. Aplique el método de la secante para a=1, b=4 y f(x)=(x-2).*(x-7).*(x-2) (Sol. aprox. 1.999989, solución exacta x=2).
Ej 3. ¿Puede aplicarse el método de la secante para a=1, b=3 y f=(x-2).^2 ?
Ej 4. Aplique el método de la secante para a=1, b=(3+sqrt(5))/2 y f(x)=-4*x^3+ 3*(3+sqrt(5))*x^2-3*(1+sqrt(5))*x-3-sqrt(5). (Sol. aprox x=1.309017)
El comando fprintf
Para mostrar en pantalla resultados con texto y valores numéricos podemos usar el comando fprintf.
fprintf('La raíz cuadrada de dos vale %d aproximadamente.\n',sqrt(2))
La raíz cuadrada de dos vale 1.414214e+00 aproximadamente.
Vemos que en el ejemplo hemos puesto una cadena de caracteres que indica el formato en el que vamos a presentar el resultado, y después de la coma el valor que vamos a mostrar.
Al mostrar el resultado, donde pone %d se sustituye el valor de sqrt(2) en formato decimal.
En el siguiente ejemplo se observa como al no poner \n el siguiente fprintf muestra la salida a continuación de la anterior. Para pasar a la siguiente línea empleamos \n. También vemos que al poner %f se sustituye el correspondiente valor en formato de punto fijo.
fprintf('Mostramos tres números.') fprintf('Uno es %d mientras que otro es %f y el tercero es %f',sqrt(2),sqrt(3),sqrt(4)) fprintf('Sin hacer cambio de linea')
Mostramos tres números.Uno es 1.414214e+00 mientras que otro es 1.732051 y el tercero es 2.000000Sin hacer cambio de linea
Ej 5. En el código de la secante, elimine el símbolo de comentario antes de fprintf para poder ver el proceso del método. Modifique el método de la secante para que se vayan contando las iteraciones efectuadas y se muestren en la instrucción fprintf.
Método de Newton-Raphson
Recordamos el pseudocódigo del método
1. Entrar f , x0, tol 2. Hacer x1 = x0 - f(x0)/f'(x0) 3. Si |x0 - x1| < tol entonces escribir x1. Fin 4. En caso contrario hacer x0 = x1 e ir a 2.
Observamos que para aplicar este método necesitamos evaluar la derivada de la función f. En la implementación newton.m utilizaremos una función fprima.m que permite evaluar la derivada de la función f.m. (Modifique el archivo f.m para que contenga la función e^x+x, y compruebe que el archivo fprima contiene la derivada de f)
Aunque en el pseudocódigo se incluye un salto a la linea 2, en la implementación hay un bucle que se repetirá mientras la variable seguir valga true. Dicha variable valdrá false cuando dos iteraciones sucesivas estén lo suficientemente próximas, lo que hará que el bucle termine. Esta implementación sólo incluye este criterio de parada.
% Método de Newton-Raphson para una función f definida. % Deben estar definidas la función f.m y su derivada fprima.m x0=2; tol=0.00001; seguir=true; while seguir x1=x0-f(x0)/fprima(x0); % fprintf('x1 = %d \n',x1); if abs(x0-x1)<tol seguir=false; end x0=x1; end
Ej 6. Represente la función f(x)=cosh(x).*cos(x)-1 en el intervalo [4,5]. Se pretende encontrar la raíz de f que está en ese intervalo mediante el método de Newton-Raphson. (La derivada de la función cosh(x), coseno hiperbólico, es la función senh(x), seno hiperbólico) ¿Se consigue aproximar dicha raíz tomando como aproximación inicial x0=4? ¿y x0=4.1? (Sol. Si x0=4, se aproxima otra raíz, y si x0=4.1 se obtiene la aproximación x1=4.730041)
Ej 7. Aplique el método de Newton-Raphson a la ecuación f(x)=x^3+94*x^2-389*x+294=0 tomando como aproximación inicial a) x0=1.9, b) x0=2.1, c) x0=2. (Sol. Para x0=1.9 se aproxima a la raíz x=1, para x0=2.1 se aproxima a la raíz x=3, para x0=2 se aproxima a la raíz x=-98).
Ej 8. Modifique el método de Newton-Raphson de forma que se detenga el programa (y muestre el mensaje 'Derivada pequeña') antes de calcular x1 si el valor de f'(x0) en valor absoluto es menor que delta=0.001. (Puede probar el programa con f(x)=(x-2).*(x-7).*(x-2) y x0=5.33333333333)
Derivación simbólica
En el método de Newton-Raphson, para aproximar la solución de f(x)=0, es necesario tener definida la derivada de la función f. En Matlab (no así en Octave) se pueden hallar derivadas simbólicas mediante el comando diff. Observe en los ejemplos cómo se incluye la función a derivar entre comillas simples.
diff(sym('sin(x^2)')) diff(sym('(sin(x)^2)'))
ans = 2*x*cos(x^2) ans = 2*cos(x)*sin(x)
El resultado del comando diff lo podemos copiar y pegar en la definición de la función fprima.
Al usar el comando diff emplearemos los operadores de multiplicación, división y potenciación sin el punto que indica operaciones con vectores. Si los usamos, por ejemplo si ejecutamos la instrucción diff(sym('x.^2')), obtendremos un error.
Ej 9. Halle de forma simbólica la derivada de la función f(x)=cos(x)/sen(x). (Sol. f '(x)= -1 -cos(x)^2/sin(x)^2).
Iteración funcional
En este apartado nos planteamos el problema de hallar un punto fijo de la función g, es decir, hallar x tal que g(x)=x. Para hallar una aproximación al punto fijo de g, implementaremos un método de iteración funcional en el que x(k+1) se obtiene a partir de x(k) de la siguiente manera : x(k+1)=g(x(k)).
Una posible implementación del método puede encontrarse en el programa iteracion.m que requiere que se haya definido la función g.m .
% Iteración funcional para encontrar un punto fijo g(x)=x % Debe proporcionarse la función g.m % Implementa tres condiciones para salir del bucle x0=2; tol1=0.00001; tol2=0.0001; Nmax=100; seguir=true; N=0; while seguir x1=g(x0); N=N+1; fprintf('N = %d, x1 = %d, g(x1)= %d \n',N,x1,g(x1)); % Primera condición: Si dos aproximaciones sucesivas están lo bastante % cerca, se sale del bucle. if abs(x1-x0)<tol1 seguir=false; end % Segunda condición: Si g(x1) es lo bastante cercano a x1, se sale % del bucle. if abs(g(x1)-x1)<tol2 seguir=false; end % Tercera condición. Si se han efectuado más de Nmax iteraciones, se % sale del bucle. if N>Nmax seguir=false; end % Actualizamos el valor de x0 x0=x1; end % Al salir del bucle mostramos el valor de la última iteración. fprintf('Fin del bucle: x0 = %d, g(x0)= %d \n',x0,g(x0));
N = 1, x1 = -4.161468e-01, g(x1)= 9.146533e-01 N = 2, x1 = 9.146533e-01, g(x1)= 6.100653e-01 N = 3, x1 = 6.100653e-01, g(x1)= 8.196106e-01 N = 4, x1 = 8.196106e-01, g(x1)= 6.825059e-01 N = 5, x1 = 6.825059e-01, g(x1)= 7.759946e-01 N = 6, x1 = 7.759946e-01, g(x1)= 7.137247e-01 N = 7, x1 = 7.137247e-01, g(x1)= 7.559287e-01 N = 8, x1 = 7.559287e-01, g(x1)= 7.276348e-01 N = 9, x1 = 7.276348e-01, g(x1)= 7.467496e-01 N = 10, x1 = 7.467496e-01, g(x1)= 7.339006e-01 N = 11, x1 = 7.339006e-01, g(x1)= 7.425676e-01 N = 12, x1 = 7.425676e-01, g(x1)= 7.367349e-01 N = 13, x1 = 7.367349e-01, g(x1)= 7.406663e-01 N = 14, x1 = 7.406663e-01, g(x1)= 7.380191e-01 N = 15, x1 = 7.380191e-01, g(x1)= 7.398028e-01 N = 16, x1 = 7.398028e-01, g(x1)= 7.386015e-01 N = 17, x1 = 7.386015e-01, g(x1)= 7.394108e-01 N = 18, x1 = 7.394108e-01, g(x1)= 7.388657e-01 N = 19, x1 = 7.388657e-01, g(x1)= 7.392329e-01 N = 20, x1 = 7.392329e-01, g(x1)= 7.389856e-01 N = 21, x1 = 7.389856e-01, g(x1)= 7.391522e-01 N = 22, x1 = 7.391522e-01, g(x1)= 7.390400e-01 N = 23, x1 = 7.390400e-01, g(x1)= 7.391156e-01 Fin del bucle: x0 = 7.390400e-01, g(x0)= 7.391156e-01
Ej 10. Represente conjuntamente las gráficas y=x e y=cos(x) en el intervalo [0,3]. Ejecute el programa iteración.m para obtener una aproximación a una raíz de cos(x)=x con x0=0.25. (Sol. aprox. x=7.390261e-01, que se obtiene realizando 22 iteraciones)
Ej 11. Modifique la ecuación x^3-2.1*x^2-x+2.1=0 para que quede en la forma g(x)=x. Utilize el programa iteracion.m para intentar encontrar una raíz cercana a x0=2. (Sol. Si se plantea la ecuación x^3-2.1x^2+2.1=x no se obtiene convergencia. Si se plantea la ecuación x=-2.1/(x^2-2.1x-1) se obtiene una aproximación de otra raíz x=1. Si se plantea x=(2.1x^2+x-2.1)/x^2) se obtiene x0=2.099932 aproximación de la raíz exacta x=2.1).
Ej 12. Modifique el programa para que muestre un mensaje explicativo según el criterio de parada que se verifique.
Ej 13. Para hallar la raíz cuadrada de un número R se propone un método iterativo x(k+1)=g(x(k)) con g(x)=x*(x^2+3*R)/(3x^2+R). Programe el método de forma que la ejecución termine si la iteración cumple abs(x(n)^2-R)<delta=0.000001. Aplique el método para x0=1, R=2. (Sol. aprox. x=1.414214e+00, que se obtiene realizando tres iteraciones)
% Matemáticas II. Grado en Ingeniería Química. (A. Palomares)
disp(date)
10-Mar-2017