%% Práctica 4. Resolución de ecuaciones (II)
% Esta práctica tratará el método de la secante, los métodos de iteración 
% funcional y el método de Newton-Raphson.

%% 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 
% esos dos valores de la función. La siguiente iteración es el corte con el
% eje OX de la recta.

%%
% El siguiente pseudocódigo permite implementar el método de la secante:
%
%  1. Entrar f , a, b, tol
%  2. Mientras |b - a| >= tol
%  3. Hacer c = b - (b-a)*f(b)/(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 - (b-a)*f(b)/(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)

%%
% 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))

%%
% 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 usando notación exponencial.

%%
% 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 (con seis dígitos tras el 
% punto decimal).
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')

%%
% 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.

%% 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));

%%
% Ej 6. 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. (Sol. aprox. x=7.390367e-001) 
%
% Ej 7. 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 8. Modifique el programa para que muestre un mensaje explicativo según
% el criterio de parada que se verifique.
%
% Ej 9. 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.

%% Método de Newton-Raphson
% Mostramos un 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 10. 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 sinh(x), seno hiperbólico) ¿Se consigue 
% aproximar dicha raíz tomando como aproximación inicial x0=4? ¿y x0=4.1? 
% (Sol. Con Matlab R2007b, si x0=4, se aproxima otra raíz, y si x0=4.1 se 
% obtiene la aproximación x1=4.730041) 
%
% Ej 11. 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 12. 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 se pueden 
% hallar derivadas simbólicas mediante el comando diff. En Octave 
% se requiere instalar un paquete externo.
% En Wolfram|Alpha se pueden hacer derivadas simbólicas mediante D,
% por ejemplo D[ Sin[x], x] .
% 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)'))

%%
% 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 13. 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).

% Ecuaciones Diferenciales y Cálculo Numérico. (A. Palomares)
% Grado en Ingeniería de Tecnologías de Telecomunicación.
disp(date)