Métodos Numéricos
Ejercicios prácticos
Autores:
Elena E. Álvarez Saiz
Código JavaScript para el libro:
Joel Espinosa Longi,
IMATE, UNAM.
Recursos interactivos:
DescartesJS, WebSim
Videos: Math-GPT
Fuentes: Lato y
UbuntuMono
Imagen de portada: ilustración generada por
Pollinations AI
Red Educativa Digital Descartes
Córdoba (España)
descartes@proyectodescartes.org
https://proyectodescartes.org
Proyecto iCartesiLibri
https://proyectodescartes.org/iCartesiLibri/index.htRes
ISBN:
Esta obra está bajo una licencia Creative Commons 4.0 internacional: Reconocimiento-No
Comercial-Compartir Igual.
Los métodos numéricos son una herramienta esencial para los ingenieros, ya que permiten resolver problemas complejos que a menudo no tienen una solución analítica simple. Estas técnicas posibilitan abordar situaciones reales mediante aproximaciones, convirtiéndose en un recurso indispensable en la práctica profesional.
Este libro está diseñado para quienes deseen adquirir una comprensión de los métodos numéricos para el cálculo de raíces de funciones no lineales, la aproximación de funciones mediante polinomios, la integración numérica, la resolución de ecuaciones diferenciales ordinarias y la optimización.
Este libro ofrece:
Esta orientación del libro busca no solo facilitar el aprendizaje, sino también despertar el interés por profundizar en este campo. Se busca así acompañar al lector en el aprendizaje de los métodos numéricos mediante una exposición clara y organizada. Se busca que cada
contenido contribuya a construir una base sólida, permitiendo conectar conceptos teóricos con su aplicación.
Este libro está inspirado en los apuntes del profesor Eduardo Casas Rentería de la Universidad
de Cantabria,
La estructura del libro comprende siete temas principales, que abarcan desde conceptos básicos hasta aplicaciones avanzadas. El primer tema introduce los métodos numéricos y el uso de Matlab/Octave, estableciendo una base para el aprendizaje y la implementación computacional. El segundo tema cubre la aritmética computacional básica, destacando los principios para trabajar con números en un entorno digital. En el tercer tema se aborda la resolución aproximada de ecuaciones no lineales, explorando métodos para encontrar raíces de funciones. El cuarto tema se dedica a la aproximación de funciones mediante polinomios, una herramienta clave en el análisis numérico. El quinto capítulo está enfocado en la integración numérica, presentando técnicas para calcular integrales de manera aproximada. Posteriormente, se introduce la resolución de ecuaciones diferenciales ordinarias, ofreciendo métodos aproximados para su solución. Finalmente, el libro concluye con el tema de optimización, abordando estrategias para encontrar soluciones óptimas a problemas matemáticos y aplicados.
Consciente de que puedan existir erratas o inconsistencias en este libro, animamos a los lectores a compartir cualquier observación o corrección que consideren necesaria. Estas contribuciones serán de
gran ayuda para mejorar futuras ediciones y garantizar que este recurso sea de utilidad.
Quiero agradecer sinceramente a la Red Educativa Digital Descartes y, en especial, a Juan Guillermo Rivera Berrio por compartir su conocimiento de manera desinteresada. Su experiencia y su compromiso con el aprendizaje han sido clave para enriquecer este trabajo.
El análisis numérico es una rama de las matemáticas que estudia métodos y algoritmos para obtener soluciones aproximadas a problemas que no pueden resolverse de forma exacta mediante métodos analíticos convencionales. Entre estos problemas se ncluyen ecuaciones algebraicas, ecuaciones diferenciales, integrales, sistemas de ecuaciones lineales y no lineales, entre otros.
El estudio de los métodos numéricos requiere abordar las siguientes cuestiones:
Para comprender la importancia de los métodos numéricos y su aplicación en problemas reales, consideremos un sistema de baterías conectado a un panel solar que genera electricidad a lo largo del día. Se quiere modelar cómo la batería se carga con el tiempo, sabiendo que:
La ecuación que podría modelizar cómo la energía se almacena en la batería podría ser la siguiente ecuación diferencial ordinaria (EDO):
$$ dE/dt = P_{solar}(t) - kE(t) $$donde:
Vamos a recurrir a técnicas numéricas para obtener una solución aproximada considerando los siguientes datos y parámetros:
La gráfica de la Figura 1.1 muestra la evolucion de la energía almacenada en la batería en un día (24 horas) en función de la potencia solar disponible utilizando como método de resolución el método de Euler que se analizará en detalle en el capítulo 6 de este libro.
Se observa cómo la batería se carga durante las horas de sol y comienza a descargarse lentamente debido a las ineficiencias cuando no hay suficiente energía solar disponible.
El ejemplo responde a un problema de valor inicial (PVI) que consiste en encontrar la curva solución de una ecuación diferencial
$\frac{dy}{dx} = f(x, y)$ que verifica una condición inicial $y(x_0)=y_0$, es decir, que pasa por el punto $(x_0,y_0)$. En estos problemas se conoce la pendiente de la curva \( y(x) \) en cualquier punto y se busca encontrar dicha función.
El método de Euler propone aproximar esta curva avanzando paso a paso a través de pequeños saltos en el eje \( x \) moviéndose suavemente por la recta tangente en cada instante ya que $f$ permite calcular la pendiente a la curva en cada punto. En este proceso, en cada paso se recalcula la dirección.
A continuación se detalla el proceso a seguir en este método numérico, que se inicia desde el punto $(x_0.y_0)$.
En el ejemplo que estamos viendo, el problema de valor inicial es
$$ dE/dt = P_{solar}(t) - kE(t) \,\,\,\, E(0)=0$$ El método de Euler permitiría encontrar la solución con la siguiente fórmula: $$ E_{n+1} = E_n + Δt \cdot (P_{solar}(t_n) - k \cdot E_n) $$donde $Δt$ es el paso de tiempo (1 hora) y $E_n$ es la energía almacenada en la batería en el instante $n$.
El método permite calcular la energía almacenada en la batería a lo largo de 24 horas utilizando la potencia solar que varía con el tiempo. Mostramos el cálculo de los primeros pasos.
Condición inicial:El proceso continuaría completando de forma similar el cálculo para las 24 horas del día.
El 25 de febrero de 1991, durante la Guerra del Golfo, una batería de misiles Patriot
estadounidense en Dharan, Arabia Saudita, no pudo interceptar un misil Scud iraquí entrante. El
Scud atacó un cuartel del ejército estadounidense y mató a 28 soldados. Un informe de la
entonces Oficina General de Contabilidad titulado "Patriot Missile Defense: Problema de software
conducido a una falla del sistema en Dhahran, Arabia Saudita" informó sobre la causa del fallo:
un cálculo inexacto del tiempo desde el inicio debido a errores aritméticos del computador
Específicamente, el tiempo en décimas de segundo, medido por el reloj interno del sistema, se
multiplicó por 10 para producir el tiempo en segundos. Este cálculo se realizó utilizando un
registro de punto fijo de 24 bits. En particular, el valor 1/10, que tiene una expansión binaria
infinita, se cortó a 24 bits después del punto de base. El pequeño error de corte, cuando se
multiplica por el gran número que da el tiempo en décimas de segundo, conduce a un error
significativo
El Programa Ariane, es un programa emprendido por Europa en 1973 para dotarse de un lanzador que
le permitiera un acceso independiente al espacio. El desarrollo del lanzador Ariane se efectúa
bajo la dirección de la Agencia Espacial Europea (ESA, por sus siglas en inglés); el Centro
Nacional de Estudios Espaciales (CNES) francés fue el contratista principal hasta mayo de 2003,
fecha en la que pasó a actuar como tal el consorcio europeo EADS (European Aeronautic Defenceand
Space Company)
La explotación comercial es gestionada por la sociedad Arianespace, creada en 1980. En total, unas 40 compañías europeas están comprometidas en el desarrollo y construcción del lanzador Ariane.
Todos los lanzamientos se efectúan desde el centro espacial de Kourou, en la Guayana Francesa, que cuenta con varias rampas de lanzamiento, y en el que trabajan de forma permanente varios cientos de personas.
El 4 de junio de 1996, el cohete Ariane 5 Flight 501,de la Agencia Europea del Espacio (ESA) explotó 40 segundos después de su despegue a una altura de 3.7km. tras desviarse de la trayectoria prevista.
Era su primer viaje tras una década de investigación que costó más de 7000 millones de euros. El cohete y su carga estaban valorados en más de 500 millones de euros. La causa del error fue un fallo en el sistema de guiado de la trayectoria provocado 37 segundos después del despegue. Este error se produjo en el software que controlaba el sistema de referencia inercial (SRI).
En concreto, se produjo una excepción debido al intento de convertir un número en punto flotante de 64 bits, relacionado con la velocidad horizontal del cohete respecto de la plataforma de lanzamiento, en un entero con signo de 16 bits. El número mas grande que se puede representar de esta forma es 32767.
Más información se puede encontrar en
El 23 de agosto de 1991, la plataforma petrolífera Sleipner A propiedad de la empresa noruega
Statoil situada en el mar del Norte a 82 metros de profundidad se hundió. La causa del error fue
un fallo en el modelado numérico de la plataforma utilizando elementos finitos. En concreto se
produjo una fuga de agua en una de las paredes de uno de los 24 tanques de aire de 12 metros de
diámetro que permitían la flotación de la plataforma de 57000 toneladas de peso que además
soportaba a mas de 20 0personas y a equipamiento de extracción con un peso adicional de unas
40000 toneladas. Las bombas de extracción de agua no fueron capaces de evacuar toda el agua. El
fallo tuvo un coste económico total de 700 millones de euros
Para modelar los tanques de la plataforma se utilizó el programa NASTRAN de elementos finitos y una aproximación mediante un modelo elástico lineal. Esta aproximación no era correcta y subestimó en un 47% los esfuerzos que debían soportar las paredes de los tanques. En particular, ciertas paredes fueron diseñadas con un grosor insuficiente.
Un análisis con elementos finitos correcto, pero posterior al accidente, demostró que el diseño de la plataforma provocaría fugas en algunas de los tanques cuando ésta estuviese sobre agua a 62 metros de profundidad. La fuga real se produjo cuando ésta estaba sobre 65 metros de agua, lo que ratifica la supuesta causa del fallo.
En este libro, se utiliza Matlab como la herramienta principal para desarrollar y resolver los ejercicios y ejemplos propuestos. Su interfaz intuitiva y su amplia biblioteca de funciones lo convierten en una herramienta ideal para ilustrar conceptos matemáticos, realizar simulaciones y aplicar métodos computacionales en un entorno práctico.
Matlab es un acrónimo de "Matrix Laboratory", lo que refleja su origen como una herramienta diseñada para trabajar principalmente con matrices. Con el tiempo, el programa ha evolucionado para convertirse en un entorno de programación de alto nivel que permite realizar cálculos numéricos, análisis de datos, simulaciones, desarrollo de algoritmos y visualización gráfica.
Aunque Matlab es un software de pago, se puede utilizar también el software libre y gratuito Octave. Para facilitar ejecutar el código propuesto en los ejemplos y ejercicios se pueden utilizar los enlaces siguientes sin realizar ninguna instalación:
A continuación se incluye como recordatorio las principales características de Matlab/Octave necesarias para seguir y comprender los ejemplos presentados en los próximos capítulos. Para una explicación más detallada, se recomienda consultar la documentación oficial o la bibliografía
La ventana de Matlab muestra un escritorio dividido en varias partes:
Cuando guardamos código en un archivo y queremos utilizarlo más adelante, es necesario indicar al programa su ubicación. Matlab tiene preconfigurada la ruta donde almacena sus funciones nativas y las de los paquetes incluidos. Sin embargo, si trabajamos con archivos propios, debemos asegurarnos de que el programa pueda acceder a ellos correctamente.
Matlab busca funciones y scripts en los directorios especificados por el comando path. El
primero de ellos es siempre el especificado en el diálogo
Current Directory.
En la figura aparece marcado en rojo el directorio actual. Tecleando el comando
path en la ventana de comandos puede ver el listado de directorios
activos.
Desde la ventana de comandos, se puede utilizar el comando help
help nombreDelComando >>help plot
Además, desde el menú help se puede acceder a la ayuda de MathWorks en el caso de Matlab.
En Matlab/Octave, crear variables es sencillo, basta asignar un valor al nombre de la variable. Las variables se almacenan en la memoria y pueden contener números, matrices, cadenas de texto, entre otros. El contenido de una variable se puede recuperar y modificar a lo largo de una sesión de trabajo.
Se detallan a continuación algunas reglas a tener en cuenta a la hora de definir los nombres de las variables:
modulo2 es un nombre
válido, pero no lo es 2modulo.
Modulo es distinta de la variable
modulo.
modulo1 es un nombre de variable válido, pero no
modulo 1.
ans, pi, Inf, i, j, ... . Para obtener una lista de palabras
reservadas teclea en la venta de comando iskeyword().
Un ejemplo de declaración de variable:
x=5; % x es una variable escalar
x=x+6; % El valor de x es 11
Las variables creadas desde la línea de comandos pertenecen al espacio de trabajo (Workspace).
Matlab ignora cualquier texto que vaya precedido por el símbolo %, por
lo tanto, este símbolo sirve para incluir comentarios.
>> x=5; 2*x; y=x^2; x=y/x;
El comando format en Matlab/Octave se utiliza para controlar la forma en
que Matlab/Octave muestra los números en la ventana de comandos. No cambia la precisión interna de
los cálculos, solo afecta a la presentación de los resultados. Aquí tienes algunos de los modos más
comunes de este comando al escribir el nùmero $\pi$.
format short
| ||
format long |
15 dígitos después del punto decimal | 3.141592653589793 |
format shortEng |
Notacion decimal con 4 dígitos | 3.1416e+00 |
format longEng |
Notacion técnica larga | 3.14159265358979e+000 |
Matlab/Octave permite realizar operaciones aritméticas básicas de manera directa. Las operaciones estándar incluyen suma (+), resta (-), multiplicación (*), división (/) y exponenciación (^). Por ejemplo,
x = 4; y = x + 3; % Suma
z = x * 2; % Multiplicación
A= 2*(x + z); % El valor de A es 20
Las expresiones se evalúan de izquierda a derecha, la potencia tiene el orden de prioridad mayor, seguido del producto y la división (ambas
tienen la misma prioridad) y por último la suma y la resta (con igual prioridad entre ellas). Para alterar este orden se deben introducir adecuadamente paréntesis. Por ejemplo,
x= 2+3^4/2; % x toma el valor 2+81/2
B =2+3^(4/2); % B toma el valor 11=2+3^2
Un script es un archivo que contiene un conjunto de comandos de Matlab/Octave que se ejecutan en orden. Para crear un script, simplemente se deben escribir los comandos en un archivo con extensión .m de la forma nombre.m siendo nombre una secuencia de caracteres que no admite caracteres blancos.
También se puede crear un guión desde el editor de Matlab/Octave eligiendo el menú New Script.
Algunas funciones sencillas, se pueden definir mediante una sola instrucción en la ventana de comandos o en un script. Basta establecer su nombre, los argumentos de los que dependen y los comandos que permiten definirla.
nombre_Funcion=@(argumentos) expresion_FuncionPor ejemplo, para definir la función que calcula el área de un rectángulo a partir de su base y altura
>>area_rectangulo=@(b,h) b*h
>>valor=area_rectangulo(4,3)
valor=
12punto_inicial:paso:punto_finalPor ejemplo, para generar un vector con punto inicial 2 hasta 7 con un salto 2.
>>v=4:2:9
v=[4 6 8]linspace(valor_inicial,valor_final,num_puntos)Por ejemplo, para generar un vector fila de 7 elementos equiespaciados siendo el primer valor 2 y el último 3
>>v=linspace(2,3,7)
v=[2.0000 2.1667 2.3333 2.5000 2.6667 2.8333 3.0000]
Una matriz es una colección de valores en filas y columnas. Cada fila se separa con un punto y coma (;) y los elementos de una fila se separan con espacios o comas. Por ejemplo,
M = [1 2 3; 4 5 6; 7 8 9]; % Es una matriz 3x3
M = [1 2 3; 4 5 6]; % Es una matriz 2x3
Veamos cómo generar algunas matrices especiales.
Matriz de ceros
zeros(numFilas,numColumnas)
% Ejemplo: una matriz 4x3 de ceros
>>zeros(4,3)
Matriz de unos
ones(numFilas,numColumnas)
%Ejemplo: Matriz 4x3 con todos sus elementos unos
>>ones(4,3)
Matriz identidad
eye(n)
% Ejemplo: Matriz identidad 3x3
>>eye(3)
Matriz diagonal
diag(vector)
%Ejemplo: Una matriz diagonal con el vector [2 3 4] en la diagonal
>>diag([2 3 4])
Matriz aleatoria
rand(m,n) randi([imin imax],m,n)
%Ejemplo: Crea una matriz de valores aleatorios de dimension 2x3 y otra del mismo tamaño con números enteros entre -5 y 5
>>rand(2,3)
>>randi([-5 5],2,3)
A = [1 2 3; 4 5 6; 7 8 9];
% Accede al elemento en la fila 2, columna 3
elemento = A(2, 3);
% Resultado: elemento = 6
A = [1 2 3; 4 5 6; 7 8 9]; fila2 = A(2, :);
% Accede a todos los elementos de la fila 2
% Resultado: fila2 = [4 5 6]
A = [1 2 3; 4 5 6; 7 8 9]; columna3 = A(:,3);
% Accede a todos los elementos de la columna3
% Resultado: columna3 = [3;6;9]
A = [1 2 3; 4 5 6; 7 8 9]; % Matriz 3x3
ultima_fila = A(end, :)
ultimas_columnas=A(:,2:end)
A = [1 2 3; 4 5 6; 7 8 9];
% Submatriz con filas 1 y 2, y columnas 2 y 3
submatriz = A(1:2,2:3);% Resultado: [2 3;5 6]
A = [1 2 3; 4 5 6; 7 8 9];
elementos_mayores_a_5 = A(A>5);
% Resultado: elementos_mayores_a_5 = [6 7 8 9]
Matlab/Octave está diseñado para trabajar de manera óptima con matrices, permitiendo realizar operaciones como la transposición, productos matriciales y otras manipulaciones comunes.
A = [1 2; 3 4];B = [5 6; 7 8];
C = A+B; % Suma de las matrices A y B
D = A*B; % Producto de las matrices A y B
X = A\B; % X es la matriz que cumple A*X=B
X = A/B; % X es la matriz que cumple X*A=B
X = A^2; % X es A*A
A = [1 2 3; 3 4 5];B = [6 7 8; 9 10 11];
C = A.*B; % Multiplicación elemento a elemento
D = A./B; % División elemento a elemento
D = A.^3; % Potencia elemento a elemento
A = [1 2 3; 3 4 5]; % A es 2x3
A' % Matriz transpuesta de A de orden 3x2
p = [3 -2 5 0 1];
% Representa 3x^4 - 2x^3 + 5x^2 + 0x + 1
Para evaluar un polinomio en un punto se utiliza el comando polyval. Para multiplicar, el comando conv y para dividir deconv.
x = 2; p1=[3 -2 5 0 1];p2=[2 3];
valor = polyval(p1, x); disp(valor) % Evalúa p(2)
pmul = conv(p1, p2);
[cociente, resto] = deconv(p1, p2);
Los programas Matlab/Octave ofrecen una amplia variedad de funciones integradas para trabajar con números, vectores, matrices, gráficos y procesamiento de datos. A continuación se muestran algunas de las más comunes, organizadas en categorías.
v = [3 4 -1]; w=[3 4 5] % Vector fila
abs(v); % Valor absoluto de cada elemento
sqrt(w); % Raíz cuadrada de cada elemento
exp(v); % Exponencial de cada elemento
log(w); % Exponencial de cada elemento<
round(w); % Redondeo al entero más cercano
floor(w); % Redondeo hacia abajo
round(w); % Redondeo hacía arriba
mod(11,4) % Residuo de la división 11 entre 4
v = [3 4 -1]; % Vector fila
% Seno, coseno, tangente de cada elemento
sin(v), cos(v), tan(v)
% Arcoseno, arcoseno, y arcotangente
asin(v), acos(v), atan(v)
% Seno, coseno, tangente hiperbólicas
asinh(v), acosh(v), atanh(v)
% Arcoseno, arcocoseno, y arcotangente
asin(v), acos(v), atan(v)
A = [1 7 3;4 8 6]; % Matriz 2x3
size(A) % Tamaño de A
length(A) % Longitud de A, mayor dimensión
trace(A) % Traza de A
sum(A) % Suma de los elementos de A
prod(A) % Producto de los elementos de A
mean(A) % Promedio de los elementos de A
max(A) % Máximo de los elementos de A
min(A) % Mínimo de los elementos de A
sort(A) % Ordena elementos de A
size(A) % Devuelve un vector con el número
% de columnas y filas de la matriz A
A = [1 2; 3 4]; % Matriz 2x2
det(A) % Determinante de A
inv(A) % Inversa de A
eig(A) % Autovalores de A
rank(A) % Rango de A
Una de las características más útiles de Matlab es su capacidad para crear gráficos de forma
rápida y sencilla. Uno de los comandos a utilizar es plot.
plot(vectorx,vectory,opciones_gráfica)
% Ejemplo: Representar $x^2$ en $[-10,10]$
>> x = -10:0.1:10; % Vector de valores x
>> y = x.^2; % Cuadrado de cada valor de x
>> plot(x, y); % Crear la gráfica
En Matlab/Octave, puedes personalizar las gráficas usando opciones como color de línea (por ejemplo, 'r' para rojo), estilo de línea ('--' para discontinua), y símbolo para representar los puntos ('o' para círculos). También puedes combinar estas opciones ('b--o' para una línea azul discontinua con círculos) y ajustar el ancho de línea ('LineWidth', 2).
Otros comandos gráficos:
figure.xlabel,
ylabel.
hold off.
grid on,
grid off.
title.subplot(m,n,p) divide la
ventana de gráficos en una cuadrícula de mxn gráficas y selecciona la subgráfica p.
% Datos de ejemplo de temperaturas a lo largo del día
horas = 0:1:23; % Horas de 0 a 23
temp_ciudad1 = [15, 14, 13, 13, 12, 12, 13, 15, 18, 22, 25, 27, 29, 30, 28, 26, 24, 22, 20, 18,
17, 16, 15, 14];
temp_ciudad2 = [10, 9, 9, 8, 8, 7, 8, 10, 13, 18, 20, 23, 25, 26, 24, 23, 21, 19, 17, 15, 13,
12, 11, 10];
% Gráfico de temperatura para la Ciudad 1
% Gráfico en azul
plot(horas, temp_ciudad1, 'b', 'LineWidth', 1.5);
hold on; % Mantiene la gráfica para añadir la segunda
% Gráfico de temperatura para la Ciudad 2
% Gráfico en rojo con línea discontinua
plot(horas, temp_ciudad2, 'r--', 'LineWidth', 1.5);
% Etiquetas y leyenda
xlabel('Hora del día');
ylabel('Temperatura (°C)');
title('Temperatura promedio durante el día en dos ciudades');
legend('Ciudad 1', 'Ciudad 2');
% Mostrar la cuadrícula para facilitar la lectura
grid on;
% Desactivar hold para futuras gráficas
hold off;
Matlab/Octave no solo son programas para realizar cálculos, también tienen un lenguaje de programación que permite crear scripts y funciones para automatizar tareas, realizar iteraciones, condicionales, y mucho más. A continuación, se presentan algunos conceptos básicos de programación con estos programas.
Funciones
Las funciones en Matlab/Octave permiten encapsular un conjunto de
operaciones dentro de un bloque de código que puede ser reutilizado. Se definen en un archivo
separado, también con la extensión .m. y con el mismo nombre que la
función. En el siguiente ejemplo se muestra la función cuadrado que
calcula, a partir de un valor de entrada, su cuadrado.
function resultado = cuadrado(x)
resultado = x^2;
end
El fichero que tiene esta función deberá llamarse cuadrado.m. La
función se invoca de la misma forma que las funciones predefinidas propias de Matlab/Octave.
a=cuadrado(2)+1
% Resultado: a=5;.
Condicionales
Matlab/Octave admite la estructura de control condicional if-else, que se usa para
ejecutar bloques de código en función de si se cumple o no una condición booleana.
En el ejemplo siguiente, si el valor de x es positivo se escribe el texto "x es positivo", en caso contrario, cuando es menor que cero, se escribe por pantalla "x es negativo" y, en otro caso, se escribe "x es cero".
if x > 0
disp('x es positivo');
elseif x <0
disp('x es negativo');
else
disp('x es cero');
end
Bucles
Matlab/Octave admite dos tipos de bucles principales: for y while.
Permiten repetir un conjunto de instrucciones varias veces.
Para salir de un bucle, se puede utilizar la instrucción break. Para
omitir el resto de las instrucciones del bucle y comenzar la siguiente iteración, se puede
utilizar la instrucción continue.
for. Repite un conjunto de instrucciones cuando varía el
índice en un conjunto de valores.
En el siguiente ejemplo se imprime por pantalla el valor de i cuando va tomando los valores 1, 2, 3,...,10.
for i = 1:10
disp(i);
end
En el siguiente ejemplo, vamos a sumar los cuadrados de los primeros 20 números naturales.
suma = 0; % Inicialización de la suma
for i=1:20
% Agrega el siguiente número a la suma
suma = suma + i^2;
end
Este cálculo se puede hacer también de la forma siguiente:
vec=1:20;
sum(vec.^2)
while. Repite un conjunto de instrucciones mientras una
condición es verdadera.
En el siguiente ejemplo, vamos a sumar los primeros números naturales hasta que la suma sea inferior a 300.
suma = 0; % Inicialización de la suma
n = 1; % Número
while suma < 300.
% Agrega el siguiente número a la suma
suma = suma + n;
% Incrementa el contador
n = n + 1;
end
En el siguiente ejemplo se calcula la suma de los números pares desde 2 a 10.
n = 2; % Número inicial
suma = 0; % Variable para almacenar la suma
% Bucle while
while n <= 10
%Sumar el valor de n a la variable suma
suma = suma +n;
n = n + 2; % Incrementar n en 2
end
suma % Mostrar el resultado
Este cálculo se puede hacer también de la forma siguiente:
sum(2:10)
fopen,
fscanf
% Lee todos los números de archivo.txt
fileID = fopen('datos.txt', 'r');
data = fscanf(fileID, '%f');
fclose(fileID);
disp(data);
fprintf
% Guarda el vector en un archivo de texto
data = [1.1, 2.2, 3.3];
fileID = fopen('output.txt', 'w');
fprintf(fileID, '%f\n', data);
fclose(fileID);
save, load
% Guarda las variables a y b en archivo.mat
a = 5;
b = [1, 2, 3];
save('variables.mat', 'a', 'b');
load('variables.mat');
Se presenta a continuación un cuestionario básico sobre Matlab/Octave.
Con el siguiente interactivo puedes realizar una autoevaluación por niveles.
Multiplo que, dados dos números \( k \) y \( h
\), devuelva:
ContarPares que reciba un vector como argumento y
devuelva tantos sus componentes pares como el número de ellas.
SumaCifras(325) debería devolver 10.
MCD(24, 15) debería devolver 3. Nota: El MCD(a,b) con a>b es b si a es múltiplo de b. En otro caso, MCD(a,b)=MCD(b,r) siendo r el resto de dividir a entre b.
SumaCuadrados(5, 3, 4) debería
devolver 41.
SuperaMedia que reciba un vector y devuelva otro
vector con los elementos que sean mayores o iguales a la media del vector de entrada.
Los aspectos clave en el estudio de los métodos numéricos son:
El objetivo de este tema es destacar la importancia de la aritmética computacional en el contexto de los métodos numéricos. Se analizará cómo los computadores representan los números y realizan cálculos, destacando que, debido a limitaciones inherentes, como el almacenamiento en binario y la precisión finita, los resultados que se realicen con el ordenador pueden contener errores.
En particular, se observará que:
A continuación, presentaremos algunos ejemplos numéricos que ilustran este tipo de problemas y los errores que pueden generar.
suma=0;n=1000;
for k=1:n
suma=suma+0.1;
end
format longEng;suma
El resultado debería ser $100$ pero se obtiene un número próximo: $99.9999999999986e+000$
f=@(x) (1-x).^6
%Consideramos el polinomio resultado de calcular la potencia
%syms x ;expand((1-x).^6)
g=@(x) x.^6-6*x.^5+15* x.^4-20*x.^3+15*x.^2-6*x+1
x1=0.996:0.0001:1.005; format longEng;
%Comparamos resultados
[x1' f(x1)' g(x1)']
plot(x1,f(x1),x1,g(x1))
En la figura se muestra la gráfica de la función $f(x)$ y la de la función $g(x)$ del código anterior.
¿Por qué ocurre esto? La razón es que los ordenadores usan un formato (punto flotante binario) que no permite representar de forma precisa números como por ejemplo 0.1, 0.2 o 0.3. Dado que los números están, en general, representados de manera inexacta, las operaciones aritméticas que se realicen también serán inexactas.
Un sistema de numeración consiste en un conjunto de símbolos que permiten representar cantidades numéricas, así como reglas para realizar operaciones con dichas cantidades. En este apartado, nos centraremos en el sistema decimal (base 10) y en el binario (base 2), fundamentales en el ámbito del cálculo numérico. Mientras que el sistema decimal es el que utilizamos en la vida cotidiana, el binario es esencial en computación y en el procesamiento numérico, ya que los ordenadores representan y manipulan los datos a través de este sistema.
En el sistema decimal los números se representan utilizando 10 dígitos: 0,1,2,3,4,5,6,7,8,9. Además, este sistema de numeración es posicional en el que las cantidades se representan utilizando como base aritmética las potencias del número diez. Así, si $x \in \mathbb{R}^{+}$
$$x = \sum\limits_{ - \infty }^n {{a_k} \cdot {{10}^k}} \,\,\,\,\, 0 \le {a_k} \le 9$$Si consideramos una base $b$, un número real $x$ positivo se puede escribir de la forma
$$x = \sum\limits_{ - \infty }^n {{a_k} \cdot {b^k}} \,\,\,\,\, 0 \le {a_k} \le b - 1$$ o también $$x = {a_n}{a_{n - 1}}...{a_0}.{a_{ - 1}}{a_{ - 2}}{a_{ - 3}}...$$El Standard de IEEE (Institute of Electrical and Electronics Engineers) propone el uso de la base $b=2$ en el almacenamiento y la aritmética sobre el computador. En este sistema de numeración los números se representan utilizando solamente las cifras cero y uno (0 y 1).
$$x = \sum\limits_{ - \infty }^n {{a_k} \cdot {2^k}} \,\,\,\,\, 0 \le {a_k} \le 1$$Para convertir un número decimal a binario, basta dividir el número y sus cocientes entre 2, acumulando los restos que conformarán el número en la nueva base. La lista de ceros y unos leídos de abajo a arriba darán la expresión en binario.
Se tiene que,
$$ \begin{aligned}{1 \over {10}} &= {1 \over {{2^4}}} + {1 \over {{2^5}}} + {0 \over {{2^6}}} +
{0 \over {{2^7}}} + {1 \over {{2^8}}} + {1 \over {{2^9}}} + {0 \over {{2^6}}} + {0 \over
{{2^7}}} + ... \\ &= 0.0001{\widehat {1001}_2} \end{aligned}$$
Para convertir un número binario a decimal, se multiplica cada dígito binario por la potencia de 2 correspondiente a su posición y luego se suman todos los valores resultantes.
En el siguiente ejemplo se muestra un número binario con parte fraccionaria y su conversión a decimal.
Se tiene que: $$\begin{aligned} x = 1·2^1 + 0·2^0 + 1·2^{-1}+ 1·2^{-2} + \\ +1·2^{-3}+ 0·2^{-4} + 0·2^{-5} + 1·2^{-6} \\ = 1+1/2+1/4+1/8+1/64=1.890625\end{aligned}$$
Un número real en base 10 se puede escribir de forma que la primera cifra no nula aparecezca inmediatamente a la izquierda del punto decimal.
$$x = {\left( { - 1} \right)^s}× \,{10^e} × \,a.m$$donde
Se denomina expresión normalizada en
Como los ordenadores utilizan aritmética binaria, en base 2, cualquier número $x$ se puede escribir normalizado en coma flotante de la forma
$$x = {\left( { - 1} \right)^s} × \,{2^e} × \,1.m$$Por lo tanto, el número $x$ se puede almacenar guardando el signo $s$, el exponente $e$ y la mantisa $m$. Como en base 2, el primer dígito que precede a la mantisa es siempre 1, no es necesario almacenarlo. De esta forma se optimiza la representación mediante el llamado bit oculto.
En el ordenador los números reales se representan según la norma IEEE 754, establecida en 1985
por el Institute of Electrical and Electronics Engineers Esta norma establece dos formatos,
Para representar un número real en precisión simple en el ordenador se utiliza:
Teniendo en cuenta que se usan 8 bits para el exponente, sus valores posibles variarían desde 0 hasta 255:
$$\underbrace {00000000}_{8\, bits}\,\,\,\,\,\,\,\,\text{hasta}\,\,\,\,\,\,\underbrace {11111111}_{8 \,bits} \equiv 1 + 2 + ... + {2^{7}} = {2^{8}} - 1 = 255$$Con objeto de poder considerar las potencias de 2 negativas se hace un desplazamiento del exponente para que represente valores desde -127 a 128. Así, en precisión simple, un exponente real de 0 se almacena
como 127, un exponente de -1 como 126, y un exponente de +1 como 128.
En general, la relación entre el exponente real y el almacenado es la siguiente:
$$e_{almacenado}=e_{real}+127$$
Los valores extremos del exponente, esto es 0 y 255, se reservan para números especiales (cero e infinito).
Se tiene que
$$3.125_{10}=11.001_{2}=1.1001 × 2^{1}$$Por lo tanto, $m=1001$ que ajustada a 23 bits es $1001\underbrace{0...0}_{19 \,\,bits}$.
Para el exponente, considerando el sesgo, se codificará $e=127+1=128=2^{7}$. Este número en binario es $10000000$.
En definitiva, la representacion en precisión simple del número decimal 3.125 es la siguiente:
$$0 10000000 1001\underbrace{0...0}_{19 \,\,bits}$$
En Matlab la precisión que usaremos para los cálculos será
Teniendo en cuenta que el exponente se almacena en 11 bits, los valores que se pueden considerar para el exponente son desde 0 hasta 2047:
$$\underbrace {0....0}_{11 \,bits}\,\,\,\,\,\,\,\,\text{hasta}\,\,\,\,\,\,\underbrace {1....1}_{11 \,bits} \equiv 1 + 2 + ... + {2^{10}} = {2^{11}} - 1 = 2047$$
Para poder considerar las potencias de 2 negativas se toma un desplazamiento o
Se consideran las siguientes representaciones especiales:
Respecto al exponente, el valor a almacenar es $e=1023+1=1024=2^{10}$, que en binario es $10000000000$.
El número se almacenará de la forma: $$ \underbrace{0}_{signo} \underbrace{10000000000}_{exponente}\,\,\underbrace{10010...0}_{mantisa, \,48\,\,últimos \,\,bits\,\,0}$$
Dado que la representación de números reales bajo estos formatos es aproximada, hay dos conceptos importantes en la aritmética en punto flotante:
realmax en Matlab
realmin en Matlab
Si consideramos que el bit escondido es cero, podemos obtener números menores que el mínimo número normalizado mediante la siguiente representación: $$x = \left( { - 1} \right)^s × 2^{- 1022} × \,\,0.f$$
Estos números se dicen que están
El mínimo número desnormalizado positivo ocurre cuando la mantisa tiene un 1 en el bit menos significativo $m=2^{-52}$. El exponente real será
$$ {2^{ - 1022}} \cdot {2^{ - 52}} = {2^{ - 1074}}$$Todos los números menores que este valor se consideran 0.
>>2^(-1074)
ans=
4.9407e-324
>>2^(-1075}
ans=
0
El número más pequeño representable después del 1 utilizando representación binaria en punto flotante con doble precisión es:
$$1 = \left( { + 1} \right)× {{\rm{2}}^0}× 1.\underbrace {00...0}_{52}$$ $$1 + \varepsilon = \left( { + 1} \right)× {{\rm{2}}^0}× 1.\underbrace {00...01}_{52}$$
El valor de $\epsilon$ se llama eps: $\varepsilon_m = {2^{ - 52}}$.
El número más pequeño $\varepsilon$ de forma que $1-\varepsilon \ne 1$ es ${2^{ - 53}}$. Por tanto, las distancias entre el número que antecede a 1 y el que le precede son respectivamente $\varepsilon_m/2$ y $\varepsilon_m$.
eps(x).
Para medir los errores, se introducen los conceptos de error absoluto y error relativo. Si $\widetilde p$ es una aproximación de $p$, se define:
El error absoluto indica cuánto se ha desviado una medición del valor real, en las mismas unidades. El error relativo expresa el error en términos fraccionarios respecto al valor real y multiplicado por 100 se obtendría en porcentaje, lo que permite comparar errores en diferentes escalas.
En general, en los métodos numéricos no se conoce el valor verdadero, por lo tanto tampoco se conocerá el error real. Lo frecuente es tener una estimación del error a partir de una cota positiva del error máximo.
En este ejemplo se muestra que si bien el error absoluto es distinto, el error relativo es el mismo. El error relativo de 3% indica que la medición tiene un error de aproximadamente el 3% respecto al valor real.
En la tabla siguiente se recogen los errores absoluto y relativo para cada estimación.
| Medición | Valor estimado | Valor real | Error absoluto | Error relativo |
|---|---|---|---|---|
| Velocidad de la luz (Foucault) | 298000 | 299776 | 1776 | 5.920661e-03 |
| Constante de Planck | 6.41e-27 | 6.626e-27 | 2.13e-28 | 3.216246e-02 |
Para determinar cuál es la medición más precisa, observamos el error relativo:
En el sistema decimal, existen dos formas de redondear un número a $t$ decimales. En el proceso
denominado
El número 3.246 se representaría como 3.24 si se usa truncamiento o bien 3.25 si se usa redondeo. El truncamiento elimina todas las cifras que no pueden ser representadas por la máquina, en este caso a partir del tercer decimal. El redondeo, toma la primera cifra no representable, en este ejemplo, sería el 6 y le suma una unidad al último dígito si es mayor o igual a 5 o lo trunca si es menor que 5.
Si consideramos el número
$$x = {\left( { - 1} \right)^s}× {2^{{e_x}}} × 1.f\,\,\,\,\,\,\,\,\,\,\,\,\,\,-1023 < {e_x} < 1024$$
y denotamos por
entonces el épsilon de la máquina es una cota superior del error de redondeo en el sistema punto flotante.
Si $\otimes$ es una operáción aritmética como una suma, resta, multiplicación o división, el cálculo sigue el siguiente proceso:
La representación de la operación es entonces
$$fl\left( {x \otimes y} \right) = fl\left( {fl\left( x \right) \otimes fl\left( y \right)} \right)$$ En consecuencia, el orden en el que se realizan las operaciones puede modificar el resultado final. Con el siguiente ejemplo se muestra que puede no cumplirse la propiedad asociativa.
x=10^20;y=-10^20;z=1;
x+(y+z)
(x+y)+z
x=1- 3*(4/3 - 1)
x =
2.2204e-16
El valor exacto es 112.7. El calculo sería de la siguiente manera
$$fl(1025)=1.02 \cdot 10^3 \,\,\,\,\,\,fl(-912.3)=-9.12\cdot 10^2$$ $$fl(1025-912.3)=fl(112.7)=1.12 \cdot 10^2$$ $$ fl(fl(x)-fl(y))=fl(90)=9.0 \cdot 10^1$$Además, cuando se realiza una cadena de operaciones, los errores de redondeo pueden acumularse de modo que resten validez al resultado. Para controlar la propagación de errores de redondeo, es necesario buscar algoritmos con el número mínimo de operaciones convenientemente dispuestas para evitar estos problemas y a la vez el tiempo de ejecucion sea aceptable.
Para un análisis detallado sobre los errores en cálculos numéricos, como son los errores de redondeo,
truncamiento y su propagación, se recomienda consultar los libros
Este error se produce cuando se calcula la diferencia de dos números bastante próximos, y el redondeo del ordenador iguala esa diferencia a cero.
Nuestro ordenador establecerá en el primer caso que $x-y=0$ cuando en realidad $x - y \approx 5,5 \cdot {10^{ - 17}}$.
% Apartado a)
x=1+2^-54;y=1;x-y
% Apartado b)
x=2^60+2^6;y=2^60;x-y
f=@(x) (exp(x)-1).^2+(1./sqrt(1+x.^2)-1).^2; k=-1:-1:-17;h=10.^(k);
v=(f(1+h)-f(1))./h;
format longEng
disp([h' v'])
Si se considera el cálculo de la segunda forma, tomando el inverso de $e^{40}$, se evita este problema de cancelación ya que todos los términos son positivos.
x=40; %Inverso del número a calcular e^-40=1/e^40
suma=1;p=1;
for k=1:99
p=p*(x/k);
suma=suma+p;
end
1/suma
x=-40; %Sin considerar el inverso
suma=1;p=1;
for k=1:99
p=p*(x/k);
suma=suma+p;
end
suma
%Valor que da Matlab
exp(-40) % 4.248354255291589e-018
En los siguientes ejemplos, se muestra como una reformulación del cálculo evita en ciertas situaciones el problema de cancelación.
%Se considera la primera expresión como resta de fracciones
f=@(x) 1./(1+2*x)-(1-x)./(1+x);
%Se realiza la suma para evitar la cancelación
g=@(x) 2*x.^2./((1+2*x).*(1+x));
format longEng
%Se considera un vector de valores próximos a cero
k=-1:-1:-17;
h=10.^k;
%Se evalúan ambas expresiones
f1=f(h);
g1=g(h);
%Se imprimen los valores de h y los valores
%obtenidos para comparar
disp([h' f1' g1'])
f=@(x) (1-cos(x))./x;
%Se multiplica y divide por (1+cos(x))
g=@(x) sin(x).^2./(x.*(1+cos(x)));
k=-1:-1:-17;
h=10.^k;
f1=f(h);
g1=g(h);
format longEng
disp([h' f1' g1'])
f=@(x) sqrt(x+1./x)-sqrt(x-1./x);
%Se multiplica y divide por el conjugado
g=@(x) 2*x./((sqrt(x+1./x)+sqrt(x-1./x)));
k=-1:-1:-25;
h=10.^k;
f1=f(1+h);
g1=g(1+h);
format longEng
disp([h' f1' g1'])
Con frecuencia, una operación aritmética con dos números válidos da como resultado un número tan
grande o tan pequeño que la computadora no puede manejarlo; como consecuencia, se produce un
realmax, por dos para ilustrar el concepto de
desbordamiento. Divide después el valor mínimo positivo normalizado, realmin, por 2 para observar que se
representa por cero.
% Obtenemos el valor máximo representable
maxVal = realmax;
% Intentamos exceder el valor máximo
overflow = maxVal * 2;
disp(['Desbordamiento resulta en: ', num2str(overflow)]);
% Obtenemos el valor mínimo positivo normalizado representable
minVal = realmin;
% Intentamos obtener un valor menor al mínimo
underflow = minVal * 1e-100;
disp(['Subdesbordamiento resulta en: ', num2str(underflow)]);
En el siguiente ejemplo veremos cómo los errores de redondeo pueden conducir a resultados incorrectos si el algoritmo no es el adecuado. Este fenómeno, conocido como inestabilidad numérica, puede evitarse seleccionando algoritmos más apropiados.
Se calcula directamente el valor de $I_{30}$ e $I_{31}$ de forma iterativa y se obtienen valores negativos. Esto es imposible porque el integrando es positivo y por tanto las integrales deben ser positivas.
El valor de $I_{30}$ que se otiene es -0.357484694711417 y para el valor de $I_{31}$ se obtiene -2.312387998625434
Despejando $I_n$ en función de $I_{n+2}$ se tiene:
\[\frac{{\left( {n + 2} \right)\left( {n + 1} \right)}}{\pi }{I_n} = 1 - \pi {I_{n + 2}}\] \[{I_n} = \frac{\pi }{{\left( {n + 2} \right)\left( {n + 1} \right)}}\left( {1 - \pi {I_{n + 2}}} \right)\] \[{I_n} = \frac{{{\pi ^2}}}{{\left( {n + 2} \right)\left( {n + 1} \right)}}\left( {\frac{1}{\pi } - {I_{n + 2}}} \right)\] Si se toma límites cuando $n$ tiene a infinito, se puede ver que la sucesión $I_n$ tiende a cero. $$\left| {\int\limits_0^1 {{x^n}sen\left( {\pi x} \right)dx\,} } \right| \le \int\limits_0^1 {{x^n}dx\,} = \left. {{{{x^{n + 1}}} \over {n + 1}}} \right|_{x = 0}^{x?1} = {1 \over {n + 1}}$$Por ello, para utilizar esta expresión hacia atrás, se debe partir de dos valores $I_{n+2}$ e $I_{n+1}$ próximos a cero para un valor de $n$ alto. En el código se ha considerado $I_{56}$ e $I_{55}$ como 0.
El código para el cálculo directo es el siguiente:
a=2/pi; v(1)=1/pi;v(2)=1/pi-2*a/pi^2;
format longEng
for k=3:31;
v(k)=1/pi-k*(k-1)*v(k-2)/pi^2;
end
disp([[1:31]' v'])
Para el cálculo inverso, las instrucciones son las siguientes.
a=56;
w(56)=0;w(55)=0;
for k=a-2:-1:1;
w(k)=pi^2/((k+2)*(k+1))*(1/pi-w(k+2));
end
disp([[1:56]' w'])
El valor de $I_{30}$, utilizando este cálculo, es $0.003139287076703$, y el valor de $I_{31}$ es $0.002950500704051$.
Se pide calcular de forma recurrente $y_{20}$ con la expresión anterior. Observa que se obtiene un valor negativo, lo cual es absurdo porque el integrando es positivo.
Para evitar este problema, considera después la fórmula de recurrencia en sentido descendente para obtener el valor de $y_{20}$
\[{y_{n - 1}} = \frac{1}{{5n}} - \frac{{{y_n}}}{5}\]Considerando la recurrencia ${y_n} = \frac{1}{n} - 5{y_{n - 1}}$, calculamos $y_{20}$.
clear all
y(1)=1-5*(log(6)-log(5));
format longEng
for k=2:20;
y(k)=1/k-5*y(k-1);
end
disp([[1:20]' y'])
Tomando límites en ${y_n} = \frac{1}{n} - 5{y_{n - 1}}$, vemos que la sucesión $y_n$ converge a 0. Por lo tanto, considerando $y_{60}=0$ y tomando la ley de recurrencia ${y_{n - 1}} = \frac{1}{{5n}} - \frac{{{y_n}}}{5}$, calculamos de nuevo el valor de $y_{20}$.
y(60)=0;
format longEng
for k=60:-1:2;
y(k-1)=1/(5*k)-y(k)/5;
end
disp([(1:60)' y'])
%y(20) es 7.997523028232164e-03
Se clasifican y analizan errores como los de redondeo y truncamiento, destacando la importancia de desarrollar algoritmos numéricos confiables para minimizar estos problemas en aplicaciones prácticas.
Los aspectos clave en el estudio de los métodos numéricos son:
Entender estas limitaciones es esencial para abordar problemas numéricos con precisión.
El propósito de los métodos numéricos para resolver ecuaciones no lineales es encontrar raíces aproximadas, ya que la mayoría de
las ecuaciones no tienen soluciones exactas o estas no son computacionalmente viables. Mediante
métodos iterativos se va consiguiendo aproxima la solución dentro de un margen de error considerado
Entre las aplicaciones más importantes de los métodos numéricos está la resolución de ecuaciones no lineales, una tarea esencial en disciplinas como física, ingeniería, economía y muchas otras áreas de la ciencia. Estas ecuaciones, donde la variable aparece de forma no lineal son clave para modelar fenómenos reales, desde el diseño de estructuras hasta la dinámica de sistemas físicos. Sin embargo, muchas veces estas ecuaciones no tienen soluciones exactas en términos algebraicos, por lo que se buscan soluciones aproximadas.
Imaginemos por ejemplo la siguiente ecuación obtenida a partir de la segunda ley de Newton para calcular la velocidad de un paracaidista,
$$v = {{g \cdot m} \over c}\left( {1 - {e^{ - {c \over m}t}}} \right)$$donde la velocidad $v$ depende de la variable independiente tiempo, $t$, $g$ es la constante de gravitación, $c$ el coeficiente de resistencia y $m$ es la masa del paracaidista.
Si quisieramos obtener el coeficiente de resistencia del paracaidista con una masa dada, para alcanzar una velocidad determinada en un periodo establecido. ¿cómo se obtendría este valor? Este valor sería el valor $c$ que hace cero a la siguiente función
$$f\left( c \right) = {{g \cdot m} \over c}\left( {1 - {e^{ - {c \over m}t}}} \right) - v$$esto es, calcular el valor de $c$ de forma que $f(c)=0$.
Cuando como en este caso, no es posible encontrar el cero de la función o pudiendo encontrarlo su cálculo es complicado, se utilizan algoritmos que determinan, en un número finito de pasos, un valor aproximado de la raíz buscada. Estos métodos de aproximación generalmente son iterativos y nos conducen a una solución a través de una sucesión $x_0, x_1, ..., x_k,...$ que converge al valor $\alpha$ que cumple $f(\alpha)=0$.
A lo largo de la historia, la resolución de ecuaciones no lineales ha evolucionado desde métodos manuales hasta complejos algoritmos implementados en computadoras modernas. En la antigüedad, los babilonios y griegos ya aplicaban principios numéricos para resolver problemas matemáticos. Durante el Renacimiento y la Ilustración, métodos como la Regla Falsa o la técnica de Newton-Raphson marcaron grandes avances.
Así, el método de bisección, que tiene sus raíces en la antigüedad, utiliza un enfoque simple: dividir un intervalo y localizar la raíz mediante evaluaciones intermedias. Otros
métodos, como el de Newton-Raphson, desarrollado en el siglo XVII por Isaac Newton
Isaac
Newton (1642-1727) fue un físico, matemático y astrónomo inglés, considerado uno de los
científicos más influyentes de la historia. Desarrolló las leyes del movimiento y la
gravitación universal, así como importantes contribuciones al cálculo, óptica y mecánica. Su
obra más destacada, Philosophiæ Naturalis Principia Mathematica, sentó las bases de la
física clásica.
y Joseph Raphson
Joseph
Raphson (1648-1715) fue un matemático inglés conocido por su trabajo en análisis numérico.
Publicó el método de Newton-Raphson en su obra Analysis Aequationum Universalis (1690), una
técnica iterativa para encontrar raíces de ecuaciones no lineales. Su versión mejorada del
método de Newton fue más algebraica y general, contribuyendo al desarrollo de la matemática
aplicada., aprovechan herramientas del cálculo diferencial
utilizando la pendiente de las tangentes de la función.
En el siglo XX, con el desarrollo de la computación, matemáticos como Alan Turing y John von
Neumann
John
von Neumann (1903-1957) fue un matemático, físico y pionero de la computación
húngaro-estadounidense. Contribuyó al desarrollo de la teoría de juegos, la mecánica
cuántica y los métodos numéricos en computación científica. Fue clave en la creación de las
primeras computadoras y en la formulación de la arquitectura de Von Neumann, base de los
sistemas computacionales modernos.
sentaron las bases para implementar estos métodos en máquinas, haciendo posible resolver
problemas que antes eran inabordables.
En el siguiente video se da una visión rápida de los métodos más importantes que se abordan en este capítulo.
En este capítulo abordaremos los principales métodos numéricos para la resolución de ecuaciones no lineales. Comenzaremos con el método de bisección, una técnica basada en el teorema de Bolzano, que garantiza la existencia de raíces en un intervalo y ofrece una convergencia segura aunque lenta. Luego, estudiaremos el método del punto fijo, que transforma la ecuación en una forma iterativa y permite analizar la convergencia local y global de la solución. Posteriormente, exploraremos el método de Newton-Raphson, que emplea la derivada de la función para obtener una convergencia rápida, aunque con ciertas limitaciones respecto a su punto de inicio. Finalmente, analizaremos el método de la secante, que mejora la eficiencia de Newton al evitar el cálculo de derivadas. A lo largo del capítulo, presentaremos ejemplos prácticos para ilustrar el comportamiento de cada método, comparando sus ventajas y limitaciones en diferentes escenarios.
DEFINICIÓN: A cualquier solución de la ecuación $f(x)=0$ se le llama
Se puede demostrar que si $f$ es $C^m(I)$, esto es, derivable hasta el orden $m$ con derivadas continuas, y $\alpha$ es un punto de $I$, la función tiene un cero de multiplicidad $m$ en $\alpha$ si
$$f(\alpha)=f'(\alpha)=...=f^{(m-1}(\alpha)=0, f^{(m}(\alpha) \ne 0$$En un primer paso para obtener las raíces de una función, se busca un intervalo donde se encuentra una o varias soluciones de la ecuación. Posteriormente, se buscan intervalos que contengan una sola raíz. El siguiente teorema ayuda en esta localización y separación de raíces.
TEOREMA DE BOLZANO: Si $f$ es una función continua en $[a,b]$ de manera que $f(a)f(b) < 0$, entonces existe al menos un cero de la función en el intervalo $(a,b)$.
El hecho de que la función $f$ tenga un cero, no significa que sea único. Si se añaden condiciones adicionales, como por ejemplo la monotonía estrica de $f$ en un intervalo $I$, sí se podría garantizar que la solución de $f(x)=0$ es única en dicho intervalo.
Se recuerda que en el caso de que la función sea derivable, si $f'(x)> 0$ en todo punto de $I$ la función es estrictamente creciente en $I$. Análogamente, si $f'(x)< 0$ en $I$ entonces es estrictamente decreciente. Si la función es estrictamente creciente o estrictamente decreciente en un intervalo donde hay cambio de signo en los extremos, la función tiene un único cero en dicho intervalo.
Este proceso depende tanto de la función $G$ como del punto inicial $x_0$. Los algorítmos a utilizar deberán dar respuesta a cuestiones como las siguientes:
Se dice que un algoritmo
Cuando la solución converge hacia la solución únicamente si el punto $x_0$ pertenece a un cierto
entorno, la convergencia se dice que es
No siempre es mejor considerar un algoritmo que converja globalmente, a veces, un algoritmo localmente convergente, puede tener una "velocidad de convergencia" mayor si bien este tipo de algoritmos plantean la dificultad de cómo elegir el intervalo adecuado donde considerar $x_0$.
El método de aproximaciones sucesivas o punto fijo, que consiste en partir de un punto $x_0$ y calcular los demás términos de la sucesión $x_n$ de la forma $x_{n+1}=g(x_n)$, para calcular la raíz de $x=g(x)$, tiene una convergencia más lenta que el método de Newton que converge localmente.
Se distinguen distintas velocidades de convergencia que pasamos a definir.
Dada una sucesión de números reales, $\left\{ {x _k} \right\}_{k = 1}^\infty $, convergente
hacia un punto $\overline x$, se dice que el
$$\left| {{x_{k + 1}} - \overline x }
\right| \le {C_p}{\left| {{x_k} - \overline x \,} \right|^p}$$ para todo $ k \ge k_0$ siendo
$k_0$ entero.
En particular,
Se tiene que
Se dice que la sucesión
Observa que si la sucesión tiene una velocidad de convergencia $p$ mayor que 1 entonces converge superlinealmente. Bastaría tomar $${\varepsilon _k} = {C_p}{\left| {{x_k} - \bar x{\mkern 1mu} } \right|^{p - 1}}\,\,{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} \forall k \ge 0$$
El siguiente resultado permite determinar un test de parada de las iteraciones cuando la convergencia es superlineal.
Si $x_k$ converge superlinealmente a ${\bar x}$ entonces $$\mathop {\lim }\limits_{k \to \infty } {{\left| {{x_{k + 1}} - {x_k}{\mkern 1mu} } \right|} \over {\left| {{x_k} - \bar x{\mkern 1mu} } \right|}} = 1$$
En este caso, una vez establecida la precisión que se quiere alcanzar,$\varepsilon > 0$, puede considerarse como test de parada:
$$\left| {{x_{k + 1}} - {x_k}} \right| < \varepsilon$$se compara así dos iteraciones consecutivas y se mira la coincidencia de dígitos que se considere (fijado por el valor de $\epsilon$).
El límite de la sucesión se obtiene resolviendo la ecuación
\[x = \frac{x}{{1.1}} + 3 \Rightarrow x\left( {1 - \frac{1}{{1.1}}} \right) = 3 \Rightarrow x = 33\] Analizamos en la siguiente tabla cómo los valores de la sucesión están próximos y, sin embargo, están lejos del valor límite.| k | xk | |xk+1 - xk| |
|---|---|---|
| 0 | 100.0000 | - |
| 1 | 93.0909 | 6.9091 |
| .. | .. | .. |
| 7 | 69.5798 | 2.1570 |
| 8 | 68.9844 | 0.5954 |
| 9 | 68.0767 | 0.9077 |
| 10 | 67.7970 | 0.2797 |
| 11 | 67.6337 | 0.1633 |
| 12 | 67.5398 | 0.0939 |
| 13 | 67.4907 | 0.0491 |
| 14 | 67.4633 | 0.0274 |
| 15 | 67.4476 | 0.0157 |
| 16 | 67.4387 | 0.0089 |
| 17 | 67.4330 | 0.0057 |
| 18 | 67.4295 | 0.0035 |
Aplicando el teorema de Bolzano, una forma sencilla de aproximar el cero de una función continua es dividir el intervalo en dos y elegir el subintervalo en el que se sigue produciendo el cambio de signo. Repitiendo este proceso se puede llegar a conseguir una aproximación de la solución de $f(x) = 0$ tan precisa como se desee.
| $f(x)$ | $a$ | $b$ |
|---|---|---|
| $cos(x)-x$ | $0$ | $2$ |
| $sen(x)$ | $-0.5$ | $\pi/2$ |
| $x^2-sen(x)$ | 0.5 | 2 |
En primer lugar, veremos que cada función tiene un único cero en el intervalo indicado.
Para calcular los valores de las tres primeras iteraciones del método de la bisección, se puede utilizar la herramienta interactiva siguiente introduciendo la función y los valores de $a$ y $b$ para cada caso. En el interactivo debes escribir el seno de x como sin(x).
Biseccion(a,b,tol,maxiter).
Calcula con esta función los ceros de la función $f(x)=x^3-6x^2+11x-6$ en el intervalo $[0, 1.5]$ tomando como tolerancia el valor $0.001$ y el valor de $30$ como máximo de iteraciones. En este ejemplo, las raíces se pueden obtener de forma exacta, son 1, 2 y 3, por lo que con este ejemplo únicamente se busca mostrar el método.
Biseccion2(a,b,tol,maxiter). Para ello, se debe pulsar en el
botón Para calcular los ceros en cada intervalo bastaría escribir:
Biseccion(0,1.5,10^-3,30) %Solución: 1.0000
Biseccion(1.5,2.5,10^-3,30) %Solución: 2.0000
Biseccion(2.5,3.5,10^-3,30) %Solución: 3.0000
Biseccion(fun,a,b,tol,maxiter).
Calcula con esta función el cero de $f(x)=x^3-6x^2+11x-6$ en el intervalo $[0,1.5]$, con una tolerancia $10^{-6}$ y un máximo de 25.
La función pedida podría ser:
function raiz = Biseccion2(f, a, b, tol, maxIter)
raiz='';
for i = 1:maxIter
c = (a + b) / 2;
fc = f(c);
fa = f(a);
if fa * fc < 0
b = c;
else
a = c;
end
if abs(f(c)) < tol
raiz = c;
break;
end
end
end
Para calcular el cero de $f(x)=x^3-6x^2+11x-6$ en el intervalo $[0, 1.5]$, con una tolerancia $10^{-6}$ y un máximo de 25 iteraciones, se deberá escribir el código siguiente.
f=@(x) x.^3-6*x.^2+11*x-6
Biseccion2(f,0,1.5,10^-6,25)
%Solución: 1.0000
Los siguientes resultados teóricos dan condiciones suficientes para que haya un punto fijo en un intervalo y para que la convergencia lo sea a un único punto fijo.
Sea $g$ continua en $[a,b]$ que aplica $[a,b]$ en $[a,b]$, entonces $g$ tiene un punto fijo en $[a,b]$.
Sea $g$ continua en $[a,b]$ con $g'$ en $(a,b)$ cumpliendo $$\left| {g'\left( x \right)} \right| \le L < 1\,\,\,\,\,\,\,\,para\,\,\,todo\,\,x \in \left( {a,b} \right)$$ Si $x_0$ es un punto cualquiera en $(a,b)$ entonces la solución $${x_{k + 1}} = g\left( {{x_k}} \right)\,\,\,\,\,\,k \ge 0$$ converge al único punto fijo en $[a,b]$.
Además, se verifica
$$\left| {{x_{k + 1}} - \overline x} \right| \le {{{L^k}\left| {{x_1} - {x_o}} \right|} \over {1 - L}}\,\,\,\,\,\,k = 1,2,...$$Si al utilizar el método del punto fijo $x_{k+1}$ y $x_k$ coinciden dentro de la exactitud especificada $\epsilon$, no significa que $\overline x \approx {x_{k+1}}$ con la misma exactitud. En general, esta afirmación no es correcta. Si $g'(x)$ está próxima a la unidad, entonces la cantidad $\left| {{x_k} - \overline x} \right|$ puede ser grande, aún cuando $\left| {{x_{k+1}} - x_k} \right|$ sea extremadamente pequeña.
g=@(x) x-x.^2/10; h=@(x) x;
% Gráfica de la función
fplot(g, [-2, 2]); hold on;
fplot(h, [-2, 2], 'r--'); % Línea y = x
xlabel('x');
ylabel('g(x)');
title('Método del Punto Fijo para g(x) = x - x^2/10');
grid on;
legend('g(x)', 'y=x');
Iterando, se ve que la convergencia es muy lenta, después de 100 iteraciones, el valor es $0.0828$.
x(1)=0.5;maxIter=100;
for k=2:maxIter
x(k)=g(x(k-1));
end
x(maxIter)
plot(1:maxIter,x,'o');
g=@(x) cos(x); xk=0.5;
format long
error=1;iter=0;
while error>10^-16
xk1=g(xk);
error=abs(xk-xk1);
xk=xk1;
iter=iter+1;
end
disp('Iteraciones')
iter
disp('Raiz')
xk
Con el siguiente código se pueden realizar un máximo de iteraciones indicado en maxIter para ver si se consigue obtener la precisión dada.
% Se representa la función g y la recta y=x
g=@(x) exp(-x.^2);x=0:0.01:1;
plot(x,g1(x),x,x)
% Número fijo de iteraciones
xk=0.5;tol=10^-10;maxIter=100;format long
for k=1:maxIter
xk1=g(xk);
if abs(xk1-xk)<tol
xk=xk1; break
end
xk=xk1;
end
[xk g(xk) k]
Se puede ver que la función cumple la condición suficiente para asegurar la convergencia, basta calcular el máximo de $g'$ para ver que si $x \in \left[ {0,1} \right]$,
$$\left| {g'\left( x \right)} \right| = \left| { - 2x\,{e^{ - {x^2}}}} \right| < 1$$
Con el siguiente código se pueden obtener las 5 primeras iteraciones. El valor de la derivada en $0.5$ no es menor que 1 en valor absoluto por lo que no se tendría garantizado la convergencia.
h=@(x) (-log(x)).^0.5;
%Derivada
h1=@(x) -0.5*((-log(x)).^-0.5)./x;
h1(0.5) %-1.201122408786450
%Iteraciones
xk=0.5;
for k=1:5
xk=h(xk)
end
Acotamos en primer lugar los ceros de la función y vemos qué funciones se pueden considerar para aplicar el punto fijo.
Se tiene que $p'(x)=3x^2+6x=3x(x+2)$ tiene dos ceros en el punto 0 y en -2. En el punto $(0,-1)$ hay un mínimo y en el punto $(-2,3)$ un máximo. Luego $p$ tiene tres raíces reales. Aislando estas tres raíces se tiene que están situadas en los intervalos [-3,-2], [-1,0], [0,1].
En la escena interactiva siguiente, se muestra paso a paso cómo generar el código de una
función: puntoFijo(g,x0, tol, maxiter) que tenga como argumentos la función, el punto inicial, la
tolerancia y el número máximo de iteraciones. Para calcular los tres ceros se puede escribir por ejemplo el siguiente código con una tolerancia de 10^-6 y considerando 100 iteraciones como máximo.
%Raiz en el intervalo [-3,-2]
g=@(x) x.^3-3;tol=10^-6;
puntoFijo(g,-2.5, tol, 100)
%Raíz en el intervalo [-1,0]
g=@(x) -sqrt((x+3)^-1); puntoFijo(g,-1.5, tol, 100)
%Raíz en el intervalo [0,1]
g=@(x) -sqrt((x+3)^-1); puntoFijo(g,0.5, tol, 100)
La idea de este método se basa en la aproximación lineal que proporciona la recta tangente. Una vez obtenido $x_k$, se calcula $x_{k+1}$ utilizando la aproximación lineal de $f(x)$ mediante la recta tangente que pasa por el punto $(x_k,f(x_k))$:
$$f(x) ≈ f(x_k )+f'(x_k)(x−x_k)$$ Se toma $x_{k+1}$ como la intersección de esta recta con el eje de abscisas, es decir, la solución de $$f(x_k ) + f'(x _k )(x − x_k ) = 0$$
Tomamos el intervalo $[0,1]$ como el intervalo $[a,b]$ donde $f(a)f(b)< 0$. Además, en ese intervalo la derivada es estrictamente positiva por lo que la función es creciente y solo hay un cero. Como punto $x_0$, consideramos el punto medio del intervalo $[0,1]$.
$x_0=0.5$
$x_1=x_0-\frac{x_0-cos(x_0)}{1+sen(x_0)}$ 0.755222417105636
$x_2=x_1-\frac{x_1-cos(x_1)}{1+sen(x_1)}$ 0.739141666149879
$x_3=x_2-\frac{x_2-cos(x_2)}{1+sen(x_2)}$ 0.739085133920807
$x_4=x_3-\frac{x_3-cos(x_3)}{1+sen(x_3)}$ 0.739085133215161
$x_5=x_4-\frac{x_4-cos(x_4)}{1+sen(x_4)}$ 0.739085133215161
El código Matlab a utilizar para calcular las cinco iteraciones es el siguiente:
x=0.5
for k=1:5
x=x-(x-cos(x))/(1+sin(x))
end
%Se puede ver que a partir de la cuarta iteración se
%estabiliza el resultado
Tomando el punto inicial $x_0$ cercano a la solución, bajo ciertas condiciones sobre la función $f$, el método de Newton es convergente localmente con orden de convergencia cuadrático (de orden 2).
f=@(x) exp(x)+2.^-x+2*cos(x)-6
df=@(x) exp(x)-log(2)*2.^-x-2*sin(x)
%se busca un intervalo [a, b] con f(a)f(b) <0
%Tomamos [1,2]
%Aplicamos el Método de Newton
x=1.5
x=x-f(x)/df(x)
%Se repetiría la instrucción para iterar
%hasta conseguir 16 decimales
En la siguiente tabla se muestran las primeras seis iteraciones viendo que se estabiliza a partir de la sexta iteración.
| $x_0$ | 1.500000000000000 | ||
| $x_1$ | 1.956489721124210 | $x_4$ | 1.829383614494166 |
| $x_2$ | 1.841533061042061 | $x_5$ | 1.829383601933849 |
| $x_3$ | 1.829506013203651 | $x_6$ | 1.829383601933849 |
Se puede observar, que en este caso, se duplica, aproximadamentemente, la cantidad de decimales de cada iteracion (convergencia cuadrática).
Para calcular la diferencia entre dos iteraciones se podría utilizar el siguiente código Matlab:
f=@(x) exp(x)+2^-x+2*cos(x)-6
df=@(x) exp(x)-2^-x*log(2)-2*sin(x)
x=1.5
x1=x-f(x)/df(x),error=abs(x-x1),x=x1;
%Se repetiría el código para realizar
%una nueva iteración
f=@(x) exp(x)+2^-x+2*cos(x)-6
df=@(x) exp(x)+2^-x*log(2)-2*sin(x)
k=1; x(k)=1.5
x(k+1)=x(k)-f(x(k))/df(x(k))
error=abs(x(k)-x(k+1))
k=k+1;
%Se repetirían las tres últimas líneas de código
polyder y polyval específicos para polinomios.
p=[1 3 0 2] %polyder deriva el polinomio
dp=polyder(p)
%polyval(p,a) evalúa el polinomio en a
%Vemos que hay cambio de signo en [-4,-2]
polyval(p,-4)*polyval(p,-2) %Cambio de signo
x=-3
%Aplicamos el método de Newton
x=x-polyval(p,x)/polyval(dp,x)
%repetimos este último paso para iterar, sol=-3.195823345445647
Si se repitiera el proceso, por ejemplo con $x_0=1$, se puede ver que no converge, las iteraciones oscilan entre negativo y positivo. Esto ocurre porque el método de Newton converge localmente.
Como condición de parada para no seguir iterando, se pueden considerar los siguientes tests:
Newton(x0, tol, maxiter) que utilice como método de parada que la
diferencia entre dos iteraciones consecutivas sea menor que la tolerancia fijada.
Se parte de un intervalo $[a,b]$ en el que $f(a)f(b) < 0$ y se toma como punto inicial $x_0=\frac{a+b}{2}$.
En cada iteración se busca una nueva aproximación $x_{k+1}$ y se actualizan $a$ y $b$ como se indica a continuación:
No se aplicará el método de Newton cuando el denominador $f'(x_k)$ sea muy pequeño
(overflow), o mejor, cuando el cociente $\frac{f(x_k)}{f'(x_k)}$ sea muy grande. Es decir,
no utilizaremos el método de Newton cuando
$$|\frac{f(x_k)}{f'(x_k)}| >
\frac{1}{\epsilon_M}$$ Esta última condición se implementa de la forma siguiente:
$$|{f'(x_k)}| < {\epsilon_M} |{f(x_k)}|$$
Mostramos a continuación el algoritmo aplicado a un ejemplo.
f=@(x) exp(x)+2.^-x+2*cos(x)-6
df=@(x) exp(x)-log(2)*2.^-x-2*sin(x)
f(1)*f(2)
em=eps/2
a=1;b=2;fa=f(a);fb=f(b);x=(a+b)/2;fx=f(x);
[a x b],[fa fx fb]
| Punto inicial | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $1$ | $1.5$ | $2$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-1.701113559804675$ | $-1.023283135733256$ | $0.806762425836365$ |
Si $f(a) \cdot f(x) < 0$ se considera $[a,x]$, en caso contrario se toma $[x,b]$. En este caso, el nuevo intervalo $[a,b]$ es $[1.5,2]$. Se comprueba antes de aplicar el método de Newton, si la derivada toma valores pequeños.
a=x;fa=fx;m=df(x);
if abs(m)>em*abs(fx),x=x-fx/m;[a x b],end
fx=f(x); [fa fx fb]
| Iteración 1 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $ 1.500000000000000$ | $1.956489721124211$ | $2.000000000000000$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-1.023283135733256$ | $0.579701373274924$ | $0.806762425836365$ |
El cambio de signo se produce en $[a,x]$.
b=x;fb=fx;m=df(x);
if abs(m)>em*abs(fx),x=x-fx/m,end
fx=f(x);[fa fx fb]
| Iteración 2 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $ 1.500000000000000$ | $1.841533061042061$ | $1.956489721124211$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-1.023283135733256$ | $0.050340951614865$ | $0.579701373274924$ |
El cambio de signo es en $[a,x]$.
b=x;fb=fx;m=df(x);
if abs(m)>em*abs(fx),x=x-fx/m,end
fx=f(x);[a x b],[fa fx fb]
| Iteración 3 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $ 1.500000000000000$ | $1.829506013203651$ | $1.841533061042061$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-1.023283135733256$ | $0.0005021213225915$ | $0.050340951614865$ |
De nuevo, el cambio de signo es en $[a,x]$.
b=x;fb=fx;m=df(x);
if abs(m)>em*abs(fx),x=x-fx/m,end
fx=f(x);[a x b],[fa fx fb]
| Iteración 4 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $ 1.500000000000000$ | $1.829383614494166$ | $1.829506013203651$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-1.023283135733256$ | $0.000000051516139$ | $0.000502121322591$ |
El cambio de signo vuelve a ocurrir en $[a,x]$.
b=x;fb=fx;m=df(x);
if abs(m)>em*abs(fx),x=x-fx/m,end
fx=f(x);[a x b],[fa fx fb]
| Iteración 5 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $ 1.500000000000000$ | $\mathbf{ \color{blue}1.829383601933849}$ | $1.829383614494166$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-1.023283135733256$ | $0.000000000000001$ | $0.0000000515161391$ |
La solución es $x=1.829383601933849$. En este ejemplo, siempre el valor de $x$ obtenido ha estado en el intervalo por lo que en cada iteración se ha aplicado Newton.
f=@(x) exp(x)+2^(-x)+2*cos(x)-6
df=@(x) exp(x)-log(2)*2^(-x)-2*sin(x)
f(-3)*f(-2)
a=-3;fa=f(a);b=-2;fb=f(b);x=(a+b)/2;fx=f(x);em=eps/2;
[a x b],[f(a) f(x) f(b)]
| $a$ | $x$ | $b$ |
| $-3.000000000000000$ | $-2.500000000000000$ | $-2.000000000000000$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $0.069802075166974$ | $-1.863347982977588$ | $-2.696958389857672$ |
El cambio de signo es en $[a,x]$.
b=x;fb=fx;dfx=df(x);
if abs(dfx)>em*abs(fx)
x=x-fx/dfx;
if a<=x && b>=x
disp('Se aplica Newton')
fx=f(x);
else
disp('Se aplica Biseccion')
x=(a+b)/2;fx=f(x);
end
else
x=(a+b)/2;fx=f(x);
end
[a x b],[fa fx fb]
| Iteración 1 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-3.000000000000000$ | $-2.750000000000000$ | $-2.500000000000000$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $0.069802075166974$ | $-1.057505574028503$ | $-1.863347982977588$ |
Repitiendo el código modificado el intervalo $[a,b]$ se obtiene la solución en 6 iteraciones.
| Iteración 2: Bisección | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-3.000000000000000$ | $-2.875000000000000$ | $-2.750000000000000$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $0.069802075166974$ | $-0.536899807501486$ | $-1.057505574028503$ |
| Iteración 3: Newton | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-3.000000000000000$ | $-2.994267548648236$ | $-2.875000000000000$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $0.069802075166974$ | $0.040014356224199$ | $-0.536899807501486$ |
| Iteración 4: Newton | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-2.994267548648236$ | $-2.986542066999646$ | $-2.875000000000000$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $0.040014356224199$ | $0.000174552940576$ | $-0.536899807501486$ |
| Iteración 5: Newton | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-2.994267548648236$ | $2.986508070038639$ | $-2.875000000000000$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $0.000174552940576$ | $0.000000003371666$ | $-0.536899807501486$ |
| Iteración 6: Newton | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-2.986508070038639$ | $\mathbf{ \color{blue} -2.986508069381928}$ | $-2.875000000000000$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $0.000000003371666$ | $0$ | $-0.536899807501486$ |
El cero en el intervalo $[-3,-2]$ es $-2.986508069381928$
NewtonBisección(x0, tol, maxiter)
Con la herramienta siguiente se muestra cómo generar paso a paso el código Matlab/Octave de
la función NewtonBiseccion(x0,tol, maxiter) que utiliza como método de
parada que la diferencia entre dos iteraciones consecutivas sea menor que la tolerancia
fijada.
El método de Newton tiene el problema de que podría darse “overflow” si se divide por un número muy pequeño, o converger a otra raíz (no a la raíz que se busca) o incluso no converger. Además, al requerir la derivada de la función podría ser difícil de obtener en el caso de ciertas funciones como las definidas por integrales o por series infinitas por ejemplo.
Este método se basa en considerar, en lugar de la recta tangente, la aproximación que da la secante en intervalos pequeños.
Para llegar a esta expresión, hay que tener en cuenta que la ecuación de la recta secante es
$$y = f\left( {{x_k}} \right) + {m_k}\left( {x - {x_k}} \right)\,\,\,\,\,\,\,\,\,\,\,\,{m_k} = {{f\left( {{x_k}} \right) - f\left( {{x_{k - 1}}} \right)} \over {{x_k} - {x_{k - 1}}}}$$El método proporciona $x_{k+1}$ como solución de
$$0 = f\left( {{x_k}} \right) + {m_k}\left( {x - {x_k}} \right)$$
Pulsando en el botón solución guiada del interactivo siguiente, se muestra la
programación de la función Secante(x0, x1, tol, maxiter).
De la misma forma que el método de Newton, resulta conveniente programar este método combinado con el método de la bisección.
Se parte de un intervalo $[a,b]$ donde $f(a)f(b)< 0$. Se toma $x_0=a$ y $x_1=b$.
En cada iteración, se busca una nueva aproximación $x_{k+1}$ y se actualiza $a$ y $b$ como se indica a continuación:
Tampoco se aplicará el método de la secante cuando el denominador $m_k$ sea muy pequeño. Esta última condición se implementa de la forma siguiente: $$|{m_k}| > {\epsilon_M} |{f(x_k)}|$$
Los test de parada a utilizar son los mismos que en el caso del método Newton-Bisección.
em=eps/2;format long
f=@(x) erf(x)-0.5
f(0)*f(1)
a=0;b=1;fa=f(a);fb=f(b);x0=a;fx0=fa;x1=b;fx1=fb;
%Calculamos siguiente punto por el mét. de la secante
m=(fx1-fx0)/(x1-x0);
if abs(m)>em*abs(fx1)
x=x1-fx1/m;fx=f(x);
end
%Se comprueba que éstá en [a,b], no necesario bisección
[a x b],[f(a) fx f(b)]
| Iteración 1 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $0$ | $0.593330401707401$ | $1.000000000000000$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-0.500000000000000$ | $0.098584503929305$ | $0.342700792949715$ |
El cambio de signo se produce en $[a,x]$. Tomamos este intervalo como el nuevo $[a,b]$. Actualizamos $x_0$ y $x_1$ que serán respectivamente $x_1$ y $x$.
b=x;x0=x1;x1=x;fx0=f(x0);fx1=f(x1);
% Calculo del siguiente punto
m=(fx1-fx0)/(x1-x0);
if abs(m)>em*abs(fx1)
x=x1-fx1/m;fx=f(x);
end
[a x b],[f(a) fx f(b)] % Se comprueba que está en [a,b]
| Iteración 2 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $0$ | $0.429099981968989$ | $0.593330401707401$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-0.500000000000000$ | $-0.043957753989566$ | $0.342700792949715$ |
El cambio de signo se produce en $[x,b]$. Tomamos este intervalo como el nuevo $[a,b]$. Actualizamos $x_0$ y $x_1$ que serán respectivamente $x_1$ y $x$.
a=x;x0=x1;x1=x;fx0=f(x0);fx1=f(x1);
% Calculo del siguiente punto
m=(fx1-fx0)/(x1-x0);
if abs(m)>em*abs(fx1)
x=x1-fx1/m;fx=f(x);
end
[a x b],[f(a) fx f(b)]
%Se comprueba que está en [a,b]
| Iteración 3 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $0.429099981968989$ | $0.479746018406641$ | $0.593330401707401$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-0.043957753989566$ | $0.002522030582461$ | $0.098584503929305$ |
El cambio de signo se produce ahora en $[a,x]$. Tomamos este intervalo como el nuevo $[a,b]$. Actualizamos $x_0$ y $x_1$ que serán respectivamente $x_1$ y $x$.
b=x;x0=x1;x1=x;fx0=f(x0);fx1=f(x1);
m=(fx1-fx0)/(x1-x0);
if abs(m)>em*abs(fx1)
x=x1-fx1/m;fx=f(x);
end
[a x b],[f(a) fx f(b)]
%Se comprueba que está en [a,b]
| Iteración 4 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $0.429099981968989$ | $0.476997923639157$ | $0.479746018406641$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-0.043957753989566$ | $0.000055407570768$ | $0.002522030582461$ |
De nuevo, el cambio de signo se produce en $[a,x]$. Tomamos este intervalo como el nuevo $[a,b]$. Actualizamos $x_0$ y $x_1$ que serán respectivamente $x_1$ y $x$.
b=x;x0=x1;x1=x;fx0=f(x0);fx1=f(x1);
m=(fx1-fx0)/(x1-x0);
if abs(m)>em*abs(fx1)
x=x1-fx1/m;fx=f(x);
end
[a x b],[f(a) fx f(b)]
%Se comprueba que está en [a,b]
| Iteración 5 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $0.429099981968989$ | $0.476936276206905$ | $0.476997923639157$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-0.043957753989566$ | $-0.000000074435110$ | $0.000055407570768$ |
El cambio de signo de $f$ se produce en $[x,b]$. Tomamos este intervalo como el nuevo $[a,b]$. Actualizamos $x_0$ y $x_1$ que serán respectivamente $x_1$ y $x$.
a=x;x0=x1;x1=x;fx0=f(x0);fx1=f(x1);
m=(fx1-fx0)/(x1-x0);
if abs(m)>em*abs(fx1)
x=x1-fx1/m;fx=f(x);
end
[a x b],[f(a) fx f(b)] % Se comprueba que está en [a,b]
| Iteración 6 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $0.476936193389100$ | $0.476936276206905$ | $0.476997923639157$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $1.0e-04 *$-0.000744351095761$ | $1.0e-04 *$$0.000000021886937$ | $1.0e-04 *$$0.554075707681623$ |
El cambio se produce en $[a,x]$, que será el nuevo $[a,b]$.
b=x;x0=x1;x1=x;fx0=f(x0);fx1=f(x1);
m=(fx1-fx0)/(x1-x0);
if abs(m)>em*abs(fx1)
x=x1-fx1/m;fx=f(x);
end
[a x b],[f(a) fx f(b)]
%Se comprueba que está en [a,b]
| Iteración 7 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $0.476936193389100$ | $\mathbf{ \color{blue} 0.476936276204470}$ | $0.476936276206905$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $1.0e-07 *$-0.744351095760543$ | $1.0e-07 *$$0$ | $1.0e-07 *$$0.000021886936707$ |
Se obtiene la solucion en el punto $0.476936276204470$.
Definimos la función y comprobamos que en en el intervalo $[1,2]$ hay cambio de signo.
%det(A-landa*I) landa=-2 valor propio
%f(x)=det(A(x)+2I)
f=@(x) det([x+2 2 -3 5 0 x^2;1 3 1 1 1 1;
0 x x+2 0 5 -2; 1 -1 1 1 2*x 0;
1 0 1 0 (x^2)+2 -2;0 0 1 0 cos(x) -x+2])
%Buscamos cambios de signo
% f(0)<0
% f(1)<0
% f(2)>0
% Hay otros intervalos donde la función cambia de signo
% pero solo piden calcular un valor de x
Definimos el intervalo inicial $[a,b]=[1,2]$ y tomamos $x_0=1$, $x_1=2$. Calculamos el punto $x_2$ bien aplicando secante si se encuentra en $[a,b]$ o utilizando bisección.
em=eps/2;format long
a=1;b=2;x0=a;x1=b;fx0=f(x0);fx1=f(x1);
m=(fx1-fx0)/(x1-x0);
if abs(m)>em*abs(fx1)
x=x1-fx1/m
[a x b]
end
%En este caso x está dentro del intervalo [a,b]
| Iteración 1 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $1.000000000000000$ | $1.264336373802500$ | $2.000000000000000$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $1.0e-02 *$-0.389735559513635$ | $1.0e-02 *$$0.168432357189391$ | $1.0e-02 *$$1.084656912121114$ |
Hay cambio de signo en $[a,x]$ por lo que este intervalo será el nuevo intervalo $[a,b]$.
b=x;fb=f(b);x1=x;fx1=fx;
m=(fx1-fx0)/(x1-x0);
if abs(m)>em*abs(fx1)
x=x1-fx1/m
[a x b]
end
%En este caso x está dentro del intervalo [a,b]
[a x b],fx=f(x);[fa fx fb]
| Iteración 2 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $1.000000000000000$ | $1.184570415928323$ | $1.264336373802500$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-38.973555951363515$ | $0.507113188549775$ | $16.843235718939074$ |
Vuelve a haber cambio de signo en $[a,x]$ por lo que ejecutando el mismo código anterior se tiene los siguientes valores.
| Iteración 3 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $1.000000000000000$ | $1.182094285551603$ | $1.184570415928323$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-38.973555951363515$ | $-0.008881043100172$ | $0.507113188549775$ |
En este caso hay cambio de signo en $[x,b]$ por lo que el nuevo intervalo $[a,b]$ es $[x,b]$.
a=x;fa=f(a);x0=x1;fx0=fx1;x1=x;fx1=fx;
m=(fx1-fx0)/(x1-x0);
if abs(m)>em*abs(fx1)
x=x1-fx1/m
[a x b]
end
%En este caso x está dentro del intervalo [a,b]
[a x b],fx=f(x);[fa fx fb]
| Iteración 4 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $1.182094285551603$ | $1.182136903509806$ | $1.184570415928323$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-0.008881043100172$ | $0.000004066033751$ | $0.507113188549775$ |
En este caso hay cambio de signo en $[a,x]$ por lo que el nuevo intervalo $[a,b]=[a,x]$.
b=x;fb=f(b);x1=x;fx1=fx;
m=(fx1-fx0)/(x1-x0);
if abs(m)>em*abs(fx1)
x=x1-fx1/m
[a x b]
end
%En este caso x está dentro del intervalo [a,b]
[a x b],fx=f(x);[fa fx fb]
| Iteración 5 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $1.182094285551603$ | $1.182136884006832$ | $1.182136903509806$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-0.008881043100172$ | $0.000000000032507$ | $0.000004066033751$ |
Como vuelve a haber cambio de signo en el intervalo $[a,x]$ se ejecuta el mismo código anterior.
| Iteración 7 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $1.182094285551603$ | $1.182136884006832$ | $1.182136903509806$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $-0.008881043100172$ | $0.000000000000003$ | $0.000000000032507$ |
Nuevamente vuelve a haber cambio de signo en $[a,x]$ por lo que se ejecuta el mismo código otra vez.
| Iteración 7 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $1.182136884006676$ | $\mathbf{ \color{blue} 1.182136884006676}$ | $1.182136884006832$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $1.0e-10 *$$0.000025575661434$ | $1.0e-10 *$$0.000025575661434$ | $1.0e-10 *$$0.325066656820725$ |
En cada iteración el valor de $x$ calculado por el método de la secante ha cumplido que está en el intervalo $[a,b]$ por lo que no se ha aplicado bisección en ninguna iteración.
f=@(x) exp(x)+2.^-x+2*cos(x)-5;
em=eps/2;
format long
a=-5;b=0;fa=f(a);fb=f(b);
x0=a;x1=b;fx0=fa;fx1=fb;
Se calcula el siguiente término de la sucesión, $x$, aplicando secante si el valor obtenido se encuentra en $[a,b]$ o, en caso contrario, utilizando bisección.
%Código cálculo siguiente término
m=(fx1-fx0)/(x1-x0);
if abs(m)>em*abs(fx1)
x=x1-fx1/m;
if a<=x & b>=x
disp('Aplicamos secante')
else
disp('Aplicamos biseccion')
x=(a+b)/2;
end
fx=f(x);
[a x b],[f(a) fx f(b)]
else
disp('División por cero')
end
| Iteración 1 | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-5.000000000000000$ | $-0.174983869789607$ | $0$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $27.574062317925538$ | $-1.062118961781644$ | $-1.000000000000000$ |
Como hay cambio de signo en $[a,x]$ este será el nuevo intervalo $[a,b]$. Se modifica el intervalo, se cambia el valor de $x_0$ por $x_1$ y $x_1$ por $x_0$ y se ejecuta después el código que está en negrita en el paso previo.
b=x;
x0=x1;fx0=fx1;
x1=x;fx1=fx;
%Ejecutar después código para calcular siguiente término
| Iteración 2. Se aplica bisección | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-5.000000000000000$ | $-2.587491934894803$ | $-0.174983869789607$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $27.574062317925538$ | $-0.615010653972961$ | $-1.062118961781644$ |
De nuevo el cambio de signo de $f$ es en $[a,x]$, por lo que este será el nuevo intervalo $[a,b]$. Se modifica el intervalo, se cambia el valor de $x_0$ por $x_1$ y $x_1$ por $x$ y se ejecuta después el código para calcular el siguiente término.
b=x;x0=x1;fx0=fx1;x1=x;fx1=fx;
%Ejecutar código cálculo siguiente término
| Iteración 3. Se aplica bisección | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-5.000000000000000$ | $-3.793745967447402$ | $-2.587491934894803$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $27.574062317925538$ | $7.301512320864664$ | $-0.615010653972961$ |
Ahora el cambio de signo de $f$ es en $[x,b]$, se cambia el intervalo $[a,b]$ y los valores de $x_0$ y $x_1$.
a=x;x0=x1;x1=x;fx0=fx1;fx1=fx;
%Repetir código para calcular el siguiente término
| Iteración 4. Se aplica secante | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-3.793745967447402$ | $-2.681202151333973$ | $-2.587491934894803$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $7.301512320864664$ | $-0.309376060449428$ | $-0.615010653972961$ |
El intervalo donde cambia de signo $f$ es $[a,x]$, cambiamos el intervalo $[a,b]$ y actualizamos $x_0$ y $x_1$ para la próxima iteración.
b=x;x0=x1;x1=x;fx0=fx1;fx1=fx;
%Repetir código para calcular el siguiente término
| Iteración 5. Se aplica secante | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-3.793745967447402$ | $-2.726426099664577$ | $-2.681202151333973$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $7.301512320864664$ | $-0.146504190440662$ | $-0.309376060449428$ |
De nuevo el intervalo donde cambia de signo $f$ es $[a,x]$, cambiamos el intervalo $[a,b]$ y actualizamos $x_0$ y $x_1$ para la próxima iteración.
b=x;x0=x1;x1=x;fx0=fx1;fx1=fx;
%Repetir código para calcular el siguiente término
| Iteración 6. Se aplica secante | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-3.793745967447402$ | $-2.767105303128926$ | $-2.726426099664577$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $7.301512320864664$ | $0.008859811280272$ | $-0.146504190440662$ |
El cambio de signo $f$ es en $[x,b]$, cambiamos el intervalo $[a,b]$ y actualizamos $x_0$ y $x_1$ para la próxima iteración.
a=x;x0=x1;x1=x;fx0=fx1;fx1=fx;
%Repetir código para calcular el siguiente término
| Iteración 7. Se aplica secante | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-2.767105303128926$ | $-2.764785524662011$ | $-2.726426099664577$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $0.008859811280272$ | $-0.000229260174994$ | $-0.146504190440662$ |
El cambio de signo de $f$ es $[a,x]$, cambiamos el intervalo $[a,b]$ y actualizamos $x_0$ y $x_1$ para la próxima iteración.
b=x;x0=x1;x1=x;fx0=fx1;fx1=fx;
%Repetir código para calcular el siguiente término
| Iteración 8. Se aplica secante | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-2.767105303128926$ | $ -2.764844038099813$ | $-2.764785524662011$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $0.008859811280272$ | $-0.000000343381642$ | $-0.000229260174994$ |
De nuevo, el cambio de signo de $f$ es en $[a,x]$, cambiamos el intervalo $[a,b]$ y actualizamos $x_0$ y $x_1$ para la próxima iteración.
b=x;x0=x1;x1=x;fx0=fx1;fx1=fx;
%Repetir código para calcular el siguiente término
| Iteración 9. Se aplica secante | ||
|---|---|---|
| $a$ | $x$ | $b$ |
| $-2.767105303128926$ | $\mathbf{ \color{blue} -2.764844125871619}$ | $-2.764844038099813$ |
| $f(a)$ | $f(x)$ | $f(b)$ |
| $0.008859811280272$ | $0.000000000013341$ | $-0.000000343381642$ |
La solución con la tolerancia indicada es $-2.764844125871619$.
SecanteBiseccion(a, b, tol, maxiter)
Con la herramienta siguiente se muestra cómo generar paso a paso el código Matlab/Octave de
la función SecanteBiseccion(a, b, tol, maxiter) que utiliza como método
de parada que la diferencia entre dos iteraciones consecutivas sea menor que la tolerancia
fijada.
En este apartado, nos ocuparemos de encontrar todas las raíces (reales y complejas) de un polinomio, es decir, calcularemos todas las soluciones de $${a_o} + {a_1}x + {a_2}{x^2} + ... +{a_n}{x^n} = 0$$
Recordemos los siguientes resultados que aplican a los polinomios:
Aunque inicialmente se ha descrito el método de Newton para calcular las raíces de una ecuación $f(x)=0$ siendo $f$ una función real de variable real, el método también es válido para encontrar raíces de funciones complejas. Además, como la derivada de un polinomio es fácilmente calculable, el método de Newton resulta adecuado para calcular las raíces de un polinomio.
Una primera idea sería
Otra posibilidad, podría ser considerar el polinomio resultado de
Dado que el hecho de que dos polinomiostengan coeficientes muy próximos no implica que tengan
las mismas raíces (por ejemplo, basta considerar $p(x)=x^2+1.99x+1.01$ y $q(x)=x^2+2x+1$). El
método alternativo a utilizar es el
Se puede comprobar que
$${{{f_m}\left( x \right)} \over {f{'_m}\left( x \right)}} = {{p\left( x \right)} \over {p'\left( x \right) - p\left( x \right)\sum\limits_{j = 1}^m {{1 \over {x - {\alpha _j}}}} }}$$Para aplicar este método, se debe tener en cuenta:
El método de Newton converge cuadráticamente hacia una raíz $\overline x$ si
$p'(\overline x)\ne0$.
Si la convergencia siguiera siendo lenta, entonces se aplicaría el método de Newton a la segunda derivada para ver si la raíz es doble o a derivadas sucesivas si la multiplicidad es mayor.
Si calculamos una raíz $\alpha$ de $p'(x)$ (o de derivadas sucesivas de p) y observamos que $$\left| {p\left( \alpha \right)} \right| \gg \left| {p'\left( \alpha \right)} \right|$$ entonces debemos concluir que $\alpha$ no es raíz de $p$, lo que ocurre es que $p$ posee dos o más raíces muy próximas. En el caso de que haya raíces próximas entonces la convergencia será lenta y, en general, no podremos obtener una precisión alta en la solución.
En los siguientes ejemplos veremos cómo aplicar el método de Maehly para encontrar todas las raíces de un polinomio utilizando el fichero itraiz.m que realiza una iteración de este método. Los parámetros de los que depende son los coeficientes del polinomio cuya raíz se quiere calcular, los coeficiente del polinomio derivada, el vector con las raíces ya calculadas y el punto a partir del cual se va a iterar.
p=[8 26 54 27 -128 -576 -864 -432]
dp=polyder(p)
v=[],x=i
format long
x=itraiz(p,dp,v,x)
%Se hacen iteraciones hasta encontrar las raíces
v=[x;conj(x)]
x=i;x=itraiz (p,dp,v,x)
v=[x;conj(x)]
x=i;x=itraiz(p,dp,v,x)
Se obtienen así las 7 raíces
v =
-1.056847723228969 + 0.262522264027478i
-1.056847723228969 - 0.262522264027478i
-1.540661467764130 + 1.474052474547216i
-1.540661467764130 - 1.474052474547216i
-0.059269229873907 + 2.202319797999819i
-0.059269229873907 - 2.202319797999819i
2.063556841734012 + 0.000000000000000i
Encuentra una raíz de \(x = \frac{1}{2} sen(x) + \frac{1}{4}\) en el intervalo \([0, \pi/2]\) utilizando el método del punto fijo. Determina previamente que la raíz es única en dicho intervalo.
Resolver \(f(x) = x^2 - e^{-x} = 0\), en el intervalo \([0,1]\) usando el método del punto fijo considerando el problema \(x = g(x)\) mediante la elección de \(g(x)\) siguientes:
Empezando la iteración con \(x_0 = 0.5\). Realizar 100 iteraciones y analizar los resultados. Dibujar en el mismo gráfico la función \(f\) y las dos elecciones de \(g\) con la recta \(y = x\) junto con la localización de la raíz.
Encuentra una raíz de $x^3 - x - 2 = 0$ en el intervalo \([1, 2]\) utilizando el método de bisección.
Implementa el método de la bisección para resolver la ecuación no lineal \(f(x) = \frac{1}{\pi} \int_0^x e^{t^2} dt - \frac{3}{4} = 0\) con un error absoluto \( \epsilon_a = 10^{-10}\).
Dada la ecuación \(3x^3 + 2x = 0\), aplicar el método de Newton comenzando las iteraciones en el punto \(x_0 = 1\). Explicar los resultados. Elegir un punto más conveniente para iniciar las iteraciones y encontrar una raíz real.
Un puente colgante sigue la ecuación de la catenaria:
\[ y = a \cosh\left(\frac{x}{a}\right) \]El punto más bajo del cable está en \( x = 0 \), y queremos encontrar la altura \( H \) de las torres en \( x = 100 \).
Si la altura mínima del cable en el centro es de 20 m, y las torres tienen una altura de 100m, determina el valor de "a" utilizando métodos numéricos.
Una empresa tiene costos y demanda modelados por:
$$ C(x) = 500 + 20x + 0.1x^2, \quad P(x) = 100 - 0.5x $$La empresa desea maximizar sus beneficios. Sabiendo que el beneficio viene dado por Beneficio = Ingresos - Costos, utiliza un método numérico para encontrar el valor de x que maximiza el beneficio
Consideremos la función \(f(x) = 0.2 sen(16x) - x + 1.75\). Esta función tiene un cero en \([1, 2]\). Inicia el método de Newton y observa que diverge. Comprueba que el método de bisección aproxima la solución.
Supongamos que la oscilación de una estructura, dotada de un sistema de amortiguación, ante un movimiento oscilatorio, viene dada por \(y(t) = e^{-2t}(10\cos(2t))\). ¿En qué instante \(t\) la posición de la estructura es \(y(t) = 4\)?
Se puede calcular la raíz n-ésima de un número \(a\) resolviendo la ecuación \(x^n - a = 0\). Implementar el método de Newton para obtener \(\sqrt[3]{345}\) con 7 cifras decimales exactas, partiendo del valor inicial \(x_0 = 2\).
La temperatura en el interior de una pared de espesor \( L \) sigue:
$$ T(x) = T_0 \cosh\left(\frac{x}{L}\right) + T_s \sinh\left(\frac{x}{L}\right) $$Encuentra \( x \) cuando \( T(x) = 50 \)°C, con \(T_0 = 100\)°C, \(T_s = 20\)°C, y \(L = 0.5\).
La función de error \(erf(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt\) está definida para \(x \geq 0\). Determinar el valor de \(x\) tal que \(erf(x) = 0.5\). Nota: Utilizar la función de Matlab `erf` para evaluar las integrales.
La velocidad \(v\) de un paracaidista que cae está dada por \(v = \frac{mg}{c}(1 - e^{-ct/m})\). Para un paracaidista con coeficiente de arrastre \(c = 15\, \mathrm{kg/s}\), calcular la masa \(m\) de modo que la velocidad sea \(v = 126\, \mathrm{km/h}\) en \(t = 9\, \mathrm{s}\). Considera \( g = 9.81 m/s^2\).
En un intercambiador de calor, la temperatura de salida \(T_s\) satisface la ecuación:
$$ T_s = T_e + (T_i - T_e)e^{-\frac{UA}{mc_p}} $$donde \(T_e = 20°C\) es la temperatura exterior, \(T_i = 80°C\) es la temperatura inicial, \(U = 200\, W/(m^2\cdot K)\) es el coeficiente global de transferencia de calor, \(m = 2\, kg/s\) es el flujo másico, y \(c_p = 4186\, J/(kg\cdot K)\) es el calor específico. Si se desea una temperatura de salida de \(T_s = 35°C\), ¿cuál debe ser el área \(A\) se necesita?
Dada la ecuación \(2e^x - x^2 - 1 = 0\):
Se necesita diseñar un tanque esférico que almacene 1000 m³ de líquido. Si la altura del líquido es de 8 metros, ¿cuál debe ser el radio de la esfera? La fórmula para el volumen de un segmento esférico es:
$$ V = \frac{\pi h^2}{3}(3R - h) $$donde \(V\) es el volumen, \(h\) es la altura del líquido y \(R\) es el radio de la esfera.
La ecuación de van der Waals para gases reales es:
$$ (P + \frac{an^2}{V^2})(V - nb) = nRT $$donde \(P\) es la presión, \(V\) es el volumen, \(n\) es el número de moles, \(R\) es la constante de los gases, \(T\) es la temperatura, y \(a\) y \(b\) son constantes específicas del gas. Para el CO₂, \(a = 3.592\, L^2\cdot atm/mol^2\) y \(b = 0.0427\, L/mol\). Si tenemos 1 mol de CO₂ a 300K en un volumen de 1L, ¿cuál es la presión?
En un canal trapezoidal, la profundidad crítica \(y_c\) satisface la ecuación:
$$ \frac{Q^2}{g} = \frac{(by + my^2)^3}{b + 2y\sqrt{1+m^2}} $$donde \(Q = 2\, m^3/s\) es el caudal, \(g = 9.81\, m/s^2\) es la gravedad, \(b = 1\, m\) es el ancho del fondo, y \(m = 1\) es la pendiente lateral. Encuentra \(y_c\).
Para la función \(f(x) = xe^x - 1\):
Se considera la ecuación \(\det(A(x)) = 0\), donde \(x\) es un número real y \(A(x)\) es la matriz:
\[ A(x) = \begin{pmatrix} \cos(x) & 1 & 2 & 3 \\ 0 & 1 & 0 & 1 \\ x & 2 & 3 & 4 \\ \sin(x) & x & 1 & \cos(x) \end{pmatrix} \]
Encuentra una raíz de \(f(x) = \det(A(x)) = 0\) en el intervalo \([0, 1]\).
Utiliza el método de la secante para resolver \(f(x) = 100x^5 - 3995x^4 + 79700x^3 - 794004x^2 + 3160075x - 0\) en el intervalo \([17, 22.2]\) con un error absoluto y relativo de \(10^{-7}\). Realiza los cálculos considerando como puntos iniciales \(x_0 = 17\) y \(x_1 = 22.2\).
Utiliza el método de Maehly para hallar todas las soluciones de las siguientes ecuaciones:
¿Cuántas raíces tiene cada ecuación? Comprueba los resultados con el comando
roots de Matlab.
Se aborda la convergencia de estos métodos, la velocidad con la que se aproximan a la solución y las condiciones necesarias para poder alicarse. Además, se discuten estrategias para la localización y separación de raíces, haciendo hincapié en la importancia de una elección adecuada del intervalo inicial y de las aproximaciones iniciales para garantizar la precisión en la obtención de soluciones numéricas.
Los aspectos clave en el estudio de los métodos numéricos de este tema son:
El ajuste por polinomios en métodos numéricos consiste en encontrar un polinomio de cierto grado que aproxime un conjunto de datos que pueden haberse obtenido de un muestreo o una experimentación o bien de una función dada. Se utiliza cuando se quiere representar una función o una serie de datos mediante una expresión polinómica para facilitar cálculos, interpolaciones o predicciones.
Imaginemos por ejemplo que se tiene una tabla que relaciona la viscosidad dinámica del agua con la temperatura
| T ºC | 0 | 5 | 10 | 20 |
| $\mu_o$ | 1.787 | 1.519 | 1.307 | 1.002 |
¿Cómo estimar la viscosidad a una temperatura de 7.5º? Una opción es encontrar un polinomio que pase por los puntos de la tabla y estimar después el valor de 7.5 con él.
En este capítulo mostraremos dos técnicas para aproximar la función por un polinomio: la interpolación y la aproximación por mínimos cuadrados. Dentro de la interpolación veremos la interpolación de Lagrange y la de Hermite.
El método de interpolación de Lagrange es una técnica que permite encontrar un polinomio que pasa exactamente por un conjunto de puntos dados.
Joseph
Louis Lagrange (1736-1813) fue uno de los más grandes matemáticos de su tiempo. Nació en
Italia pero se nacionalizó Francés. Hizo grandes contribuciones en todos los campos de la
matemática y también en mecánica. Su obra principal es la “Mécanique analytique”(1788). En
esta obra de cuatro volúmenes, se ofrece el tratamiento más completo de la mecánica clásica
desde Newton y sirvió de base para el desarrollo de la física matemática en el siglo
XIX.
de una función $f$ definida en $[a,b]$ es el polinomio de grado menor o igual a $n$ que cumple
$p(x_i)=f(x_i)$ para cada $i$ entre $0$ y $n$. A los puntos ${x_0,x_1,...,x_n}$ se les llama
El video siguiente muestra un ejemplo práctico que calcula los polinomios base de Lagrange con tres nodos.
Los datos son:
$$\color{red}{x_o = - 1}\,\,\,\,\,\,\,\,\,\,\color{blue}{x_1 = 1}\,\,\,\,\,\,\,\,\,\,\,\color{green}{x_2 = 2}$$ $$f\left( {{x_o}} \right) = 2\,\,\,\,\,\,\,\,\,\,f\left( {{x_1}} \right) = 1\,\,\,\,\,\,\,\,\,\,f\left( {{x_2}} \right) = 1$$Los polinomios base de Lagrange son:
$$\begin{array}{lll} {l_o}\left( x \right) & = \frac{{x - {\color{blue}{1}}}}{{{\color{red}{-1}} - {\color{blue}{1}}}} \cdot \,\frac{{x - {\color{green}{2}}}}{{{\color{red}{-1}} - {\color{green}{2}}}} \\ \\ & = \frac{1}{6}\left( {{x^2} - 3x + 2} \right) = \frac{1}{6}{x^2} - \frac{1}{2}x + \frac{1}{3} \end{array}$$ $$\begin{array}{lll} {l_1}\left( x \right) & = \frac{{x - \left( { {\color{red}{-1}}} \right)}}{{{\color{blue}{1}} - \left( { {\color{red}{-1}}} \right)}} \cdot \frac{{x - {\color{green}{2}}}}{{{\color{blue}{1}} - {\color{green}{2}}}} \\ \\ & = - \frac{1}{2}\left( {x + 1} \right)\left( {x - 2} \right) = - \frac{1}{2}{x^2} + \frac{1}{2}x + 1 \end{array}$$ $$\begin{array}{lll} {l_2}\left( x \right) & = \frac{{x - \left( {{\color{red}{-1}}} \right)}}{{{\color{green}{2}} - \left( { {\color{red}{-1}}} \right)}} \cdot \frac{{x - {\color{blue}{1}}}}{{{\color{green}{2}} - {\color{blue}{1}}}} \\ \\ & = \frac{1}{3}\left( {x + 1} \right)\left( {x - {\color{blue}{1}}} \right) = \frac{1}{3}{x^2} - \frac{1}{3} \end{array}$$El polinomio interpolador es:
$$p\left( x \right) = f(x_o){l_0}\left( x \right) + f(x_1){l_1}\left( x \right) + f(x_2){l_2}\left( x \right) =$$ $$p\left( x \right) = 2\,{l_0}\left( x \right) + {l_1}\left( x \right) + {l_2}\left( x \right) = $$ $$={1 \over 6}{x^2} - {1 \over 2}x + {4 \over 3}$$Se puede comprobar que:
$$p\left( { - 1} \right) = 2\,\,\,\,\,\,\,\,\,\,p\left( 1 \right) = 1\,\,\,\,\,\,\,\,\,\,p\left( 2 \right) = 1$$
%Añadir código
Aunque la fórmula del teorema anterior posee un interés teórico, existen otras formas de obtener el polinomio interpolador. Uno de estos ellos, se muestra a continuación.
Como se quiere calcular $p(x) = {a_o} + {a_1}x + ... + {a_n}{x^n}$ de forma que $p({x_i}) = f\left( {{x_i}} \right)$, se tendrá que cumplir que:
$$\begin{pmatrix} {x_{o}^n}&{x_{o}^{n-1}}&{...}&{x_0}&{1}\\ {{x_{1}^n}}&{{x_{1}^{n-1}}}&{...}&{x_{1}}&{1}\\ {...}&{...}&{...}&{...}&{...}\\ {{x_{n-1}^n}}&{{x_{n-1}^{n-1}}}&{...}&{x_{n-1}}&{1}\\ {{x_{n}^n}}&{{x_{n}^{n-1}}}&{...}&{x_{n}}&{1}\\ \end{pmatrix} \begin{pmatrix} {a_n} \\ {a_{n-1}} \\ {...}\\ {a_1}\\ {a_0} \end{pmatrix} = \begin{pmatrix} {f\left( {{x_o}} \right)}\\ {f\left( {{x_1}} \right)}\\ {...}\\ {f\left( {{x_{n-1}}} \right)}\\ {f\left( {{x_n}} \right)} \end{pmatrix}$$El método consiste en construir la matriz $A$ de orden $(n+1)$, el vector b de $\mathbb{R}^n$ y resolver el sistema $Aa=b$, donde $a$ es el vector de los coeficientes. La matriz $A$ es invertible La matriz A es de Vandermonde ya que presenta una progresión geométrica en cada fila. Esta matriz es invertible si y solo si todos los $x_i$ son distintos entre sí. y el problema de interpolación de Lagrange tiene una única solución.
x=[-1 1 2]';b=[2 1 1]';n=2;
A=ones(n+1,1);t=A;
for k=1:n
t=t.*x;
A=[t,A];
end
a=A\b
%Comprobación
polyval(a,x)
vander para hallar la matriz de Vandermonde.
x=[-1 1 2]';
A=vander(x)
Con los cálculos que se muestran en el documento siguiente se observa que el máximo de la diferencia entre la función y el polinomio para $n=5$ es $2.654830101089800e-06$, para $n=10$ es $1.794120407794253e-13$ y para $n=15$ es $ 1.110223024625157e-14$.
Evidentemente, estos valores se pueden obtener dado que se conoce la función.
El error de interpolación mide la diferencia entre la función real y el polinomio interpolante, proporcionando una estimación cuantitativa de la precisión de la aproximación.
Si llamamos $M$ a la cota de la derivada tercera de $f$ en el intervalo, se tendrá que para cada $x$ en $[-1,2]$ se cumple,
$$\left| {f\left( x \right) - p\left( x \right)} \right| \le {M \over {3!}}\left| {\left( {x - {x_o}} \right)\left( {x - {x_1}} \right)\left( {x - {x_2}} \right)} \right|$$ $$\left| {f\left( x \right) - p\left( x \right)} \right| \le {M \over {3!}}\left| {\left( {x + 1} \right)\left( {x - 1} \right)\left( {x - 2} \right)} \right|$$Para obtener el máximo de $q(x)=(x+1)(x-1)(x-2)$ se calculan las raíces de $q'(x)$, que serán los puntos donde tiene un máximo o mínimo, y se evalúa $q$ en esos valores obteniendo la cota como el valor más grande de todos ellos en valor absoluto.
format long
q=poly([-1 1 2]);dq=polyder(q);
r=roots(dq);v=polyval(q,r)
max(abs(v))
%2.112611790922380
La cota del error en funcion de $M$ será entonces:
$$M\cdot 2.112611790922380/6$$Según la expresión del error, el error en un punto $x$ entre $0$ y $pi$ sería
\[\left| {sen\left( x \right) - p\left( x \right)} \right| \leqslant \frac{{\mathop {\max }\limits_{c \in \left[ {0,\pi } \right]} \left| {{f^{(n + 1}}\left( c \right)} \right|}}{{3!}}\,\left| {x - 0} \right|\left| {x - \frac{\pi }{2}} \right|\left| {x - \pi } \right|\] $$ \le \frac{1}{{3!}}\left| {x - 0} \right|\left| {x - \frac{\pi }{2}} \right|\left| {x - \pi } \right|$$Para encontrar una cota en todo el intervalo, se considera el polinomio $q(x)=x(x-\pi/2)(x-\pi)$ y se calculan sus máximos y mínimos, esto es, las raíces de la derivada. El valor absoluto máximo de $q$ en estos valores será una cota de $\left| {q\left( x \right)} \right|$.
format long
q=poly([0 pi/2 pi]);dq=polyder(q);
r=roots(dq);v=polyval(q,r)
max(abs(v))
%1.491790182328260
La cota del error sería entonces,
$$\left| {sen\left( x \right) - p\left( x \right)} \right| \le \frac{1}{{3!}}\left| {x - 0} \right|\left| {x - \frac{\pi }{2}} \right|\left| {x - \pi } \right| \le \frac{{1.492}}{6} \approx {\rm{0.248632}}$$Dado que en este ejemplo se conoce la función, vamos a calcular una estimación del error de la aproximación entre la función y el polinomio interpolante evaluando en varios puntos y calculando el máximo de la diferencia entre $f$ y $p$ en todos ellos.
format long
%a es el polinomio interpolante
v=linspace(0,pi,300);
max(abs(polyval(a,v)-sin(v)))
%0.248631697054710
%Nodos: 0.9, 1, 1.1
%Función log(x)
format long
f=@(x) log(x);
n=2;x=[0.9 1 1.1]';fx=f(x);
A=[x ones(n+1,1)];t=x;
for k=1:n-1
t=t.*x;
A=[t A];
end
p=A\fx
%Comprobación
polyval(p,x)-f(x)
%cota de error
q=poly(x);dq=polyder(q);r=roots(dq);
v=polyval(q,r);M=max(abs(v))
%Derivada tercera de f: 2/x^3 el máximo es 2/(0.9^3)
%El error es M*2/(6*0.9^2)
error=M/(3*0.9^3)
% 1.759945950892248e-04
La elección de nodos que permite minimizar el valor de la expresión siguiente,
\[\mathop {\max }\limits_{x \in \left[ {a,b} \right]} \left| {\prod\limits_{j = 0}^n {\left( {x - {x_j}} \right)} } \right|\]que aparece en la
cota de error de interpolación de Lagrange, son los 
Pafnuti
Lvóvich Tchebyschev (1821 - 1894). El más prominente miembro de la escuela de matemáticas de St.
Petersburg. Hizo investigaciones en Mecanismos, Teoría de la Aproximación de Funciones, Teoría
de los Números, Teoría de Probabilidades y Teoría de Integración. Sin embargo, escribió acerca
de muchos otros temas: formas cuadráticas, construcción de mapas, cálculo geométrico de
volúmenes, etc..
Estos nodos en [0, 1] son las raíces del polinomio del $n+1$-ésimo
Se puede demostrar que
$${T_{n + 1}}(x) = \cos \left( {\left( {n + 1} \right)arc\cos \left( x \right)} \right)\,$$
n=8;a=-3;b=3;
%Para obtener los nodos en un vector
x=(a+b)/2+(b-a)/2*cos((2*(0:n)+1)*pi/(2*n+2));
plot(x,zeros(n+1),'o')
En la siguiente figura se muestran la distribución de los nodos de Chebyshev y se observa que no están igualmente distribuidos en el intervalo $[-3,3]$, se concentran más en los extremos.
Representa la función y los dos polinomios obtenidos.
n=20;% Número de nodos
f=@(x) 1./(1+25*x.^2);
u=linspace(-1,1);
plot(u,f(u))
%polinomio
x=linspace(-1,1,n)'
b=f(x)
m=n-1;
A=x.^(m:-1:0)
p=A\b
hold on
plot(u,polyval(p,u))
hold off
n=20;%Grado polinomio
f=@(x) 1./(1+25*x.^2);
x=linspace(-1,1);
plot(x,f(x))
%Nodos
x=cos((2*(0:n)+1)*pi/(2*n+2))';%21 nodos
%polinomio
b=f(x)
A=x.^(n:-1:0)
p=A\b
hold on
plot(x,polyval(p,x))
hold off
El método de diferencias divididas de Newton es una técnica de interpolación que permite construir un polinomio de forma incremental, evitando recalcular desde cero cuando se añade un nuevo nodo. Gracias a su estructura recursiva, solo es necesario actualizar los términos nuevos, lo que facilita su aplicación en problemas donde se agregan puntos progresivamente.
*** Pendiente de hacer ***Los polinomios de Hermite extienden la interpolación de Lagrange al incluir no solo los valores de la función en ciertos puntos, sino también sus derivadas. Aseguran que la interpolación no solo pase por los puntos dados, sino que también respete las derivadas en esos puntos.
En el siguiente video se muestra un ejemplo de interpolación de Hermite.
Determina el polinomio sin utilizar Matlab/Octave.
Se muestra a continuacion la expresión del error que nos permite estimar la precisión de la interpolación de Hermite.
Se define la función, su derivada y los nodos del ejemplo.
format long
f=@(x) exp(x/2);
df=@(x) exp(x/2)/2
x=[0;pi/2;pi];
El grado del polinomio a calcular es 2+1+1+1=5. Consideremos el término independiente del sistema a resolver.
dx=[f(x);df(x)]
Calculamos la matriz del sistema,
A=[x ones(3,1)];t=x;
for k=1:4
t=t.*x
A=[t A];
end
A
B=[2*x ones(3,1) zeros(3,1)];t=x;
for k=3:5
t=t.*x;
B=[k*t B];
end
B
Los coeficientes del polinomio se obtienen resolviendo un sistema de ecuaciones.
p=[A;B]\dx
%Comprobación
polyval(p,x)-f(x)
dp=polyder(p);polyval(dp,x)-df(x)
Calculamos una cota del error de la aproximación.
q=poly(x);q=conv(q,q);
dq=polyder(q);r=roots(dq);
v=polyval(q,r);
M=max(abs(v));
error=M*exp(pi/2)/((2^6)*factorial(6))
Da una cota del error de interpolación.
En el caso de los métodos de interpolación, no siempre se consigue una mejor aproximación aumentado el número de puntos de interpolación.
Además, cuando se considera un número grande de nodos la manipulación y evaluación el cálculo resulta muy costoso computacionalmente y puede llevar también problemas de redondeo importantes.
Con el siguiente método veremos que en el caso de que se tenga mucha información sobre la función en ciertos puntos, en lugar de aumentar de forma considerablemente el grado del polinomio en los métodos de interpolación, es más ventajoso utilizar el método de mínimos cuadrados.
En el siguiente video y con la siguienta herramienta interactiva se puede ver cómo hacer el ajuste a una recta.
Se escribe a continuación la formulación del método de mínimos cuadrados.
El problema que tratamos de resolver es encontrar el polinomio $p$ de grado menor o igual a $k$ que minimiza la siguiente expresión $$\sum\limits_{j = 0}^n {{\omega _j}} {\left( {p\left( {{x_j}} \right) - f\left( {{x_j}} \right)} \right)^2}$$
Para ver cómo se puede formular algebraicamente este problema, consideremos un polinomio $p$ de grado $k$
$$p\left( x \right) = {a_o} + {a_1}x + {a_2}{x^2} + .. + {a_k}{x^k}$$y las siguientes matrices
Teniendo en cuenta la norma euclidea $${\left\| {Aa - b} \right\|^2} = \sum\limits_{j = 0}^n {{\omega _j}} {\left( {p\left( {{x_j}} \right) - f\left( {{x_j}} \right)} \right)^2}$$ los coeficientes del polinomio $p$ que se busca serían la solución del problema: $$\mathop { (P2) \,\,\,\,{\rm{minimizar}}}\limits_{x \in {{\mathbb R}^{k + 1}}} {\left\| {Ax - b} \right\|^2}$$
La siguiente figura muestra el ajuste a una recta, en el sentido de mínimos cuadrados, a los puntos $(-1,3)$, $(2,-1)$, $(3,3)$, $(4,2)$ y $(5,0)$.
En la siguiente figura se muestra el ajuste por mínimos cuadrados a un polinomio de grado 2 de los mismos puntos.
Imaginemos un muelle colgado verticalmente dentro de un cilindro abierto por los extremos a altura 0, por encima de la apertura superior. La longitud del muelle es desconocida porque su extremo inferior está centro del cilindro y es inalcanzable.
| $h$ | 6 | 7 | 10 |
| $F$ | 3 | 8 | 12 |
Se desea determinar L, y de paso la constante del muelle usando la ley de Hooke: $F=k(h-L)=kh-kL$ donde $F$ denota la fuerza aplicada al muelle y $h$ es la distancia desde el extremo superior del muelle. Se disponen de los siguientes datos
El problema de mínimos cuadrados consiste en encontrar el vector $x$ de ${\mathbb R}^2$ de forma que sea mínima $\left\| {Ax - b} \right\|$.
Veamos que para resolver el problema $$\mathop {{\rm{minimizar}}}\limits_{x \in {^{k + 1}}} {\rm{ }}{\left\| {Ax - b} \right\|^2}$$ en general es conveniente utilizar la descomposición QR de la matriz $A$. El método que se propone es el siguiente
Veamos la justificación del punto 3. Como Q es ortogonal se cumple
$${\left\| {Qv} \right\|^2} = {\left( {Qv} \right)^T}\left( {Qv} \right) = {v^T}\left( {{Q^T}Q} \right)v = {v^T}Iv = {v^T}v = {\left\| v \right\|^2}$$Por lo tanto, $${\left\| {Ax - b} \right\|^2} = {\left\| {QRx - b} \right\|^2} = {\left\| {Q\left( {Rx - {Q^T}b} \right)} \right\|^2} = $$ $${\left\| {Rx - {Q^T}b} \right\|^2} = {\left\| {\,\left( {_O^{\widehat R}} \right)x - \left( {_{{Z^T}}^{{Y^T}}} \right)b\,} \right\|^2} = {\left\| {\,\left( {_{ - {Z^T}b}^{\widehat Rx - {Y^T}b}} \right)\,} \right\|^2} = $$ $${\left\| {\,\widehat Rx - {Y^T}b\,} \right\|^2} + {\left\| {\,{Z^T}b\,} \right\|^2}$$
El mínimo valor del problema se obtiene para el valor de $x$ que hace cero el término $${\left\| {\,\widehat Rx - {Y^T}b\,} \right\|^2}$$
De la expresión anterior se deduce que, equivalentemente, la solución del problema $(P_2)$ se obtiene calculando el valor $x$ solución del sistema $${\widehat Rx = {Y^T}b}$$
Si $k$ es el grado del polinomio buscado y $QR$ es la factorización de $A$, la matriz $\widehat R$ es la submatriz de $R$ con las primeras $k+1$ filas y la matriz $Y$ es la submatriz de $Q$ considerando las primeras $k+1$ columnas.| $h$ | 6 | 7 | 10 |
| $F$ | 3 | 8 | 12 |
[Q,R]=qr([6 1;7 1;10 1])
$$ Q= \begin{pmatrix} -0.4411 & -0.6777 & -0.5883\\ -0.5147 & -0.3460 & 0.7845\\ -0.7352 &
0.6488 & -0.1961\\ \end{pmatrix} $$ $$ R= \begin{pmatrix} -13.6015 & -1.6910\\ 0 & -0.3749\\ 0 &
0 \\ \end{pmatrix}$$
Como $k$ es 1 al ser una recta, las matrices $Y$, $\widehat R$ son
$$ Y= \begin{pmatrix} -0.4411 & -0.6777 \\ -0.5147 & -0.3460 \\ -0.7352 & 0.6488 \\ \end{pmatrix} \,\,\,\,\, \widehat R= \begin{pmatrix} -13.6015 & -1.6910\\ 0 & -0.3749\\ \end{pmatrix}$$y el sistema a resolver es entonces, $\widehat Rx=Y^Tb$.
Es decir,
$$ \begin{pmatrix} -13.6015 & -1.6910\\ 0 & -0.3749\\ \end{pmatrix} \begin{pmatrix} a_1\\ a_0\\ \end{pmatrix} $$ $$ \,\,\,\,\, = \begin{pmatrix} -0.4411 & -0.5147 & -0.7352 \\ -0.6777 & -0.3460 & 0.6488 \\ \end{pmatrix} \begin{pmatrix} 3\\ 8\\ 12\\ \end{pmatrix}$$ Operando, $$ \begin{pmatrix} -13.6015 & -1.6910\\ 0 & -0.3749\\ \end{pmatrix} \begin{pmatrix} a_1\\ a_0\\ \end{pmatrix}=\begin{pmatrix} -14.2632\\ 2.9847\\ \end{pmatrix} $$ La solución de este sistema es $a_1=2.0385$ y $a_o=-7.9615$
f=@(x) 1./(1+x.^2);
n=100;k=20;a=-5;b=5;
x=linspace(a,b,n+1)';fx=f(x);
A=[ x ones(n+1,1)];t=x;
for iter=1:k-1
t=t.*x;
A=[t A];
end;
[Q,R]=qr(A);
p=R(1:k+1,:)\Q(:,1:k+1)'*fx;
%Representamos la función
s=linspace(a,b,1001)';
fs=f(s);ps=polyval(p,s);
plot(s,[ps fs])
| T (K) | 8 | 10 | 12 | 14 | 16 | 18 |
| C(T) | 0.0236 | 0.0475 | 0.083 | 0.1736 | 0.202 | 0.25 |
Halla el polinomio interpolador y da un valor aproximado del calor específico a 15 K. Representa gráficamente los datos y el polinomio interpolador.
Da una cota del error de interpolación y compara el resultado con el máximo en valor absoluto de la diferencia entre la función y el polinomio tomando 500 puntos del intervalo $[0,1]$.
Dada la función \(f(x) = \sin(\pi x)\), calcula el polinomio de interpolación \(p(x)\) que resuelve los siguientes problemas:
Representa gráficamente ambos polinomios junto con la función \(f(x)\).
| x | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| y | 0.5 | 2.5 | 2 | 4 | 3.5 | 6 | 5.5 |
Determina el polinomio interpolador que se ajusta a estos valores y grafica las soluciones junto con el conjunto \((x, y)\) en la misma figura.
El capítulo también aborda la elección óptima de nodos, como los nodos de Chebyshev, para mejorar la precisión de la interpolación y evitar problemas como el fenómeno de Runge. Además, se discuten las limitaciones de los métodos de interpolación en casos de funciones con comportamientos complejos,
Los aspectos clave en el estudio de los métodos numéricos son:
Los métodos numéricos de integración sirven para aproximar el valor de una integral definida $\int\limits_a^b {f\left( x \right)dx} $ con una precisión dada cuando no es posible resolverla de forma exacta o analítica. El interés de estas técnicas viene motivado por distintas razones:
Estas razones obligan a encontrar métodos alternativos para calcular de forma aproximada el valor de una integral
definida en un intervalo. Por razones históricas, el cálculo aproximado de
integrales se denomina
Veamos un ejemplo de aplicación. La gestión sostenible del agua, objetivo de Desarrollo sostenible número 6, requiere monitorizar el caudal de los ríos para garantizar su uso sostenible. Sin embargo, el caudal no siempre se mide directamente en todo el curso del río, sino que se calcula a partir de datos como la velocidad del agua en diferentes puntos y las características de la sección transversal del río.
Para calcular el volumen de agua que pasa por una sección transversal del río en un día, se debe integrar la velocidad del agua a lo largo de la profundidad del río. La fórmula para ello es la siguiente,
$$Q = \int\limits_a^b {v\left( x \right)dx} $$donde $Q$ es el caudal total, $v(x)$ es la velocidad del agua en el punto $x$ de la seccion transversal y $a$ y $b$ son los límites de la sección transversal. En la siguiente tabla se recogen unos datos para calcular el valor de $Q$ de forma aproximada.
| Profundidad (m) | Velocidad (m/s) |
|---|---|
| 0 | 0.0 |
| 1 | 1.2 |
| 2 | 2.1 |
| 3 | 2.4 |
| 4 | 2.0 |
| 5 | 1.5 |
| 6 | 0.0 |
En la siguiente figura, se muestra la velocidad del agua en función de la profundidad. El área bajo la curva es el valor a calcular.
Si utilizamos una de las aproximaciones que se verá en este tema, la fórmula de los trapecios, se podría obtener un valor aproximado mediante la siguiente expresión :
\[ Q \approx \frac{h}{2} \left[v(x_0) + 2 \sum_{i=1}^{n-1} v(x_i) + v(x_n)\right] \]
donde:
Así, si $p(x)$ es un polinomio interpolador de $f$ en $n+1$ nodos,
$$\left\{ {{x_i}} \right\}_{i = 0}^n \subset \left[ {a,b} \right]\,\,\,\,\,\,\,\,\,\,\,\,p\left( {{x_i}} \right) = f\left( {{x_i}} \right)\,\,i = 0,1,...,n$$una aproximación del valor de la integral de $f$ en $[a,b]$ será
$$\int\limits_a^b {f\left( x \right)dx} \approx \int\limits_a^b {p\left( x \right)dx} $$
se obtendría, $$\int\limits_a^b {p\left( x \right)dx} = \sum\limits_{i = 0}^n {f\left( {{x_i}} \right)\int\limits_a^b {{l_i}\left( x \right)dx} } = $$ $$ = \left( {b - a} \right)\sum\limits_{i = 0}^n {f\left( {{x_i}} \right)\underbrace {\left( {{1 \over {b - a}}\int\limits_a^b {{l_i}\left( x \right)dx} } \right)}_{{w_i}}} $$
Para calcular la calidad de una fórmula de cuadratura se define su grado de precisión o exactitud.
Para ello se divide el intervalo $[a,b]$ en $k$ subintervalos $$a = {\alpha _o} < {\alpha _1} < ... < {\alpha _k} = b$$ obteniendo $$\int\limits_a^b {f\left( x \right)dx} = \sum\limits_{i = 1}^k {\int\limits_{{\alpha _{i - 1}}}^{{\alpha _i}} {f\left( x \right)dx} } $$
Posteriormente, se aplica la fórmula de cuadratura a cada subintervalo. Este tipo de fórmulas
se denominan de
Analizaremos en este apartado las
En el siguiente interactivo se puede obtener un número variable de nodos equidistantes en un intervalo.
Dado que esta fórmula es interpolatoria, será exacta con orden de exactitud $n$. Por lo tanto, los nodos se podrán calcular resolviendo un sistema de ecuaciones.
De cara a obtener los pesos, resulta de utilidad el resultado siguiente.
Este resultado significa que si se cambia de intervalo $[a,b]$ a $[c,d]$ y se consideran los nodos $x_i$ en $[a,b]$ y $t_i$ en $[c,d]$ relacionados de la forma siguiente: $${t_i} = c + {{d - c} \over {b - a}}\left( {{x_i} - a} \right)$$ se tendrá que los pesos de la fórmula interpolatoria en ambos casos coinciden.
En el siguiente interactivo se muestra dos intervalos, $[a,b]$ y $[c,d]$ y cómo está vinculado el punto $x$ en $[a,b]$, en rojo, con el correspondiente punto $t$ en $[c,d]$ que se muestra en naranja teniendo en cuenta el mapeo anterior.
Calculamos a continuación las primeras fórmulas de Newton-Cotes tomando diferentes número de nodos $x_0$, $x_1$, ..., $x_n$. Para ello, utilizaremos que el cálculo de los pesos no depende del intervalo siempre que se consideren equidistantes en dicho intervalo.
En este caso se considerarán dos nodos en el intervalo $[a,b]$, esto es, $x_0=a$ y $x_1=b$. La fórmula de cuadratura será de la forma, $$\int\limits_a^b {f\left( x \right)dx} \approx \left( {b - a} \right)\left[ {{\color{red}{{\omega _o}}}f\left( a \right) + {\color{red}{{\omega _1}}}f\left( b \right)} \right]$$
Teniendo en cuenta que el cálculo de los pesos, no depende del intervalo considerado, tomaremos $[a,b]=[0,1]$. Para obtener $\omega_0$ y $\omega_1$, aplicamos que esta fórmula de cuadratura es exacta de grado $n=1$, por lo que los valores buscados cumplirán el siguiente sistema:
$$\begin{aligned} & \text{Exacta para }\,\,f(x)=1\,\,\,\,\,\,\,\,\,& 1 = \int\limits_0^1 {dx} = {\omega _o} + {\omega _1} \\ & \text{Exacta para }\,\,f(x)=x\,\,\,\,\,\,\,\,\,&{1 \over 2} = \int\limits_0^1 {xdx} = {\omega _1} \end{aligned}$$Resolviendo el sistema ${\omega _o} = {\omega _1} = {1 \over 2}$. Por lo tanto, la fórmula de cuadratura es
En este caso, el número de nodos es 3 nodos equidistantes en $[a,b]$: $x_o=a$, $x_1=\frac{a+b}{2}$ y $x_2=b$ y la fórmula de cuadratura es $$\int\limits_a^b {f\left( x \right)dx} \approx \left( {b - a} \right)\left[ {{\color {red}{\omega _o}}f\left( a \right) +{\color{red} {\omega _1}}f\left( {{{a + b} \over 2}} \right) + {\color{red}{\omega _2}}f\left( b \right)} \right]$$ Para determinar los nodos consideramos el intervalo $[a,b]=[-1,1]$. Resolviendo el sistema de 3 ecuaciones con 3 incógnicas, $${\text{Exacta para} \,\,f(x)=1 \,\,\,\,\,\,\,\,\, 2 = \int\limits_{ - 1}^1 {dx} = 2\left[ {{\omega _o} + {\omega _1} + {\omega _2}} \right]}$$ $${\text{Exacta para} \,\,f(x)=x\,\,\,\,\,\,\,\,\, 0 = \int\limits_{ - 1}^1 {x\,dx} = 2\left[ { - {\omega _o} + 0 \cdot {\omega _1} + {\omega _2}} \right]}$$ $${\text{Exacta para}} \,\,f(x)=x^2,\,\,\,\,\,\,\,\,{2 \over 3} = \int\limits_{ - 1}^1 {{x^2}\,dx} = 2\left[ {{\omega _o} + 0 \cdot {\omega _1} + {\omega _2}} \right]$$
Resolviendo el sistema, obtenemos ${\omega _o} = {\omega _2} = {1 \over 6}$, ${\omega _1} = {2 \over 3}$. Por lo tanto, la fórmula de cuadratura será:
Esta regla se llama
En la siguiente escena interactiva se puede ver las aproximaciones que proporciona esta cuadratura eligiendo la función y el intervalo a considerar. Puede observarse cómo mejora la aproximación cuando la longitud del intervalo es más pequeño.
En este caso el número de nodos equidistantes a considerar en el intervalo $[a,b]$ es 4: $x_0=a$, $x_1=a+\frac{b-a}{3}$, $x_2=a+2 \frac{b-a}{3}$ y $x_3=b$ y la fórmula de cuadratura es $$\int\limits_a^b {f\left( x \right)dx} \approx \left( {b - a} \right)\left[ {{\omega _o}f\left( a \right) + {\omega _1}f\left( { {{2a+b} \over 3}} \right) + } \right.{\rm{ }}$$ $$\left. +{{\omega _2}f\left( {{{a+2b} \over 3}} \right) + {\omega _3}f\left( b \right)} \right]$$ Para determinar los nodos se considera el intervalo $[0,3]$ y se resuelve el sistema de 4 ecuaciones con 4 incógnitas asociado, obteniendo ${\omega _o} = {\omega _3} = {1 \over 8}$, ${\omega _1} = {\omega _2} = {3 \over 8}$ Por lo tanto, la fórmula de cuadratura es:
Procediendo de la misma forma, se pueden obtener las fórmulas de Newton-Cotes que se recogen a continuación.
En este apartado se analizará error de la aproximación, es decir, la diferencia entre el valor exacto dado por la integral y el valor aproximado que se obtiene con la fórmula de Newton-Cotes $$E = \int\limits_a^b {f\left( x \right)dx} - \left( {b - a} \right)\sum\limits_{j = 0}^n {{\omega _j}f\left( {{x_j}} \right)} $$ El siguiente resultado nos da una expresión para este error.
Este teorema se puede escribir de la siguiente manera:
$$E = \left\{ \begin{aligned} &{C_n}{h^{n + 2}}{f^{(n + 1}}\left( \xi \right)\,\,\,\,\,\text{si n impar} \\ &{C_n}{h^{n + 3}}{f^{(n + 2}}\left( \xi \right)\,\,\,\,\,\text{si\ n par} \end{aligned} \right. $$ donde la constante $C_n$ tiene el siguiente valor: $$C_n = \left\{ \begin{aligned} &{1 \over {\left( {n + 1} \right)!}}\int\limits_0^n {t \left( {t - 1} \right)\left( {t - 2} \right)} ...\left( {t - n} \right)dt \,\,\,\,\,\text{si n impar} \\ & {1 \over {\left( {n + 2} \right)!}}\int\limits_0^n {t^2 \left( {t - 1} \right)\left( {t - 2} \right)} ...\left( {t - n} \right)dt \,\,\,\,\,\text{si n par} \end{aligned} \right. $$Consideremos la fórmula de Newton-Cotes para $n=1$, que es impar. Tomamos $m=n+1=2$. Se tiene que $$\begin{aligned} E=&\color{red}{\int\limits_a^b {f\left( x \right)dx}} \color{black}{-} \color{green}{{{b - a} \over 2}\left( {f\left( a \right) + f\left( b \right)} \right)} \\ E=&{C_1} h^3 f''\left( \xi \right) \end{aligned}$$
Para determinar $C_1$, consideramos la función $f(x)=(x-a)^m=(x-a)^2$. Para esta función se tiene que por un lado se cumple
$$\color{red}{\int\limits_a^b {f(x)dx}} \color{black}{=\int\limits_a^b {{(x-a)^2}dx} = {1 \over 3}{\left( {b - a} \right)^3} =} \color{red}{{ {h^3} \over 3}}$$y por otro
$$\color{green}{{{b - a} \over 2}\left( {f\left( a \right) + f\left( b \right)} \right)} \color{black}{= {{b - a} \over 2}\left( {0 + {{\left( {b - a} \right)}^2}} \right) =} {{{{\left( {b - a} \right)}^3}} \over 2} = \color{green}{ \frac {h^3} {2}}$$Por lo tanto,
$${C_1}{h^3}{f^{''}}\left( \xi \right) = {{{h^3}} \over 3} - {{h^3} \over 2}$$es decir,
$$2 {C_1} = {1 \over 3} - {1 \over 2} = - {1 \over 6}$$ $${C_1} = - {1 \over 12}$$Consideremos la fórmula de Newton-Cotes para $n=2$ que es par, consideramos $m=n+2=4$. Se tiene que $$\begin{aligned} E=&\int\limits_a^b {f\left( x \right)dx} - {{{b - a} \over 6}\left[ {f\left( a \right) + 4f\left( {{{a + b} \over 2}} \right) + f\left( b \right)} \right]} \\ E=&{C_2} h^5 f''\left( \xi \right) \end{aligned}$$ Para determinar $C_2$ consideramos la función $f(x)=(x-a)^m=(x-a)^4 $. Se tiene que por un lado se cumple $$\int\limits_a^b {{(x-a)^4}dx} = {1 \over 5}{\left( {b - a} \right)^5} = { {(2h)^5} \over 5}$$ y por otro $${{b - a} \over 6}\left[ {0 + 4\, \cdot {h^4} + {(2h)^4}} \right] = {{2h^5} \over 6} \cdot 20$$
Consideremos la fórmula de Newton-Cotes para $n=3$, que es impar. Tomamos $m=n+1=4$. Se tiene que $$E=\int\limits_a^b {f\left( x \right)dx} - $$ $${{\left( {b - a} \right)} \over 8}\left[ {f\left( a \right) + 3f\left( {{{2a + b} \over 3}} \right) + 3f\left( {{{a + 2b} \over 3}} \right) + f\left( b \right)} \right] $$
Para determinar $C_3$ se considera la función $f(x)=(x-a)^m=(x-a)^4$. Se cumplirá $$\int\limits_a^b {{(x-a)^4}dx} = {1 \over 5}{\left( {b - a} \right)^5} = { {(3h)^5} \over 5}$$ $${{3h} \over 8}\left[ 0 + 3\, \cdot {h^4} + 3\, \cdot (2h)^4 + (3h)^4 \right] = {{99h^5} \over 2} $$ Por lo tanto, $$E={C_3} h^5 f^{(4}\left( \xi \right)=\frac{-9}{10}h^5$$ $$4!{C_3} = {{ - 9} \over {10}}\,\,\,\,\,\, \Rightarrow \,\,\,{C_3} = - {3 \over {80}}$$
En la siguiente imagen se muestra la expresión del resto para cada una de las fórmulas de Newton-Cotes vistas anteriormente.
El error de la fórmula de los trapecios
$$E = \color{red}{\int\limits_a^b {f\left( x \right)dx} } - \color{green}{{\rm{ }}{{\left( {b - a} \right)} \over 2}\left[ {f\left( a \right) + f\left( b \right)} \right]}\color{black}{ = - {{{{\left( {b - a} \right)}^3}} \over {12}}{f^{''}}\left( \xi \right)}$$con $\xi \in \left[ {a,b} \right]$. En este caso
$$E = \color{red}{\int\limits_0^\pi {{e^{sen\left( x \right)}}dx}} -\color{green}{ {{\pi} \over 2}\left[ {f\left( 0 \right) + f\left( {\pi} \right)} \right] }\color{black}{= - {{{\pi ^3}} \over {12}}{f^{''}}\left( \xi \right){\rm{ }}\,\,\,\,\xi \in \left[ {0,\pi } \right]}$$Dado que
$$f'\left( x \right) = \cos \left( x \right){e^{sen\left( x \right)}}{\rm{ }}$$$$ f''\left( x \right) = {\cos ^2}\left( x \right){e^{sen\left( x \right)}} - sen\left( x \right){e^{sen\left( x \right)}}$$se tiene que
$$\left| {f''\left( x \right)} \right| \le 2e$$Una cota del error es, por tanto,
$$\left| E \right| = \left| { - {{{\pi ^3}} \over {12}}{f^{''}}\left( \xi \right)} \right| \le {{{\pi ^3}e} \over 6}$$Observamos que, en el caso de la fórmula de los trapecios, el error que se comete en la aproximación depende del cubo de la longitud del intervalo. Si se disminuye esta longitud, el error disminuye rápidamente. Eso se puede lograr dividiendo el intervalo en subintervalos, y aplicando la cuadratura a cada uno de ellos y sumando después los resultados obtenidos.
Cuando la distancia entre los nodos $h$ es de longitud grande, mayor que 1, el error que se comete puede ser grande. Por conseguir una mejor aproximación, se divide el intervalo $[a,b]$ en $N$ subintervalos $$a = {\alpha _o} < {\alpha _1} < ... < {\alpha _N} = b$$ calculando la integral como $$\int\limits_a^b {f\left( x \right)dx} = \sum\limits_{i = 1}^N {\int\limits_{{\alpha _{i - 1}}}^{{\alpha _i}} {f\left( x \right)dx} } $$ aplicando la fórmula de cuadratura a cada uno de los subintervalos $[{\alpha _{i-1}},{\alpha _i}]$.
La
El error en esta aproximación se corresponde con
$$-\frac{(b-a)^3}{12N^2}f''(\xi ) \,\,\,\,\, \xi\epsilon (a,b)$$siendo $N$ el número de subintervalos.
La
El error en esta aproximación es
$$ - {\left( {{{b - a} \over {2N}}} \right)^5}{N \over {90}}{f^{(4}}(\xi ) = - {\rm{ }}{{{{\left( {b - a} \right)}^5}} \over {180 \cdot {{(2N)}^4}}}{f^{(4}}(\xi ) \,\,\,\,\,\,\,\,\,\,\xi \in (a,b)$$En Los métodos compuestos se considera la longitud de cada subintervalo igual, por lo que no tienen en cuenta el comportamiento de la función.
En primer lugar, se calcula la aproximación $I_1$ de la integral de $f$ en $[a,b]$ mediante la cuadratura de los trapecios en el intervalo $[a,b]$.
Además, se considera $I_2$ la aproximación de la integral utilizando la regla del trapecio compuesta con $N=2$ siendo $ \color{red}{E}$ el error de esta aproximación.
Se tendrá:
$$\int\limits_a^b {f\left( x \right)dx} = I_1 - {\color{blue}{\frac{(b-a)^3}{12} f''(\xi _1)}}$$ $$\int\limits_a^b {f\left( x \right)dx} = I_2 {\color{red}{-\frac{(b-a)^3}{12 \cdot 4} f''(\xi _2)}}=I_2+ {\color{red}{E}}$$Si suponemos que $f''(\xi _1)\approx f''(\xi _2)$, restando las dos expresiones anteriores, podemos considerar que una estimación del error utilizando el método de los trapecios en dos subintervalos
$$\left|I_2-I_1\right| \approx \left | \color{red}{E}-\color{blue}{4E} \right|= \left |3E\right|$$ $$\color{red}{\left |E\right|=\left |\int\limits_a^b {f\left( x \right)dx} - I_2\right|\approx \left |\frac{1}{3}(I_2-I_1)\right|}$$Teniendo en cuenta la estimación anterior del error, el método recursivo para calcular la integral consiste en aplicar lo siguiente,
Llamando $I_1$ a la fórmula de cuadratura de Simpson en $[a,b]$ y llamando $I_2$ a la fórmula de cuadratura de Simpson compuesta con $N=2$, se tendrá,
$$\int\limits_a^b {f\left( x \right)dx} = I_1 -{\color{blue}{ \frac{h^5}{90 } f''(\xi _1)}}$$ $$\int\limits_a^b {f\left( x \right)dx} = I_2 - {\color{red}{\frac{h^5}{90 \cdot 2^4 } f''(\xi _2)}}=I_2+\color{red}{E}$$llamando $h=(b-a)/2$.
Si suponemos que $f''(\xi _1)\approx f''(\xi _2)$, restando las dos expresiones anteriores, podemos considerar que una estimación del error utilizando el método de Simpson en dos subintervalos
$$\left|I_2-I_1\right| \approx \left |E-16E \right|= \left |15E\right|$$ $$\color{red}{\left |E\right|=\left |\int\limits_a^b {f\left( x \right)dx} - I_2\right|\approx \left |\frac{1}{15}(I_2-I_1)\right|}$$Teniendo en cuenta la estimación anterior del error, el método recursivo para calcular la integral consiste en aplicar lo siguiente,
Veamos cómo encontrar una estimación del error cuando se utiliza una fórmula de Newton-Cotes considerando dos pasos diferentes $h$ y $h/2$ y una estrategia para aplicar este método adaptativo en general.
Si llamamos $I$ a la aproximación de la integral $$\int\limits_a^b {f\left( x \right)dx} = I + E \,\,\, \text{con} \,\,E = {C_n}{h^{m + 1}}{f^{(m}}\left( \xi \right)$$ siendo $$m = \left\{ \begin{aligned} & n + 1\,\,si\,\,n\,\,impar \\ & n + 2\,\,si\,\,n\,\,par \\ \end{aligned} \right.$$
Aplicando la misma fórmula a los intervalos $[a,\frac{a+b}{2}]$ y $[\frac{a+b}{2},b]$ se tendrá, $$\begin{aligned} \int\limits_a^{\frac {a+b}{2}} {f\left( x \right)dx} &= {I_1} + {E_1} = {I_1} + {C_n}{\left( {{h \over 2}} \right)^{m + 1}}{f^{(m}}\left( {{\xi _1}} \right)\\ & {\xi _1} \in \left[ {a,{{a + b} \over 2}} \right] \end{aligned}$$ $$\begin{aligned} \int\limits_{{{a + b} \over 2}}^b {f\left( x \right)dx} &= {I_2} + {E_2} = {I_2} + {C_n}{\left( {{h \over 2}} \right)^{m + 1}}{f^{(m}}\left( {{\xi _2}} \right)\\ &{\xi _2} \in \left[ {{{a + b} \over 2},b} \right] \end{aligned}$$
Si $[a,b]$ es pequeño, $$I + E{\mkern 1mu} {\mkern 1mu} = {I_1} + {I_2} + {C_n}{{{h^{m + 1}}} \over {{2^m}}}\left[ {{{{f^{(m}}\left( {{\xi _1}} \right) + {f^{(m}}\left( {{\xi _2}} \right)} \over 2}} \right]$$ dado que el intervalo es pequeño, $$I + E{\mkern 1mu} {\mkern 1mu} \approx {I_1} + {I_2} + \underbrace {{C_n}{{{h^{m + 1}}} \over {{2^m}}}{f^{(m}}\left( \xi \right)}_{E/{2^m}}$$ $$I + E \approx {I_1} + {I_2} + {E \over {{2^m}}}$$ $${{{2^m} - 1} \over {{2^m}}}E \approx {I_1} + {I_2} - I$$ $$E \approx {{{2^m}} \over {{2^{m}-1}}}\left( {{I_1} + {I_2} - I} \right)$$
Por lo tanto, fijada una precisión $\epsilon$ deseada de la integral, se describe la estrategia para obtener la aproximación con la precisión deseada.
[itg error]=iehardy(a,b,f)
%Devuelve tanto la estimación de la integral como
%el error cometido. Requiere que la
%función f se programe vectorialmente.
a=0;b=pi;
format long
f=@(x) exp(sin(x).*cos(x));
error=10^-12;
x=[a (a+b)/2 b];
[i1 e1]=iehardy(x(1),x(2),f);
[i2 e2]=iehardy(x(2),x(3),f);
E=[e1 e2];I=[i1 i2];
error=sum(E)
itg=sum(I)
[s k]=max(E)
el valor que se obtiene es $k=1$, Se continúa repitiendo el proceso en el intervalo
$[x(k),x(k+1)]$ hasta conseguir la precisión considerada.
xn=(x(k)+x(k+1))/2;
x=[x(1:k) xn x(k+1:end)];
[i1 e1]=iehardy(x(k),xn,f);
[i2 e2]=iehardy(xn,x(k+2),f);
E=[E(1:k-1) e1 e2 E(k+1:end)]
I=[I(1:k-1) i1 i2 I(k+1:end)];
itg=sum(I)
error=sum(E)
integral(f,0,pi)
Usualmente, si la integral es un número grande, el algoritmo se detiene cuando el error relativo es inferior a $\epsilon_r$, mientras que si $|I|$ es pequeño, el algoritmo finaliza cuando el error absoluto es inferior a $\epsilon_a$.
En general, se tendrá que el código Matlab a ejecutar es
%ea error absoluto
%er error relativo
x=[a (a+b)/2 b]
[i1 e1]=iehardy(x(1),x(2),f)
[i2 e2]=iehardy(x(2),x(3),f)
I=[i1 i2];E=[e1 e2];
error=sum(E);
if error<ea
itg=sum(I);
[itg error]
else
[s,k]=max(E)
end
parar=0
while parar==0
xn=(x(k)+x(k+1))/2
x=[x(1:k) xn x(k+1:end)]
[i1 e1]=iehardy(x(k),xn,f);
[i2 e2]=iehardy(xn,x(k+2),f);
E=[E(1:k-1) e1 e2 E(k+1:end)];
I=[I(1:k-1) i1 i2 I(k+1:end)];
error=sum(E);
if error < ea+er*abs(itg)
itg=sum(I);
[itg error]
parar=1;
else
[s,k]=max(E);
end
end
En este apartado, veremos cómo proceder cuando el intervalo $[a,b]$ no es acotado. En el primer caso, se considerará que $a$ o $b$ o ambos no son finitos (integrales impropias de primera especie) y en el segundo caso cuando $f$ tenga una singularidad en un punto $x_o$ en $[a,b]$ de forma que se cumple $\mathop {\lim }\limits_{x \to {x_o}} \left| {f\left( x \right)} \right| = + \infty $ (integrales impropias de segunda especie).
En este tipo de integrales la función $f$ es continua en $[a,b]$ salvo en un punto $x_0$ donde $$\mathop {\lim }\limits_{x \to {x_0}} \left| {f\left( x \right)} \right| = + \infty $$
Se procede de la siguiente forma:
En este capítulo, se aborda la integración numérica, una herramienta fundamental para calcular integrales cuando las soluciones analíticas no son posibles o prácticas. Se centra en los métodos numéricos para aproximar el valor de una integral definida, describiendo varias técnicas y evaluando su precisión y eficiencia. Además, se estudian métodos específicos para tratar integrales impropias y la forma de gestionar errores asociados a estas aproximaciones.
Los principales aspectos tratados en este tema son los siguientes.
Describen cómo aproximar una integral utilizando funciones polinómicas interpolantes, mostrando su utilidad en casos donde se conoce el valor de la función en puntos específicos.
Estas son técnicas basadas en interpolación polinómica que incluyen reglas conocidas como la del trapecio y Simpson. Se analizan sus derivaciones, aplicaciones y errores asociados.
Extienden las fórmulas de Newton-Cotes a subintervalos más pequeños, mejorando la precisión para integrales de funciones complejas. Además, se introducen métodos adaptativos que ajustan dinámicamente los intervalos según la función.
Se presentan métodos para evaluar integrales que tienen singularidades o límites infinitos, diferenciando entre aquellas de primera y segunda especie.
Se analizan las fuentes de error en los métodos numéricos, como el error de truncamiento y los errores asociados a las aproximaciones polinómicas, con énfasis en cómo minimizarlos.
El capítulo incluye ejemplos prácticos y ejercicios diseñados para aplicar y comprender los conceptos, así como autoevaluaciones para reforzar el aprendizaje.
En Ingeniería, el comportamiento de un sistema puede describirse mediante ecuaciones que relacionan la tasa de variación de una magnitud con la magnitud misma y otras variables asociadas. Igualmente, las leyes físicas se expresan con frecuencia como una razón de cambio, por ejemplo, la caída de un cuerpo, el cambio de temperatura, etc. Estas relaciones suelen expresarse a través de ecuaciones diferenciales ordinarias (EDO). Las ecuaciones diferenciales tienen gran importancia para modelar y explicar el comportamiento de muchos fenómenos.
En este capítulo vamos a estudiar la resolución numérica de dos problemas asociados a las ecuaciones diferenciales ordinarias:
En la siguiente escena se muestra un ejemplo de aplicación de las ecuaciones diferenciales a un sistema muelle-resorte.
En el siguiente ejemplo, las ecuaciones diferenciales modelizan al vaciado de un tanque.
En su forma general, se considera
$$y\left( t \right) = \left( \begin{aligned} &{{y_1}\left( t \right)} \\ & {{y_2}\left( t \right)} \\ &\vdots \\ & {{y_n}\left( t \right)} \\ \end{aligned} \right)\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,f\left( {t,y\left( t \right)} \right) = \left( \begin{aligned} & {{f_1}\left( {t,{y_1},{y_2},...,{y_n}} \right)} \\ &{{f_2}\left( {t,{y_1},{y_2},...,{y_n}} \right)} \\ & \vdots \\ &{{f_n}\left( {t,{y_1},{y_2},...,{y_n}} \right)} \\ \end{aligned} \right)$$ $$\begin{aligned} & y'_1={{f_1}\left( {t,{y_1},{y_2},...,{y_n}} \right)} \\ & y'_2={{f_2}\left( {t,{y_1},{y_2},...,{y_n}} \right)} \\ & \vdots \\ & y'_n={{f_n}\left( {t,{y_1},{y_2},...,{y_n}} \right)} \\ \end{aligned}$$Esta forma permite escribir muchos problemas. Veamos algún ejemplo.
Llamando
$${y_1} = y \,\,\,\,\,\,\,\,{y_2} = y'$$ se puede escribir $$\left\{ \begin{aligned} & y{'_1} = {y_2} \\ & y{'_2} = 1 - {y_1} \\ & y_1(0)=1, \,\,\,\,y_2(0)=1 \end{aligned} \right.$$Considerando
$$y\left( t \right) = \left( \begin{aligned} {{y_1}\left( t \right)} \\ {{y_2}\left( t \right)} \\ \end{aligned} \right)\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,f\left( {t,y\left( t \right)} \right) = \left( \begin{aligned} {{y_2}} \\ {1- {y_1}} \\ \end{aligned} \right)$$la ecuación diferencial puede escribirse de la forma
$$ \left( \begin{aligned} {{y'_1}\left( t \right)} \\ {{y'_2}\left( t \right)} \\ \end{aligned} \right)=\left( \begin{aligned} {{y_2}} \\ {1- {y_1}} \\ \end{aligned} \right)$$con las siguientes condiciones en el punto inicial $t_0=0$
$$y\left( 0 \right) = \left( \begin{aligned} 1 \\ 1 \\ \end{aligned} \right)$$Llamando
$${y_1} = y \,\,\,\,\,\,\,\,{y_2} = y'$$ se puede escribir $$\left\{ \begin{aligned} & y{'_1} = {y_2} \\ & y{'_2} = 2t{y_1} +t^2\\ & y_1(0)=1, \,\,\,\,y_2(0)=0 \end{aligned} \right.$$Considerando
$$y\left( t \right) = \left( \begin{aligned} &{{y_1}\left( t \right)} \\ & {{y_2}\left( t \right)} \\ \end{aligned} \right)\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,f\left( {t,y\left( t \right)} \right) = \left( \begin{aligned} {{y_2}} \\ {2t{y_1} + {t^2}} \\ \end{aligned} \right)$$la ecuación diferencial puede escribirse de la forma
$$ \left( \begin{aligned} {{y'_1}\left( t \right)} \\ {{y'_2}\left( t \right)} \\ \end{aligned} \right)=\left( \begin{aligned} {{y_2}} \\ {2t {y_1}+t^2} \\ \end{aligned} \right)$$con las siguientes condiciones en el punto inicial
$$y\left( 0 \right) = \left( \begin{aligned} 1 \\ 0 \\ \end{aligned} \right)$$
Para asegurar que un problema de valor inicial (PVI) tenga una única solución, se considera el
siguiente
Para que la solución sea única se puede exigir que la función cumpla la
En este apartado se estudiarán los métodos de Runge-Kutta de un solo paso, donde la solución aproximada se obtiene progresivamente de un punto al siguiente a partir de la información de la derivada en puntos intermedios dentro de cada intervalo, sin emplear los valores de pasos previos. Así, dado el problema
$$y'\left( t \right) = f\left( {t,y\left( t \right)} \right)\,\,\,\,\,\,\,\,\,t \in \left[ {{t_o},{t_f}} \right] $$la integración numérica del problema consiste en determinar unos instantes
$${t_o} < {t_1} < {t_2} < ... < {t_N} = {t_f} \,\,\,\,\,\,h_n=t_{n+1}-t_n$$
y valores $$\left\{ {{y_o},\,\,{y_1},\,\,{y_2}...,\,\,{y_N}} \right\}$$ de forma que a partir de
$$y\left( {{t_n}} \right) \approx {y_n}$$ se pueda obtener la aproximación de $y_{n+1}$. En este
apartado veremos distintos métodos de Runge-Kutta
Carl
David Tolmé Runge (30 de agosto de 1856 - 3 de enero de 1927) fue un matemático y físico
alemán conocido por su trabajo en análisis numérico y ecuaciones diferenciales. Nació en
Bremen, Alemania, y estudió matemáticas y física en las universidades de Múnich y Berlín,
obteniendo su doctorado en 1880. Runge es reconocido por desarrollar, junto con Martin
Kutta, los métodos de Runge-Kutta para resolver ecuaciones diferenciales ordinarias de
manera numérica. También realizó importantes contribuciones a la espectroscopía y la teoría
de aproximación. Fue profesor en la Universidad de Hannover y, más tarde, en la Universidad
de Gotinga.
Martin
Wilhelm Kutta (3 de noviembre de 1867 - 25 de diciembre de 1944) fue un matemático alemán
conocido principalmente por su trabajo conjunto con Carl Runge en el desarrollo de los
métodos de Runge-Kutta para la solución numérica de ecuaciones diferenciales ordinarias.
Nacido en Pitschen (entonces parte de Prusia), estudió en las universidades de Breslavia y
Múnich, donde trabajó con destacados matemáticos de la época. Además de su contribución a
los métodos de Runge-Kutta, Kutta realizó estudios en mecánica y aerodinámica, desarrollando
la teoría de sustentación conocida como el "teorema de Kutta-Joukowski". Fue profesor en la
Universidad Técnica de Stuttgart.. Su idea principal es aproximar la solución mediante una combinación ponderada de evaluaciones
de la función derivada en puntos intermedios dentro de un intervalo. Estas evaluaciones permiten
calcular el siguiente valor de la solución con más exactitud, sin necesidad de derivadas
superiores, lo que los hace altamente versátiles y aplicables en problemas científicos e
ingenieriles.
En la imagen siguiente se muestra, para una ecuación diferenciable de primer orden, la solución de un problema de valor inicial y la solución aproximada con 7 pasos considerando el intervalo $[0,4]$. El método numérico deberá indicar cómo obtener cada punto de la solución numérica a partir del anterior. En el caso del gráfico, el avance al siguiente punto utiliza la pendiente en el punto anterior.
El método de Euler es el procedimiento más sencillo para aproximar la solución de un problema de valor inicial. Veamos cómo obtener a partir de $y(t_n) \approx y_n$ el valor $y_{n+1}$.
Dado que, $${y(t_{n + 1})} = {y(t_n)} +\color{blue}{ \int\limits_{{t_n}}^{{t_{n + 1}}} {f\left( {t,y\left( t \right)} \right)dt}}$$
se puede usar la fórmula de cuadratura con un nodo para obtener un valor aproximado de la integral de la forma siguiente
$${y(t_{n + 1})} \approx {y_{n+1}}={y_n} + \color{blue}{{h_n}f\left( {{t_n},{y_n}} \right)}$$En la siguiente figura se muestra gráficamente cómo obtener $y_1$ a partir de $y_0$, esto es, a partir del valor de la solución en $t_0$ se muestra cómo obtener el valor en $t_1$.
Con esta aproximación, vemos que el método usa la pendiente dada por la ecuación diferencial en el extremo inicial de cada paso para “mover” la solución al extremo siguiente.
En este caso, $f(t,y)=0.7y-t^2+1$, el intervalo es $[1,2]$ y la condición inicial es $y(1)=1$. Como se pide que la solución pase por 11 puntos, se considera ($N = 10$). En este caso el valor $h_n$ es $(2-1)/10=0.1$.
Calculamos los dos primeros pasos de la aproximación de la curva solución del PVI
$${t_1} = 1.1 \,\,\,\,\,\,\,\, {y_1} = y\left( {1.1} \right) = y\left( 1 \right) + 0.1 \cdot \,f\left( {1,\,\,y\left( 1 \right)} \right) = 1.07$$ $${t_2} = 1.2 \,\,\,\,\,\,\,\, {y_2} = y\left( {1.2} \right) = y\left( {1.1} \right) + 0.1\, \cdot f\left( {1.1,\,\,y\left( {1.1} \right)} \right) = 1.1239$$ $${t_3} = 1.3 \,\,\,\,\,\,\,\, {y_3} = y\left( {1.3} \right) = y\left( 1.2 \right) + 0.1 \cdot \,f\left( {1.2,\,\,y\left( 1.2 \right)} \right) = 1.15857$$ $${t_4} = 1.4 \,\,\,\,\,\,\,\, {y_4} = y\left( {1.4} \right) = y\left( {1.3} \right) + 0.1\, \cdot f\left( {1.3,\,\,y\left( {1.3} \right)} \right) = 1.17067$$ $$...$$y así sucesivamente con el resto de puntos. Puedes ver estos valores en el siguiente interactivo. Observa que aumentando el valor del número de nodos, la aproximación de la solución (marcada en azul) es mejor. En la parte inferior del interactivo, verás la comparativa con otros métodos que se irán viendo en este capítulo (Heun y Runge-Kutta de 3 y 4 etapas).
Se muestra seguidamente el código Matlab para obtener los cálculos de la solución aproximada del PVI en los 11 nodos.
clear all
f=@(t,y) 0.7*y-t^2+1;%Función
a=1;b=2;y0=1;%Intervalo
N=10;%Pasos (N+1 puntos)
h=(b-a)/N; %Longitud del paso
t(1)=a;
y(1)=y0;
for n=1:N
t(n+1)=t(n)+h;
%dy/dt=f(t,y)
y(n+1)=y(n)+h*f(t(n),y(n));
end
[t y]
Para dibujar en Matlab en una misma gráfica tanto los valores obtenidos como la solución exacta, se añadiría al código anterior el siguiente.
plot(t,y,'o')
hold on
%Solución exacta
syms s(u)
eqn=diff(s)==f(u,s)
cond=s(a)==y0
sol(u)=dsolve(eqn,cond)
t1=linspace(a,b,50);
plot(t1,sol(t1))
hold off
xlabel(['Metodo de Euler con N=',num2str(N),' pasos'])
Nótese que en cada paso se necesita el valor $y(t_n)$, por lo que cuanto más cercano sea el valor $y_n$ al real, $y(t_n)$, más preciso será el método. Así, para calcular $y_1$ se utiliza el valor $y_0$, para calcular $y_2$ sustituimos el valor exacto $y(t_1)$ desconocido por su valor aproximado $y_1$, para calcular $y_3$, sustituimos el valor $y(t_2)$ por su valor aproximado $y_2$, y así sucesivamente.
Como en el caso anterior, para encontrar el valor de $y_{n+1}$ a partir de $y_n$, tendremos en cuenta que
$${y(t_{n + 1})} = {y(t_n)} +\color{blue}{ \int\limits_{{t_n}}^{{t_{n + 1}}} {f\left( {t,y\left( t \right)} \right)dt}}$$El valor aproximado de $y(t_{n+1})$ se obtendrá utilizando la fórmula de cuadratura del trapecio para la integral, esto es,
$$\int\limits_{{t_n}}^{{t_n} + {h_n}} {g\left( t \right)dt} \approx {{{h_n}} \over 2}\left( {g\left( {{t_{_n}}} \right) + g\left( {{t_{n + 1}}} \right)} \right)$$De esta manera teniendo en cuenta que ${y_n} \approx y\left( {{t_n}} \right)$, se tendrá que el siguiente paso se podrá caclular como $y_{n+1}$ de la forma siguiente
$${y_{n + 1}} = {y_n} + {{{h_n}} \over 2}\left[ {f\left( {{t_n},{y_n}} \right) + f\left( {{t_{n + 1}},\color{blue}{{y_{n + 1}}}} \right)} \right]$$ y considerando la fórmula de Euler para aproximar $\color{blue}{{y_{n + 1}}}$ $${y_{n + 1}} = {y_n} + {{{h_n}} \over 2}\left[ {f\left( {{t_n},{y_n}} \right) + f\left( {{t_n} + {h_n},\color{blue}{\underbrace {{y_n} + {h_n}f\left( {{t_n},{y_n}} \right)}_{Euler}}} \right)} \right]$$Este método podría escribirse de la forma siguiente:
$$\begin{aligned} & \,\,\,\,\,\,\,\,\,\,{K_{n,1}} = f\left( t_n,y_n \right) \\ & \,\,\,\,\,\,\,\,\,\,{K_{n,2}} = f\left( {t_n}+{\color{blue}{1}}\,\, h_n, y_n + {\color{magenta}{1}} h_n K_{n,1} \right) \\ & {y_{n+1}} = {y_n} + {h_n} \left[ {\color{red}{\frac{1}{2}}}\,\, K_{n,1} + {\color{red}{\frac{1}{2}}}\,\, {K_{n,2}} \right] \\ \end{aligned}$$Este algoritmo puede asociarse al siguiente tablero
$$\begin{array} {r|r r } 0 & 0 & 0 \\ \color{blue}{1} & \color{magenta}{1} & 0 \\ \hline & \color{red} {1/2} & \color{red} {1/2} \\ \end{array}$$
A continuación se muestra el código Matlab para calcular los pasos del método de Heun.
f=@(t,y) (t-y)/2%Función
t0=0;tf=3;%Intervalo
N=30;%Pasos (N+1 puntos)
h=(tf-t0)/N; %Longitud del paso
tn=t0;t=[tn];yn=1;y=[yn];
for n=1:N %
K1=f(tn,yn);
K2=f(tn+h,yn+h*K1);
yn1=yn+(K1+K2)*h/2;
yn=yn1; tn=tn+h;
y=[y yn];t=[t tn];
end
plot(t,y,'*')
Llamando $y_1=y$, $y_2=y'$, se tiene
$$\begin{aligned} &{y'_1} = {y_2} \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, {y'_2} = - 7seny - 0.1\cos t \\ &{y_1}\left( 0 \right) = 0\,\,\,\,\,\,\,\,\,\,\,\,\,\,{y_2}\left( 0 \right) = 1 \end{aligned} $$En este caso,
$$ \left( \begin{aligned} {{y'_1}\left( t \right)} \\ {{y'_2}\left( t \right)} \\ \end{aligned} \right)=\left( \begin{aligned} {{y_2}} \\ {-7sen({y_1})-0.1cos(t)} \\ \end{aligned} \right)\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \left( \begin{aligned} {{y_1}\left( 0 \right)} \\ {{y_2}\left( 0 \right)} \\ \end{aligned} \right)=\left( \begin{aligned} {{0}} \\ {1} \\ \end{aligned} \right)$$El código para resolver la ode con matlab considerando $N=30$ pasos se muestra a continuación.
Dada el PVI
$$\begin{aligned} & y'\left( t \right) = f\left( {t,y\left( t \right)} \right)\,\,\,\,\,\,\,\,\,t \in \left[ {{t_o},{t_f}} \right] \\ & y\left( {{t_o}} \right) = {y_0} \end{aligned} $$ veamos como calcular, $y(t_{n+1})$ a partir de $y(t_n)$. Para ello, se considera la siguiente igualdad $$y\left( {{t_{n + 1}}} \right) = y\left( {{t_n}} \right)\, + \int\limits_{{t_n}}^{{t_n+h_n}} {y'\left( t \right)dt} $$ es decir,En la imagen se representa un ejemplo con $m=6$ donde se muestran los valores $c_i$ y $t_{h,i}$.
En la expresión anterior observamos que se precisa también calcular una aproximación de los valores $\color{red}{y(t_n+c_ih_n)}$. Considerando
$${y_{n,i}} = y\left( {{t_n} + {c_i}{h_n}} \right) = y\left( {{t_n}} \right) + \int\limits_{{t_n}}^{{t_n} + {c_i}{h_n}} {f\left( {t,y\left( t \right)} \right)dt} $$ se puede utilizar alguna fórmula de cuadratura de $m$ pasos obteniendo $$\color{red}{y\left( {{t_n} + {c_i}{h_n}} \right)} \approx {y_n} + {h_n}\sum\limits_{j = 1}^m {{a_{ij}} {\underbrace {f({t_n} + {c_j}{h_n},{y_{n,j}})}_{{K_{n,j}}} } } $$El método puede escribirse también la forma siguiente:
$$ \begin{aligned} & {t_{n,i}} = {t_n} + {\color{red} {c_i}} {h_n}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,,\,\,\,\,\,\,\,\,\,\,\,1 \le i \le m\\ &\,\,\,\,\,\,\,\,\,{K_{n,1}} = f\left( {{t_{n,1}},{y_{n,1}}} \right)\\ &\,\,\,\,\,\,\,\,\, {K_{n,i}} = f\left( {t_n} + {\color{red} {c_i}} {h_n},y_n+{h_n} \sum\limits_{j = 1}^{m} {\color{green}{a_{ij}}} {K_{n,j}} \right) \,\,\,\,\,\,2 \le i \le m\\ &{y_{n + 1}} = {y_n} + {h_n}\sum\limits_{j = 1}^m {\color{blue}{b_j}}{K_{n,j}}\ \end{aligned} $$
Los métodos de Runge-Kutta pueden describirse a trávés del llamado
John
Charles Butcher (nacido en 1933 en Auckland, Nueva Zelanda) es un matemático reconocido por
sus contribuciones a los métodos numéricos para ecuaciones diferenciales, en particular los
métodos de Runge-Kutta y el "Butcher tableau". Estudió en el Auckland University College y
la Universidad de Sídney, donde obtuvo su Ph.D. en 1961. En 2013, fue nombrado Oficial de la
Orden del Mérito de Nueva Zelanda por sus servicios a las matemáticas.
donde se pueden disponer los parámetros $A$, $b$ y $c$ de la forma siguiente:
| $c$ | $A$ |
| $b^T$ |
Cuando los métodos son explícitos la matriz $A$ es triangular inferior.
Vamos a ver algunos métodos de Runge -Kutta con su tablero de Butcher.
En este caso, se considerará la fórmula de cuadratura de tres puntos
$$\int\limits_{{t_n}}^{{t_n} + {h_n}} {g\left( t \right)dt} \approx {{{h_n}} \over 6}\left[ {g\left( {{t_n}} \right) + 4g\left( {{t_n} + \\ {1 \over 2}{h_n}} \right) + g\left( {{t_n} + {h_n}} \right)} \right]$$Se tendrá, $m=3$, $c_1=0$, $c_2=1/2$, $c_3=1$, $b_1=1/6$, $b_2=3/2$, $b_3=1/6$ obteniendo
$${y_{n + 1}} = {y_n} + {{{h_n}} \over 6}\left[ {f\left( {{t_n},{y_n}} \right) + 4f\left( {{t_n} + {1 \over 2}{h_n},y\left( {{t_n} + {1 \over 2}{h_n}} \right)} \right)} \right.$$ $${ + f\left( {{t_{n + 1}},y\left( {{t_{n + 1}}} \right)} \right)}$$siendo el tablero de Butcher
$$\begin{array} {r|r r r } 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 \\ 1 & -1 & 2 & 0 \\ \hline & 1/6 & 4/6 & 1/6 \\ \end{array}$$Este método se podría escribir:
$$\begin{aligned} & {t_{n,1}} = {t_n}\,\,\,,\,\,\,{t_{n,2}} = {t_n}+h_n/2\,\,\,,\,\,\,{t_{n,3}} = {t_n}+h_n \\ & \,\,\,\,\,\,\,\,\,\,{y_{n,1}} = {y_n} \\ & \,\,\,\,\,\,\,\,\,\,{y_{n,2}} = {y_n} + \frac{{h_n}}{2} f\left( {{t_{n,1}},{y_{n,1}}} \right) \\ & \,\,\,\,\,\,\,\,\,\,{y_{n,3}} = {y_n} + {h_n}\left[ {-f\left( {{t_{n,1}},{y_{n,1}}} \right) + 2f\left( {{t_{n,2}},{y_{n,2}}} \right)} \right] \\ & {y_n} = {y_n} + {{{h_n}} \over 6}\left[ {f\left( {{t_{n,1}},{y_{n,1}}} \right) + 4f\left( {{t_{n,2}},{y_{n,2}}} \right) + f\left( {{t_{n,3}},{y_{n,3}}} \right)} \right] \end{aligned}$$También podría escribirse de la forma siguiente.
Se incluyen los correspondientes tableros de Butcher de algunos métodos de Runge-Kutta.
Método de Euler
$$\begin{array} {r|r } 0 & 0 \\ \hline & 1 \\ \end{array}$$Método de Heun
$$\begin{array} {r|r r } 0 & 0 & 0 \\ 1 & 1 & 0 \\ \hline & 1/2 & 2/2 \\ \end{array}$$Método de Kutta de Tres Etapas
$$\begin{array} {r|r r } 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 \\ 1 & -1 & 2 & 0 \\ \hline & 1/6 & 4/6 & 1/6 \\ \end{array}$$El método dado por la tabla se corresponde con el siguiente algoritmo
$$\begin{aligned} & K_1=f(t_n,y_n)\\ & K_2=f(t_n+\frac{1}{2}h,y_n+\frac{1}{2} h K_1)\\ & K_3=f(t_n+\frac{3}{4}h,y_n+\frac{3}{4} h K_2)\\ \end{aligned}$$ $$y_{n+1}=y_n+\frac{h}{9}(2K_1+3K_2+4K_3)$$Para aproximar la solución del PVI del enunciado en $t=0.4$, partiendo de $t_0=0$ e $y_0=0$ con paso $h=0.2$, siendo $f(t,y)=t+e^y$, hemos de realizar dos pasos.
Con el siguiente código Matlab/Octave se pueden realizar los cálculos.
%Definición de la función
f=@(t,y) t+exp(y)
%Valores h=0.2;t0=0;y0=0;
%Paso 1
K1=f(t0,y0)
K2=f(t0+h/2,y0+h/2*K1)
K3=f(t0+3*h/4,y0+3*h/4*K2)
y1=y0+h/9*(2*K1+3*K2+4*K3)
%Paso 2
t1=t0+h;
K1=f(t1,y1)
K2=f(t1+h/2,y1+h/2*K1)
k3=f(t1+3*h/4,y1+3*h/4*K2)
y2=y1+h/9*(2*K1+3*K2+4*k3)
En el siguiente interactivo se calcula los dos pasos con otros métodos de Runge-Kutta.
Escribimos el código Matlab para resolver el problema definiendo los valores $K_i$ correspondientes al tablero.
Modificamos el código para definir los valores $K_i$ a través de las matrices A,B y C que definen el tablero.
Utiliza el método de Runge-Kutta clásico para encontrar la solución con las condiciones iniciales siguientes
$$y_1(0) = 1,\,\,y_2(0) =1,\,\,\, y_3(0) = 0,\,\,\,y_4(0) = 1$$Se muestra a continuación el código Matlab/Octave para resolver el problema.El tablero de Butcher para el método de Runge-Kutta clásico de 4 pasos es
$$\begin{array} {r|r r r r } 0 & 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 & 0 \\ 1/2 & 0 & 1/2 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ \hline & 1/6 & 2/6 & 2/6 & 1/6\\ \end{array}$$En la siguiente imagen se muestra la trayectoria del satélite una vez resuelto el problema.
Utilizándolas, el código Matlab/Octave para encontrar la solución se muestra a continuación.
%Definición de la función
f=@(t,y) [y(3);y(4);-y(1)/(y(1)^2+y(2)^2)^(3/2);-y(2)/(y(1)^2+y(2)^2)^(3/2)];
%Datos PVI
t0=0;tf=8*pi;
y0=[1;1;0;1];N=100;hn=tf/N;
%Método Runge Kutta
[a,b,c]=mirkf4c;
for iter=1:N
[yn1]=pasorkf(f,tn,hn,yn,a,b,c);
tn=tn+hn;yn=yn1;
t=[t tn];y=[y yn];
end
%Representación de la solución
plot(y(1,:),y(2,:),'r')
axis equal
Veamos a continuación cómo analizar el error en los métodos de Runge-Kutta. Observemos en primer lugar que los métodos pueden escribirse de la forma $${y_{n + 1}} = {y_n} + {h_n}{\Phi _f}\left( {{t_n},{y_n},{h_n}} \right)\,\,\,\,\,\,\,0 \le n \le N - 1$$ Se definen a continuación dos tipos de errores.
Puede demostrarse que si el método es de orden $p$, entonces el error global es estimado por la siguiente expresión:
$$\left\| {{e _n}} \right\| \le {C \over A}\left[ {{e^{\left( {{t_n} - {t_o}} \right)}} - 1} \right]{h^p}$$donde $h = \max {h_n}$ y $A$ es la constante positiva que existe, al ser $f$ derivable, cumpliendo
$$\left\| {{\phi _f}\left( {t,y,h} \right) - {\phi _f}\left( {t,z,h} \right)} \right\| \le A\left\| {y - z} \right\|$$Para obtener el orden de un método en términos de los coeficientes del mismo, Butcher estableció las condiciones que se deberían cumplir. Tomando el vector de $m$ componentes $e=(1,1,...,1)^T$ se debe exigir a todos los métodos que $Ae=c$ y además, a continuación, se indican las condiciones que se deben cumplir para los tres primeros órdenes.
Los métodos vistos hasta este momento tienen una dependencia directa con el paso $h$
considerado, por lo que su eleccíon puede condicionar enormemente el funcionamiento del propio
ḿetodo. Nos planteamos en este apartado ćomo determinar cúal es el paso de tiempo $h$ mas
indicado analizando los llamados
Estos ḿetodos consisten en una variación del paso $h$ a partir del control de los errores de aproximacion obtenidos a medida que se calculan las diferentes iteraciones del ḿetodo.
Así, si llegado el momento, el error de aproximación supera la tolerancia establecida para ́este, el método varía el paso para tener una mejor aproximacion. Por el contrario, si el error es menor que la tolerancia, el ḿetodo amplía el paso con el fin de acelerar el rendimiento.
Estos métodos emplean dos fórmulas de diferente orden, que comparten evaluaciones intermedias de la función derivada, para calcular dos aproximaciones del siguiente paso: una más precisa y otra menos precisa. La diferencia entre ambas proporciona una estimación del error, que puede usarse para ajustar automáticamente el tamaño del paso en los métodos de paso variable, optimizando el equilibrio entre precisión y eficiencia en la solución.
donde las matrices $A$ y $C$ del método de orden $p$ están encajados en $\widehat A$ y $\widehat c$ de orden $q > p$.
Fehlberg construyó pares de métodos encajados de orden 4 y 5. En este método los parametros $A$ y $c$ del ḿetodo de orden 4 son coincidentes con las 5 primeras filas de los paŕametros $\widehat A$ y $\widehat c$ del método de orden 5. Por lo tanto, los valores $K_{n,i}$ para $1 \le i \le 5$ son coincidentes para ambos metodos. Esto simplifica notablemente el numero de cálculos, particularmente el ńumero de evaluaciones de la función $f$.
Para pasar de $t_n$ a $t_{n+1}$ aplicamos el método de orden 5 para obtener $\widehat y_{n+1}$. Para hallar una cota del error local de aproximacion de $y(t_{n+1})$ se aplica el método de orden 4 para obtener $y_{n+1}$. Una estimación del error cometido es $$\left\| {{{\widehat y}_{n + 1}} - {y_{n+1}}} \right\|$$
Además, si $${y_{n + 1}} = {y_n} + {h_n}\sum\limits_{j = 1}^m {{b_j}{K_{n,j}}} $$ $${\widehat y_{n + 1}} = {y_n} + {h_n}\sum\limits_{j = 1}^{\widehat m} {{{\widehat b}_j}{K_{n,j}}} \,\,\,\,\,\,\,\,\,\,\,\widehat m > m$$se tiene que
por lo que basta conocer los parámetros del método de orden $q$ y
La información que se precisa se incluye en el siguiente tablero.
$$ {\begin{array} {r|r r r r r} 0 & & & & & \\ 2/9 & 2/9 & & & & \\ 1/3 & 1/12 & 1/4 & & & \\ 3/4 & 69/128 & -243/128 & 135/64 & & \\ 1 & -17/12 & 27/4 & -27/5 & 16/15 & \\ 5/6 & 65/432 & -5/16 & 13/16 & 4/27 & 5/144 & \\ \hline orden 5 & 47/450 & 0 & 12/25 & 32/225 & 1/30 & 6/25\\ estimadores & 1/150 & 0 & -3/100 & 16/75 & 1/20 & -6/25\\ \end{array}}$$La función rkf45a.m desarrollada por Eduardo Casas (ficheros Matlab) devuelve los valores $A$, $b$, $c$, $est$ para el par de métodos encajados de orden 4 y 5 del método de Runge-Kutta-Fehlberg #1.
Runge-Kutta-Felhberg 4-5 #2
El siguiente tablero describe otro método Runge-Kutta-Fehlberg donde también se encajan también un método de orden 4 y otro de orden 5.
La tolerancia deseada se mide entre el error absoluto $e_a$ y el error relativo $e_r$. Dado que el error de $q$ es menor que el de $p$, se considera para la tolerancia el del método de orden
$$TOL=e_a+e_r||{\hat y_{n + 1}}||$$Con todo lo anterior, a partir de dos métodos encajados de órdenes $p$ y $q$ con $q <q$, utilizando los estimadores, $est$, el algoritmo que daría un paso de método es el siguiente.
La función rkf.m desarrollada por Eduardo Casas (ficheros Matlab) implementa el algoritmo. Para ello necesita:
Esta función devuelve
Antes de aceptar el valor obtenido, se debe comprobar si el error es menor que la tolerancia. Si es menor se acepta el paso y si es mayor se rechaza. En cualquier paso se debe proceder con la estrategia de control de paso siguiente:
Elección de $h_0$
$${h_o} = \max \left\{ {e_a + e_r{{\left\| {{y_o}} \right\|}^{{1 \over {p + 1}}}},{\rm{h_{min}}}} \right\}$$ donde podemos fijar por ejemplo $${h_{\min }} = {{{t_f} - {t_o}} \over {{{10}^6}}}$$
$$e_n = \left\| {{{\hat y}_{n + 1}} - {y_{n + 1}}} \right\|{\rm{ < }}\,\,TOL = {e_a} + {e_r}||{\hat y_{n + 1}}||$$
Previsión del paso siguiente. Si ha habido un fallo del paso en la determinación de una de las dos últimas aproximaciones de la solución, entonces se toma $r_{max}=1$, en caso contrario se fija $r_{max}=5$.
Se calcula $r_n$ de la forma indicada en el punto 3 y se toma $${h_{n + 1}} = \left\{ \begin{aligned} & {h_n}\,\,si\,\,\,{r_n} \in \left[ {1,\,\,1.1} \right]\\ & \min \left\{ {{r_n},{r_{\max }}} \right\}\,\text{en otro caso} \end{aligned} \right.$$ En el caso de que el paso proporcione un error inferior a la tolerancia, tenemos que analizar la posibilidad de incrementar este paso para acelerar la integración. El factor de incremento del paso nunca será mayor que $r_{max}=5$, esto se hace para evitar fallos.
En este apartado vamos a estudiar el problema siguiente:
$$\begin{aligned} & y''\left( t \right) = f\left( {t,y\left( t \right),y'\left( t \right)} \right)\,\,\,\,\,\,\,\,t \in \left[ {{t_o},{t_f}} \right] \cr & y\left( {{t_o}} \right) = {y_o}\,\,\,\,\,,\,\,\,\,y\left( {{t_f}} \right) = {y_f} \end{aligned}$$En este caso no se puede aplicar el método de Runge-Kutta para integrar esta ecuación porque no tenemos $y'(t_0)$. Lo que haremos será resolver el siguiente problema:
$$(P) \,\,\,\,\left\{ \begin{aligned} & y{'_1}\left( t \right) = {y_2}\left( t \right) \\ & y{'_2}\left( t \right) = f\left( {t,{y_1}\left( t \right),{y_2}\left( t \right)} \right)\,\,\,\,\,\,\,t \in \left[ {{t_o},{t_f}} \right] \\ & y_1\left( {{t_o}} \right) = {y_o}\,\,\,\,\,,\,\,\,\,y_2\left( {{t_o}} \right) = x \end{aligned} \right.$$ y buscar el valor de $x$ de forma que la solución $y_x=(y_{x,1},y_{x_2})^T$ cumpla que $y_{x,1}$ para el valor de $t=t_f$ sea $y_f$.Para resolver este problema utilizaremos la función oderkf.m desarrollada por Eduardo Casas (ficheros Matlab que resuelve el problema de valor inicial (P) para un valor de $x$ dado.
Esta función tiene la siguiente sintáxis:
$$[t,y]=oderkf(f,t0,y0,tf,ea,er)$$ donde:Para encontrar el valor $x$ se debe resolver, por el método de la secante, el valor de $x$ que cumple que la solución del problema correspondiente tiene en $t_f$ el valor dado $y_f$. Se construye la función g a la cual aplicar el método de la secante.
function [gx,tx,yx]=g(x)
f=....;
to=. . .;
yo=[y0;x];
tf=. . .;
ea=10^−14;
er=10^−12;
[tx,yx]=oderkf(f,to,yo,tf,ea,er);
%la condición es y(tf)=yf;
gx=yx(1,end)-yf;
end
Veamos el procedimiento en un ejemplo.
Paso 1. Definición de la ecuación diferencial en la función en el fichero odeEj.m
function z=odeEj(t,y)
z=[y(2);1/t^3*(y(1)-t*y(2))^2];
Paso 2. Definición de la función $fe$ a la que aplicar el método de la secante.
function [fx,tx,yx]=fe(x)
f='odeEj';
to=1;
tf=2;
yo=[0;x];
ea=10^-14;
er=10^-12;
[tx yx]=oderkf(f,to,yo,tf,ea,er);
fx=yx(1,end)-2;
end
Paso 3. Aplicamos método de la secante a la función $fe(x)$ para hallar $x$ de forma que $fe(x)=0$.
El método puede adaptarse fácilmente al siguiente ejemplo.
Paso 1. Definición de la ecuación diferencial en la función en el fichero odeEj.m
function z=odeEj(t,y)
z=[y(2);1/t^3*(y(1)-t*y(2))^2];
Paso 2. Definición de la función $fe$ a la que aplicar el método de la secante.
function [fx,tx,yx]=fe(x)
f='odeEj';
to=1;
tf=2;
yo=[0;x];
ea=10^-14;
er=10^-12;
[tx yx]=oderkf(f,to,yo,tf,ea,er);
fx=yx(2,end)-2;
end
Paso 3. Aplicamos método de la secante a la función $fe(x)$ para hallar $x$ de forma que $fe(x)=0$.
Paso 1. Definición de la ecuación diferencial en la función en el fichero ode.m
function z=ode(t,y)
z=[y(2);t*y(1)^3-1-sin(t)];
Paso 2. Definición de la función $g$ a la que aplicar el método de la secante.
function [fx,tx,yx]=g(x)
f='ode';
to=0;
tf=pi/2;
yo=[0;x];
ea=10^-14;
er=10^-12;
[tx yx]=oderkf(f,to,yo,tf,ea,er);
fx=yx(1,end)-1;
end
Paso 3. Aplicamos método de la secante a la función $g(x)$ para hallar $x$ de forma que $g(x)=0$.
function z=ode(t,x)
z=[y(2);sin(y(1))+0.01*t*y(1)+0.1*y(3);
sin(y(1)+y(3))+y(4);cos(y(1)+y(3)+y(4)]
end
En este capítulo se trata la integración numérica de ecuaciones diferenciales, una herramienta fundamental en la resolución de problemas matemáticos y científicos que no tienen solución analítica. Este capítulo presenta diversas técnicas para resolver problemas de valor inicial mediante métodos numéricos y analiza sus características, errores y aplicabilidad.
Los contenidos principales del capítulo son los siguientes.
Se introduce el concepto de ecuaciones diferenciales ordinarias (EDO) y su resolución a través de condiciones iniciales dadas, planteando el problema como una base para aplicar métodos numéricos.
Explicación de las condiciones bajo las cuales las soluciones de las EDO están garantizadas y son únicas, un punto clave para entender la validez de los métodos numéricos.
Este es el método numérico más básico para resolver EDO. Se describe su implementación, ventajas y desventajas, además de discutir el error global acumulado.
Una mejora del método de Euler que utiliza una corrección iterativa para aumentar la precisión. También conocido como método del trapecio, es un ejemplo de método de orden 2.
Se analizan los métodos de mayor precisión basados en evaluaciones de la pendiente intermedia, destacando el método de cuarto orden como estándar debido a su equilibrio entre precisión y esfuerzo computacional.
Se detalla cómo se evalúan y minimizan los errores locales y globales en los métodos Runge-Kutta, incluyendo estrategias para ajustar el tamaño del paso.
Introducción a los métodos adaptativos que ajustan dinámicamente el tamaño del paso en función de una estimación del error, optimizando la eficiencia y precisión del cálculo.
Este apartado aborda problemas de contorno, en los cuales las condiciones iniciales y finales están especificadas. Se introduce el método de tiro como técnica para aproximar soluciones.
La optimización es una disciplina matemática y computacional que busca identificar la mejor solución posible para un problema, maximizando o minimizando una función objetivo dentro de un conjunto de condiciones o restricciones. Su importancia radica en su amplia aplicabilidad en diversos campos como la ingeniería, la economía, la logística y la inteligencia artificial, donde es esencial mejorar procesos, reducir costos, etc. Desde problemas simples, como encontrar la ruta más corta entre dos puntos, hasta desafíos complejos en el diseño de sistemas, la optimización proporciona herramientas y métodos que permiten tomar decisiones óptimas en escenarios con recursos limitados.
La optimización y los métodos numéricos están estrechamente relacionados, ya que muchos problemas de optimización requieren soluciones que no pueden encontrarse de manera analítica y, por lo tanto, dependen de enfoques numéricos. Los métodos numéricos proporcionan algoritmos para aproximar soluciones óptimas en problemas donde la función objetivo es compleja, no lineal o involucra grandes conjuntos de datos.
La optimización es una rama de las matemáticas y la ingeniería que se centra en encontrar la mejor solución posible para un problema, de acuerdo con un conjunto de criterios o restricciones. Su objetivo principal es maximizar o minimizar una función objetivo, como el costo, el tiempo o el rendimiento.
