%% Práctica 1. Ecuaciones Diferenciales Ordinarias.
% Veremos en esta práctica cómo resolver simbólicamente ecuaciones
% diferenciales usando Wolfram|Alpha y cómo obtener soluciones 
% numéricas mediante el método de Euler.

%% Solución simbólica de una ecuación diferencial
% Para calcular una solución exacta (analítica) de una ecuación diferencial
% basta escribirla adecuadamente en la página de Wolfram|Alpha

%% Solución general de una ecuación diferencial
% Para este ejemplo abrimos la página www.wolframalpha.com
% y escribimos  DSolve[ y'(x)-2*x*y(x)=0 ]
% En los resultados puede ver la clasificación de la ecuación, la solución
% general en términos de una constante c_1, el campo de pendientes e 
% incluso una muestra de la familia de soluciones en términos del valor de y(0).

%%
% Observe que hemos escrito la función incógnita con su variable y(x) .
% Wolfram|Alpha es flexible con la notación que usamos, a cambio debemos
% prestar atención a cómo ha interpretado nuestra petición y,  en  ocasiones, 
% nos puede pedir información adicional.

%%
% Ej 1 . Halle la solución general de x*y'(x)-y^2+9=0 .
% (solución : y(x) = -(3 (e^(6 c_1) x^6 - 1))/(e^(6 c_1) x^6 + 1)
% )

%% Solución particular de un PVI
% De forma similar, incluyendo las condiciones iniciales, se obtiene la
% solución particular de una problema de valores iniciales.
% Introduzca la instrucción DSolve[ y'(x)-2*x*y(x)==0, y(0)==2 ]

%%
% Vemos que en este caso la solución no incluye la constante c_1 y
% obtenemos la gráfica de y(x) .
%
% Ej 2. Halle la solución de x*y'(x)-y^2+9=0 que cumple que y(1)=7 .
% (solución : y(x) = (6 x^6 + 15)/(5 - 2 x^6) ).

%% Circuito RLC
% También es posible resolver ecuaciones de orden superior.
% Escriba la instrucción
% DSolve [ 0.2*i''(t)+330*i'(t)+(1/(20*10^-6))*i(t)=0, i(0)=0, i'(0)=30 ]
% donde la función i(t) es la corriente, variable del problema, no es la
% unidad de los números complejos.
% La ecuación está tomada de https://youtu.be/w9Z4Ly98Lcw?t=143
% donde se resuelve el problema paso a paso. 

%% Solución numérica un PVI. Método de Euler.
% Nos planteamos aplicar el método de Euler al siguiente problema
% de valor inicial
% y'(x)=f(x,y) para x en el intervalo (x0,xn),
% y(x0)=y0 .

%%
% Comenzamos definiendo dos variables para la condición inicial 
% primero x0 y a continuación y0 .
x0=0.0
y0=3.0

%%
% Para implementar el método de Euler aproximamos el valor de la solución
% en unos puntos x0, x1, ..., xn también llamados nodos.
% Definimos ahora el valor de último punto xn, el número de intervalos que 
% definen los puntos n y la distancia entre los puntos h. 
xn = 1
n = 10
h = (xn - x0)/n

%%
% Ahora necesitamos definir la función f. Veremos más adelante cómo definir
% funciones en un script, pero ahora lo haremos mediante una función anónima
% con la siguiente sintaxis con lo que tendremos f(x,y)=2*x*y .
f=@(x,y) 2*x*y

%%
% Definimos los puntos x1, x2, ..., xn en un vector, indicando el primer
% punto, x0, luego la longitud de los intervalos o paso, h, y luego el 
% último punto.
% Inicializamos Y como un vector de ceros (una fila, n+1 columnas).
X=[x0 : h : xn]
Y=zeros(1,n+1)

%%
% Ahora definimos los valores de los puntos yi . Como ilustración hemos
% sacado la definición de Y(2) fuera del bucle.
Y(1)=y0;
Y(2)=Y(1) + f(X(1), Y(1)) * (X(2) - X(1));

for i= 3:n+1
Y(i)=Y(i-1) + f(X(i-1), Y(i-1)) * (X(i) - X(i-1));
end

%% Tabla de resultados
% Para ver los valores hallados, trasponemos el vector con los puntos 
% de X para que quede como una columna, hacemos lo mismo con Y
% y mostramos los dos vectores uno junto a otro.
[X' , Y']

%% Gráfica de la poligonal
% Podemos ver la gráfica de la poligonal formada por los puntos (xi, yi) .
% (El código del método está en el script MetEuler.m )
plot(X,Y)

%% Comparación con los valores exactos
% En este caso podemos comparar los resultados obtenidos por el método
% numérico con los valores exactos porque podemos definir la función exacta
% (analítica) y evaluarla en los puntos del vector X, para lo cual hemos puesto
% un punto antes del signo de elevado al definir la función
yexacta=@(x) 3*exp(x.^2)
Ye=yexacta(X)

%% Tabla de errores
% Para ver mejor los resultados los mostramos en columnas. La primera
% columna son los puntos x0, x1, ..., xn. En la segunda columna se ven
% los valores hallados por el método de Euler. En la tercera columna están
% los valores exactos y en la cuarta la diferencia entre los valores
% exactos y los aproximados, con lo que podemos ver los errores que se
% cometen.
[X', Y', Ye', (Ye-Y)']

%%
% Vemos que en los primeros nodos los errores son relativamente pequeños
% 0.0302, 0.0624,... pero que van creciendo conforme nos alejamos de la
% condición inicial

%% Gráfica
% Para comparar gráficamente los resultados con la solución exacta del problema,
% definimos un vector XE con más puntos, definimos una función yexacta
% con la solución exacta del problema 3*e^(x^2), y evaluamos la función
% en todos los puntos de XE, para lo cual hemos puesto un punto antes
% del signo de elevado al definir la función.
XE=[x0:h/4:xn];
YE=yexacta(XE);
plot(X,Y,XE,YE)

%%
% Ej 3. En la gráfica se representan dos funciones ¿cuál es la exacta y cuál
% la aproximada? (Indicación: mire la tabla de errores.)

%% Ejercicios finales
% Ej 4. Se quiere emplear el método de Euler con h = 0.05 para hallar una solución
% aproximada de y(2) en la misma ecuación diferencial del ejemplo 
% anterior (y'(x)=f(x,y)=x+y, y(0)=3).
% 
% a) Halle el número n de partes en las que se debe dividir el intervalo [0, 2] 
% y obtenga dicha solución aproximada.
%
% b) Use el mismo paso h=0.05 para hallar una aproximación a y(3).

%%
% Ej 5. Se considera el problema de valor inicial y'(x) = 3y^2 -x^2,
% y(0)=0.15
%
% a) Halle una aproximación de y(3) usando h=0.1.
%
% b) Halle una aproximación de y(3) usando h=0.05.
%
% c) Suponiendo que y(3) vale -1.673381304457, ¿cuál de los dos valores 
% de h permite obtener la mejor aproximación de y(3)?
% (solución :  Exacto - aproximado=0.0001, 
% Exacto - aproximado=5.418e-05, más pequeño )

%%
% Ej 6. Dentro del bucle del método de Euler, ¿se puede sustituir (X(i) - X(i-1)) por h?
% Justifique la respuesta.

%%
% Ej 7. Halle una aproximación de la solución de 
% y'(x)=2y-(3/10)y^2,  y(0)=1/10
% en x=4 usando el método de Euler en el intervalo (0,4) para n=25.
% Represente la poligonal del Euler.
%
% Sabiendo que la solución exacta es y(x)= (20 e^(2x)) / (197 + 3e^(2x)) ¿Qué error
% absoluto se está cometiendo? (Valor exacto menos valor aproximado).
% (Sol. y(4)=~ 6.5164. Error=0.00653 .)

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