Merge project of MATLAB utilities to this repo

.gitattributes

@@ -0,0 +1,2 @@
+# Auto detect text files and perform LF normalization
+* text=auto

Numerical Methods - MATLAB&Octave/Diff eq/Tests_euler-euler_mod-heun.m

@@ -0,0 +1,63 @@
+% TESTS FOR FUNCTIONS:
+%   euler.m
+%   euler_mod.m
+
+% USING: TEMA 5 - EJERCICIO 9
+%{
+*** ENUNCIADO ***
+Sea el problema del valor inicial (PVI) { y' + 30y = 0, 0<=t<=1, y(0) = 1 }
+cuya solución exacta es e^(-30t),
+aproximar por el método de Euler de pasos 0.1, 0.05 y 0.01
+y representar gráficamente junto con la solución exacta.
+Comentar los resultados
+%}
+
+% Clear workspace and screen, set up format configuration
+clear, clc;
+format long g;
+
+% Runtime config variables
+%   Visualización de datos
+nx_plots = 1;
+ny_plots = 3;
+
+% Función y dato inicial
+f = @(t, y) - 30 .* y;
+y0 = 1;
+
+% Solución exacta
+y_sol = @(t) exp(-30 .* t);
+
+% Intervalo para aproximar
+a = 0; b = 1;
+
+% Tamaños de paso y pasos
+h10 = 0.01; M10 = round( ( b - a ) / h10 )
+h20 = 0.01; M20 = round( ( b - a ) / h20 )
+h30 = 0.01; M30 = round( ( b - a ) / h30 )
+
+% Cálculo de las soluciones por Euler
+E10 = euler(f, a, b, y0, M10); disp(E10);
+E20 = euler_mod(f, a, b, y0, M20); disp(E20);
+E30 = heun(f, a, b, y0, M30); disp(E30);
+
+% Representación gráfica de las soluciones
+% La exacta
+t = a:0.01:b;
+y = y_sol (t);
+
+% Aproximaciones
+subplot (nx_plots, ny_plots, 1);
+plot(t, y, 'r', E10(:,2), E10(:,3), 'b- ' , "LineWidth", 2 );
+title("Metodo de Euler");
+grid on;
+
+subplot (nx_plots, ny_plots, 2);
+plot(t, y, 'r', E20(:,2), E20(:,3), 'b- ' , "LineWidth", 2 );
+title("Metodo de Euler modificado");
+grid on;
+
+subplot (nx_plots, ny_plots, 3);
+plot(t, y, 'r', E30(:,2), E30(:,3), 'b- ' , "LineWidth", 2 );
+title("Metodo de Heun");
+grid on;

Numerical Methods - MATLAB&Octave/Diff eq/euler.m

@@ -0,0 +1,34 @@
+% Author: Echedey Luis Álvarez
+% Date: 21/05/2021
+% 
+% Abstract: Aplicación general de aproximación de ec. diferenciales por el método de Euler
+
+function E = euler( f, a, b, y_a, M )
+%{
+    Args:
+      f: 2-vars anonym function which represents explicitly the derivative of y
+      a & b: begin and end of interval
+      y_a: known value of f, f(a)
+      M: desired number of iterations. M gets rounded to the nearest int32
+    Output:
+      E: (M+1) x 3 matrix
+        E(:,1): index column (first one is 0)
+        E(:,2): independent value (t_k, x_k)
+        E(:,3): dependent value (y_k)
+%}
+
+    if ( a >= b )
+      error("Cannot aproximate within given interval");
+    end
+    M = double ( uint32 (M) );
+    if (M == 0)
+      error("Number of iterations must be a natural number");
+    end
+
+    h = (b - a) / M;
+    E = NaN( [M+1, 3] );
+    E(1, :) = [0, a, y_a];
+    for j=1:M
+       E(j+1, :) = [ j, a + j .* h, E(j, 3) + h * f( E(j, 2), E(j, 3) ) ];
+    end % !for
+end % !if

Numerical Methods - MATLAB&Octave/Diff eq/euler_mod.m

@@ -0,0 +1,37 @@
+% Author: Echedey Luis Álvarez
+% Date: 24/05/2021
+%
+% Abstract: Aplicación general de aproximación de ec. diferenciales por el método de Euler modificado
+%   Se utiliza la pendiente en el punto medio de del paso, es decir, en t_k + h/2
+
+function E = euler_mod( f, a, b, y_a, M )
+%{
+    Args:
+      f: 2-vars anonym function which represents explicitly the derivative of y
+      a & b: begin and end of interval
+      y_a: known value of f, f(a)
+      M: desired number of iterations. M gets rounded to the nearest int32
+    Output:
+      E: (M+1) x 3 matrix
+        E(:,1): index column (first one is 0)
+        E(:,2): independent value (t_k, x_k)
+        E(:,3): dependent value (y_k)
+%}
+
+    if ( a >= b )
+      error("Cannot aproximate within given interval");
+    end
+    M = double ( uint32 (M) );
+    if (M == 0)
+      error("Number of iterations must be a natural number");
+    end
+
+    h = (b - a) / M;
+    E = NaN( [M+1, 3] );
+    E(1, :) = [0, a, y_a];
+    for j=1:M
+      f1 = f( E(j, 2), E(j, 3) );
+      f2 = f( E(j, 2) + h/2, E(j, 3) + h/2 * f1 );
+      E(j+1, :) = [ j, a + j .* h, E(j,3) + h * f2 ];
+    end % !for
+end % !if

Numerical Methods - MATLAB&Octave/Diff eq/heun.m

@@ -0,0 +1,38 @@
+% Author: Echedey Luis Álvarez
+% Date: 25/05/2021
+%
+% Abstract: Aplicación general de aproximación de ec. diferenciales por el método de Heun
+%   Se utiliza la media de las pendientes en  t_k y t_(k+1)
+
+function E = euler_mod( f, a, b, y_a, M )
+%{
+    Args:
+      f: 2-vars anonym function which represents explicitly the derivative of y
+      a & b: begin and end of interval
+      y_a: known value of f, f(a)
+      M: desired number of iterations. M gets rounded to the nearest int32
+    Output:
+      E: (M+1) x 3 matrix
+        E(:,1): index column (first one is 0)
+        E(:,2): independent value (t_k, x_k)
+        E(:,3): dependent value (y_k)
+%}
+
+    if ( a >= b )
+      error("Cannot aproximate within given interval");
+    end
+    M = double ( uint32 (M) );
+    if (M == 0)
+      error("Number of iterations must be a natural number");
+    end
+
+    h = (b - a) / M;
+    E = NaN( [M+1, 3] );
+    f1 = f2 = 0;
+    E(1, :) = [0, a, y_a];
+    for j=1:M
+      f1 = f( E(j, 2), E(j, 3) );
+      f2 = f( E(j, 2) + h, E(j, 3) + h * f1 );
+      E(j+1, :) = [ j, a + j .* h, E(j,3) + h/2 * ( f1 + f2 ) ];
+    end % !for
+end % !if

Numerical Methods - MATLAB&Octave/Equations - Non linear/Test_funcionesNRySecante.m

@@ -0,0 +1,25 @@
+% El único propósito de este programa es comprobar que las funciones
+% funcionan, valga la redundancia.
+% El documento y el nombre de las variables son self-explanatory
+
+clear, clc
+format long
+
+f = @(x) x.^4 - 3*x - 5; % Esto me lo saco de la manga
+df = @(x) 4*x.^3 - 3; % Y esto es la manga diferenciada
+
+pto0 = 1; % Con el Geogebra se grafica f(x) y tomo un valor cualquiera
+pto1 = 1.1; % Para el met. de la secante
+err_rel = 0.001; % Porque es un buen número
+
+[NR_p, NR_deltaRel] = metodo_newton_raphson(f, df, pto0, 20, err_rel);
+
+[SEC_p, SEC_deltaRel] = metodo_secante(f, pto0, pto1, 20, err_rel);
+
+% Debe converger a 1.795176530
+disp("Por el método de Newton-Raphson");
+porNewtonRaphson = table(NR_p, NR_deltaRel)
+disp("Por el método de la secante");
+porSecante = table(SEC_p, SEC_deltaRel)
+
+% Solo quiero decir que esto es maravilloso

Numerical Methods - MATLAB&Octave/Equations - Non linear/metodo_newton_raphson.m

@@ -0,0 +1,53 @@
+% Autor: Echedey Luis Álvarez
+% Fecha: 07/03/2021
+% 
+% Aplicación general del método de Newton-Raphson.
+
+function [p, err_rel] = metodo_newton_raphson(f, df, p0, max_iter, deltaRel)
+% ----Aproxima la raíz de una función por el método de Newton-Raphson------
+% Argumentos:
+%   f: función a calcular raíz
+%   df: primera derivada de f
+%   p0: valor de la abscisa inicial
+%   max_iter: número máximo de iteraciones
+%   deltaRel: error relativo entre iteraciones consecutivas. Si es menor,
+%       el programa finaliza aunque no haya alcanzado max_iter. Indique 0
+%       para que se hagan max_iter iteraciones
+% Salida:
+%   p: vector columna de tamaño i con la sucesión de puntos p que convergen
+%       a alguna solución de f(p) = 0
+%   err_rel: vector columna de tamaño i. Tiene el error relativo entre
+%       p(n) y p(n-1) para todos los p excepto el primero. El primer p
+%       tiene error NaN porque no existe número anterior a él que permita
+%       calcularlo.
+%           El tamaño i es el número de iteraciones que se han hecho hasta
+%           que se ha finalizado el algoritmo.
+
+% Allocamos la memoria máxima que usa o puede usar el programa
+% Inicializar en NaN para luego eliminar lo que sobre.
+p = NaN(max_iter, 1);
+err_rel = p; % El primero no se calculará porque no se puede sin otro dato inicial
+
+% Inicialización variables
+p(1)=p0; % Asignación dato inicial para la convergencia
+
+j=1; % Inicialización contador
+g = @(x) x - f(x)./df(x); % Función del método iterativo
+
+% Calcula la sucesión
+while ( (err_rel(j)>=deltaRel) || isnan(err_rel(j)) ) && (j <= max_iter)
+    p(j+1)=g(p(j));
+    err_rel(j+1)=abs(p(j+1)-p(j))/abs(p(j+1) + eps);
+    j=j+1;
+end % Del bucle de iteraciones
+
+% Elimina los NaNs finales (recorta los vectores columna)
+invID = max_iter;
+while isnan(p(invID)) % Así no obtengo un NaN. Véase función.
+    invID = invID - 1;
+end % Del bucle para cortar los arrays
+
+p = p(1:invID);
+err_rel = err_rel(1:invID);
+
+end % De la función

Numerical Methods - MATLAB&Octave/Equations - Non linear/metodo_secante.m

@@ -0,0 +1,58 @@
+% Autor: Echedey Luis Álvarez
+% Fecha: 07/03/2021
+% 
+% Aplicación general del método de la secante. Es una copia barata del de Newton-Raphson, pero con un par de cambios.
+
+
+function [p, err_rel] = metodo_secante(f, p0, p1, max_iter, deltaRel)
+% ------Aproxima la raíz de una función por el método de Newton-Raphson--------
+% Argumentos:
+%   f: función a calcular raíz
+%   p0: primer valor de la abscisa inicial
+%   p1: segundo valor de la abscisa inicial
+%   max_iter: número máximo de iteraciones
+%   deltaRel: error relativo entre iteraciones consecutivas. Si es menor,
+%       el programa finaliza aunque no haya alcanzado max_iter. Indique 0
+%       para que se hagan max_iter iteraciones
+% Salida:
+%   p: vector columna de tamaño i con la sucesión de puntos p que convergen a alguna
+%       solución de f(p) = 0
+%   err_rel: vector columna de tamaño i. Tiene el error relativo entre
+%       p(n) y p(n-1) para todos los p excepto el primero. El primer p
+%       tiene error NaN porque no existe número anterior a él que permita
+%       calcularlo.
+%           El tamaño i es el número de iteraciones que se han hecho hasta que se
+%           ha finalizado el algoritmo.
+
+% Allocamos la memoria máxima que usa o puede usar el programa
+% Inicializar en NaN para luego eliminar lo que sobre.
+p = NaN(max_iter, 1);
+err_rel = p;
+
+% Inicialización variables
+% Asignación datos iniciales para la convergencia
+p(1)=p0;
+p(2)=p1;
+% Asignación inicial para el error. Imposible calcular
+err_rel(1) = NaN;
+err_rel(2) = abs(p(2)-p(1))/(abs(p(2)) + eps); % Nótese el uso de eps para no dividir por cero
+
+j=2; % Inicialización contador
+
+% Calcula la sucesión
+while ( (err_rel(j)>=deltaRel) || isnan(err_rel(j)) ) && (j <= max_iter)
+    p(j+1)=p(j) - f(p(j)).*(p(j)-p(j-1))./(f(p(j))-f(p(j-1)));
+    err_rel(j+1)=abs(p(j+1)-p(j))/(abs(p(j+1)) + eps); % Nótese el uso de eps para no dividir por cero
+    j=j+1;
+end % Del bucle de iteraciones
+
+% Elimina los NaNs finales (recorta los vectores columna)
+invID = max_iter;
+while isnan(p(invID)) % Así no obtengo un NaN. Véase función.
+    invID = invID - 1;
+end % Del bucle para cortar los arrays
+
+p = p(1:invID);
+err_rel = err_rel(1:invID);
+
+end % De la función

Numerical Methods - MATLAB&Octave/Exercises/Cuestionario-5-integracion/formula_simpson_composite.m

@@ -0,0 +1,21 @@
+% Autor: Echedey Luis Álvarez
+% Fecha: 11/05/2021
+% 
+% Aplicación general de aproximación de integrales por la Fórmula de Simpson Compuesta
+
+function aprox = formula_simpson_composite(h, y)
+%{ 
+Aproxima la integral de una función que pasa por (xi, yi)
+Argumentos:
+  y: array puntos x de los vectores
+  h: distancia homogénea entre cada par de puntos
+Salida:
+  aprox: valor aproximado de la integral
+%}
+if ( mod(size( y ), 2) == 0)
+  error("formula_simpson_composite: Number of nodes must be odd (to match an even number of subintervals).");
+end
+
+aprox = ( h / 3 ) * ( y(1) + y(end) + 4 * sum( y( 2:2:end-1 ) ) + 2 * sum( ( y(3:2:end-2 ) ) ) );
+
+end % De la función

Numerical Methods - MATLAB&Octave/Exercises/Cuestionario-5-integracion/formula_trapezium_composite.m

@@ -0,0 +1,18 @@
+% Autor: Echedey Luis Álvarez
+% Fecha: 11/05/2021
+% 
+% Aplicación general de aproximación de integrales por la Fórmula del Trapecio Compuesta
+
+function aprox = formula_trapezium_composite(h, y)
+%{ 
+Aproxima la integral de una función que pasa por (xi, yi)
+Argumentos:
+  y: array puntos x de los vectores
+  h: distancia homogénea entre cada par de puntos
+Salida:
+  aprox: valor aproximado de la integral
+%}
+
+aprox = h * ( (y(1) + y(end) ) / 2 + sum( y( 2:end-1 ) ) );
+
+end % De la función