Esta obra está bajo una licencia Creative Commons 4.0 internacional: Reconocimiento-No
Comercial-Compartir Igual.
Tabla de contenido
Prefacio
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:
Ejercicios resueltos paso a paso, que ilustran cómo aplicar las técnicas numéricas.
Explicaciones detalladas que conectan la teoría con la práctica, facilitando la comprensión
de cómo y por qué funcionan los algoritmos.
Una introducción gradual a herramientas como Matlab/Octave, permitiendo a los lectores
implementar las soluciones y comprobar sus resultados.
Herramientas interactivas gráficas y visuales que refuerzan las explicaciones, ayudando a
que los conceptos abstractos sean más comprensibles
Este enfoque 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, , cuyo enfoque teórico ha sido una referencia clave
para la estructura conceptual de este contenido. Este material aporta un enfoque práctico,
combinando teoría, ejemplos y ejercicios diseñados para facilitar la comprensión y el dominio de
los conceptos fundamentales. Además, se utiliza la programación en Matlab/Octave para facilitar
los cálculos en el desarrollo y resolución de los ejemplos.
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,
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.
Capítulo
Introducción a los Métodos Numéricos y Matlab
Los métodos numéricos
El análisis numérico es una rama de las matemáticas que se centra en el estudio de métodos y
algoritmos para obtener soluciones aproximadas a problemas matemáticos que no pueden resolverse
de manera exacta mediante métodos analíticos convencionales. Estos problemas incluyen ecuaciones
algebraicas, ecuaciones diferenciales, integrales, sistemas de ecuaciones lineales y no
lineales, entre otros.
El estudio de los métodos numéricos exige abordar las siguientes cuestiones:
Desarrollo de métodos. Diseñar y analizar algoritmos que permitan
obtener soluciones aproximadas que sean lo suficientemente precisas y rápidas para aplicarse
a problemas reales.
Estimación y control del error. Dado que las soluciones numéricas
son aproximadas, se debe estudiar y cuantificar el error asociado a las soluciones
obtenidas. Esto incluye tanto el error debido a la naturaleza aproximada de los métodos como
el error de redondeo y truncamiento producido por las limitaciones en la representación
numérica de las computadoras.
Convergencia de los métodos. Esto es, analizar si un algoritmo
converge a una solución exacta cuando el número de iteraciones o cálculos aumenta. La
convergencia es crucial para garantizar que el método numérico mejorará la aproximación con
más cálculos.
Estabilidad y eficiencia. Los algoritmos numéricos deben ser
estables, es decir, pequeñas perturbaciones en los datos de entrada no deben provocar
grandes cambios en los resultados.
Además, los métodos deben ser eficientes en términos del tiempo de cálculo y los recursos
computacionales que utilizan, especialmente cuando se aplican a problemas grandes y
complejos.
Ejemplo de introducción
Como aplicación de un método numérico, consideremos que tenemos 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 potencia solar $P_{solar}(t)$ varía a lo largo del día, dependiendo de la luz solar.
La batería tiene una capacidad máxima de almacenamiento.
La velocidad de carga disminuye a medida que la batería se acerca a su capacidad máxima
(comportamiento típico de una batería).
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:
$E(t)$: Energía almacenada en la batería (kWh).
$P_{solar}(t)$: Potencia solar disponible en función del tiempo (kW).
$k$: Coeficiente de ineficiencia de la batería.
$t$: Tiempo en horas.
Se consideran los siguientes datos y parámetros
Capacidad máxima de la batería: 100 kWh.
Coeficiente de ineficiencia (k): 0.05
Tiempo de simulación: 24 horas.
La potencia solar varía a lo largo del día como la función
senoidal siguiente (máximo de 50 kW) considerando $t=0$ el momento en el que sale el sol
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.
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.
Simulación de almacenamiento de energía en batería solar
La línea azul muestra cómo la energía almacenada en la batería varía con el tiempo,
dependiendo de la energía solar recibida.
La línea discontinua naranja representa la potencia solar disponible a lo largo del día, que
sigue una curva senoidal para simular el ciclo diurno.
¿Cómo se aplica un método numérico para obtener $E(t)$?
Supongamos que tenemos una ecuación diferencial: $\frac{dy}{dx} = f(x, y)$
La ecuación indica que conocemos la pendiente de la función \( y(x) \) en cualquier punto, pero
no cómo obtener $y(x)$ solución de la EDO. El Método de Euler permite
aproximar \( y(x) \) avanzando a través de pequeños pasos en el eje \( x \).
Imagina que intentas seguir una curva en un gráfico, pero en lugar de moverte suavemente a lo
largo de la curva, solo conoces la dirección (pendiente) en el punto donde estás parado. La idea
es dar un paso pequeño en esa dirección, recalculando la dirección en el nuevo punto, y seguir
avanzando de la misma forma. Así es como el Método de Euler aproxima la curva.
Pasos del Método de Euler
Punto inicial: Necesitamos un valor inicial \( y(x_0) = y_0 \) en un punto
conocido \( x_0 \). Es decir, queremos que la solución pase por el punto \( (x_0 , y_0
)\),
Aproximación por pasos: A partir del punto inicial, avanzamos con pasos pequeños
\( \Delta x \), calculando el valor aproximado de \( y \) en el siguiente punto usando
la pendiente. Se obtendría así una lista de puntos \( (x_0 , y_0) \), \( (x_1 , y_1) \),
... \( (x_n , y_n) \) donde $x_n=x_{n-1}+\Delta x$
Fórmula de actualización: En cada paso, el valor de \( y \) se actualiza según la
fórmula: $ y_{n+1} = y_n + \Delta x \cdot f(x_n, y_n) $
En el ejemplo que estábamos viendo, la fórmula de Euler que permite encontrar la solución
aproximada de
Pero, como la capacidad máxima de la batería es 100 kWh, limitamos el valor a \(E(5) = 100
\, \text{kWh}\).
...
Se continuaría de esta manera completando el cálculo para las 24 horas del día.
Errores en cálculos numéricos
Misiles Patriot
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
.
Misil Patriot
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 .
La explosión del cohete Ariane
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.
Cohete Ariane
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 en puede encontrarse en .
El hundimiento de la plataforma petrolífera Sleipner A
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.
Plataforma petrolífera
Introducción a Matlab/Octave
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.
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.
Aunque Matlab es un software de pago, se puede utilizar para realizar los ejercicios el
programa Octave, que es software libre y gratuito. Para ejecutar el código de los ejemplos
sin realizar ninguna instalación puede utilizar los siguientes enlaces:
La información que se recoge seguidamente pretende ser un recordatorio de las principales
características de Matlab/Octave que se utilizarán en los ejemplos propuestos en los siguientes
capítulos. Para una mayor información se puede consultar .
Entorno de trabajo
La ventana de Matlab muestra un escritorio dividido en varias partes:
Las órdenes se escriben en la ventana de comandos, Command Window.
La ventana Workspace proporciona información sobre las variables
utilizadas y permite explorar datos que cree o importe de archivos.
Current Folder: carpeta actual donde se ejecutan los archivos.
Ventana Matlab
Directorio de trabajo
Cuando se guarda código en un fichero y queremos utilizarlo en un momento dado, debemos decir al
programa dónde se encuentran. Matlab ya sabe, porque así viene configurado, dónde tiene
guardadas las funciones propias del programa y de los distintos paquetes pero no sabe dónde
están las que hemos escrito.
Matlab busca funciones y scripts en los directorios especificados por la función 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.
Directorio de trabajo
Ayuda de Matlab
Desde la ventana de comandos, se puede utilizar el comando help
help nombreDelComando Ejemplo: >>help plot
Además, desde el menú help se puede acceder a la ayuda de MathWorks en el caso de Matlab.
Creación de variables
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.
Reglas para los nombres de las variables
El nombre de una variable puede tener letras, números y el guion de subrayar.
El primer carácter tiene que ser una letra, modulo2 es un nombre
válido, pero no lo es 2modulo
Las mayúsculas y las minúsculas tienen valor distintivo. La variable
Modulo es distinta de la variable
modulo.
Dentro de un nombre de variable no puede haber espacios en blanco,
modulo1 es un nombre de variable válido, pero no
modulo 1.
Existen nombres que deben evitarse porque tienen significado propio en Matlab:
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.
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 format al escribir el nùmero $\pi$.
format short
4 dígitos después del punto decimal
3.1416
format long
15 dígitos después del punto decimal
3.141592653589793
format short e
Notacion decimal con 4 dígitos
3.1416e+00
format long e
Notacion decimal con 15 dígitos
3.141592653589793e+00
Operaciones aritméticas
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 (^):
x = 4; y = x + 3; % Suma
z = x * 2; % Multiplicación
A= 2*(x + z); % El valor de A es 20
Precedencia de operadores
Las expresiones se evalúan de izquierda a derecha, la potencia tiene el orden de prioridad más , 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.
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
Utilizar Matlab/Octave como calculadora para calcular \[\frac{{\left( {\frac{\pi }{2} + 4}
\right){2^{3 + {2^2}}}}}{{1 + {3^2}}}\]
Un proyecto cuesta \(C_{\text{inst}} = 10,000\) euros y genera un beneficio anual neto de
\(B_{\text{anual}} = 2,000\) euros durante \(N = 5\) años. Calcula el retorno de inversión
(\(ROI\)) con una tasa de descuento de \(r = 4\%\). Beneficios totales y retorno de la inversión
(RPI): \[ B_{\text{total}} = \sum_{t=1}^{N} \frac{B_{\text{anual}}}{(1 + r)^t} \,\,\, \,\,\,
\,\,\, ROI = \frac{B_{\text{total}} - C_{\text{inst}}}{C_{\text{inst}}} \cdot 100 \]
Un script es un archivo que contiene un conjunto de comandos de
Matlab/Octave que se ejecutan en orden. Para crear un script, simplemente escribe tus 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.
Menú New Script
Funciones anónimas
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_Funcion Ejemplo: Definir la función que calcular 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= 12
Escribe una función que permita pasar de grados Celsius a grados Fahrenheit.
Vector fila: se define separando los elementos con espacios o
comas.
Vector columna: se define separando los elementos con un punto y
coma (;).
Vector equiespaciado en un intervalo
punto_inicial:paso:punto_final Ejemplo: Genera un vector con punto inicial 2 hasta llegar a 7 con un salto 2 >>v=4:2:9 v es el vector [4 6 8]
Comando linspace. Genera un vector fila con números igualmente
espaciados considerando un punto inicial, un punto final y el número de elementos.
linspace(valor_inicial,valor_final,num_puntos) Ejemplo: Genera un vector fila de 7 elementos equiespaciados siendo el primer valor
2 y el último 3 >>linspace(2,3,7) Resultado v=[2.0000 2.1667 2.3333 2.5000 2.6667 2.8333 3.0000]
Genera los siguientes vectores eligiendo un valor cualquiera para $n$
El vector cuyas componentes son $n,\,2n,3n,...,10n$
El vector cuyas componentes son $n - 1,\,n - 3,n - 5,...,n - 41$
El vector cuyas componentes son ${n^2} + 1,\,{n^2} + 3,{n^2} + 5,...,{n^2} + 2n - 1$
Un vector con $10+n$ componentes regularmente espaciadas entre 0 y $\pi$
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.
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(m,n) Ejemplo: Matriz 4x3 de ceros >>zeros(4,3)
Matriz de unos
ones(m,n) 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: Crea una matriz diagonal con el vector [2 3 4] en la diagonal >>diag([2 3 4])
Matriz aleatoria
rand(m,n) Ejemplo: Crea una matriz de valores aleatorios de dimension 2x3 >>rand(2,3)
Acceso a los elementos de una matriz
Elemento específico indicando su fila y columna
A = [1 2 3; 4 5 6; 7 8 9]; elemento = A(2, 3); % Accede al elemento en la fila 2, columna 3 % Resultado: elemento = 6
Fila
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]
Columna
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]
Submatrices especificando rango de filas y columnas
A = [1 2 3; 4 5 6; 7 8 9]; submatriz = A(1:2,2:3); % Submatriz con el rango de filas 1 y 2 % y columnas 2 y 3 de A % Resultado: columna3 = [2 3;5 6]
Submatrices generadas a partir de condiciones
A = [1 2 3; 4 5 6; 7 8 9]; elementos_mayores_a_5 = A(A>5); % Submatriz con el rango de filas 1 y 2 % y columnas 2 y 3 de A % Resultado: elementos_mayores_a_5 = [6 7 8 9]
Operaciones con vectores y matrices
Matlab/Octave está diseñado para trabajar de manera eficiente con matrices, permitiendo realizar
operaciones como transposiciones, productos matriciales y otras manipulaciones comunes.
Operaciones matriciales básicas
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
Multiplicación, división, exponenciación elemento a elemento
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
Traspuesta de una matriz
A = [1 2 3; 3 4 5]; % A es 2x3
A' % Matriz traspuesta de A de orden 3x2
Dados los vectores fila \[x = \left( {4,6,2} \right),\,y = \left( {3, - 2,4} \right)\] y el
número $k=3$, calcular
Crea de un vector que va desde 1 a 22 con salto de 5 puntos. Crea otro vector t con el mismo
número de elementos que el anterior de forma que la suma de los dos vectores tenga todos sus
componentes iguales.
Matlab/Octave ofrecen una amplia variedad de funciones integradas en el programa 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:
Funciones matemáticas básicas
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 de 11 entre 4
Funciones trigonométricas
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)
A = [1 7 3;4 8 6]; % Matriz 2x3
size(A) % Tamaño de A
length(A) % Longitud de A, mayor dimensión
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
Funciones para álgebra lineal
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
Ejercicios sobre funciones en Matlab/Octave
Gráficas 2D
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: Para representar la función $x^2$ en el intervalo $[-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:
Abrir una ventana gráfica: figure.
Añadir etiquetas en los ejes: xlabel,
ylabel.
Superponer gráficos dentro de una misma figura: hold on,
hold off.
Poner o quitar una rejilla al gráfico: grid on,
grid off.
Título de un gráfico: title.
Matríz de gráficos en una misma gráfica: 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;
Una compañía eléctrica tiene la siguiente tarifa. Los primeros 100Kwh se pagarán a 2€ el Kwh,
para los siguientes 200 Kwh costará 3 € y 6 de allí en adelante. Expresa el valor de la factura
como una función de la cantidad de Kwh consumida al mes y representa la función.
Matlab/Octave no solo son programas para realizar cálculos, también es 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(4)+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.
En el ejemplo siguiente, si el valor de x es positivo se escribe el texto "x es positivo", en
caso contrario se escribe por pantalla "x es negativo" si es menor que cero y "x es cero" en
otro caso.
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 .
Ciclo 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 i 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 50 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 podría hacerse con la orden:
vec=1:20;
sum(vec.^2)
Bucle 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: sum(2:10)
Manejo de archivos y datos
Abrir y leer dados de un archivo de texto: fopen,
fscanf
% Lee todos los números de archivo.txt
fileID = fopen('datos.txt', 'r');
data = fscanf(fileID, '%f');
fclose(fileID);
disp(data);
Escritura de archivos: 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);
Carga y guardado de variables: 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');
En el siguiente interactivo puedes practicar con ciclos y bucles.
Ejercicios sobre ciclos y bucles
Autoevaluación
Se presenta a continuación un cuestionario básico sobre Matlab/Octave.
Test básico de Matlab/Octave
Con el siguiente interactivo puedes realizar una autoevaluación por niveles.
Test con distintos niveles
Ejercicios
Define las tres variables variables a=1.5, b=4, c=3.5 y calcula el valor de \(d =
\frac{a}{{\frac{b}{c} + \frac{{{c^2}}}{a}}}\).
Consideremos un gas que escapa de un tanque presurizado bajo un proceso adiabático. Su flujo másico se puede aproximar con la función:
\[
f(P_{\text{ext}}, P_{\text{int}}, K) =
P_{\text{int}} \sqrt{\frac{2K}{K-1} \left( 1 - \left( \frac{P_{\text{ext}}}{P_{\text{int}}} \right)^{\frac{K-1}{K}} \right)}
\]
donde \( P_{\text{int}} \) es la presión interna en el tanque, \( P_{\text{ext}} \) es la presión externa fuera del tanque t \( K \) es el coeficiente adiabático del gas. Escribe esta expresión en Octave y calcula su valor para los valores de \(K=1.4\), \( P_{\text{ext}} =1\) y \( P_{\text{int}} =2\).
Genera un vector \( x \) de 30 componentes regularmente espaciadas entre 0 y \( \pi \).
Evalúa \( x \) en cada una de las funciones siguientes:
<\[ f(x) = \frac{\log(x+2)}{x} f(x) = x^2 + x - e^x f(x) = e^{x^2} \sin(x) \]
Diseña un script en Matlab/Octave que genere una matriz cuadrada \( A \) de tamaño \( n
\times n \), con \( n \) proporcionado por el usuario. Llena la matriz con valores
aleatorios entre 1 y \( n \), duplica los valores de la diagonal principal y reduce a la
mitad los valores de la diagonal secundaria. Imprime la matriz original y la modificada.
Escribe un programa en Matlab/Octave que cree dos vectores \( X \) e \( Y \) de longitud 10
con números aleatorios entre 1 y 20. Compara elemento a elemento:
Si \( X[i] > Y[i] \), almacena la suma en un nuevo vector \( Z \).
Si \( X[i] \leq Y[i] \), almacena el valor absoluto de su diferencia en \( Z \).
Escribe un script que clasifique un conjunto de 24 mediciones de temperaturas en tres
categorías: 'Frío' (< 10°C), 'Templado' (10-20°C), y 'Caluroso' (> 20°C). Muestra los
conteos de cada categoría.
Para \( n = 10, 20, 40, 80, \ldots, 10240 \), calcula la suma de los \( n \) primeros
números naturales aplicando y sin aplicar la fórmula \( S(n) = \frac{n(n+1)}{2} \). Muestra
los resultados.
Construye un vector con los diez primeros términos de la sucesión: \[ x_{n+1} = \frac{1}{2 -
x_n} \quad \text{con } x_1 = 1 \quad \text{para } n \geq 1 \]
Utilizando el algoritmo de Euclides, escribe una M-función que calcule el máximo común
divisor de dos números. Ejemplo: MCD(24, 15) debería devolver 3.
Escribe una M-función que, a partir de tres números dados, calcule la suma de los cuadrados
de los dos números mayores. Por ejemplo, SumaCuadrados(5, 3, 4) debería
devolver 41.
Escribe una M-función llamada SuperaMedia que reciba un vector y devuelva otro
vector con los elementos que sean mayores o iguales a la media del vector de entrada.
En este capítulo, se ha visto la necesidad de utilizar los métodos numéricos para resolver
problemas matemáticos complejos que carecen de soluciones analíticas exactas. Los campos de
aplicación son múltiples, por ejemplo, resolución de ecuaciones no lineales, ecuaciones
diferenciales, cálculo de integrales y aproximación por polinomos.
Los aspectos clave en el estudio de los métodos numéricos son:
Desarrollo de métodos. Diseño y análisis de algoritmos que
proporcionen soluciones aproximadas con precisión y rapidez adecuadas para problemas
reales.
Estimación y control del error. Cuantificación de errores
asociados a soluciones numéricas, incluyendo errores de redondeo y truncamiento
debido a las limitaciones computacionales.
Convergencia de los métodos. Evaluación de si un
algoritmo se aproxima a la solución exacta al aumentar el número de iteraciones o
cálculos.
Estabilidad y eficiencia. Asegurar que pequeñas
perturbaciones en los datos de entrada no causen grandes variaciones en los
resultados, y que los métodos sean eficientes en términos de tiempo de cálculo y
recursos computacionales.
Para facilitar la resolución de problemas numéricos y la visualización de resultados, se ha
introducio la herramienta computacional Matlab/Octave. Los contenidos que se han repasado son
los siguientes:
Entorno de trabajo, directorio de trabajo y ayuda de matlab.
El programa incluye una interfaz intuitiva para gestionar variables, explorar el
directorio de trabajo y acceder a ayuda interactiva para funciones y comandos
Creación de variables y operadores numéricos. Las variables
se crean automáticamente al asignarles un valor, y se pueden manipular mediante
operadores aritméticos, de comparación y lógicos.
Scripts y guiones. Son conjuntos de comandos almacenados en
archivos .m que permiten automatizar cálculos y organizar programas estructuradamente.
Vectores y matrices.El programa está diseñado para trabajar
con vectores y matrices como estructuras de datos principales, permitiendo cálculos
algebraicos y vectoriales de forma eficiente.
Funciones nativas. Incluye un conjunto extenso de funciones
nativas para realizar cálculos, optimizando las tareas numéricas comunes.
Gráficas 2D. Permite crear gráficas 2D para representar datos,
incluyendo opciones para personalizar títulos, etiquetas y leyenda
Introducción a la programacón. Su programación se basa en el
uso de scripts, funciones definidas por el usuario y control de flujo (condicionales y
bucles).
Capítulo
Cuestiones básicas sobre aritmética computacional
Introduccion
El objetivo de este apartado es introducir la importancia de la aritmética computacional en el
contexto de los métodos numéricos. Se aborda cómo los computadores representan números y
realizan cálculos, destacando que, debido a limitaciones inherentes (como la precisión finita y
el almacenamiento en binario), los resultados pueden contener errores.
Se verá que:
Los computadores no pueden representar de forma exacta todos los números reales,
especialmente números irracionales o números muy pequeños o muy grandes.
Las operaciones aritméticas pueden introducir errores acumulativos debido a limitaciones en
la precisión.
Veamos a continuación algunos ejemplos numéricos que comete el ordenador.
Analiza el siguiente código Matlab
.
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$
En la figura se muestra la gráfica de la función $f(x)$ y su desarrollo como polinomio $g(x)$.
¿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 entonces las operaciones
aritméticas también se hacen de manera inexacta.
Gráficas de las funciones del ejemplo 2.1.
Sistemas de numeración
Un sistema de numeración consiste en reglas y símbolos que permiten representar cantidades
numéricas de forma clara y organizada. Sirve como fundamento para contar, calcular y realizar
operaciones. En este apartado, nos centraremos en el sistema decimal y el binario (base 2).
Sistema decimal
En este sistema los números se representan utilizando 10 dígitos, es un sistema de numeración
posicional en el que las cantidades se representan utilizando como base aritmética las potencias
del número diez. Así,
Si consideramos una base $b$, un número real $x$ 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}}...$$
Sistema binario
El Standard de IEEE propone el uso de la base 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)
El procedimiento es sencillo, basta dividir el número en base decimal 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.
En el caso de que el número $x$ sea entero, basta dividir $x$ entre $2$ obteniendo de cociente
$x_1$ y resto $r_1$ (0 o 1) pudiendo escribir $x=2\cdot {\color {red} x_1} + {\color {blue} r_1}
$. Se repite el proceso con $x_1$ y el resto de cocientes hasta que el cociente sea menor que 2
$$x=2\cdot {\color {red} (2\cdot x_2+r_2)}+{\color {blue} r_1} =2^2 {\color {red} x_2}+2{\color
{blue} r_2}+ {\color {blue} r_1} ... $$ Los números $r_1$, $r_2$... nos darán las cifras de la
expresión binaria de $x$.
Conversor decimal a binario
Para la parte fraccionaria $0.\,{f_1}{f_2}{f_3}...$, se deben buscar $a_1,a_2,...$ de forma que
$$0.\,{f_1}{f_2}{f_3}... = {{{a_1}} \over 2} + {{{a_3}} \over {{2^2}}} + {{{a_3}} \over {{2^3}}}
+ ...$$ Estos valores se pueden obtener multiplicando por 2 y teniendo en cuenta que
$\sum\limits_{k = 1}^\infty {{1 \over {{2^k}}}} = 1$.
Conversión a binario. Escena creada por Héctor A. Tabares Ospina y Juan Guillermo Rivera Berrío
El número 0.1 tiene infinitas cifras distintas de cero cuando se escribe en base 2.
En la escena que se presenta a continuación, se muestra cómo hacer la conversión paso a paso a
base 2 de ciertos números decimales. En la escena se puede elegir números enteros y números con
parte fraccionaria, observa el procedimiento para este último caso.
Escena creada por Héctor A. Tabares Ospina y Juan Guillermo Rivera Berrío
Conversión de binario a decimal
Puedes practicar con la siguiente herramienta interactiva introduciendo un número binario sin
parte fraccionaria para convertirlo a decimal.
Conversión a binario. Escena creada por Héctor A. Tabares Ospina y Juan Guillermo Rivera Berrío
En el siguiente ejemplo se muestra un número binario con parte fraccionaria y su conversión a
decimal.
Convierte el número $x=10.111001_2$ 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}$$
Representación de los números reales en el ordenador
Un número real en base 10 puede escribirse de la forma siguiente en el que la primera cifra no
nula aparece inmediatamente a la izquierda del punto decimal.
Observa que para almacenar $x$ basta almacenar el signo $s$, el exponente $e$ y la mantisa $m$.
No es necesario almacenar el 1 (bit oculto).
En el ordenador los números reales se representan según la norma IEEE 754, establecida en 1985
por el Instituto de Ingenieros Eléctricos y Electrónicos. Esta norma establece dos formatos,
simple y doble precision. Para los dos casos se utiliza una expresión
normalizada del número que se denomina coma flotante.
En un computador, la mantisa y el exponente se almacenan en palabras de memoria con un número
limitado de bits. La cantidad de bits asignados a cada uno define el rango y la precisión de los
números representados.
Simple precisión
En este caso se utiliza:
1 bit para el signo
8 bits para el exponente
23 bits para la mantisa
Teniendo en cuenta que se utilizan 8 bits para el exponente, sus valores posibles variarían
desde 0 hasta 255:
Para poder considerar las potencias de 2 negativas se toma un desplazamiento o
sesgo del exponente de 1023. Así, es posible representar los siguientes números:
$$x = \left( { - 1} \right)^s × 2^{e - 1023}× \,\,1.m \,\,\,\,\,\, \,\,\,\, 0 < e < 2047 $$
Almacenamiento del número $3.125_{10}$
Se tiene $$3.125_{10}=11.001_{2}=1.1001 × 2^{1}$$ donde $e=1023+1=1024=2^{10}$, $m=1001$.
Convirtiendo el exponente en binario, se tendrá $10000000000$. La mantisa ajustada a
52 bits es $1001\underbrace{0...0}_{48 \,\,bits}$.
El número se almacenará de la forma: $$ \underbrace{0}_{signo}
\underbrace{10000000000}_{exponente}\,\,\underbrace{10010...0}_{mantisa, \,48\,\,últimos
\,\,bits\,\,0}$$
Se consideran las siguientes representaciones especiales:
Cero: $e=0$, $m=0$.
Infinito: $e=2047$, $m=0$.
NAN: $e=2047$, $m \ne 0$.
Conversión a IEEE 64
Rango y precisión
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:
Rango: Nos da el conjunto de intervalos donde existen números representables, este
valor depende del exponente.
Precisión: Nos da la diferencia entre dos números representables consecutivos,
depende del número de bits de la mantisa.
El número más pequeño que se puede representar 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 épsilon de la máquina y en Matlab se identifica con
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$.
Para cualquier otro número real, la distancia entre x y el número máquina más próximo en Matlab
se utiliza eps(x)
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 desnormalizados.
Los números representados en la máquina no están igualmente distribuidos en la recta real, hay
muchos más cerca del cero.
Todos los números menores que este valor se consideran 0.
>>2^(-1074}
ans= 4.9407e-324
>>2^(-1075}
ans=
0
Errores de redondeo y aritmética computacional
Error absoluto y relativo
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:
Error absoluto: $\left| {\widetilde p - p} \right|$
Error relativo:${{\left| {\widetilde p - p} \right|} \over {\left| p \right|}}$
En general, en los métodos numéricos no se conoce el valor verdadero, por lo tanto tampoco el
error real. Lo frecuente es tener una estimación del error dada por una cota positiva al tamaño
máximo del error.
Calcular el error absoluto y relativo siendo $p=3$ y $\widetilde p=3.1$. Repetir el cálculo para
$p=3000$ y $\widetilde p=3001$.
En 1862 el físico Foucault, utilizando un espejo giratorio, calculó en $298 000$ km/s la
velocidad de la luz. Aceptando como exacta la velocidad de $299 776$ km/s, estimar el error
absoluto y el error relativo cometido por Foucault. Por otra parte, la determinación de la
constante universal h realizada por Planck en 1913 dio el valor de $6.41×10−27$ erg seg.
Adoptando el valor de $6.623×10−27$ erg seg, estimar el error absoluto y relativo cometido por
Planck. ¿Cuál de las dos medidas es más precisa?
Errores de redondeo
En el sistema decimal, existen dos formas de redondear un número a t de decimales. En el proceso
denominado truncamiento (o redondeo hacia el cero) donde simplemente se eliminan todos los
números ubicados a la derecha del t-ésimo decimal.La magnitud del error puede ser tan grande
como $10^{-t}$.El redondeo al más cercano (tambien llamado redondeo óptimo) debe elegirse el
número con $s$ decimales que sea más cercano a x.
En una máquina hipotética que utilizara un sistema decimal con dos dígitos decimales, escribe el
número 3.246 usando truncamiento o redondeo.
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
respresentable, 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.
y denotamos por $fl(x)$ a su representación en punto flotante, es decir, al número que
mejor representa a $x$ dentro de las limitaciones del sistema de representación numérica, se
tendrá
$$\left| {x - fl\left( x \right)} \right| \le {2^{{e_x}}} \cdot {2^{ - 53}} \Rightarrow {{\left| {x
- fl\left( x \right)} \right|} \over {\left| x \right|}} \le {2^{ - 53}}$$
entonces el épsilon de la máquina es una cota superior del error de redondeo en el sistema a
punto flotante.
Recuerda que $fl(x)$
Es un número representable en el sistema.
Es el número más cercano a $x$ dentro del rango representable.
Si $x$ está exactamente entre dos números representables, el valor de $fl(x)$ dependerá
de la política de redondeo (por ejemplo, redondeo al más cercano, truncamiento, etc.).
Propagación de errores
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:
Almacena los números $x$ e $y$ en coma flotante.
Después realiza la operación con los números obtenidos.
Finalmente da el resultado en punto flotante.
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.
Consideremos $x=10^{20}$, $y=-10^{20}$, $z=1$ y calcula con Matlab/Octave $x+(y+z)$ y $(x+y)+z$.
Observa que se obtiene respectivamente los valores 0 y 1.
Calcula con Matlab la siguiente operación $e=1 - 3*(4/3 - 1)$.
e =
2.2204e-16
Imaginemos una máquina hipotética con base decimal y 2 dígitos de mantisa. Calcular la
diferencia de 1025 y de 912.3.
El valor exacto es 112.7. El calculo sería de la siguiente manera
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 un número mínimo de operaciones convenientemente dispuestas,
manteniendo a la vez que el tiempo de ejecucion sea aceptable.
Para un análisis detallado sobre los errores en cálculos numéricos, como errores de redondeo,
truncamiento y su propagación, se recomienda consultar los libros
y .
Error por cancelación
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 0.
Dados $x=1+2^{-54}$ e $y=1$, calcula con Matlab $x-y$
Nuestro ordenador establecerá que $x-y=0$ cuando en realidad $x - y \approx 5,5 \cdot {10^{ -
17}}$.
Dados $x=2^{60}+2^6$ e $y=2^{60}$, calcula con Matlab $x-y$.
Matlab establecerá que $x-y=0$ cuando en realidad $x - y =64$.
Considera la función \[f\left( x \right) = \left( {{e^x} - 1} \right) + {\left( {\frac{1}{{\sqrt
{1 + {x^2}} }} - 1} \right)^2}\]Se quiere calcular $f'\left( 1 \right) \approx \frac{{f\left( {1
+ h} \right) - f\left( 1 \right)}}{h}$ para valores de h pequeños. Genera un vector de valores
de $h$ desde $10^{-1}$ hasta $10^{-17}$ observa que valor se obtendría para $f'(1)$. Realiza el
cálculo derivando la función y sustituyendo en el punto 1.
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 long
disp([h' v'])
A partir del valor $k=-16$ se obtiene 0. Realizando la derivada y sustituyendo se obtiene el
valor $ 9.54865532212975587$.
Se sabe que ${e^x} = \sum\limits_{n = 0}^\infty {{{{x^n}} \over {n!}}} $. Calcula una
aproximación de $e^{-40}$ con los 100 primeros términos de la serie.
Calcula el mismo valor teniendo en cuenta que ${e^{ - x}} = {1 \over {{e^x}}}$. ¿Qué
observas?
Observa que como $x< 0$ los términos de la serie cambian de signo y si $|x|$ es grande, el
valor de los sumandos puede devolver un valor pequeño. En algún momento estaremos con un error
de cancelación cuando sumemos "el siguiente" término.
Si se considera el cálculo de la segunda forma, tomando el inverso de $e^{40}$, se evita el
problema de cancelación.
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
En los siguientes ejemplos, se muestra como una reformulación del cálculo puede evitar el
problema de cancelacion en ciertas situaciones.
Considera la función \[f\left( x \right) = \frac{1}{{1 + 2x}} - \frac{{1 - x}}{{1 + x}}\]
Observa que para valores próximos a 0 se puede producir cancelacion. Evalúa esta función en
puntos próximos a 0 directamente y considerando la siguiente expresión equivalente \[f\left( x
\right) = \frac{{2x^2}}{{\left( {1 + 2x} \right)\left( {1 + x} \right)}}\]
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));
k=-1:-1:-17;
h=10.^k;
f1=f(h);
g1=g(h);
format long
disp([h' f1' g1'])
Repite el ejemplo anterior considerando la función $f\left( x \right) = \frac{{1 - \cos x}}{x}$
en los puntos próximos a cero y la expresión equivalente resultado de multiplicar y dividir por
$1+\cos x$.
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 long
disp([h' f1' g1'])
Repite el ejemplo anterior considerando la función \[f\left( x \right) = \sqrt {x + \frac{1}{x}}
- \sqrt {x - \frac{1}{x}} \] en los puntos próximos a uno. Considera como expresión equivalente
el resultado de multiplicar y dividir por el conjugado.
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 long
disp([h' f1' g1'])
Error por desbordamiento
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
overflow o un underflow, respectivamente.
Multiplica el valor máximo representable en punto flotante por dos para ilustrar el concepto de
desbordamiento. Divide el valor mínimo positivo normalizado 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)]);
Inestabilidad numérica
Vamos a ver en el siguiente ejemplo, como los errores de redondeo pueden llevar a resultados
incorrectos si se elige mál el algoritmo. Este problema se llama inestabilidad numérica que
puede evitarse eligiendo algoritmos más adecuados.
Para cada $n \ge 0$ se considera la integral $${I_n} = \int\limits_0^1 {{x^n}sen\pi x\,dx} $$
Sabiendo que
$$\begin{aligned} & {I_o} = {2 \over \pi }\,\,\,,\,\,\,\,{I_1} = {1 \over \pi }\\ &{I_{n + 2}} =
{1 \over \pi }\left[ {1 - {{\left( {n + 2} \right)\left( {n + 1} \right)} \over \pi }{I_n}}
\right]\,\,\,\,\,\,\,\,n \ge 0 \end{aligned}$$ usar esta recurrencia para hallar $I_{30}$ e
$I_{31}$. ¿Son los resultados aproximaciones satisfactorias de la integral? Explica los
resultados y propón una alternativa eficiente para calcular $I_{30}$ e $I_{31}$. Obtén $I_{30}$
e $I_{31}$ por este nuevo procedimiento.
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. 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 long
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}$ es 0.003139287076703. El valor de $I_{31}$ es 0.002950500704051
Observa que el cálculo recurrente llega a obtener un valor de $y_{20}$ negativo, que es absurdo
porque el integrando es positivo. Para evitar este problema, considera después la fórmula de
recurrencia en sentido descendente \[{y_{n - 1}} = \frac{1}{{5n}} - \frac{{{y_n}}}{5}\] Calcula
de esta manera $y_5$.
En el cálculo recurrente de ${y_n} = \frac{1}{n} - 5{y_{n - 1}}$ el error se produce porque se
restan dos números del mismo signo de magnitud similar.
Considerando la recurrencia ${y_n} = \frac{1}{n} - 5{y_{n - 1}}$
clear all
y(1)=1-5*(log(6)-log(5));
format longEng
for k=2:21;
y(k)=1/k-5*y(k-1);
end
disp([[1:40]' y'])
Tomando límites en ${y_n} = \frac{1}{n} - 5{y_{n - 1}}$ se observa que la sucesión $y_n$
converge a 0. Por lo tanto, tomando $y_{60}=0$ y considerando la ley de recurrencia ${y_{n - 1}}
= \frac{1}{{5n}} - \frac{{{y_n}}}{5}$ se puede calcular por ejemplo $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'])
Autoevaluación
**** Incluir preguntas del tema ****
Actividad de autoevaluación
Ejercicios
Calcular los números cuya representación en punto flotante de precisión doble es
Resuelve el siguiente sistema con Matlab \[\left\{ \begin{array}{l} 17{x_1} + 5{x_2} = 22\\
1.7{x_1} + 0.5{x_2} = 2.2 \end{array} \right.\] Una solución obvia es $x_1=x_2=1$ y el
sistema tiene infinitas soluciones (la segunda ecuacion es la primera dividida por 10). ¿Qué
sucede?
Se quiere calcular $p_n$ como $(1/3)^n$ para $n>0$ considerando $p_0=1$ definiendo
$p_n=(1/3)p_{n-1}$ para $n>1$. Otra forma de generar la sucesión es definiendo $p_0=1$ y
$p_1=1/3$ y calculando para $n>2$ \[{p_n} = \frac{{10}}{3}{p_{n - 1}} - {p_{n - 2}}\]
¿Alguno de los dos métodos es inestable?
Evaluar $\sum\limits_{n = 1}^{10000000} {\frac{1}{n}} $ primero en orden usual y luego en
orden opuesto. Explicar las diferencias obtenidas e indique cuál es el resultado más
preciso.
Justifica qué ocurre cuando se evalúa la expresión \[\sqrt {{x^2} + 1} - 1\] para valores
próximos a cero. ¿Cómo se podría hacer más estable? Evalúa ambas expresiones para 20 puntos
equidistantes entre 0 y 0.2.
Indica cómo evaluar con Matlab la expresión \[\frac{1}{{x + h}} - \frac{1}{x}\] para los
valores de $h$ siguientes ${10^{ - 2}},\,{10^{ - 4}}{,...,10^{ - 22}}$ y el valor de $x$ que
se considere. Escribe una reformulación de esta expresión de forma que no se produzca
cancelación. Utiliza el formato longEng.
Calcular una fórmula alternativa para la expresión \[ - \log \left( {x - \sqrt {{x^2} - 1} }
\right)\] que no sea susceptible de pérdida de cifras significativas cuando x es grande y
positiva.
En este capítulo se han visto los fundamentos sobre los errores en cálculo numérico,
destacando cómo los sistemas computacionales representan y manipulan números, y cómo esto
puede originar errores que afectan los resultados. A través de ejemplos históricos, como el
fallo de los misiles Patriot en 1991 y la explosión del cohete Ariane 5 en 1996, se ilustra
la relevancia de comprender la aritmética computacional.
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:
Introducción a los errores en cálculo numérico. Los
errores en los sistemas computacionales son inevitables debido a la naturaleza
finita con la que las computadoras representan números reales y realizan operaciones
matemáticas. Esta limitación se debe al uso de representaciones en punto flotante,
que introducen imprecisiones al redondear o truncar los valores.
Entender estas limitaciones es esencial para abordar problemas numéricos con
precisión.
Errores históricos relevantes. Casos como el fallo de los
misiles Patriot en 1991, donde un error en la representación del tiempo provocó un
cálculo incorrecto, y la explosión del cohete Ariane 5 en 1996, causada por un fallo en
la conversión de datos numéricos, son ejemplos clave que muestran cómo errores
aparentemente pequeños pueden llevar a consecuencias catastróficas en sistemas críticos.
Clasificación de los errores Los errores numéricos se
clasifican principalmente en dos tipos: errores de redondeo y errores de truncamiento.
Los errores de redondeo ocurren al ajustar un número real a la precisión limitada de la
computadora, mientras que los errores de truncamiento surgen al aproximar funciones o
procesos matemáticos, como al usar series finitas para representar funciones infinitas.
Además, los errores algorítmicos aparecen cuando los métodos numéricos no convergen
correctamente hacia la solución deseada.
Minimización de errores Reducir el impacto de los errores
numéricos requiere el desarrollo de algoritmos robustos y estrategias específicas, como
el control de precisión, el análisis de estabilidad de los métodos y la elección de
representaciones numéricas adecuadas. Estas herramientas permiten minimizar la
propagación de errores a lo largo de los cálculos y garantizar la confiabilidad de los
resultados, especialmente en aplicaciones críticas como la ingeniería, la física y las
finanzas.
Capítulo
Resolución aproximada de ecuaciones no lineales
Introduccion
El propósito de los métodos numéricos es encontrar soluciones aproximadas, ya que la mayoría de
las ecuaciones no tienen soluciones exactas o estas no son computacionalmente viables. Este
enfoque iterativo nos acerca a la solución deseada dentro de un margen de error aceptable
.
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, lo que motiva el uso de métodos numéricos para encontrar 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,
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$ la masa del paracaidista.
Simulador de la velocidad de un paracaidista
Imaginemos que queremos 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? La solucion del problema sería obtener el valor de $c$ que hace cero a la siguiente
función
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
pero efectivo: 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 para encontrar soluciones con gran rapidez,
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
NeumannJohn
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.
El estudio de estos métodos no solo es relevante por su aplicabilidad, sino también por su
historia, que refleja el ingenio humano para buscar soluciones a problemas complejos.
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.
Métodos de resolución de ecuaciones no lineales
Raíz o cero de una función
A cualquier solución de la ecuación $f(x)=0$ se le llama cero o raíz de la funcion.
Sea $f$ una función continua en un intervalo abierto $I$, $\alpha$ un punto de $I$ y $m \ge
1$ un número entero, se dice que $\alpha$ es un cero de $f$ de multiplicidad $m$ si
$$f\left( x \right) = {\left( {x - \alpha } \right)^m}q\left( x \right)\,\,\,\,\,\,\,\,x \ne
\alpha $$
Se puede demostrar que si f es $C^m(I)$ y $\alpha$ es un punto de $I$, la función tiene un cero
de multiplicidad $m$ en $\alpha$ si
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.
Algoritmos iterativos.
Para la aproximación al cero de la función $f$ se utilizan métodos iterativos que en el caso de un
paso consiste en:
Partir de un valor inicial $x_0$ que se supone suficientemente próximo a la solución
buscada.
Iterar, es decir, obtener mediante la fórmula $x_{k+1}=G(x_k)$ con $k \ge 0$.
Este proceso depende tanto de la función $G$ como del punto inicial $x_0$. Un algoritmo se debe
dar respuesta a las siguientes cuestiones:
¿es la sucesión convergente? (convergencia del método)
en caso de ser convergente, ¿es el límite de la sucesión la solución de $f(x)=0$?
¿cuál es el número de iteraciones que hay que realizar para conseguir que el error de
aproximación sea menor que una cantidad prefijada? (velocidad de convergencia).
¿cómo evoluciona el error al realizar las iteraciones?
Convergencia de un algoritmo
Se dice que un algoritmo converge globalmente, si la sucesión generada es convergente
hacia una solución de la ecuación cualquiera que sea el punto inicial $x_0$ que tomemos.
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 local.
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.
Velocidad de convergencia
Se distinguen distintos tipos 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 orden o velocidad de convergencia es $p$,
con $p \ge 1$ si existe una constante $C_p \ge 0$ tal que: $$\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,
si $p$=1 y $C_p< 1$, la convergencia es lineal. Si $C_p \ge 1$, la convergencia es
sublineal.
si $p$=2, la convergencia es cuadrática.
si $p$=3, la convergencia es cúbica.
Consideremos la sucesión ${x_n} = {1 \over 6} - {1 \over 3}{\left( { - {1 \over 2}} \right)^n}$
que converge a $1/6$. Veamos que la convergencia es lineal.
Supongamos que $C_p=0.5$ y que $\left| {{x_k} - \overline x} \right| < 0.01$. Considerando
$p=1$ y $p=2$, analizar cómo se aproxima al valor límite considerando el término $k+1$ y $k+2$.
Se dice que la sucesión converge superlinealmente, si existe una sucesión $\left\{
{{\varepsilon _k}} \right\}_{k = 1}^\infty $ convergente a 0 tal que $$\left| {{x_{k + 1}} -
\bar x} \right| \le {\varepsilon _k}\left| {{x_k} - \bar x{\mkern 1mu} } \right|\,\,{\mkern 1mu}
{\mkern 1mu} {\mkern 1mu} \forall k \ge 0$$
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.
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$).
Cuando la convergencia es lineal puede ocurrir que la diferencia entre los valores obtenidos
en dos iteraciones, $|x_{k+1}-x_k|$, sea pequeño y, sin embargo, $x_{k+1}$ esté alejado de
la solución. Considerar la sucesión $x_{k+1}=\frac{x_k} {1.1}+3$ con $x_0=100$ que converge
a 33 y analizar la diferencia entre dos términos consecutivos para $k$ desde 1 hasta 18.
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
xn
|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
Método de la bisección
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.
Método de la bisección
Los pasos del algoritmo son:
Elige $[a, b]$ donde la función $f$ cambia de signo.
Calcula el punto medio de $[a, b]$: $c = (a + b) / 2$.
Si $f(a)$ y $f(c)$ tienen signos opuestos, la raíz está en el intervalo $[a, c]$; de lo
contrario, está en $[c, b]$.
Redefine el intervalo $[a, b]$ como $[a, c]$ si $f(a)$ y $f(c)$ tienen signos opuestos,
o como $[c, b]$ si $f(c)$ y $f(b)$ tienen signos opuestos.
Repite los pasos 2-4 hasta que $[a, b]$ sea lo suficientemente pequeño o se alcance una
precisión deseada.
Este resultado asegura la convergencia del método de la Bisección.
Considerando las funciones y los intervalos de la tabla , calcular las tres primeras iteraciones
del método de la bisección. Comprueba previamente que la función en cada intervalo tiene
únicamente un cero.
$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.
La función $f(x)=\cos(x)-x$ es continua y cumple además que $f(0)\cdot f(2) <0$. Esto
significa que al menos tiene un cero en $[0,2]$. Como $f'(x)=-sen(x)-1 <0$ la función
es decreciente en ese intervalo.
La función $f(x)=sen(x)$ es continua y cumple que $f(-0.5)\cdot f(\pi/2) <0$. Esto
significa que al menos tiene un cero en $[-0.5,\pi/2]$. Como $f'(x)=\cos(x) >0$ la
función es creciente en ese intervalo.
La función $f(x)=x^2-sen(x)$ es continua y cumple además que $f(0.5)\cdot f(2) <0$.
Esto significa que al menos tiene un cero en $[0.5,2]$. Como $f'(x)=2x+\cos(x) >0$ la
función es creciente en ese intervalo.
Se puede comprobar los valores de las tres primeras iteraciones, utilizando la herramienta
interactiva que se muestra a continuación introduciendo la función y los valores de $a$ y
$b$ en cada caso.
Ejemplos del método de la bisección
Dada la función $f(x)=x^3-6x^2+11x-6$, calcular el cero de la función en el intervalo $[3,5]$
tomando como tolerancia el valor $0.001$ y $30$ como número máximo de iteraciones.
Define una función Matlab/Octave que tenga como parámetros los extremos del intervalo $a$ y
$b$, la tolerancia y el máximo de iteraciones. Utilizar esta función para calcular los ceros
de la función $f(x)$ en los intervalos $[-1,0]$ y $[0,1]$.
La herramienta interactiva siguiente construye paso a paso el código de la función que se
solicita, bisectionMethod(a,b,tol,maxiter). Para ello, se debe pulsar en el
botón solucion guiada.
Programación método de la bisección
Para calcular los ceros en cada intervalo bastaría escribir
Modifica el código para que la función cuyo cero se quiere calcular sea también un parámetro,
esto es, programar una función Matlab/Octave:
bisectionMethod2(fun,a,b,tol,maxiter).
Definida bisectionMethod2(fun,a,b,tol,maxiter), para cacular el cero de
$f(x)=x^3-6x^2+11x-6$ en el intervalo $[3,5]$, con una tolerancia $10^{-6}$ y un máximo de 25
iteraciones, se deberá escribir el código siguiente.
Dada la ecuación $f(x)=0$, los pasos a seguir son:
Se despeja $x$ como $x=g(x)$ buscando un punto fijo $\overline x$ de $g$, esto es un
$\overline x$ cumpliendo $\overline x=g(\overline x)$.
Partiendo de $x_0$ se busca la raíz haciendo iteraciones $$x_{k+1}=g(x_k),
k=0,1,...$$ de manera que ${x_k}$ converja a $\overline x$
Método del punto fijo
Veamos algunas condiciones suficientes para que haya un punto fijo en un intervalo y condiciones
suficientes de convergencia a un único punto fijo.
Sea $g$ continua en $[a,b]$ con $g(x)\in [a,b]$ para todo $x \in [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
cualquier punto 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]$.
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.
Considera la función $g(x)=x-x^2/10$ que tiene como punto fijo el 0. Comprueba que la derivada
en las proximidades del 0 es próxima a 1 y analizar los valores que se obtienen con el método
del punto fijo.
A partir de la iteración 89 se estabiliza la sucesión obteniendo 0.739085133215161
Gráficas del ejemplo Ejemplos del método del punto fijo
Calcular los ceros de $f(x)=x-e^{-x^2}$ en $[0,1]$ utilizando el método del punto fijo
considerando las funciones $g(x)=e^{-x^2}$ y $g(x)={-log(x)}^{-1/2}$. Empieza en el punto $0.5$
realizando 100 iteraciones y analiza los resultados.
Con ayuda de Matlab/Octave se han representando las dos funciones comprobando que, mientras que
con la primera función se cumple la condición para asegurar la convergencia al punto fijo,
mientras que en la segunda no. Para la primera función el valor obtenido como punto fijo de $g$
y cero de $f$ es $ 0.652918623454627$.
Acotar las raíces de $p(x)=x^3+3x^2-1=0$. Encontrar las raíces programando el método del punto
fijo con una función Matlab/Octave que tenga como argumentos la función, el punto inicial, la
tolerancia y el número máximo de iteraciones fixedPointMethod(g,x0, tol, maxiter).
Analizamos en primer lugar, los ceros de la función y qué función 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 el intervalo [-3,-2], dividiendo por $x^2$ al ser $x$ distinto de 0, habría que
calcular las raíces de $x+3-1/x^2=0$. Se puede considerar $g(x)=1/x^3-3$. En este caso
la sucesión es convergente ya que $$\left| {g'\left( x \right)} \right| = \left| { -
2{x^{ - 3}}} \right| \le {1 \over 4} < 1$$
En el intervalo [-1,0], la raíz a calcular será de la ecuación $x^2(x+3)=1$. Se puede
considerar $g(x)=-\sqrt{(x+3)^{-1}}$ que aplica el intervalo $[-1,0]$ dentro del
intervalo $[-1,0]$. En este caso $$\left| {g'\left( x \right)} \right| = \left| {{1
\over 2}{{\left( {x + 3} \right)}^{ - 3/2}}} \right| \le {1 \over {2 \cdot 2^{3/2}}}
< 1$$
En el intervalo [0,1], la raíz a calcular será de la ecuación $x^2(x+3)=1$. Se puede
considerar $g(x)=\sqrt{(x+3)^{-1}}$ que aplica el intervalo $[0,1]$ dentro del intervalo
$[0,1]$. En este caso $$\left| {g'\left( x \right)} \right| = \left| {{1 \over
2}{{\left( {x + 3} \right)}^{ - 3/2}}} \right| \le {1 \over {2 \cdot 3^{3/2}}} < 1$$
En la escena interactiva siguiente, se muestra paso a paso cómo generar el código de esta
función: fixedPointMethod(g,x0, tol, maxiter).
Este método aproxima la solución de la ecuación $f(x)=0$ mediante una sucesión $x_k$ contruida,
a partir de un valor inicial $x_0$, utilizando la expresión iterativa: $$x_{k+1} = x_k −
\frac{f(x_k )}{ f'(x_k )} \,\,\,\,\,k ≥ 0$$
La idea es, una vez obtenido $x_k$, calcular $x_{k+1}$ utilizando la aproximación lineal de $f(x)$
que proporciona su recta tangente en 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$$
Método de Newton
Representación del método de Newton
Obtener cinco iteraciones del Método de Newton para encontrar la solución de la ecuación no
lineal $$x-cos(x)=0$$
Consideramos $f(x)=x-cos(x)$. Se cumple que $f'(x)=1+sen(x)$.
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]$.
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
%A partir de la cuarta iteración se
%estabiliza el resultado
Con ayuda de la siguiente herramienta puedes encontrar las primeras iteraciones del Método de
Newton para distintas funciones introduciendo el punto inicial.
Ejemplos del método de Newton
Convergencia
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).
Convergencia local. Sea $f$ un función de $\mathbb R$ en $\mathbb R$ derivable y tal que
que $f'$ es una función continua. Supongamos que $f\left( {\overline x } \right) = 0$ y
$f'\left( {\overline x } \right) \ne 0$. Entonces existe $\varepsilon > 0$ tal que si $${x_o}
\in \left( {\overline x - \varepsilon ,\,\overline x + \varepsilon } \right)$$ la sucesión
generada por el método de Newton, $\left\{ {{x_k}} \right\}_{k = 0}^\infty $ está bien definida
(esto es, está en el intervalo anterior) y $$\mathop {\lim }\limits_{k \to \infty } {x_k} =
{\overline x } $$ Convergencia cuadrática. Además, si $f$ es dos veces derivable con
$f''$ continua, existe una constante $C > 0$ tal que: $$\left| {{x_{k + 1}} - {\overline x }
} \right| \le C{\left| {{x_k} - {\overline x }} \right|^2}\,\,\,\,\,\,\forall k \ge 0$$
Calcular las raíces de $e^x+2^{-x}+2cos(x)=6$ en $[0,1]$
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)
%repetimos la instrucción para iterar hasta conseguir
%16 decimales
En la siguiente tabla se muestran las primeras seis iteraciones.
$x_0$
1.500000000000000
$x_1$
1.956489721124210
$x_5$
1.829383601933849
$x_2$
1.841533061042061
$x_6$
1.829383601933849
$x_3$
1.829506013203651
$x_4$
1.829383614494166
Observad 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:
Obtener la raíz del polinomio $p\left( x \right) = {x^3} + 3{x^2} + 2$ que se encuentra en el
intervalo $[-4 -2]$.
Podemos considerar como punto de inicio el punto medio del intervalo $x_0=-3.25$.
p=[1 3 0 2]
%polyder deriva el polinomio
dp=polyder(p)
%polyval(p,a) evalúa el polinomio en a
%método de Newton
%intervalo [-4,-2]
polyval(p,-4);polyval(p,-2)
x=-3.25
x=x-polyval(p,x)/polyval(dp,x)
%repetimos este último paso para iterar, sol=-3.195823345445647
Si se repite el proceso con $x_0=1$, no converge, las iteraciones oscilan entre negativo y
positivo. Esto es porque el método de Newton converge localmente.
Test de parada
Se pueden considerar los siguientes tests:
que la diferencia entre dos iteraciones consecutivas sea menor que un cierto valor
$\epsilon$ $$\left| {{x_{k + 1}} - {x_k}} \right| < \varepsilon \,\,\,\,\text{error
absoluto}$$ $$\left| {{x_{k + 1}} - {x_k}} \right| < \varepsilon |{x_{k+1}}|
\,\,\,\,\text{error relativo}$$ Combinando estos dos errores, sería $$\left| {{x_{k +
1}} - {x_k}} \right| < \varepsilon \max \left\{ {1,\left| {{x_{k + 1}}} \right|}
\right\}$$
dado que se busca $f(x)=0$ se parará cuando $f(x_{k+1})$ sea pequeño, lo que
comprobaremos con el test $$|f(x_{k+1})| < \epsilon_f$$ Para funciones sencillas un
valor habitual es $\varepsilon_f = \varepsilon _M^{3/4}$.
Construye una función Matlab/Octave que permita aplicar el método de Newton para una función
dada considerando como parámetros el punto inicial, el máximo de iteraciones y la tolerancia.
Con la herramienta siguiente se muestra cómo generar paso a paso el código Matlab/Octave de la
función newtonMethod(x0, tol, maxiter) que utiliza como método de parada que la
diferencia entre dos iteraciones consecutivas sea menor que la tolerancia fijada.
Programación guiada del método de Newton
Método de Newton combinado con bisección
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:
Si $x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k}$ está en $[a,b]$ se acepta, en caso contrario se
utiliza bisección, es decir, $x_{k+1}=\frac{a_k+b_k}{2}$
Se considera $[a_{k+1}, b_{k+1}]$ como el intervalo $[a_{k}, x_{k+1}]$ o $[x_{k+1},
b_{k}]$ según haya cambio de signo de $f$ en los extremos.
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:
Encontrar la raiz de $e^x+2^{-x}+2cos(x)-6=0$ en $[1,2]$ utilizando el método de Newton
combinado con Bisección.
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]
$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]$
Iteración 1
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$
Iteración 2
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$
Iteración 3
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$
Iteración 4
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$
Iteración 5
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 siempre Newton.
Encontrar la raiz de $e^x+2^{-x}+2cos(x)-6=0$ que se encuentra entre $[-3,-2]$ utilizando el
método de Newton combinado con bisección.
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$
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.
Escribir el código Matlab/Octave para programar el método de la secante combinado con bisección
a partir de dos puntos imiciales donde la función cambie de signo, la tolerancia y el máximo de
iteraciones permitidas: NewtonBisecciónMethod(x0, tol, maxiter)
Con la herramienta siguiente se muestra cómo generar paso a paso el código Matlab/Octave de
la función NewtonBisecciónMethod(x0,tol, maxiter) que utiliza como método de
parada que la diferencia entre dos iteraciones consecutivas sea menor que la tolerancia
fijada.
Programación guiada del método de Newton combinado con bisección
Método de la secante
Este método se basa en considerar que la recta secante aproxima a la recta tangente en
intervalos pequeños.
Este método parte de dos aproximaciones iniciales $x_o$ y $x_1$. eEn el paso $k+1$, se calcula
$x_{k+1}$ usando $x_k$ y $x_{k-1}$ como la intersección con el eje X de la recta secante, que
pasa por los puntos $(x_{k-1}, f(x_{k-1})$ y $(x_{k}, f(x_{k})$. Es decir, $${x_{k + 1}} = {x_k}
- {{f\left( {{x_k}} \right)} \over {{m_k}}}\,\,\,\,\,\,\,\,\,\,\,\,{m_k} = {{f\left( {{x_k}}
\right) - f\left( {{x_{k - 1}}} \right)} \over {{x_k} - {x_{k - 1}}}}$$
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)$$
Método de la secante
Representación del método de la secante
Calcular, utilizando el método de la secante, las 5 primeras iteraciones para calcular el cero
de la función $f(x)=cos(x)$ considerando como $x_0=0$ y $x_1=1$.
Con la ayuda del siguiente interactivo puede calcular las tres primeras iteraciones del método
de la secante comenzando en los puntos $0$ y $2$.
Ejemplos del método de la secante
Escribir una función Matlab/Octave para aplicar el método de la secante que tenga como
argumentos los puntos $x_0$, $x_1$, la tolerancia y el número máximo de iteraciones permitidas.
Pulsando en el botón solución guiada del interactivo siguiente, se muestra la
programación de la función secantMethod(x0, x1, tol, maxiter).
Programación guiada del método de la secante
Convergencia
El orden de convergencia de este método, en un punto cercano a la solución, es el número áureo,
$\phi = \frac{{1 + \sqrt 5 }}{2} \approx 1.618$, por lo que se trata de una convergencia superlineal
inferior a la del método de Newton-Raphson. En caso de que la aproximación inicial sea demasiado
lejana o la raíz no sea simple, este método no asegura la convergencia y tiene un comportamiento
similar al de Newton-Raphson.
Método de la secante combinado con bisección
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:
Si $x_{k+1}=x_{k}-\frac{f(x_{k})}{m_k}$ está en $[a_k,b_k]$ se acepta, en caso contrario
se utiliza bisección, es decir, $x_{k+1}=\frac{a_k+b_k}{2}$
Se modifica el intervalo $[a_{k+1},b_{k+1}]$ tomando o bien $[a_k,x_{k+1}]$ o bien
$[x_{k+1},b_k]$ según que haya cambio de signo de $f$ en los extremos del intervalo.
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.
La función \(erf(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt\) está definida para \(x \geq
0\). Determina el valor de \(x\) tal que \(erf(x) = 0.5\). Nota: Utilizar la función de Matlab
\(erf\) para evaluar las integrales.
Consideramos en primer lugar un intervalo donde la función cambie de signo, por ejemplo
$[a,b]=[0,1]$. Definimos $x_0=0$ y $x_1=0$ para empezar a iterar el método de la secante. Se
calcula la siguiente iteración $x$ utilizando el método de la secante y se comprueba si está en
$[0,1]$.
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 respectiamente $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 respectiamente $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 respectiamente $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 respectiamente $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 respectiamente $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$.
Dada la matriz $$A (x)= \begin {pmatrix} x & 2 & { - 3} & 5 & 0 & {{x^2}} \\ 1 & 1 & 1 & 1 & 1 &
1 \\ 0 & x & x & 0 & 5 & { - 2} \\ 1 & { - 1} & 1 & { - 1} & {2x} & 0 \\ 1 & 0 & 1 & 0 & {{x^2}}
& { - 2} \\ 0 & 0 & 1 & 0 & {\cos (x)} & { - x} \\ \end{pmatrix}$$ determinar un valor de \(x\)
tal que \(-2\) sea un valor propio.
Utilizando el método de la secante, dada la dificultad de calcular la derivada, se debe calcular
el valor de x de manera que $f(x)=det(A(x)+2I)=0$ siendo I la matriz identidad 6x6.
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.
Aplicando el método de la secante bisección calcular el cero de la función
$f(x)=e^x+2^{-x}+2cos(x)-5$ en el intervalo $[-5,0]$ con una tolerancia $10^{-7}$.
Se comienza viendo si hay cambio de signo de $f$ en $[-5,0]$. Como $f(-5)f(0)<0$, se
considera $[a,b]=[-5,0]$, y se toma $x_0=-5$, $x_1=0$ los puntos iniciales para aplicar el
método de la secante.
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 cálculo 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.
La solución con la tolerancia indicada es $-2.764844125871619$.
Escribir el códgo Matlab/Octave para programar el método de la secante combinado con bisección a
partir de dos puntos imiciales donde la función cambie de signo, la tolerancia y el máximo de
iteraciones permitidas: SecanteBisecciónMethod(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 SecanteBisecciónMethod(a, b, tol, maxiter) que utiliza como método
de parada que la diferencia entre dos iteraciones consecutivas sea menor que la tolerancia
fijada.
Programación guiada del método de la secante combinado con bisección
Raíces de un polinomio
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:
Teorema Fundamental del Álgebra: Un polinomio de grado $n$ tiene exactamente $n$
raíces en el conjunto de los números complejos, donde cada una de ellas debe contarse
con su multiplicidad.
Si un polinomio tiene coeficientes reales, siempre que tenga una raiz compleja no real,
también tendrá como raíz su conjugada.
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 vez calculada una raíz de un polinomio por el Método de Newton, nos planteamos cómo
encontrar otra evitando que el método de Newton vuelva a converger a la raíz ya calculada.
Una primera idea sería elegir otro punto de inicio, sin embargo, esto no garantizaría
que el método converja necesarimente a otra raíz.
Otra posibilidad, podría ser considerar el polinomio resultado de
dividir $p(x)$ entre $x-\alpha$ siendo $\alpha$ la raíz calculada previamente y aplicar
el método de Newton a este nuevo polinomio. Sin embargo, este método, conocido como
deflacción, presenta también problemas:
La operación de división conlleva errores de redondeo.
El método de Newton nos devuelve una aproximación de la raíz exacta, $\alpha *$ por lo que
el cociente $\frac{p(x)}{x-\alpha}$ no coincide con $\frac{p(x)}{x-\alpha *}$ y, en
consecuencia, no tienen por qué tener las mismas raíces.
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 Método de Maehly.
La idea de este método es, una vez calculadas las raíces, $$\left\{ {{\alpha _1},{\alpha
_2},...,{\alpha _m}} \right\}$$ del polinomio $p(x)$, aplicar el método de Newton al polinomio
$${f_m}\left( x \right) = {{p\left( x \right)} \over {\left( {x - {\alpha _1}} \right)\left( {x
- {\alpha _2}} \right)...\left( {x - {\alpha _m}} \right)}}$$ Nota: Sin hacer la división
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}}}} }}$$
El método consiste en aplicar el método de Newton aplicado a la función $f_m$ considerando:
$$\begin{aligned} & {x_o} \in \mathbb C \\ & x_{k + 1} = {x_k} - {\rm{ }}{{p\left( x_k \right)}
\over {p'\left( x_k \right) - p\left( x_k \right)\sum\limits_{j = 1}^m {{1 \over {x_k - {\alpha
_j}}}} }}\,\,\,\,\,\, k \ge 0 \end{aligned}$$
En los ejercicios prácticos usaremos un programa del profesor Eduardo Casas (ficheros Matlab) para realizar una iteraciòn del método de Maehly (ficheros Matlab)
Puedes descargar el fichero Matlab itraiz.m
.
Para aplicar este método, se debe tener en cuenta:
A la hora de
calcular la primera raiz, el algoritmo se reduce a aplicar el método de Newton
al polinomio: $$\begin{aligned} & {x_o} \in \mathbb C \\ & x_{k + 1} = {x_k} - \frac{p(x_k)}
{p'(x_k)} \,\,\,\,\,\, k \ge 0 \end{aligned}$$
Si los coeficientes del polinomio son todos reales y $x_0$ es real, la sucesion ${x_k}$ será
de números reales y no podrá converger a una raíz compleja. Por lo tanto,
deberá tomarse $x_0$ complejo si se quiere encontrar una raíz compleja.
El método de Newton converge cuadráticamente hacia una raíz $\overline x$ si
$p'(\overline x)\ne0$.
Si se observa que los valores de $x_k$ convergen pero lentamente, puede deberse a
que se está aproximando a una raíz múltiple o hay raíces muy próximas.
En este caso se aplica el método de Newton a la derivada del polinomio comenzando en el
último punto $x_k$ calculado en las iteraciones sobre $p$.
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.
Si el
polinomio tiene coeficientes reales y $\alpha=a+bi$ es una raiz compleja entonces su
conjugada $\overline \alpha = a - bi$ también es raíz del polinomio . Además, si $\alpha$ tuviera multiplicidad $m$, entonces $\overline \alpha$ tendría la
misma multiplicidad.
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.
Calcular todas las raíces del polinomio $p\left( x \right) = 25{x^8} + 85{x^7} + 181{x^6} +
238{x^5} + 226{x^4} + 142{x^3} + 61{x^2} + 13x + 1$.
Calcular todas las raíces del polinomio $p\left( x \right) = 3{x^7} + 5{x^6} + 21{x^5} + 35{x^4}
+ 48{x^3} + 80{x^2} + 36x + 60$.
Autoevaluación
**** Incluir preguntas del tema ****
Actividad de autoevaluación
Ejercicios
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:
\(g(x) = e^{-x^2}\)
\(g(x) = \sqrt{\ln(x)}\)
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.
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\)?
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.
Consideremos la función \(f(x) = \sin(0.2x) - 16 + 1.75x\). 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.
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\).
La función \(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}\).
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\).
Implementa el método de la bisección para resolver la ecuación no lineal \(f(x) =
\frac{1}{\pi} \int_0^x \frac{\sin(t)}{1+t^2} dt - \frac{3}{4} = 0\) con un error
absoluto \( arepsilon_a = 10^{-10}\).
El cápítulo se centra en la resolución aproximada de ecuaciones escalares no lineales. Se
analizan métodos iterativos para encontrar soluciones de ecuaciones de la forma $f(x)=0$,
destacando técnicas como el método de bisección, el método de Newton-Raphson y el método de
la secante.
El capítulo aborda la convergencia de estos métodos, la velocidad con la que se
aproximan a la solución y las condiciones necesarias para su aplicación efectiva.
Además, se discuten estrategias para la localización y separación de raíces, enfatizando
la importancia de una elección adecuada del intervalo inicial y de las aproximaciones
iniciales para garantizar la eficiencia y precisión en la obtención de soluciones
numéricas.
Los aspectos clave en el estudio de los métodos numéricos son:
Introducción Se plantea la necesidad de resolver
ecuaciones no lineales de la forma $f(x)=0$ que, en muchos casos, no admiten
soluciones analíticas exactas, lo que motiva el uso de métodos numéricos para
obtener soluciones aproximadas..
Generalidades Se definen conceptos fundamentales como
raíces o ceros de una función, multiplicidad de una raíz y las fases del proceso de
resolución: localización, separación y aproximación de raíces.
Método de la bisecciónSe presenta este método como una técnica
que, partiendo de un intervalo donde la función cambia de signo, divide repetidamente el
intervalo en mitades para aproximar la raíz con una precisión deseada.
Método de Newton-Raphson Se introduce este método iterativo
que utiliza la derivada de la función para aproximar la raíz, ofreciendo una
convergencia rápida bajo condiciones adecuadas, pero que requiere una buena estimación
inicial y que la derivada no sea nula en las proximidades de la raíz.
Método de la secante
Se describe este método que, similar al de Newton-Raphson, no requiere el cálculo
explícito de la derivada, sino que aproxima la pendiente mediante diferencias finitas,
siendo útil cuando la derivada es difícil de calcular.
Convergencia de los métodos Se analizan las condiciones bajo
las cuales los métodos iterativos convergen a una solución, la velocidad de convergencia
y cómo la elección de las aproximaciones iniciales influye en la eficiencia del proceso.
Capítulo
Aproximación de funciones por polinomios
Introduccion
A menudo se dispone de una tabla de datos obtenida por muestreo o experimentación y se supone
que los datos corresponden a los valores de una función $f$ desconocida o, aunque sea conocida,
se precisa cambiarla por una más sencilla de calcular. Un caso particular de ajuste de curvas
consiste en construir un polinomio $P(x)$ que pase por los puntos de esta tabla de valores.
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 estirmar 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.
Comparación interpolación vs. regresión
Interpolación de Lagrange
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.
DEFINICIÓN. Dados $n+1$ puntos distintos ${x_0,x_1,...,x_n}$ en el intervalo $[a,b]$, el
polinomio de LagrangeJoseph
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
nodos de interpolación.
TEOREMA. El polinomio de interpolación de Lagrange posee una única solución y viene dado
mediante la fórmula: $$p\left( x \right) = \sum\limits_{i = 0}^n {f\left( {{x_i}}
\right){l_i}\left( x \right)} $$ donde $\left\{ {{l_i}} \right\}_{i = 0}^n$ son los
polinomios de base de Lagrange asociados a los nodos $\left\{ {{x_i}} \right\}_{i =
0}^n$ que vienen dados por las expresiones: $${l_i}\left( x \right) = \prod_{ j = 0 \,\,\,\, j
\ne i \,\,\,\,}^n {{{x - {x_j}} \over {{x_i} - {x_j}}}} $$
El video siguiente muestra un ejemplo práctico con tres nodos calculando los polinomios base de
Lagrange.
Método de Lagrange
Calcular los polinomios de base de Lagrange asociados a los nodos $x_0=-1$, $x_1=1$ y $x_2=2$.
Calcula el polinomio que interpola en dichos nodos a la función $f$ real en $[-1,2]$ que
satisface $f(x_0)=2$, $f(x_1)=1$ y $f(x_2)=1$.
Los datos son: $${x_o} = - 1\,\,\,\,\,\,\,\,\,\,{x_1} = 1\,\,\,\,\,\,\,\,\,\,\,{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 $${l_o}\left(
x \right) = {1 \over 6}{x^2} - {1 \over 2}x + {1 \over 3}$$ $${l_1}\left( x \right) = -{1 \over
2}{x^2} + {1 \over 2}x + 1 \,\,\,\,\,{l_2}\left( x \right) = {1 \over 3}{x^2} - {1 \over 3}$$ 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_0}\left( x \right) =$$ $$={1 \over 6}{x^2} - {1 \over 2}x + {4 \over 3}$$
Con ayuda de la siguiente herramienta interactiva calcula el polinomio de Lagrange con 5 nodos
introduciendo los datos que consideres. Los puntos deben ser distintos.
Polinomios de Lagrange. Autor: William
Aunque la fórmula del teorema anterior posee un interés teórico, no es la forma eficiente de
calcular el polinomio interpolador. Un método sencillo y más rápido es el que 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 $\alpha_i$ son distintos entre sí.
y el problema de interpolación de Lagrange tiene una única solución.
Calcular con Matlab/Octave el polinomio interpolador de grado 2 de una función que verifica
$f(-1)=2$, $f(1)=1$ y $f(2)=1$.
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)
Interpolar con Matlab $f(x)=e^x$ en el intervalo $[0,1]$ con polinomios de grado 5, 10, 15.
Evaluar la función y el polinomio en 500 puntos del intervalo y ver cuál es el máximo de la
diferencia entre la función y el polinomio en esos puntos.
Con los cálculos que se muestran en el siguiente fichero 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, este valor se
puede obtener dado que se conoce la función.
Error en la interpolación
COTA DEL ERROR. Supongamos que la función $f$ es $n+1$ veces derivable en el intervalo
$[a,b]$ y que $f^{(n+1}$ es una función continua. Una cota del error de interpolación de
Lagrange viene dado por la expresión $$\left| {f\left( x \right) - p\left( x \right)} \right|
\le {{\mathop {\max }\limits_{c \in \left[ {a,b} \right]} \left| {{f^{(n + 1}}\left( c \right)}
\right|} \over {\left( {n + 1} \right)!}}\left| {\prod_{j = 0}^n {\left( {x - {x_j}} \right)} }
\right|$$
Si consideramos una función $f$ en el intervalo $[-1,2]$ cumpliendo $f(-1)=2$, $f(1)=1$ y
$f(2)=1$, calcular una cota del error en $[-1,2]$ del polinomio interpolador de Lagrange que
aproxima a la función en esos puntos en términos del máximo de $|f^{(n+1}(x)|$ siendo $n$ esel
grado del polinomio, en este caso 2.
Si llamamos M a la cota de la derivada tercera de $f$ en el intervalo, se tendrá que $$\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 es entonces: M*2.112611790922380/6.
Nodos de Chebyshev
La elección de nodos que permite minimizar el valor de la siguiente expresión que aparece en la
cota de error de interpolación de Lagrange
son los puntos de ChebyshevPafnuti
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 polinomio de Chebyshev que se
define de forma recurrente: $${T_o}\left( x \right) = 1\,\,\,\,\,\,\,\,\,\,\,{T_1}(x) = x$$
$${T_k}(x) = 2x{T_{k - 1}}\left( x \right) - {T_{k - 2}}\left( x \right)\,\,\,\,\,\,\,k \ge 2$$
Nodos de Chebyshev. Los nodos de Chebyschev en $[a,b]$ se calculan como $${x_i} = {{a +
b} \over 2} + {{b - a} \over 2}\cos \left( {{{\left( {2i + 1} \right)\pi } \over {2n + 2}}}
\right)\,\,\,\,\,\,\,0 \le i \le n$$
Considerando el intervalo $[-3,3]$ y que $n=8$ es el grado del polinomio de interpolación,
calcular con Matlab/Octave los nodos de Chebyshev (9 nodos) y dibujarlos.
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));
Representacón de 9 nodos de Chebyshev
Considerando estos nodos se cumple: $$\max \left| {\prod\limits_{j = 0}^n {\left( {x - {x_j}}
\right)} } \right| = {{{{\left( {b - a} \right)}^{n + 1}}} \over {{2^{2n + 1}}}}$$ por lo que la
cota del error de Lagrange será: $$\left| {f\left( x \right) - p\left( x \right)} \right| \le
{{\mathop {\max }\limits_{c \in \left[ {a,b} \right]} \left| {{f^{(n + 1}}\left( c \right)}
\right|} \over {\left( {n + 1} \right)!}} {{{{\left( {b - a} \right)}^{n + 1}}} \over {{2^{2n +
1}}}} $$
Calcular ...
Interpola la función $f(x)=log(x)$ en los nodos ${0.99, 1, 1.1}$ utilizando polinomios de
Lagrange y obtén aproximaciones de $log(x)$ en los puntos $0.95678$ y $1.03597$. Compara los
valores obtenidos con los exactos.
Aproxima también $f$ en los mismos puntos mediante el polinomio de grado 2 de Lagrange
considerando los nodos de Chebyshev en el intervalo $[0.9,1.1]$.
Interpolación de Hermite
DEFINICIÓN. Sea $f$ una función real definida en el intervalo $[a,b]$, ${x_0,x_1,...,x_k}$ son
$k+1$ puntos distintos del intervalo $[a,b]$ y ${\alpha_0,\alpha_1,...,\alpha_k}$ $k+1$ enteros
no negativos. El problema de interpolación de Hermite consiste en determinar un
polinomio de grado menor o igual que $n=k+\alpha_0+\alpha_1+...+\alpha_k$ satisfaciendo
$$p\left( {{x_o}} \right) = f\left( {{x_o}} \right){\mkern 1mu} \,\,\,\,\,p'\left( {{x_o}}
\right) = f'\left( {{x_o}} \right)\,\,\,...\,\,\,{p^{({\alpha _o}}}\left( {{x_o}} \right) =
\,{f^{({\alpha _o}}}\left( {{x_o}} \right)$$ $$p\left( {{x_1}} \right) = f\left( {{x_1}}
\right){\mkern 1mu} \,\,\,\,\,p'\left( {{x_1}} \right) = f'\left( {{x_1}}
\right)\,\,\,...\,\,\,{p^{({\alpha _1}}}\left( {{x_1}} \right) = \,{f^{({\alpha _1}}}\left(
{{x_1}} \right)$$ $$...$$ $$p\left( {{x_k}} \right) = f\left( {{x_k}} \right){\mkern 1mu}
\,\,\,\,\,p'\left( {{x_k}} \right) = f'\left( {{x_k}} \right)\,\,\,...\,\,\,{p^{({\alpha
_k}}}\left( {{x_k}} \right) = \,{f^{({\alpha _k}}}\left( {{x_k}} \right)$$ Nota: En el caso
particular en que $\alpha_i=1$ para cada $i$, entonces el problema se denomina
interpolación de Hermite clásica.
Es importante que en cada nodo $x_i$ que se imponga una condición sobre la derivada de orden
$\alpha_i$, se imponga tambien una condición sobre cada derivada inferior a $\alpha_i$ y el
valor de la función. En otro caso, los problemas serían de interpolación de Birkhoff que no
tienen solución en todos los casos.
TEOREMA. El problema de interpolación de Hermite posee una única solución.
En el siguiente video se muestra un ejemplo de interpolación de Hermite.
Ejemplo de interpolación de Hermite con dos nodos
Con la siguiente herramienta interactiva puedes ver cómo varía el polinomio según se modifiquen
los cuatro nodos y los valores de las derivadas en cada uno de ellos.
Método de Hermite
Error de interpolación
COTA DEL ERROR. Supongamos que la función $f$ es $n+1$ veces derivable en el intervalo
$[a,b]$ y que $f^{n+1}$ derivable en el intervalo $[a,b]$ y que $f^{(n+1}$ es una función
continua, entonces una cota del error de interpolación de Hermite viene dada por la expresión:
$$\left| {f\left( x \right) - p\left( x \right)} \right| \le {{\mathop {\max }\limits_{c \in
\left[ {a,b} \right]} \left| {{f^{(n + 1}}\left( c \right)} \right|} \over {\left( {n + 1}
\right)!}}\left| {\prod\limits_{j = 0}^{} {{{\left( {x - {x_j}} \right)}^{{\alpha _j} + 1}}} }
\right|$$
Calcular ...
Calcular ...Calcular ...
Mínimos cuadrados
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.
Función de Runge
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.
La idea consiste en encontar un polinomio de grado $k$ a partir de un conjunto de $m$ puntos de
forma que el polinomio tenga grado inferior al que se podría obtener pasando por todos ellos y
se minimice el error que existe entre el valor del polinomio y las ordenadas de los puntos.
En el siguiente video y con la siguienta herramienta interactiva se puede ver cómo hacer el
ajuste a una recta.
Ajuste a una recta por mínimos cuadrados Ajuste a una recta por mínimos cuadrados
Se escribe a continuación la formulación del método de mínimos cuadrados.
DEFINICIÓN. Dados $n+1$ puntos distintos ${x_0,x_1,...,x_n}$ en el intervalo $[a,b]$, $n+1$
números reales positivos $\left\{ {{\omega _0},{\omega _1},...,{\omega _n}} \right\}$, $f$ una
función definida en $[a,b]$ y un entero $1 \le k \le n$, el
problema de mínimos cuadrados consiste en resolver $$\mathop {{(P) \,\,\, \,
\rm{minimizar}}}\limits_{p \in {P_k}} \sum\limits_{j = 0}^n {{\omega _j}{{\left( {p\left(
{{x_j}} \right) - f\left( {{x_j}} \right)} \right)}^2}} $$ donde $P_k$ denota el espacio de los
polinomios de grado menor o igual $k$.
La razón de incluir los pesos $\omega_i$ se justifica, por ejemplo, para ajustar el polinomio a
los valores más precisos de $f$ o para obtener más precisión en unas zonas que en otras de un
intervalo.
Formulación algebraica del problema
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$
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}$$
Formular algebraicamente el problema de ajustar un polinomio de grado 2 de la forma
$y=a_o+a_1x+a_2x^2$ a un conjunto de observaciones: $(x_0,y_0)$, ($x_1,y_1)$,...,$(x_n,y_n)$ con
$n \gt 2$ en el sentido de mínimos cuadrados.
El problema de mínimos cuadrados es encontrar $a_o$, $a_1$, $a_2$ de forma que sea mínima la
siguiente expresión $$\sum\limits_{i = 0}^n {{{\left( {{y_i} - \left( {{a_o} +
{a_1}{x_i}+{a_2}{x_i^2}} \right)} \right)}^2}} $$ Se considera $$A = \begin{pmatrix} {{x_0^2}}
&{{x_0}} & 1 \\ {{x_1^2}} &{{x_1}} & 1 \\ {{...}} & {{...}}& {{...}}\\ {{x_n^2}}&{{x_n}} & 1 \\
\end{pmatrix} \,\,\, x = \begin{pmatrix} {{a_2}} \\ {{a_1}} \\ {{a_0}} \\ \end{pmatrix} \,\,\,\,
b= \begin{pmatrix} {{y_0}} \\ {{y_1}} \\ {{...}} \\ {{y_n}} \\ \end{pmatrix}$$ El problema se
reduce a buscae el vector $x$ de ${\mathbb R}^3$ solución del problema $$\min \left\| r \right\|
= \min \left\| {Ax - b} \right\|$$ La siguiente figura muestra gráficamente un ejemplo con cinco
puntos.
Ajuste mínimos cuadrados
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
Llamamos $a_1=k$ y $a_2=-kL$. El sistema resultante de escribir la ley de Hooke es $$ 6{a_1} +
{a_2} = 3 $$ $$ 7{a_1} + {a_2} = 8 $$ $$ 10{a_1} + {a_2} = 12 $$ En forma matricial, $Ax=b$
siendo $$ A= \begin{pmatrix} 6& 1 \\ 7& 1 \\ 10& 1 \\ \end{pmatrix} \,\,\,\,\, x=
\begin{pmatrix} a_1 \\ a_2\\ \end{pmatrix}\,\,\,\,\, b= \begin{pmatrix} 3 \\ 8\\ 12\\
\end{pmatrix}$$
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\|$.
Resolución del problema mediante la factorización QR
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
Realizamos la factorización $QR$ de $A$ donde $Q$ es ortogonal, es decir, ${Q^T}Q =
Q{Q^T} = I$ y tiene dimensiones $(n+1)\text{x}(n+1)$ y $R$ es una matriz triangular superior
$(n+1)\text{x}(k+1)$
Consideramos en $Q$ y $R$ las submatrices $Y$ y ${\widehat R}$:
$Q=[Y Z]$ con $Y$ matriz $(n+1)\text{x}(k+1)$ y $Z$ matriz $(n+1)\text{x}(n-k)$
$R = \begin{pmatrix} {\widehat R} \\ 0 \\ \end{pmatrix}$ siendo ${\widehat R}$ es
una matriz triangular superior$(k+1)\text{x}(k+1)$ y O es una matriz formada por
ceros de orden $(n-k)\text{x} (k+1)$
Los coeficientes del polinomio solución
$$a = {\left( {{a_0},{a_1},...,{a_k}} \right)^T}$$ se obtienen
resolviendo el sistema de ecuaciones lineales $$\widehat Rx = {Y^T}b$$
y el residuo es $${\mathop{\rm Re}\nolimits} s = \left\| {{Z^T}b} \right\|$$
Veamos la justificación del punto 3. Como Q es ortogonal se cumple
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.
Hay que hacer notar que la matriz ${\widehat R}$ es inversible debido a que los nodos $\left\{
{{x_j}} \right\}_{j = 0}^n$ son todos distintos. Por lo tanto, la solución del problema de
mínimos cuadrados (Q) existe y es única.
Considerando el ejemplo anterior de la ley de Hooke, se seguirá el método visto para considerar
la recta que resuelve el problema de mínimos cuadrados con los siguientes datos
$h$
6
7
10
$F$
3
8
12
En este caso, $$ A= \begin{pmatrix} 6& 1 \\ 7& 1 \\ 10& 1 \\ \end{pmatrix} \,\,\,\,\, x=
\begin{pmatrix} a_1 \\ a_2\\ \end{pmatrix}\,\,\,\,\, b= \begin{pmatrix} 3 \\ 8\\ 12\\
\end{pmatrix}$$ La descomposición $QR$ de la matriz $A$ se obtiene con Matlab de la forma
Calcular los polinomios de base de Lagrange asociados a los nodos \(x_0, x_1, x_2\). Calcula
el polinomio que interpola en dichos nodos a la función \(f\) definida en \([-1,2]\) y que
satisface:
Interpolar con MATLAB la función \(f(x) = e^x\) en el intervalo \([0,1]\) con polinomios de
grado 5, 10, y 15. Evaluar la función y el polinomio en 500 puntos del intervalo y encontrar
el máximo de la diferencia entre la función y el polinomio en esos puntos.
Interpola la función \(\log(x)\) en los nodos \({0.9, 1, 1.1}\) utilizando el polinomio
interpolador de Lagrange para aproximar en los puntos \(0.956678\) y \(1.03597\). Luego,
realiza la aproximación con los nodos de Chebyshev en el intervalo \([0.9, 1.1]\).
Tomando como nodos \(x_0\) y \(x_1\), obtén el polinomio de interpolación de Hermite
\(p(x)\) que interpola la función \(f(x) = \sin(x)\) en el intervalo \([0,1]\). Calcula
también una cota del error de interpolación.
Tomando como nodos \(x_0, x_1, x_2\), obtén el polinomio de interpolación de Hermite que
interpola la función \(f(x) = e^x\) en los puntos \(x = 0, 1, \pi/2\).
Consideremos un muelle colgado verticalmente dentro de un cilindro abierto por los extremos
a altura 0. Se desea determinar la longitud del muelle \(L\) y la constante del muelle \(K\)
usando la ley de Hooke \(F = k(h - L)\), a partir de los siguientes datos:
Dada la función de Runge \(f(x) = \frac{1}{1+x^2}\), halla los polinomios interpoladores de
Lagrange de grado 10 y 20 respectivamente en puntos igualmente espaciados en \([-5, 5]\).
Representa gráficamente la función y los polinomios interpoladores.
Dada la función de Runge \(f(x) = \frac{1}{1+x^2}\), halla el polinomio de grado menor o
igual a 20 que mejor aproxima la función en el sentido de mínimos cuadrados en \([-5, 5]\),
tomando como datos los valores de \(f(x)\) en nodos igualmente espaciados.
Repite el ejercicio con la función \(f(x) = e^{-x^2}\).
Este capítulo, se centra en la aproximación de funciones de una variable mediante
polinomios, presentando dos técnicas principales: la interpolación y la aproximación por
mínimos cuadrados. En la interpolación, se destacan los métodos de Lagrange y Hermite, que
permiten construir polinomios que pasan exactamente por un conjunto de puntos dados. Por
otro lado, la aproximación por mínimos cuadrados se utiliza cuando se dispone de una gran
cantidad de datos y se busca un polinomio que minimice el error global, sin necesariamente
pasar por todos los puntos.
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:
Introducción Se presenta la necesidad de aproximar
funciones desconocidas o complejas mediante polinomios, utilizando técnicas como la
interpolación y la aproximación por mínimos cuadrados.
Interpolación de Lagrange Se introduce el polinomio de
Lagrange como una herramienta para construir un polinomio de grado n que pasa por
n+1 puntos dados, garantizando una única solución y proporcionando una fórmula
explícita para su cálculo.
Método de mínimos cuadrados
Se presenta la técnica de aproximación por mínimos cuadrados, que permite ajustar un
polinomio de grado menor que el número de puntos disponibles, minimizando el error
global y proporcionando una mejor aproximación en presencia de datos con ruido o
inconsistencias.
Errores en la interpolación. Se analizan los errores asociados
a la interpolación polinómica, incluyendo el fenómeno de Runge, que ilustra cómo
aumentar el número de nodos no siempre mejora la aproximación y puede, en algunos casos,
empeorarla
Nodos de Chebyshev.Se discute la elección óptima de nodos para
la interpolación, destacando que los nodos de Chebyshev minimizan el error de
interpolación y mejoran la estabilidad del polinomio interpolador..
Aplicaciones prácticas. Se proporcionan ejemplos y ejercicios
que ilustran la aplicación de las técnicas de interpolación y aproximación por mínimos
cuadrados en problemas reales, como la estimación de la viscosidad del agua a una
temperatura específica.
Capítulo
Integración numérica
Introduccion
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. Su interés viene motivado por distintas razones:
La posibilidad de aplicar
la Regla de Barrow se restringe a una pequeña parte de funciones elementales que
tienen primitivas también elementales. Muchas funciones, como por ejemplo $e^{x^2}$ que es
continua, no tienen primitiva elemental.
Aunque existiera un procedimiento para encontrar la primitiva,
obtenerla puede resultar en muchos casos muy complicado. Pensemos por ejemplo en
las funciones racionales o en aquellas que requieren de cambios de variable más o menos
ingeniosos.
Por otro lado,
en muchas aplicaciones sólo se tiene una tabla de valores de la función $f$ o la
posibilidad de obtener su valor en un número finito de puntos.
Estas razones obligan a encontrar métodos alternativos para calcular el valor de una integral
definida en un intervalo de forma aproximada. Por razones históricas, el cálculo aproximado de
integrales se denomina cuadratura dejando la denominación de integración para la
resolución de ecuaciones diferenciales.
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
$$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. Se tiene los siguientes datos
Profundidad (m)
Velocidad (m/s)
0
0.0
1
1.2
2
2.1
3
2.4
4
2.0
Profundidad (m)
Velocidad (m/s)
5
1.5
6
0.0
El siguiente gráfico muestra la velocidad del agua en función de la profundidad, con el área
bajo la curva que representa el volumen de agua que fluye por unidad de tiempo.
Velocidad del agua en función de la profundidad
Se verá que la fórmula de los trapecios para el cálculo del caudal permitirá obtener un valor
aproximado:
\(h\): Ancho del intervalo entre los puntos de profundidad (\(h = x_{i+1} - x_i\)).
\(v(x_i)\): Velocidad del agua en el punto \(x_i\) de la profundidad.
\(x_0\) y \(x_n\): Extremos de la profundidad medida.
Como es conocido, las sumas de Riemann permiten aproximar el valor de la integral definida de
una función integrable. Veremos en este tema otras fórmulas utilizando polinomios
interpoladores.
Sumas de Riemann
Fórmulas de cuadratura interpolatorias
Para encontrar una aproximación de la integral definida de $f$ en un intervalo $[a,b]$, se
considerará la integración de un polinomio interpolador de $f$ en $n$ nodos en el intervalo.
Así, si $p(x)$ es un polinomio interpolador de $f$ en $n+1$ nodos,
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} $$
Considerando el polinomio de interpolación de Lagrange, definido de la forma siguiente:
$$p\left( x \right) = \sum\limits_{i = 0}^n {f\left( {{x_i}} \right){l_i}\left( x \right)}
\,\,\,\, , \,\,\,\,\,\,\, {l_i}\left( x \right) = \prod\limits_{j = 0\,\,\,j \ne i}^n {{{x -
{x_j}} \over {{x_i} - {x_j}}}} \,\,\,\,\,\, 0 \le i \le n$$
La aproximación será $$\int\limits_a^b {f\left( x \right)dx} \approx \left( {b - a}
\right)\sum\limits_{i = 0}^n {{w_i}f\left( {{x_i}} \right)} \,\,\,\,\,\,\,\,siendo$$ $${w_i} =
{1 \over {b - a}}\int\limits_a^b {{l_i}\left( x \right)dx} $$ Este tipo de cuadratura se dice
que es una fórmula interpolatoria. Los puntos ${w_i}$ se denominan pesos y los
puntos ${x_i}$ son los nodos de la fórmula de cuadratura.
Consideremos la función $f(x)=e^{-x^2}$, el intervalo $[0,1]$ y los 4 nodos siguientes $x_0=0$,
$x_1=1/3$, $x_2=2/3$ y $x_3=1$. Calculemos la fórmula de cuadratura.
Se verá posteriormente que el cálculo de los pesos $\omega_i$ se puede obtener de otra forma más
sencilla.
La fórmula $\int\limits_0^2 {f\left( x \right)dx} \approx 3f\left( 0 \right) + f\left( 2
\right)$ no es de tipo interpolatorio.
En efecto, se tiene $x_0=0$ y $x_1=2$ con pesos $\omega_0=3/2$ y $\omega_1=1/2$
$$\int\limits_0^2 {f\left( x \right)dx} \approx 2 \left( {\frac{3}{2}} f\left( 0 \right) +
{\frac{1}{2}} f\left( 2 \right)\right) $$ Si fuera de tipo interpolatorio, se tendría que
$w_o=3/2$ $${\frac{1}{2}}\int\limits_0^2 {{l_o}\left( x \right)dx} = {\frac{1}{2}}
\int\limits_0^2 {{{x - {x_1}} \over {{x_o} - {x_1}}}\,} dx = {\frac{1}{2}} \int\limits_0^2 {{{x
- 2} \over {0 - 2}}\,} dx = {\frac{1}{2}} \ne {\frac{3}{2}} $$
Para calcular la calidad de una fórmula de cuadratura se define su grado de precisión o
exactitud.
DEFINICIÓN: Una fórmula de cuadratura tiene grado de exactitud $n\gt0$ si calcula
exactamente la integral para cada polinomio de grado menor o igual que $n$, pero no es exacta
para el menos un polinomio de grado $n+1$. Esto significa que para cada polinomio $q(x)$ de
grado menor o igual a $n$, se cumple $$\int\limits_a^b {q\left( x \right)dx} = \left( {b - a}
\right)\sum\limits_{i = 0}^n {{\omega _i}\,q\left( {{x_i}} \right)} $$
TEOREMA. La fórmula de cuadratura $$\int\limits_a^b {f\left( x \right)dx} \approx \left( {b - a}
\right)\sum\limits_{i = 0}^n {{w_i}f\left( {{x_i}} \right)} $$ es de tipo interpolatorio si y
solo sí es exacta de grado $n$.
Comprobar que la fórmula $$\int\limits_0^4 {f\left( x \right)dx} \approx 2f\left( 1 \right) +
2f\left( 3 \right)$$ es de tipo interpolatorio.
Se puede comprobar que es exacta para $f(x)=1$, $$\int\limits_0^4 {1\,dx} = 2 + 2$$ Es exacta
para $f(x)=x$ $$\int\limits_0^4 {x\,dx} = 2 + 2 \cdot 3$$ Se puede ver además que el grado de
precisión es exactamente uno ya que la fórmula no es exacta para $f(x)=x^2$ ya que
$$\int\limits_0^4 {{x^2}\,dx} = {{64} \over 3} \ne 2 \cdot {1^2} + 2 \cdot {3^2} = 20$$
Esta tipo de fórmulas de cuadratura son eficientes en intervalos de longitud pequeña. Cuando el
intervalo es de longitud grande, se descompone en subintervalos pequeños y se aplica en cada uno
de ellos fórmulas de cuadratura.
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} } $$
y, posteriormente, se aplica la fórmula de cuadratura a cada subintervalo. Este tipo de fórmulas
se denominan de cuadratura compuesta y se analizan en el apartado 5.4 de este tema.
Fórmulas de cuadratura de Newton-Cotes
Veremos en primer lugar cómo obtener fórmulas de cuadratura elementales y el error asociado
limitándonos a las fórmulas de Newton-Cotes. Posteriormente veremos cómo construir fórmulas de
cuadratura compuesta para aproximar la integral con un error inferior a la precisión dada.
Las fórmulas de cuadratura de Newton-Cotes son de tipo interpolatorio $$\int\limits_a^b
{f\left( x \right)dx} \approx \left( {b - a} \right)\sum\limits_{i = 0}^n {{w_i}f\left( {{x_i}}
\right)} \,\,\,\,\,\,\,\,\,\,\,\,\,{w_i} = {1 \over {b - a}}\int\limits_a^b {{l_i}\left( x
\right)dx} $$ donde los nodos se toman equidistantes.
DEFINICIÓN. Se llama fórmula de Newton-Cotes (también fórmula cerrada de Newton-Cotes)
a aquella fórmula de cuadratura numérica de tipo interpolatorio, donde los nodos vienen dados
por la expresión $${x_j} = a + jh\,\,\,,\,\,\,0 \le j \le n\,\,\,\,\,,\,\,\,h = {{b - a} \over
n}$$ Es decir, se tienen $n+1$ nodos equidistantes $$a = {x_o} < {x_1} < {x_2} < ...
< {x_n} = b$$ $$h = {x_j} - {x_{j - 1}}\,\,\,\,\,\,j = 1,...,n$$
En el siguiente interactivo se puede obtener un número variable de nodos equidistantes en un
intervalo.
Intervalo [a, b] con n+1 nodos equidistantes
En la práctica resulta de utilidad el resultado siguiente a la hora de calcular los pesos.
Los pesos son independientes del intervalo siempre que se elijan adecuadamente los nodos, esto
es, estén relacionados por la transformación de $[a,b]$ en $[c,d]$ dada de la siguiente forma,
$$\begin{aligned} & \left[ {a,b} \right] &\to &\left[ c,d \right] \\ & x &\to & t = c + {{d - c}
\over {b - a}}\left( {x - a} \right) \end{aligned}$$
Mapeo entre dos intervalos [a, b] y [c, d]
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 en ambos casos coinciden.
En efecto, $$\widehat {{\omega _i}} = {1 \over {d - c}}\int\limits_c^d {\widehat {{l_i}}\left( t
\right)dt} = {1 \over {b - a}}\int\limits_a^b {{l_i}\left( x \right)dx} = {\omega _i}$$ donde
$$\widehat {{l_i}} \left( t \right) = \prod\limits_{j = 0\,\,\,j \ne i}^n {{{t - {t_j}} \over
{{t_i} - {t_j}}}} \,\,\,\,\,\, 0 \le i \le n$$ $${l_i}\left( x \right) = \prod\limits_{j =
0\,\,\,j \ne i}^n {{{x - {x_j}} \over {{x_i} - {x_j}}}} \,\,\,\,\,\, 0 \le i \le n$$
Con los resultados expuestos anteriormente vamos a ver que la construcción de las fórmulas de
Newton-Cotes no precisa calcular los polinomios base de Lagrange y luego realizar las integrales
para obtener los correspondientes pesos. Bastará resolver un sistema de $n+1$ ecuaciones con
$n+1$ incógnitas que se obtendrá aplicando que la fórmula de cuadratura es exacta para todos los
polinomios de grado menor o igual a $n$ como se muestra a continuación.
Dados los nodos $x_i$ para $i=0,1,...,n$ se pueden determinar los pesos formando un sistema de
$n+1$ ecuaciones con $n+1$ incógnitas considerando las siguientes $n+1$ igualdades:
$$\int\limits_a^b {{x^k}dx} = \left( {b - a} \right)\sum\limits_{j = 0}^n {{\omega _j}x_j^k}
\,\,\,\,\,\,\,k = 0,1,...,n$$ Este sistema tiene una única solución ya que la matriz de
coeficientes es de Vandermonde.
Cálculo de las fórmulas de cuadratura Newton-Cotes
Calculamos las primeras fórmulas de Newton-Cotes aplicando también que el cálculo de los pesos
no depende del intervalo siempre que se elijan adecuadamente en él.
Primera fórmula de Newton-Cotes (n=1)
Dados dos nodos en el intervalo $[a,b]$, esto es, $x_0=a$ y $x_1=b$, 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( b \right)} \right]$$
Para calcular los pesos se considera el intervalo $[0,1]$ y se resuelve el sistema siguiente
$$\begin{aligned} & k = 0\,\,\,\,\,\,\,\,\,& 1 = \int\limits_0^1 {dx} = {\omega _o} + {\omega
_1} \\ & k = 1\,\,\,\,\,\,\,\,\,&{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,
$$\int\limits_a^b {f\left( x \right)dx} \approx {{b - a} \over 2}\left[ {f\left( a \right) +
f\left( b \right)} \right]$$ Esta fórmula de Newton-Cotes de dos puntos se llama
fórmula del trapecio.
En este caso el número de nodos 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:
$$\int\limits_a^b {f\left( x \right)dx} \approx $$ $$ \approx {{\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] $$ Esta fórmula se llama
regla de los tres octavos.
Regla de los tres octvavos
Notar que los pesos son simétricos, es decir, el primero es igual al último, el segundo es igual
al penúltimo, etc.
De forma análoga pueden obtenerse las siguientes fórmulas de Newton-Cotes.
Obtener con Matlab la aproximación de la integral $\int\limits_0^\pi {{e^{senx\,\,\cos x}}dx} $
utilizando la regla de los tres octavos ($n=3$) y la de Hardy ($n=6$).
****Añadir solución *****
Error en la fórmula de cuadratura de Newton-Cotes.
Se analizará en este apartado el error de la aproximación, es decir, la diferencia entre el
valor exacto y el aproximado por 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.
TEOREMA. Si $n$ es impar y la función $f$ es $n+1$ veces derivable y $f^{(n+1}$ es continua
entonces existe $\xi \in \left[ {a,b} \right]$, de forma que el error de la fórmula viene dada
por $$E = {{{h^{n + 2}}} \over {\left( {n + 1} \right)!}}{f^{(n + 1}}\left( \xi
\right)\int\limits_0^n { \pi \left( x \right)} \,dx$$ donde $h=\frac{b-a}{n}$ y $\pi \left( x
\right) = \left( x \right)\left( {x - 1} \right)...\left( {x - n} \right)$ Si $n$ es par y la
función $f$ es $n+2$ veces derivable y $f^{(n+2}$ es continua entonces existe $\xi \in \left[
{a,b} \right]$, de forma que el error de la fórmula viene dada por $$E = {{{h^{n + 3}}} \over
{\left( {n + 2} \right)!}}{f^{(n + 2}}\left( \xi \right)\int\limits_0^n {x\pi \left( x \right)}
\,dx$$
Según este teorema, si $n$ es par, entonces el grado de precisión de la fórmula de Newton-Cotes
es $n+1$. Cuando $n$ es impar, el grado de precisión es $n$.
Este teorema se puede escribir de la siguiente manera:
El error de la fórmula de cuadratura de Newton Cotes se puede escribir de la forma siguiente:
$$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.$, $C_n$ una
constante que viene dada por la expresión anterior y $\xi \in \left[ {a,b} \right]$.
Vamos a determinar el valor de las constantes $C_n$ para las primeras fórmulas de Newton-Cotes.
Error en la fórmula del trapecio ($n=1$)
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
Se cumple entonces $$\int\limits_a^b {f\left( x \right)dx = {{b - a} \over 2}\left[ {f\left( a
\right) + f\left( b \right)} \right]} - {1 \over {12}}{h^3}{f^{''}}\left( \xi \right)$$ para
algún valor de $\xi$ en $[a,b]$.
Error en la fórmula de Simpson
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$$
Se cumple entonces $$\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$$ $$E = - {1
\over {90}}{h^5}{f^{(iv}}\left( \xi \right)$$ para algún valor de $\xi$ en $[a,b]$.
Error en la fórmula de cuadratura 3/8
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] $$
$$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)}}$$
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.
Fórmulas compuestas
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 cuadratura del trapecio compuesta consiste en aproximar una integral en $[a,b]$
utilizando $N$ trapecios. En primer lugar, se divide el intervalo en $N$ subintervalos de
longitud $h=\frac{b-a}{N}$. Despues de aplicar en cada uno de ellos la fórmula de cuadratura
de los trapecios se obtiene
La cuadratura de Simpson compuesta consiste en aproximar una integral en $[a,b]$
dividiendo el intervalo en $N$ subintervalos de longitud $h=\frac{b-a}{N}$ y aplicar en cada
uno de ellos la regla de Simpson. En este caso se precisarían $2N+1$ nodos $x_i=a+ih/2,
\,\,0 \le i \le 2N$ se obtiene
Calcular la siguiente integral $$\int\limits_0^1 {\left( {1 + sen\left( {{x^2}} \right)}
\right)\,dx} $$ utilizando la cuadratura de los trapecios compuesta dividiendo el intervalo en 4
subintervalos.
*** Añadir explicación ****
Métodos adaptativos
Observa que los métodos compuestos no tienen en cuenta las zonas de mayor variabilidad de la
función ya que se considera la longitud de cada subintervalo siempre igual.
Dependiendo de la función, puede ser conveniente utilizar subintervalos más pequeños donde haya
más variabilidad de la función. En estos métodos adaptativos se estima el error cometido en cada
intervalo de forma que permita decidir si hay que subdividir el intervalo.
Cuadratura adaptiva trapecios
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]$.
Por otro lado, 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.
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
Una estimación del error considerando como valor de la integral el valor de $I_1$ es decir,
utilizando la cuadratura de los trapecios en el intervalo $[a,b]$, es $\color{blue}{\left
|\frac{4}{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,
Si la estimación del error es menor que la tolerancia, devolvemos el valor de $I_2$ con la
estimación del error.
Si esta estimación del error es mayor que la tolerancia, repetimos el proceso en cada
subintervalo $[a, \frac{a+b}{2}$ y $[\frac{a+b}{2},b]$ permitiendo en cada subintervalo solo
la mitad de la tolerancia.
Calcular $$\int\limits_0^1 {\left( {1 + sen\left( {{x^2}} \right)} \right)\,dx} $$ utilizando la
cuadratura de los trapecios adaptativa con un error menor que $10^{-5}$.
Cuadratura adaptiva Simpson
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á, llamando $h=(b-a)/2$
$$\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}$$ 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|}$$
Una estimación del error utilizando Simpson en un intervalo $[a,b]$ es, por tanto,
$\color{blue}{\left |\frac{16}{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,
Si esta estimación del error es menor que la tolerancia, devolvemos el valor de $I_2$ con la
estimación del error.
Si esta estimación del error es mayor que la tolerancia, repetimos el proceso en cada
subintervalo $[a, \frac{a+b}{2}$ y $[(\frac{a+b}{2},b]$ permitiendo en cada subintervalo
solo la mitad de la tolerancia.
Calcular $$\int\limits_0^1 {\left( {1 + sen\left( {{x^2}} \right)} \right)\,dx} $$ utilizando el
método de Simpson adaptativo con un error menor que $10^{-5}$.
**** Añadir solución ****
Calcular la siguiente integral $\int\limits_0^\pi {{e^{{\rm{sen}}x\cos x}}\,dx} \,$ con un error
inferior a $10^{-12}$.
Solución: 3.341031544735790
Calcular la siguiente integral $\int\limits_0^4 {{e^{ - {x^2}}}dx} $ con un error inferior a
$10^{-12}$.
Solución: 0.886226911789732
Calcular la siguiente integral $\int\limits_0^4 {{e^{{x^2}}}dx} $ con un error inferior a
$10^{-6}$.
Solución: 1.149400634590278e+06
Cuadratura adaptiva
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.$$
Se ha obtenido que si se utiliza la fórmula de Newton-Cotes de orden $n$ para aproximar la
función en el intervalo $[a,b]$, tomando $m=n+1$ si n es impar y $m=n+2$ si n es par y
considerando las aproximaciones $I$, $I_1$ e $I_2$ respectivamente en los subintervalos $[a,b]$,
$[a,\frac{a+b}{2}]$ y $[\frac{a+b}{2},b]$, se cumple que el error $E$ de la aproximación en
$[a,b]$ es $$\left| E \right| \approx \left| {{{2^m}} \over {{2^{m}-1}}}\left( {{I_1} + {I_2} -
I} \right) \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.
Elegimos una fórmula de Newton-Cotes para aproximar la integral en $[a,b]$.
Subdividimos el intervalo $[a,b]$ en dos subintervalos de la misma longitud y calculamos las
aproximaciones de las integrales $I_1$ e $I_2$ en ambos intervalos, con las correspondientes
estimaciones de los errores $E_1$ y $E_2$. Definimos los vectores $I=[I_1, \,\,I_2]$ y
$E=[E_1, \,\,E_2]$
Si la suma de los errores, sum(E), es inferior a $\epsilon$, tomamos la suma de las
aproximaciones $I$ como valor aproximado de la integral y se habrá terminado. Si no es así,
$sum(E) \ge \epsilon$, entonces vamos al paso siguiente.
Buscamos el error $E_k$ más grande de los contenidos en el vector $E$. Subdividimos el
intervalo asociado a $E_k$ en dos mitades y calculamos las aproximaciones de la integral en
cada una de las mitades y los correspondientes errores. Eliminamos $I_k$ y $E_k$ y ponemos
en su lugar los dos valores de la integral calculados y los correspondientes errores.
Volvemos al paso 3.
Para las prácticas utlizaremos los ficheros iehardy.m que a su vez
utiliza la función hardy.m, para obtener la aproximación de la
integral utilizando la la fórmula de Hardy. Estas dos funciones han sido creadas por el profesor
Eduardo Casas (ficheros Matlab).
[itg error]=iehardy(a,b,f)
%Devuelve una estimación de la integral y
%el error cometido. Requiere que la
%función f se programe vectorialmente.
Calcular la integral $\int\limits_0^\pi {{e^{senx\,\cos x}}dx} $ con un error inferior a
$\epsilon_a=10^{-12}$ considerando la regla de Hardy.
a=0;b=pi;
format long
f=@(x) exp(sin(x).*cos(x));
error=10^-12;
Dividimos el intervalo en dos partes y aplicamos la fórmula de Newton-Cotes de Hardy. Se tendrá
que
Si el error es menor que $\epsilon_a$ habremos terminado y el valor de la integral es itg. Como
en este caso no es así, se considera el subintervalo con mayor error. Ejecutando la orden,
[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.
Los vectores $x$, $E$, $I$ van creciendo según se van tomando subdivisiones en los intervalos.
Se repite el proceso 26 veces obteniendo que el valor de la integral es $ 3.341031544735853$. El
valor que devuelve matlab es $3.341031544735852$
integral(f,0,pi)
Si se considera el error absoluto en la precisión de la solución el algoritmo se detiene si $E
< \epsilon_a$. Si se usa el error relativo, el algoritmo se detenderá si $E <
\epsilon_r|I|$. A menudo se utiliza un test de parada que combina ambos errores: $$E <
\epsilon_a+\epsilon_r |I|$$
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
Cálculo de Integrales impropias
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).
Integrales impropias de primera especie
TIPO I. Consideraremos en este caso integrales del tipo $$\int\limits_a^\infty {f\left( x
\right)dx} = \mathop {\lim }\limits_{b \to \infty } \int\limits_a^b {f\left( x \right)dx} $$ donde
$f$ es continua. Cuando el límite anterior existe, la integral impropia es convergente. En este caso
se tendrá que $$\mathop {\lim }\limits_{b \to \infty } \int\limits_b^\infty {f\left( x \right)dx} =
0$$
En primer lugar se hallará un valor de $b$ de forma que $$\left| {\int\limits_b^\infty {f\left(
x \right)dx} } \right| < {\varepsilon \over 2}$$
y luego se buscará una aproximación de la intengral $Itg$ de forma que $$\left| {\int\limits_a^b
{f\left( x \right)dx} - Itg} \right| < {\varepsilon \over 2}$$
TIPO II. De forma análoga se podrá hacer para calcular las integrales
$${\int\limits_{ - \infty }^a {f\left( x \right)dx} }$$
TIPO III. En el caso de una integral de la forma $${\int\limits_{ - \infty }^\infty
{f\left( x \right)dx} }$$ se tendrá en cuenta que $$\int\limits_{ - \infty }^\infty {f\left( x
\right)dx} = \int\limits_{ - \infty }^0 {f\left( x \right)dx} + \int\limits_0^\infty {f\left( x
\right)dx} $$ buscando una aproximación de cada una de ellas con un error menor que
$\epsilon/2$. Alternativamente se puede buscar números $a<0$ y $b>0$ de forma que $$\left|
{\int\limits_{ - \infty }^a {f\left( x \right)dx} } \right| < {\varepsilon \over
4}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left| {\int\limits_b^\infty {f\left( x \right)dx} } \right|
< {\varepsilon \over 4}\,\,\,\,\,\,\,\,\,\left| {\int\limits_a^b {f\left( x \right)dx} }
\right| < {\varepsilon \over 2}$$
Dado que $$\int\limits_0^\infty {{e^{ - x}}{{\cos }^2}\left( x \right)dx} = \underbrace
{\int\limits_0^b {{e^{ - x}}{{\cos }^2}\left( x \right)dx} }_{ = {I_1}} + \underbrace
{\int\limits_b^\infty {{e^{ - x}}{{\cos }^2}\left( x \right)dx} }_{ = {I_2}}$$ se debe tomar $b$
de forma que $$\left| {{I_2}} \right| \le {{{{10}^{ - 12}}} \over 2}$$ Puede considerarse
$${\log \left( {2 \cdot {{10}^{12}}} \right) \le b}$$ Calbulando $I_1$ con un error menor que
$10^{-12}/2$, el valor de la integral es $0.599999999999007$.
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 $$
Caso 1 Si el punto $x_0=a$ será $$\int\limits_a^b {f\left( x \right)dx} = \mathop {\lim
}\limits_{\delta \searrow 0} \int\limits_{a + \delta }^a {f\left( x \right)dx} $$
Se procede de la siguiente forma:
Se busca $\delta$ cumpliendo $$\left| {\int\limits_a^{a + \delta } {f\left( x \right)dx} }
\right| < {\varepsilon \over 2}$$ y luego una aproximación verificando $$\left|
{\int\limits_{a + \delta }^b {f\left( x \right)dx} } \right| < {\varepsilon \over 2}$$
Caso 2 Si el punto $x_0=b$ será $$\int\limits_a^b {f\left( x \right)dx} = \mathop {\lim
}\limits_{\delta \searrow 0} \int\limits_a^{b - \delta } {f\left( x \right)dx} $$
Se procede de la misma forma que en el caso anterior
Caso 3 Si $x_0$ es un punto interior se define esta integral impropia de la forma
siguiente: $$\int\limits_a^b {f\left( x \right)dx} = \mathop {\lim }\limits_{\delta \searrow 0}
\int\limits_a^{{x_o} - \delta } {f\left( x \right)dx} + \mathop {\lim }\limits_{\delta \searrow 0}
\int\limits_a^{{x_o} + \delta } {f\left( x \right)dx} $$
Se busca un valor de $\delta > 0$ verificando $$\left| {\int\limits_{{x_o} - \delta }^{{x_o}}
{f\left( x \right)dx} } \right| + \left| {\int\limits_{{x_o}}^{{x_o} + \delta } {f\left( x
\right)dx} } \right| < {\varepsilon \over 2}$$ Después se buscan aproximaciones de las
integrales cumpliendo $$\left| {\int\limits_a^{{x_o} - \delta } {f\left( x \right)dx} } \right|
< {\varepsilon \over 4}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left| {\int\limits_{{x_o}
+ \delta }^b {f\left( x \right)dx} } \right| < {\varepsilon \over 4}$$
Calcular $$\int\limits_0^1 {{{\cos \left( x \right)} \over {2\pi \,{\rm{sen}}\left( {\sqrt x }
\right)}}dx} $$ considerando $e_a=10^{-12}$.
Solución Matlab: 0.302993744656398
Autoevaluación
**** Incluir preguntas del tema ****
Actividad de autoevaluación
Ejercicios
Consideremos la función \(f(x) = e^{-2x}\), el intervalo \([0,1]\), y los nodos: \(x_0 =
0\), \(x_1 = 1/3\), \(x_2 = 2/3\), \(x_3 = 1\). Calcular la fórmula de cuadratura a partir
de los polinomios base de Lagrange.
Utiliza la fórmula de los trapecios compuesta con tres subintervalos para aproximar la
integral \(\int_0^{\pi} \sin(x) dx\). Calcular cuántos subintervalos serían necesarios para
que el error cumpla \(E \leq 0.5 \times 10^{-8}\).
Aproximar el valor de \(\pi\) mediante la integral: \[ \int_0^1 \frac{4}{1+x^2} dx, \]
utilizando las fórmulas compuestas del trapecio y de Simpson con 1024 nodos.
Calcular el área limitada por las gráficas de las funciones: \[ y = \frac{\sin(x)}{1+x},
\quad y = 0, \quad x = 1, \] con un error inferior a \(10^{-8}\). Calcular los puntos de
corte utilizando el método de Newton.
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.
Fórmulas de cuadratura interpolatorias
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.
Fórmulas de Newton-Cotes
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.
Fórmulas compuestas
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.
Cálculo de integrales impropias
Se presentan métodos para evaluar integrales que tienen singularidades o límites
infinitos, diferenciando entre aquellas de primera y segunda especie.
Errores en integración numérica
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.
Ejercicios y autoevaluaciones
El capítulo incluye ejemplos prácticos y ejercicios diseñados para aplicar y
comprender los conceptos, así como autoevaluaciones para reforzar el aprendizaje.
Capítulo
Integración numérica de ecuaciones diferenciales
Introduccion
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:
el problema de Cauchy o problema de valor inicial
el problema de contorno.
En un problema de valor inicial las condiciones se establecen sobre el instante inicial y en un
problema de contorno, las condiciones se establecen tanto en el instante inicial como en el
final.
En la siguiente escena se muestra un ejemplo de aplicación de las ecuaciones diferenciales a un
sistema muelle-resorte.
Un problema de valor inicial (PVI) se formula a partir de una ecuación diferencial ordinaria junto
con un valor específico de la función en un instante dado, de modo que la solución buscada deba
satisfacer simultáneamente la relación diferencial y dicha condición inicial.
El problema de valor inicial, PVI, puede escribirse de forma general como $$\left\{
\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} \right.$$
siendo $$ - \infty < {t_0} < {t_f} < \infty $$ $$f:\left[ {{t_o},{t_f}}
\right]\text{x}{\mathbb R^n} \to \mathbb R^n$$
Esta forma permite escribir muchos problemas. Veamos algún ejemplo.
Consideremos el siguiente problema
$$\left\{ \begin{aligned} & y''\left( t \right) + y\left( t \right) = 1\,\,\,\,\,\,\,\,\,\,t \in
\left[ {0,2\pi } \right] \\ & y\left( 0 \right) = 1,\,\,y'\left( 0 \right) = 1 \end{aligned}
\right.$$ Escribirlo en la forma dada para un problema de valor inicial.
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
Para asegurar que un problema de valor inicial (PVI) tenga una única solución, se considera el
siguiente teorema de existencia y unicidad.
Dado el problema $$\left\{ \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} \right.$$ con $$f:\left[ {{t_o},{t_f}} \right]\text{x}{\mathbb R^n} \to \mathbb
R^n$$ si $f$ es continua y acotada en un dominio que contenga a $(t_0,y_o)$ entonces el problema
tiene al menos una solución.
Para que la solución sea única se puede exigir que la función cumpla la
condición de Lipschitz para la variable $y$, es decir, exista $L> 0$ de forma que
para todo $t$ en $[t_0,t_f]$
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
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.
Representación de la solución exacta y aproximada de un PVI
Método de Euler
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
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$.
Primer paso del método de Euler. Método de Lagrange
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.
Consideremos el problema de valor inicial
$${{dy} \over {dt}} = 0.7y - {t^2} +
1\,\,\,\,\,\,\,,\,\,\,\,\,\,y\left( 1 \right) = 1\,\,\,\,\,\,\,,\,\,\,\,\,\,\,\,t \in \left[
{1,2} \right]$$ Encontrar la solución aproximada al problema utilizando por el método de Euler
con 11 puntos.
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
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).
Cálculo de la solución de la EDO por el método de Euler
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'])
Gráfica de la solución exacta y aproximada
.
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.
**** poner ejemplo de una edo de orden mayor ****
Método de Heun
Como en el caso anterior, para encontrar el valor de $y_{n+1}$ a partir de $y_n$, tendremos en
cuenta que
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
El valor de $y_{n+1}$ se obtiene a partir del punto $(t_n,y_n)$ considerando la media ponderada
de la pendiente en ese punto, y la pendiente en el punto que se obtendría aplicando el método de
Euler, $(t_n+h_n,y_n+h_nf(t_n,y_n)))$.
Método de Heun
Este método podría escribirse de la forma siguiente:
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}$$
Calculo del primer paso por el método de Heun (cambiar la imagen***)
Considera el siguiente PVI:
$$y' = (t-y)/2 \,\,\,\,\,, t \in [0, 3] \,\,\,\, y(0) = 1$$ Calcula por el método de Heun con 10
pasos la solución aproximada. La solución exacta para esta ecuación es $y(t) = 3e^{t/2} - 2 +
t$.
En el siguiente interactivo puedes ver los cálculos para obtener la solución aproximada. En este
caso $f(t,y)=\frac {t-y}{2}$, $t_0=0$, $t_f=3$, $y_0=1$ y $h=3/10$.
.
Solución aproximada por el método de Heun
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,'*')
El código para resolver la ode con matlab considerando $N=30$ pasos se muestra a
continuación.
Método general de Runge Kutta
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,
$$y\left( {{t_{n + 1}}} \right) = y\left( {{t_n}} \right)\, +\color{blue}{
\int\limits_{{t_n}}^{{t_{n + 1}}} {f\left( {t,y\left( t \right)} \right)dt}} $$ La aproximación,
se obtendrá calculando la integral de $f(t,y(t))$ utilizando una fórmula de cuadratura de $m$
nodos $$y\left( {{t_{n + 1}}} \right) \approx {y_n} + \color{blue}{{h_n}\sum\limits_{i = 1}^m
{{b_i}f\left( {{t_n} + {c_i}{h_n},\color{red}{y\left( {{t_n} + {c_i}{h_n}} \right)}} \right)}}
$$ donde $b_i$ son los pesos y los $c_i$ son los números entre 0 y 1 que definen los nodos
$$t_{n,i}=t_n+c_i h_n$$ en el intervalo $[t_n,t_{n+1}]$.
En la imagen se representa un ejemplo con $m=6$ donde se muestran los valores $c_i$ y $t_{h,i}$.
Representación de los nodos de la forma de cuadratura con $m=6$
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
Observad que el método de Runge-Kutta queda definido a partir de los valores $c_i$, los pesos
$b_i$ y la matriz $(a_{ij})$.
Se define un método Runge-Kutta de $m$ etapas como el método numérico que dada una
aproximación $y_n$ nos da una aproximación a dicha solución en el punto $t_{n+1}=t_n+h_n$
mediante las siguientes fórmulas: $$ \begin{aligned} & {t_{n,i}} = {t_n} + {\color{red} {c_i}}
{h_n} , 1 \le i \le m\\ & \left\{ \begin{aligned} &\,\,\,{y_{n,1}}
= {y_n} , {K_{n,1}} = f\left( {{t_{n,1}},{y_{n,1}}} \right)\\
&\,\,\,{y_{n,i}} = {y_n} + {h_n}\sum\limits_{j = 1}^{m} {\color{green}{a_{ij}}}
{K_{n,j}} , {K_{n,i}} = f\left( {{t_{n,i}},{y_{n,i}}} \right)\\&\,\,\,\,\,\,\,\,\,\,2
\le i \le m \end{aligned} \right.\\ &{y_{n + 1}} = {y_n} + {h_n}\sum\limits_{j = 1}^m
{\color{blue}{b_j}}{K_{n,j}}\\ \end{aligned}$$
El método puede escribirse también la forma siguiente:
Los métodos de Runge-Kutta pueden describirse a trávés del llamado
tablero de ButcherJohn
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.
Tableros de Butcher
Método de Runge de tres pasos
En este caso, se considerará la fórmula de cuadratura de tres puntos
Método de Heun de Tres Etapas
$$\begin{array} {r|r r r } 0 & 0 & 0 & 0 \\ 1/3 & 1/3 & 0 & 0 \\ 2/3 & 0 & 2/3 & 0 \\ \hline & 1/4 &
0 & 3/4 \\ \end{array}$$
Método de Runge-Kutta Clasico
$$\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}$$
Método de Runge-Kutta asociado a tableros
Escribir el algoritmo del método Runge-Kutta explícito de tres etapas cuyo tablero tiene la
forma $$\begin{array} {r|r r } 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 \\ 3/4 & 0 & 3/4 & 0 \\ \hline
& 2/9 & 3/9 & 4/9 \\ \end{array}$$ Dado el PVI, $$y'(t)=t+e^{y(t)}, y(0)=0$$ aproximar mediante
dicho método, la solución en el tiempo $t=0.4$ tomando como paso $h=0.2$:
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;
En el siguiente interactivo se calcula los dos pasos con otros métodos de Runge-Kutta.
Método de Runge-Kutta de 3 pasos
Resolver el problema $y'=yt^3-1.5y$ con la condición inicial $y(0)=1$ en el intervalo $[0,1.8]$
considerando 9 pasos en el método Runge-Kutta clásico (de cuarto orden): $$\begin{array} {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 escena se visualiza, para distintos métodos de Runge-Kutta, incluido el clásico
de cuarto orden, la integral del problema. En este caso puede obtenerse la solución exacta que
es la función $y(t)={e^{{{t\left( {{t^3} - 6} \right)} \over 4}}}$.
Método de Runge-Kutta de 4 pasos
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.
Consideremos el problema plano de dos cuerpos puntuales que se mueven atraídos según la ley de
la gravitación universal de Newton. En unidades apropiadas, las ecuaciones del movimiento
adoptan la forma siguiente: $$\begin{aligned} & y_1' = {y_3} \\ & y_2' = {y_4} \\ & y_3' = -
{{{y_1}} \over {{{\left( {{y_1}^2 + y_2^2} \right)}^{3/2}}}}\\ & y_4' = - {{{y_2}} \over
{{{\left( {{y_1}^2 + y_2^2} \right)}^{3/2}}}} \end{aligned}$$ donde $(y_1,y_2)$ son las
coordenadas de uno de los cuerpos (satélite) en un sistema con origen en el otro cuerpo central.
Los valores $(y_3,y_4)$ representan el vector velocidad del cuerpo satélite.
Utiliza el método de Runge-Kutta clásico para encontrar la solución con las condiciones
iniciales siguientes
Ségún la primera ley de Kepler: ”los planetas describen ́orbitas eĺıpticas alrededor del sol, que
ocupa uno de los focos de esa elipse”. En general, las soluciones del problema de dos cuerpos
son cónicas (pueden ser elipses, paŕabolas o hiṕerbolas), pero si se toman las condiciones
iniciales siguientes: $$y_1(0) = 1,\,\,y_2(0) = y_3(0) = 0,\,\,\,y_4(0) = 1$$ la solución es una
circunferencia de centro el origen radio 1 y con periodo $2\pi$.
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
En la siguiente imagen se muestra la trayectoria del satélite una vez resuelto el problema.
Para resolver este ejercicio de forma alternativa, vamos a utilizar dos funciones creadas por el
profesor Eduardo Casas (ficheros Matlab).
la función mirkf4c.m para definir los
parámetros del método de Runge-kutta clásico: a, b, c
la función pasorkf.m para definir un paso
del método con los parámetros: f, tn, hn, yn, a, b, c
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
Error de los métodos de Runge-Kutta
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.
El error local del método en el paso $n$ es la cantidad $${\varepsilon_n} = y\left(
{{t_n} + {h_n}} \right) - \left[ {y\left( {{t_n}} \right) + {h_n}{\Phi _f}\left( {{t_n},y\left(
{{t_n}} \right),{h_n}} \right)} \right]$$ El error global del método en el paso $n$ es
la cantidad $${e _n} = y\left( {{t_n} + {h_n}} \right) - \left[ {{y_n} + {h_n}{\Phi _f}\left(
{{t_n},{y_n},{h_n}} \right)} \right]$$ Diremos que el método es de orden $p$ si para
toda función $f$ definiendo la ecuación diferencial se tiene que $||{\varepsilon _n}||=O(h_n^{p
+ 1})$, es decir, si existe $C>0$ de forma que $||{\varepsilon_n}|| \le C h_n^{p+1}$.
Observa que:
El error local es que se comete al dar un paso partiendo del valor exacto de la solución
en el instante $t_n$.
El error global es el que se comete al dar un paso a partir de la aproximación de ${y_n}
\approx y(t_n)$.
Puede demostrarse que si el método es de orden $p$, entonces el error global es estimado por la
siguiente expresión:
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.
Orden 1. $b^Te=1$
Orden 2. La anterior más la condición $b^Tc=1/2$
Orden 3. La anterior más la condición $b^Tc^2=1/3$, $b^TAc=1/6$
Puede comprobarse que el método de Euler es de orden 1, los de Heun de 2 y 3 etapas poseen,
respectivamente, orden 2 y 3. El método de Kutta de 3 etapas es de orden 3 y el de Runge-Kutta
clásico de orden 4.
Métodos de Runge Kutta encajados
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 métodos encajados.
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.
Se considera entonces dos métodos: $${y_{n + 1}} = {y_n} + {h_n}{\Phi _f}\left(
{{t_n},{y_n},{h_n}} \right)\,\,\,\,\,\,\,\,\,\,orden\,\,p$$ $${\widehat y_{n + 1}} = {y_n} +
{h_n}{\psi _f}\left( {{t_n},{y_n},{h_n}} \right)\,\,\,\,\,\,\,\,\,\,orden\,\,q > p$$
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
los estimadores (${{{\widehat b}_j} - {b_j}})$ para calcular $$e_n=||{\widehat y_{n +
1}}-y_{n+1}||$$
La información que se precisa se incluye en el siguiente tablero.
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 función rkf45b.m desarrollada por Educardo Casas
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 #2.
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:
f: nombre de la función que define la ecuación diferencial
tn: el instante del que se parte
hn: el paso que se va a dar
yn: la aproximación de la solución conocida en el instante $t_n$
a,b,c,est: los parámetros del algorítmo A,b,c y los estimadores
Esta función devuelve
yn1: la aproximación de $y(t_{n+1})$ con $t_{n+1}=t_n+h_n$
en: una estimación del error $e_n=||y(t_{n+1})-y_{n+1}||$
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:
Reducción del paso tras un fallo. Si $e_n>TOL$ se debe reducir $h_n$ con un factor
de reduccion $r_n<1$. Se considera $${r_n} = 0.9{\left( {{TOL} \over {e_{n}} }
\right)^{{1 \over {p + 1}}}}$$ y como medida de seguridad se reducirá el paso al menos a la
mitad, por lo que se considerará $${h_n} \to \min \left\{ {{r_n},{r_{\min }}}
\right\}{h_n}$$ Se toma un valor entre 0.1 y 0.5, tomaremos $r_{min}=0.5$
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:
$e_a=10^{-14}$, $e_r=10^-12$.
$f$ es el nombre del fichero .m que define a la función diferencial
$t$ es el vector fila de instantes $t_n$ donde se ha calculado la aproximación de la
solución de (P).
$y$ es el vector solución que tiene dos filas $y_{x,1}$ para $y_1$ e $y_{x,2}$ para
$y_2$.
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.
Calcular la solución de $$t^3y''=(y-t*y')^2$$ $$y(1)=0,\,\,\,y(2)=2 \,\,\,\,\,\,\,\,t \in \left[
{1,2 } \right]$$
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.
Calcular la solución de $$t^3y''=(y-t*y')^2$$ $$y(1)=0,\,\,\,y'(2)=2 \,\,\,\,\,\,\,\,t \in
\left[ {1,2} \right]$$
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$.
Resolver el siguiente problema:
$$\begin{aligned} & y''\left( t \right) = sen\left(
{y\left( t \right)} \right) + 0.01\,t\,y\left( t \right) + 0.1\,z\left( t \right)\\ & z'\left( t
\right) = sen\left( {y\left( t \right) + z\left( t \right)} \right) + w\left( t \right)\\ &
w'\left( t \right) = \cos \left( {y\left( t \right) + z\left( t \right) + w\left( t \right)}
\right)\\ & y\left( 0 \right) = 0 = z\left( 0 \right),\,\,\,\,\,w\left( 0 \right) =
1\,\,\,\,\,\,\,\,t \in \left[ {0,\pi } \right]\\ & y\left( \pi \right) + z\left( \pi \right) +
w\left( \pi \right) = 0 \end{aligned}$$
Consideramos $${y_1} = y,\,\,\,{y_2} = y',\,\,\,\,\,{y_3} = z,\,\,\,\,{y_4} = w$$ pudiendo
escribirse $$\begin{aligned} & y'_1 = {y_2}\\ & y'_2 = sen\left( {{y_1}} \right) + 0.01\,\,
t\,\, {y_1} + 0.1\,\, {y_3}\\ & y'_3\left( t \right) = sen\left( {{y_1} + {y_3}} \right) +
{y_4}\\ & y'_4\left( t \right) = \cos \left( {{y_1} + {y_3} + {y_4}} \right)\\ & {y_1}\left( 0
\right) = 0 = {y_3}\left( 0 \right),\,\,\,\,\,{y_4}\left( 0 \right) =
1\,\,,\,\,\,\,\,\,{y_2}\left( 0 \right) = x \end{aligned}$$ y se buscará x para que
$${y_{x,1}}\left( \pi \right) + {y_{x,3}}\left( \pi \right)+ {y_{x,4}}\left( \pi \right) = 0$$
La función que definirá la ecuación diferencial será
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 la función g habrá que considerar $$yo = [0;x;0;1]$$ $$fx = {y_x}\left( {1,end} \right) +
{y_x}\left( {3,end} \right) + {y_x}\left( {4,end} \right)$$
Autoevaluación
**** Incluir preguntas del tema ****
Actividad de autoevaluación
Ejercicios
Resolver el problema de valor inicial: \[ \frac{dy}{dt} = -0.7y + t, \quad y(1) = 1, \quad t
\in [1,2]. \] Utiliza el método de Euler con 10 pasos y representa la solución con MATLAB.
Resuelve el siguiente problema de valor inicial con los métodos de Euler, Heun, y
Runge-Kutta de cuarto orden, utilizando 30 pasos: \[ y'' + 7\sin(y) + 0.1\cos(t) = 0, \quad
y(0) = 0, \, y'(0) = 1, \quad t \in [0,5]. \]
Implementa el método de Runge-Kutta explícito de tres etapas (tablero de Butcher) para
resolver: \[ \frac{dy}{dt} = t + e^t, \quad y(0) = 0, \quad t = 0.4, \quad h = 0.1. \]
Resolver el sistema de ecuaciones diferenciales que describe el movimiento de dos cuerpos
puntuales bajo la ley de gravitación universal: \[ \begin{aligned} y_1' &= y_3, \quad y_2' =
y_4, \\ y_3' &= -\frac{y_1}{(y_1^2 + y_2^2)^{3/2}}, \quad y_4' = -\frac{y_2}{(y_1^2 +
y_2^2)^{3/2}}. \end{aligned} \] Utiliza el método de Runge-Kutta clásico con 20 pasos en el
intervalo \( [0, 8\pi] \).
Resuelve el siguiente sistema que describe la competencia entre dos especies por un mismo
recurso: \[ \begin{aligned} y_1' &= y_1(4 - 0.0003y_1 - 0.0004y_2), \\ y_2' &= y_2(2 -
0.0002y_1 - 0.0001y_2), \end{aligned} \] con \( y_1(0) = y_2(0) = 10000 \) y \( t \in [0,4]
\), utilizando el método de Runge-Kutta de cuarto orden.
Resolver el siguiente sistema de ecuaciones: \[ \begin{aligned} y_1' &= y_2, \\ y_2' &=
-\sin(y_1), \end{aligned} \] utilizando el método de Runge-Kutta con \( h = 0.1 \) y \( t
\in [0, 2\pi] \).
Resuelve el sistema de ecuaciones diferenciales: \[ \begin{aligned} x' &= y + t^2, \\ y' &=
x - t, \end{aligned} \] con \( t \in [0, 5] \), utilizando un método implícito y comparando
con Runge-Kutta clásico.
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.
Formulación del problema de valor inicial
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.
Existencia y unicidad de soluciones
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.
Método de Euler
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.
Método de Heun
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.
Métodos de Runge-Kutta
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.
Error en los métodos de Runge-Kutta
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.
Métodos de Runge-Kutta encajados
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.
Problemas de contorno y método de tiro
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.
Capítulo
Programación lineal
Introduccion
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.
Definición de Optimización
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.
Clasificación de problemas de optimización
Según la Naturaleza del Objetivo
Optimización de máximo: Buscar el valor más alto de la función objetivo.
Optimización de mínimo: Buscar el valor más bajo de la función objetivo.
Según el Tipo de Variables
Continuas: Las variables pueden tomar cualquier valor dentro de un rango.
Discretas: Las variables solo pueden tomar valores específicos.
Enteras: Las variables son números enteros.
Mixtas: Combinación de variables continuas y enteras.
Según las Restricciones
Con restricciones: Incluyen limitaciones explícitas en la solución.
Sin restricciones: No tienen limitaciones explícitas en la solución.
Según la Naturaleza del Problema
Lineal: La función objetivo y las restricciones son lineales.
No lineal: La función objetivo o las restricciones no son lineales.
Convexo: La región factible y la función objetivo son convexas.
No convexo: La región factible o la función objetivo no son convexas.
Optimización Determinista vs. Estocástica
Determinista: Los parámetros y datos son conocidos con certeza.
Estocástica: Los parámetros y datos incluyen incertidumbre.
RESUMEN
En este capítulo,
**** Añadir *****
Bibliografía
Aubanell, A., Benseny, A., & Delshams, A. (1993). Útiles básicos de cálculo numérico.
Barcelona: Editorial Labor, S.A.Burden, R. L., & Faires, J. D. (2015). Numerical analysis (10ª ed.). Cengage
Learning..Casas Rentería, Eduardo (2024). Métodos Matemáticos para la Ingeniería. Cálculo Numérico.
Recuperado de:
https://personales.unican.es/casase/MMI/pdfs/Apuntes.pdfArnold, D. N. (2000). The Patriot Missile Failure. Institute for Mathematics and its
Applications. Recuperado de http://www.ima.umn.edu/~arnold/disasters/patriot.htmlArnold, D. N. (2000). The Explosion of the Ariane 5 Rocket. Institute for Mathematics and its
Applications. Recuperado de http://www.ima.umn.edu/~arnold/disasters/ariane.htmlSIAM News. (1996, octubre). Vol. 29, Number 8.Arnold, D. N. (2000). The sinking of the Sleipner A offshore platform. Institute for
Mathematics and its Applications. Recuperado de
http://www.ima.umn.edu/~arnold/disasters/sleipner.htmlInfante del Río, J. A., & Rey Cabezas, J. M. (2022). Métodos numéricos: teoría, problemas y
prácticas con MATLAB. Ediciones Pirámide..Rivera, J. G., Álvarez, E. E., Galo, J. R., & Tabares, H. A. (2017). Métodos numéricos
interactivo. Primera edición. Recuperado de
Métodos Numéricos. InteractivoFaires, J. D., & Burden, R. L. (2012). Numerical methods, 4th. Cengage Learning.
Solución
>
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 para encontrar soluciones con gran rapidez,
utilizando la pendiente de las tangentes de la función.
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.
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
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..
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.
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:
