Título de la obra
Métodos Numéricos Interactivo

Juan Guillermo Rivera Berrío
Elena Esperanza Álvarez Sáiz
José Román Galo Sánchez
Héctor Aníbal Tabares Ospina
Primera edición: 2017

Diseño de cubierta: Diana María Velásquez García
Librería turn.js: Emmanuel García
Herramienta de edición: DescartesJS
Fuente: Amaranth

Fondo Editorial Pascual Bravo
Calle 73 73A-226
PBX: (574) 4480520
Apartado 6564
Medellín, Colombia
www.pascualbravo.edu.co
ISBN: 978-958-58510-6-1

Esta obra está bajo una licencia Creative Commons 4.0 internacional: Reconocimiento-No Comercial-Compartir Igual. Todos los objetos interactivos y los contenidos de esta obra colectiva están protegidos por la Ley de Propiedad Intelectual.

iv

vi

Introducción

Los problemas de la ingeniería, matemáticos, estadísticos o, en general, todo problema que se pueda modelar es posible resolverlo por algún método analítico, gráfico o numérico. Los métodos numéricos permiten resolver muchos problemas usando operaciones aritméticas. Un ejemplo introductorio es el problema planteado por Chapra y Canale (2007, pág. 12Chapra, S., & Canale, R. (2007). Métodos numéricos para ingenieros (Quinta ed.). México: Mc Graw Hill.), que modelamos en la siguiente escena interactiva:

Al hacer clic en el botón Iniciar, observarás la solución analítica del problema y, además, una solución numérica o aproximada conocida como "método de Euler". Este método, que trabajaremos en la última parte de este libro, resuelve una ecuación diferencial utilizando una expansión de la serie de Taylor.

viii

Es de advertir que esta obra no es un libro más de métodos numéricos, pues su contenido se ha diseñado de tal forma que permita al estudiante interactuar con los objetos de conocimiento que en cada apartado se describen y explican. Al interactuar, el estudiante logra una mayor percepción de estos objetos de conocimiento.

Dicho de otra forma, esta obra es el primer libro digital interactivo de métodos numéricos.

Por otra parte, el libro se ajusta al modelo del aula invertida (flipped classroom), pues el usuario del libro puede obtener información en un tiempo y lugar que no depende de la presencia del profesor, además de poder interactuar tanto en ordenadores como en dispositivos móviles, aumentando su percepción y, en consecuencia, su velocidad de aprendizaje.

Esta primera edición la hemos dividido en seis partes. La primera parte trata sobre la teoría de errores. La segunda, presenta los diferentes métodos numéricos empleados para solucionar ecuaciones no lineales. La tercera estudia los algoritmos que solucionan sistemas de ecuaciones lineales. El tema de interpolación se explica en la cuarta parte. La penúltima parte versa sobre diferenciación e integración numérica. La última parte es acerca de la solución numérica de Ecuaciones Diferenciales, que inician con el modelo simple de Euler (utilizado en la escena introductoria) y termina con los métodos de Runge-Kutta.

Actualmente, existe una gran variedad de libros y textos, algunos de ellos reseñados en la bibliografía, que incluyen algoritmos para ser aplicados en diferentes aplicativos. El más usual de estos aplicativos es el Matlab y, en algunos casos, las macros de Visual Basic de Excel. No osbtante esta tendencia, hemos optado por usar aplicaciones libres, que permitan una mayor accesibilidad de nuestros usuarios, el libro mismo tiene esta naturaleza.

x

Introducción

El uso de los ordenadores para resolver problemas por métodos numéricos, pueden generar errores cuyas causas es necesario comprender. Un error puede surgir en la digitación del algoritmo que resuelve un problema, por una mala formulación del modelo, por la elección del modelo (el método de Euler de la introducción) o por problemas propios de la máquina que realiza los cálculos numéricos. De esta última causa nos ocuparemos en esta primera parte del libro, para ello, iniciamos con el siguiente algoritmo escrito en wxMaxima.

14

1. Sistemas numéricos

La humanidad históricamente ha utilizado diferentes sistemas numéricos para representar cantidades, actualmente usamos el sistema decimal que usa 10 dígitos (del cero al nueve). Sin embargo, los ordenadores trabajan con bits (dígitos binarios), conformando grupos denominados palabras, las cuales pueden ir desde ocho bits hasta 64. El manejo de esta gran extensión de bits se facilita al dividirla en grupos de ocho bits, denominados bytes. Internamente, entonces, el ordenador trabaja con bits, no obstante, su interpretación (interface hombre - máquina) se hace tanto en sistemas binarios, como octales y hexadecimales. En este apartado recordaremos algunos aspectos del sistema binario, para los otros sistemas se recomienda consultar la bibliografía.

1.1 Sistema digital (0,1), base 2

Es un sistema de numeración en el que los números se representan utilizando solamente las cifras cero y uno (0 y 1).

Ejemplo

Observa un número binario de varias cifras, haz clic en el botón para observar otros números.

16

Ejemplo

Haz clic en el botón "Siguiente paso" para observar cómo se obtienen los dígitos binarios correspondientes al número decimal dado.

Una vez hayas practicado y comprendido el procedimiento, puedes realizar los ejercicios que desees en la página siguiente, sólo hemos incluido la parte entera.

18

Conversión de binario a decimal

Para realizar la conversión de binario a decimal sigue los siguientes pasos:

1. Inicia por el lado derecho del número binario, cada cifra multiplícala por 2 elevado a la potencia consecutiva (comenzando por la potencia 0, 1, 2, ....).

2. Después de realizar cada una de las multiplicaciones, súmalas y el número resultante será el equivalente al sistema decimal.

20

2.1 Aritmética del punto flotante

La restricción dada por la cantidad de bits permitidos en un ordenador ha obligado a algunas técnicas de almacenamiento de los números, una de ellas es la representación del punto flotante, que para base binaria y un número x dado, sería:

x = m·2^e

En la expresión, e es el exponente o característica y m es conocida como la mantisa, la cual es un número decimal de la forma 0.d₁d₂d₃d₄.... Si los primeros decimales valen cero, esto generarían bits innecesarios, por lo que se recurre a correr (flotar) el punto al primer decimal no nulo. Veamos dos ejemplos:

En el primer ejemplo, la coma la desplazamos dos cifras a la izquierda, por ello, el exponente es 10 (dos en decimal). En el segundo ejemplo, eliminamos el cero innecesario en la primera cifra decimal desplazando la coma una cifra a la derecha, por ello, el exponente es -1. Este proceso se conoce como normalización de la mantisa.

Retornando a los años 80 del siglo pasado, supongamos que nuestro ordenador tiene una capacidad de almacenamiento de 16 bits. En la siguiente escena interactiva, observa cómo se almacena un número entero. El primer bit es el signo del número (uno si es negativo y cero si es positivo), el segundo bit es el signo del exponente, los bits 3 a 6 almacenan el valor el exponente y los restantes el valor de la mantisa.

22

2.2 Errores en un ordenador

Presentamos, a continuación, algunos tipos de errores:

2.2.1 Errores inherentes

Son errores inherentes a los datos, a la incertidumbre en las medidas, a los números que se exceden en el tamaño de bits permitidos por el ordenador, al diseño del modelo, a la digitación del algoritmo, entre otros.

2.2.2 Error absoluto

Si p^* es una aproximación de p y si p es el valor real, entonces,
Error Absoluto = |p-p^*|, o sea el valor absoluto de p menos p^*.

El valor absoluto genera sólo errores positivos, lo que se traduce en una suma o acumulación de errores.

El error absoluto se propaga en los cálculos de la siguiente manera:

Producto E_x*y = yE_x + xE_y
Cociente E_x/y = (E_x/x) + (E_y/y)
Suma E_x+y = E_x + E_y
Resta E_x-y = E_x + E_y

24

2.2.4 Error de redondeo

Este tipo de error ya lo hemos explicado en la introducción de este apartado, recordemos que las limitaciones de almacenamiento de bits generan cálculos aproximados de los números verdaderos, es decir, los números y errores están sujetos al sistema numérico de punto flotante. Similar al sistema binario, un número real positivo puede ser normalizado para que adquiera la forma:
    y = 0.d₁d₂d₃...d_kd_k+1... x10ⁿ
    1≤d_i≤9, con i = 2,3,4,...k
La forma de punto flotante fl(y) se obtiene terminando (recortando) la mantisa de y en k dígitos decimales. Existen dos métodos de terminar:
1. Cortando los dígitos d_k+1d_k+2...   fl(y) = 0.d₁d₂d₃...d_k x10ⁿ
2. Redondeando el número
Si (d_k+1 ≥ 5) Entonces
    Se agrega uno a d_k para obtener fl(y)
FinSi

Si (d_k+1 < 5) Entonces
    Se cortan todos excepto los primeros k dígitos.
FinSi

2.2.5 Error de truncamiento

Similar al error de redondeo es originado en un proceso que requiere un número infinito de pasos pero se detiene en un número finito de pasos. Sin embargo, estos errores no son inherentes al ordenador o al sistema numérico utilizado, pues surgen al usar una aproximación en lugar de un procedimiento matemático exacto (Chapra & Canale, 2007, pág. 78Chapra, S., & Canale, R. (2007). Métodos numéricos para ingenieros (Quinta ed.). México: Mc Graw Hill.). El problema del paracaidista es un ejemplo de ello

26

2.3 Algoritmos y estabilidad

"La condición de un problema matemático relaciona su sensibilidad con los cambios en los datos de entrada. Se dice que un cálculo es numéricamente inestable si la inexactitud de los valores de entrada se aumenta considerablemente por el método numérico" (Chapra & Canale, 2007, pág. 98Chapra, S., & Canale, R. (2007). Métodos numéricos para ingenieros (Quinta ed.). México: Mc Graw Hill.).

Algunas recomendaciones

Evitar restar números casi iguales.
Evitar la división por un número cercano a cero.
Minimizar la ejecución del número de operaciones en el algoritmo.
Utilizar datos del tipo doble precisión cuando sea necesario.
Realizar escalamiento de los datos si es posible.

28

En esta parte del libro estudiaremos los métodos numéricos para determinar las raíces de una Ecuación No Lineal con una variable. Analizaremos los problemas de convergencia y error que puedan presentarse (Burden, 2002, pág. 47Burden, R., & Faires, D. (2002). Análisis numérico (Séptima ed.). México: Thomson Learning.).

En la mayoría de los casos no es posible determinar las raíces de una función, es decir, los valores x^* tal que f(x^* ) = 0, aplicando un algoritmo que termine en un número finito de pasos. Para ello, tenemos que usar métodos de aproximación que, en general, son iterativos.

En el tercer capítulo observaremos cómo podemos determinar si una función no lineal tiene raíces, a través de un análisis gráfico y uno análitico. Al final, haremos una introducción a los métodos iterativos

En los dos capítulos siguientes veremos los métodos iterativos usuales: bisección, regula falsi, punto fijo, Newton y el método de la secante.

Incluiremos algunos objetos interactivos diseñados con el editor Descartes, que permitirán un mayor acercamiento a los conceptos, teoremas y métodos iterativos descritos en cada capítulo. No obstante esta ayuda, se recomienda utilizar el cuadernillo de trabajo para realizar algunos ejercicios propuestos

3. Raíces de una ecuación

Para hallar las raíces de una ecuación es necesario saber si la ecuación tiene solución. En caso de la existencia de una solución, convendría identificar en qué intervalo se encuentra, pues ello optimizaría el algoritmo utilizado para hallar la raíz de la ecuación.

32

En la escena interactiva podrás ver algunas funciones y hallar sus raíces. Puedes escribir tu propia función, pero antes de hacerlo observa cómo se escriben las funciones comparando la expresión que aparece en la parte superior (sintáxis del ordenador) con la expresión de la parte inferior del interactivo (sintáxis matemática).

Haz clic para ampliar la escena

34

3.3 Método iterativo

Los métodos numéricos para resolver ecuaciones no lineales son de tipo iterativo, los cuales son modelos de aproximaciones sucesivas. En una forma simple, el método consiste en elegir un valor inicial x₀ (semilla), a partir del cual se genera una sucesión de iteraciones x₁, x₂, x₃, ..., hasta que éstas converjan a la solución.

Si una ecuación puede ponerse en la forma f(x) = x, podemos, en principio, empezar con un punto x₀ y calcular x₁ sustituyendo en la ecuación el valor de x₀, el valor obtenido será el nuevo valor a sustituir en la ecuación, y así sucesivamente hasta encontrar un valor de convergencia. Matemáticamente, considerando x_n+1 = f(x_n) para n ≥ 1, si f es continua y la secuencia {x_n} con n ≥ 1 converge, lo hará a la solución x.

Tratemos de resolver la siguiente ecuación x² - 2x + 1 = 0, cuya solución sabemos que es uno (1), para ello podemos hacer el siguiente despeje: x = √2x - 1. Supongamos que nuestra semilla es x₀ = 3, lo que significa que x₁ = 2.23607, calculamos x₂ reemplazando el valor anterior en la ecuación y obtenemos x₂ = 1.86337, en este momento hemos realizado sólo dos iteraciones, si continúas iterando notarás que la convergencia es lenta, lo que justifica el uso de otros métodos que veremos en el siguiente capítulo.

En la escena interactiva de la siguiente página podrás encontrar el número de iteraciones necesarias para hallar la solución en la que, para facilitar el proceso, sólo hemos considerado dos cifras decimales. Observa que la columna x_n inicia con el valor semilla, en la columna x_n+1 aparecen los cálculos realizados en la ecuación, los cuales se convierten en una nueva semilla, repitiendo (iterando) el proceso hasta que hallamos la solución.

36

4. Métodos iterativos abiertos

En el último ejemplo de la escena interactiva anterior, observamos que la convergencia a la solución se realizaba con pocas iteraciones, esto depende del tipo de ecuación y del método iterativo utilizado. En cada ejemplo observamos una convergencia a un punto, el cual se suele llamar “punto fijo atractivo”. Los métodos abiertos utilizan una ecuación similar a la que hemos usado en al apartado anterior para predecir la raíz, la cual puede desarrollarse como una iteración simple de punto fijo (también llamada iteración de un punto o sustitución sucesiva o método de punto fijo). En este capítulo analizaremos los métodos de punto fijo, de Newton y el de la Secante.

4.1 Método de punto fijo

Este método lo explicaremos a través de la unidad didáctica interactiva del proyecto Un_100Métodos de punto fijo, la cual inicia con un método para calcular la raíz cuadrada de un número.

Método de Herón para calcular la raíz cuadrada

Herón de Alejandría (126-51 a.C.) afirmó que si y es una aproximación de la raíz cuadrada de un número positivo a entonces:

es una aproximación mejor. En base a ello, partiendo de una aproximación inicial y₀

38

Formulación general del método iterativo

40

Convergencia y divergencia

42

Resolviendo ecuaciones

44

Actividad con GeoGebra

Practica con la siguiente escena creada por Adrian Phelipp Villada

Haz clic para ampliar la escena

46

4.2 Método de Newton o de las tangentes

El método de Newton se describe en su libro Method of Fluxions, el cual fue explicado por J. Raphson en 1690, por lo que se suele llamar el método de Newton-Raphson. Es un método cuya convergencia es mucho mejor que la del punto fijo.

Este método lo explicaremos a través de la unidad didáctica interactiva del proyecto Un_100Resolución numérica de ecuaciones, en la cual se describe otros tres métodos que usaremos en este libro.

Considérese una ecuación no lineal de la forma f(x) = 0, donde f es continua y derivable en un intervalo [a, b] con derivada no nula.

Definimos un método iterativo en el que partiendo de un valor inicial x₀ construimos una sucesión (x_n), donde x_n es el punto de intersección de la recta tangente a f(x) en (x_n-1, f(x_n-1)) con el eje de abscisas. Es decir,

El método de Newton es siempre convergente a una raíz si el valor inicial es lo suficientemente próximo a esa raíz, pero la convergencia es dependiente de ese valor inicial. Por ejemplo, compruebe en la escena interactiva que para f(x) = x² - 2 si consideramos x₀ = 0.5 convergerá a √2 y si consideramos x₀ = -0.5 convergerá a -√2.

48

La dificultad de este método radica en la necesidad de calcular los valores de la derivada primera de la función. Para evitar ello hay variantes del método de Newton que obvian ese cálculo, por ejemplo, el método de la Secante y el de “Regula Falsi” o de la falsa posición que combina el de bisección y el de la secante, los cuales veremos más adelante.

No obstante, en programas de cálculo simbólico o de algebra computacional (CAS) como GeoGebra o Maxima, la derivada de la función no es un problema, como veremos a continuación.

En la escena de GeoGebra de la página siguiente hemos usado la misma función de la escena anterior, desplaza uno de los puntos sobre el segmento azul. Puedes observar que a medida que el valor incial esta cerca de una de las dos raíces, el número de iteraciones se reduce. La escala de la gráfica se puede modificar con la rueda del ratón (mouse)".

Ejercicios

1. Obtener la solución única de x³ + 4x² - 10 = 0 mediante el método de Newton. Compara los resultados obtenidos interactuando con la escena de Descartes.

2. Resuelve las ecuaciones planteadas en el apartado 3.3. Puedes hacerlo con Descartes, GeoGebra o Maxima

3. Usa una calculadora para resolver, paso a paso, la ecuación dada en el ejercicio 1.

50

Programación del método de Newton en wxMaxima

A continuación se presenta el código en wxMaxima para el método de Newton.

/* Definimos en la variable g la función a resolver */
g:x^2-2;
/* Definimos en la variable gd la derivada de la función */
gd:diff(g,x);
/* Para empezar, eliminamos las asignaciones a la variable newton */
kill(newton);
/* Creamos la función de iteración (newton(pn)) con x = pn */
newton(pn):=ev(x-g/gd,x=pn);
/* La primera iteración en pn = 0.2 */
pn:.2;
/* Repetimos el proceso desde 1 hasta 10 iteraciones */
i:0$
while i<10 do (
    i:i+1,
    print(i,float(pn)),
    pn:newton(pn));
/* Dibujamos la función (Paso opcional)*/
plot2d([x^2-2],[x,-2,2]);

Ejecuta el programa combinando las teclas control + R. En la siguiente página se muestran los resultados, observa que en la octava iteración se obtiene la solución. Puedes digitar el código, copiarlo de esta página o hacer clic en la imagen para descargar el archivo.

52

4.3 Método de la Secante

En el método de Newton, para obviar la dificultad de calcular la derivada de la función en cada punto se considera una aproximación

Con lo que el método iterativo quedaría expresado como

Para poder aplicar el método son necesarias dos aproximaciones iniciales que permitan iniciar el proceso. Hay diversas alternativas como puede ser considerar una aproximación inicial y obtener la segunda aplicando el método de Newton. Aquí consideraremos un intervalo inicial y sus extremos serán esas dos primeras aproximaciones.

El nombre del método es debido a que gráficamente se está sustituyendo la recta tangente del método de Newton por la recta secante.

El método de la Secante puede formularse de manera independiente del método de Newton y como mejora del método de la bisección (ver apartado 5.2) donde en lugar del punto medio del intervalo se selecciona un valor proporcional a los valores de la función en los extremos del intervalo, y de aquí que el método se denomine también “de las partes proporcionales”.

54

El método de la secante como variación del método de Newton busca evitar el cálculo de la derivada de la función, aproximando la pendiente a la recta secante que une los puntos (x₁, f(x₁)) y (x₂, f(x₂)) a la recta tangente en el punto de intersección de la función con el eje x.

En la escena de GeoGebra de la página siguiente, aproxima los puntos x₁ y x₂ a la intersección y observa cómo la solución converge. Desplaza el deslizador para que observes paso a paso la aproximación de la secante a la tangente.

Ejercicios

1. Obtener la solución única de x³ + 4x² - 10 = 0 mediante el método de la Secante. Compara los resultados obtenidos interactuando con la escena de Descartes.

2. Resuelve las ecuaciones planteadas en el apartado 3.3. Puedes hacerlo con Descartes, GeoGebra o Maxima

3. Halla las raíces de la función f(x) =−11−22x+17x² −2.5x³

56

Programación del método de la Secante en wxMaxima

A continuación se presenta el código en wxMaxima para el método de la Secante.

/* Creamos la función de iteración (secante(arg)) */
secante(arg):=block([a:arg[1],b:arg[2],f:arg[3],xsec], xsec:(a*ev(f,x=b)-b*ev(f,x=a))/(ev(f,x=b)-ev(f,x=a)), [b,xsec,f] );
/* La primera iteración en a = 1, b = 2 y f = x^2 - 2 */
ini:[1,2,x^2-2];
/* Repetimos el proceso desde 1 hasta 10 iteraciones */
for i:1 thru 10 do (print(i,float(ini),float(ini:secante(ini))));
/* Dibujamos la función (Paso opcional)*/
plot2d([x^2-2],[x,-2,2]);

En el código hemos empleado la estructura tipo bloque, que permite evaluar varias expresiones y devolver el último resultado. Por ejemplo, la estructura block([a,b],a:2,b:3,a+b); devolvería el valor de 5.

Recuerda de ejecutar el programa con la combinación de las teclas control + R. En la siguiente página se muestran los resultados para dos funciones diferentes, observa que la primera salida es el bloque donde aparece la expresión de la pendiente. Luego aparecen pares de bloque así: [a, b, f(x)] [b, f(b), f(x)].

Puedes digitar el código, copiarlo de esta página o hacer clic en la imagen para descargar el archivo.

58

5. Métodos iterativos cerrados

Los métodos cerrados están basados en el teorema de Bolzano y la convergencia está garantizada bajo las condiciones del teorema. Los métodos abiertos, anteriormente explicados, no nos garantizan la convergencia, por lo que deberemos definir un número máximo de iteraciones y elegir en el proceso iterativo un buen punto inicial. Los métodos cerrados utilizan intervalos dentro de los que se debe encontrar la raíz de la ecuación, es decir, en ese intervalo la función cambia de signo. El intervalo se determina por medio de la gráfica de la ecuación o del método de búsqueda incremental. En este capítulo trabajaremos con el método de la Bisección y el de Regula Falsi.

5.1 Teorema de Bolzano

Si una función f(x) es continua en un intervalo cerrado [a,b] y f(a) y f(b) son de distinto signo, existe por lo menos un punto entre a y b para el cual f(c)=0.

5.2 Método de la Bisección

También conocido como método de la Bipartición o de la Dicotomía. Para este método presentaremos la programación del código en las herramientas de autor GeoGebra y Descartes, que serán una base conceptual para programar los demás métodos. Igualmente, si ya descargaste este libro a tu ordenador, es posible explorar y analizar los códigos de programación de los métodos anteriores, pues estos se encuentran almacenados en la carpeta interactivos.

60

Actividad con Descartes

Escena creada por José R, Galo, cambia funciones y compara con la escena de GeoGebra.

Haz clic para ampliar la escena

62

Programación del Método de la Bisección con Descartes

En el siguiente vídeo podrás observar cómo se programa una escena de la Bisección con Descartes. Sugerimos verlo en forma ampliada y en el navegador Firefox.

Haz clic para ampliar la escena

64

Programación del método de la Bisección en wxMaxima

A continuación se presenta el código en wxMaxima para el método de la Bisección.

/* Definimos el número de decimales a mostrar */
decimales:4;
prec:10^decimales;
/* Definimos la función f(x) y el intervalo [a, b] */
f(x):=cos(x)-x;
a:0.0;
b:2.0;
/* Usamos un bucle tipo for para las iteraciones */
for i:1 thru 15 do
(
   c:(a+b)/2, /* calculamos el punto medio */
   if f(a)*f(c)>0 /* ¿Son del mismo signo en [a,c]? */
     ; ; then a:c /* elegimos [c,b] */
     ; ; else b:c, /* elegimos [a,c] */
   print(i,float(round(a*prec)/prec),float(round(b*prec)/prec)) /* escribimos los resultados con cuatro decimales */
)$

En el código hemos usado el comando round(), el cual trunca los decimales. Recuerda ejecutar el programa con la combinación de las teclas control + R. Puedes digitar el código, copiarlo de esta página o hacer clic en la imagen para descargar el archivo. Recuerda que puedes incluir la gráfica con una expresión como plot2d(f(x),[x,0,2]);

66

Actividad con Descartes

Escena creada por José R, Galo, cambia funciones y compara con la escena de GeoGebra.

Haz clic para ampliar la escena

68

Introducción

Una ecuación lineal, que denotaremos como E_1, con n variables o incógnitas x₁,x₂,...,x_n es de la forma

E₁: a₁x₁ + a₂x₂ + ... + a_nx_n = b

Tanto los coeficientes (a₁,a₂,...,a_n ) como el término independiente b son números reales. La solución de la ecuación lineal es un conjunto de valores para las variables o incógnitas que satisfacen la ecuación. Por otra parte, los sistemas de ecuaciones lineales que surgen al modelar problemas de diferentes áreas del conocimiento, en especial las disciplinas ingenieriles, son un conjunto de m ecuaciones lineales E₁,E₂,...,E_m que se pueden representar así:

E₁ : a₁₁x₁ + a₁₂x₂ + a₁₃x₃ + ... + a_1nx_n = b₁

E₂ : a₂₁x₁ + a₂₂x₂ + a₂₃x₃ + ... + a_2nx_n = b₂

E₃ : a₃₁x₁ + a₃₂x₂ + a₃₃x₃ + ... + a_3nx_n = b₃

E_m : a_m1x₁ + a_m2x₂ + a_m3x₃ + ... + a_mnx_n = b_m

A este conjunto de ecuaciones se le llama sistema lineal de dimensión m x n. El sistema lineal se puede escribir de la forma Ax = b, donde A es la matriz de coeficientes, b es el vector de términos independientes y x = (x₁,x₂,...,x_n), es la solución.

72

6.1 Representación matricial del sistema

El sistema lineal, como vimos antes, se puede escribir de la forma Ax = b. En el siguiente interactivo puedes observar la matriz representativa del sistema.

74

6.2 Tipos de matrices

Existen varios tipos de matrices de acuerdo a su forma y a su contenido. Por ejemplo, la matriz M = [3  2   1  2   4  1],  se conoce como matriz fila, por su forma. Otros tipos de matrices las puedes reconocer fácilmente, seleccionando algunas de las opciones que se presentan en el siguiente interactivo. Cuando aparezca el botón "Otro ejemplo", puedes hacer clic varias veces, para observar varias matrices.

76

Ejercicio

En este ejercicio practicaremos la suma de matrices a través de colores. El propósito es recordar el procedimiento de la suma.
La suma de A y B da como resultado otra matriz de colores. Recuerda que dos colores iguales dan el mismo color, para colores diferentes, ten en cuenta:

Amarillo + azul = verde, rojo + amarillo = naranja, rojo + blanco = rosa, azul + blanco = turquesa, rojo + azul = magenta y blanco + amarillo = amarillo claro

78

Ejercicio

En este ejercicio practicaremos la suma de matrices y el producto de un escalar por una matriz. Realiza las operaciones en tu cuaderno y luego haz clic en "Verificar".
Observa que las matrices tienen un contenedor diferente al que veníamos usando (paréntesis en lugar de corchetes), esta representación es usual en algunos textos de métodos numéricos.

80

6.6 Matriz transpuesta

La matriz traspuesta se utiliza con frecuencia en las operaciones entre matrices. Los elementos a_ij de una matriz A, se cambian por a_ji y se obtiene la transpuesta de A.

La matriz transpuesta es una de las operaciones más fáciles de realizar, basta con intercambiar las filas por las columnas y las columnas por las filas. Pruébalo con la siguiente escena:

82

6.7 Determinante de una matriz

El determinante de una matriz cuadrada A es un número real asignado a ella. Si A es una matriz de orden n, su determinante lo denotaremos por det(A) o también por |A| (no confundir con el valor absoluto). En la escena interactiva se explica, paso a paso, el cálculo del determinante de matrices de orden 2, 3 y 4. El método utilizado se conoce como el desarrollo por cofactores de primera fila. Para ello, se toma la primera fila, multiplicando cada elemento a_1i por su cofactorEs importante aclarar que en los textos de Álgebra Lineal existen dos versiones. La primera versión usa el término "cofactores" para el cálculo de determinantes y de la matriz adjunta; la segunda, usa el término "adjuntos". En este libro usaremos la primera versión, más utilizada en los libros clásicos de Álgebra Lineal., es decir, el determinante de la matriz que se obtiene eliminando la fila y columna correspondiente a dicho elemento multiplicado por (-1)¹⁺ⁱ. La suma de todos los productos es igual al determinante.

84

Cálculo matricial con wxMaxima

La sintaxis en Maxima es muy similar a la de GeoGebra, tal como se observa en la imagen, con algunas diferencias, tales como: Maxima exige el uso del punto (.) para el cálculo de productos de matrices, los comándos transpuesta y determinante están en inglés. Observa, además, que GeoGebra usa paréntesis para representar una matriz, mientras que Maxima usa corchetes.

86

Cuando un sistema lineal presenta una matriz de coeficientes de la forma triangular superior, la solución del sistema es inmediata, usando una sustitución regresiva. En el siguiente interactivo puedes generar varios ejemplos de un sistema 4x4, observa lo directa que es su solución. En el siguiente apartado, obtener un sistema triangular es el principal objetivo del método de solución directa.

88

Haz clic para ampliar la escena

90

Programación del método de Gauss en wxMaxima

En estas dos páginas presentamos el código para el método de eliminación de Gauss. Un primer bloque son las declaraciones, en el que se ingresan la matriz de coeficientes A y los vectores B y X, igualmente, hemos declarado el vector X0 en forma simbólica para la presentación de resultados. El tamaño de la matriz se obtiene con el comando length. Para la presentación de la matriz A, hemos usado cuatro decimales de aproximación, que se definen en las variables: decimales y prec. Luego de ejecutar el código, los resultados aparecen de esta forma:

92

7.3 Método de Gauss con pivoteo parcial

Cuando un elemento pivote es cero se presentan problemas, pues se origina una división entre cero, igual ocurre cuando el pivote es cercano a cero que pueden introducir errores de redondeo. En el siguiente interactivo haz clic en el botón Triangularizar, luego intenta con el botón Pivote cero y observa los mensajes de error que se presentan.

94

7.4 Método de Gauss-Jordan

Este método es una extensión del método anterior. "La principal diferencia consiste en que cuando una incógnita se elimina en el método de Gauss-Jordan, ésta es eliminada de todas las otras ecuaciones, no sólo de las subsecuentes. Además, todos los renglones se normalizan al dividirlos entre su elemento pivote. De esta forma, el paso de eliminación genera una matriz identidad en vez de una triangular" (Chapra & Canale, 2007, pág. 277Chapra, S., & Canale, R. (2007). Métodos numéricos para ingenieros (Quinta ed.). México: Mc Graw Hill.). En la siguiente escena de GeoGebra practica cambiando los valores de los coeficientes y del vector B en la hoja de cálculo.

Haz clic para ampliar la escena

96

El sistema anterior lo podemos expresar en notación matricial, así: UX - D = 0. Ahora, supongamos que existe una matriz diagonal inferior L con la diagonal compuesta de elementos de valor uno:

Se puede demostrar la siguiente igualdad: L(UX - D) = AX - B, lo que nos permite concluir que:

LU = A (Véase la escena de la siguiente página)

LD = B

Para resolver el sistema, se siguen los siguientes pasos:

a) Factorización LU. Se descompone A en las matrices L y U. Donde L, como se dijo antes, es una matriz triangular inferior con los elementos de la diagonal iguales a uno (1), y U es una matriz triangular superior que se obtiene después de la eliminación hacia adelante.

b) Determinación del vector D mediante una sustitución hacia adelante.

c) Determinación del vector X por sustitución regresiva.

98

Factorización LU de Doolittle con wxMaxima

Presentamos el algoritmo en wxMaxima para hallar las matrices L y U. Usamos las ecuaciones del primer ejemplo de la escena de Descartes, para efectos de comparación.

Los resultados obtenidos son los siguientes:

100

Método de Cholesky en Descartes

En la siguiente escena se explica, paso a paso, la solución de un sistema lineal con el método de Cholesky. El primer ejercicio presenta una matriz simétrica igual al de la imagen anterior de wxMaxima.

102

Matriz inversa desde la matriz adjunta

Al dividir la matriz adjunta (Adj[A]) de una matriz cuadrada A entre su determinante, se obtiene la matriz inversa.

Es importante aclarar que en los textos de Álgebra Lineal existen dos versiones para calcular la martiz adjunta. La primera versión corresponde a la matriz de cofactores traspuesta; la segunda, se corresponde a la matriz de cofactores. En este libro usaremos la primera versión, más utilizada en los libros clásicos de Álgebra Lineal. En la siguiente escena de Descartes, puedes observar la explicación de esta primera versión.

104

En la siguiente escena puedes observar, paso a paso, el cálculo de la matriz inversa de una matriz A de orden 3. Obviamente, el procedimiento para ecuaciones lineales de grandes dimensiones requiere del diseño de un algoritmo más complejo y la formulación de diferentes modelos, como son los métodos iterativos.

Observa que inicialmente se calcula el determinante de la matriz, pues de ser cero no sería posible calcular la inversa.

106

7.8 Solución del sistema con la matriz inversa

Con lo aprendido hasta aquí, puedes encontrar los valores de las incógnitas que satisfacen las ecuaciones cuando el sistema es compatible y determinado (o cuando la solución existe y es única). En la escena de abajo, tenemos un sistema de 3x3 (tres ecuaciones y tres incógnitas). Haz clic en el botón para que observes varios ejemplos.

108

Resolución del sistema con wxMaxima

En estas dos páginas presentamos el código para la solución de un sistema de ecuaciones lineales usando la matriz inversa. Hemos incluido la matriz adjunta para efectos de comparación con las escenas anteriores, igualmente, se incluye la matriz inversa calculada por el método de la adjunta y por el método de factorización LU. En la página siguiente se muestran los resultados de un ejemplo de un sistema de 3x3. No obstante, es posible modificar el archivo para sistemas de un orden mayor.

110

Ejercicio

Para el ejercicio puedes usar GeoGebra, wxMaxima, una hoja de cálculo como Excel o en forma manual con tu calculadora.

Haz clic para ampliar el vídeo

112

8.1 Método de Jacobi

Este método parte de la descomposición de la matriz A en la suma de tres matrices A = L + D + U, donde D es diagonal, L triangular inferior y U triangular superior, formadas con los elementos respectivos de A. Así las cosas, en el método de Jacobi:

M = D    y     N = L + U

Reemplazando en el sistema Ax = B,

Dx^(k+1) = -(L + U)x^(k) + B

Cuyos componentes se calculan, así:

Veamos un ejemplo sencillo con el siguiente sistema:

Al cual le hacemos los siguientes despejes de variables:

114

En la siguiente escena interactiva de Descartes hemos dejado como primer ejercicio el calculado en Excel, la convergencia se halla en la iteración 12 porque la aproximación utilizada es de tres cifras decimales.

Haz clic para ampliar la escena

116

Método de Jacobi con wxMaxima

En este algoritmo hemos usado el mismo sistema de ecuaciones para efectos de comparación con los resultados obtenidos en Descartes y Excel. Adrede, no hemos cerrado el bloque de instrucciones del bloque jacobi(A,B,n) con el signo $, con el fin de observar en los resultados la expresión general de Jacobi, para ello, hemos invocado el comando sum(), el cual puedes cambiar o reprogramar con un bucle tipo for, en el caso de usar este algoritmo en otra herramienta de programación.

118

8.1 Método de Gauss-Seidel

Este método es muy semejante al método de Jacobi, pero con una convergencia (si la hay) en menos iteraciones. Recuerda que en Jacobi usamos el valor de las incógnitas calculadas para determinar una nueva aproximación, en el de Gauss-Seidel se van utilizando los valores de las incógnitas recien calculados en la misma iteración, y no en la siguiente. Observa, con un ejemplo, el procedimiento empleado:

Sea el siguiente sistema de ecuaciones:

Despejamos la incognita sobre la diagonal, para cada ecuación:

Suponiendo la primera aproximación con x₂ = 0 y x₃ = 0, obtenemos:

120

En la siguiente escena interactiva de Descartes hemos dejado como primer ejercicio el calculado en Excel, la convergencia se halla en la iteración siete porque la aproximación utilizada es de tres cifras decimales.

Haz clic para ampliar la escena

122

Método de Gauss-Seidel con wxMaxima

En este algoritmo hemos usado el mismo sistema de ecuaciones para efectos de comparación con los resultados obtenidos en Descartes y Excel. Descarga el archivo y ejecuta el algoritmo.

124

Introducción

La interpolación se usa para obtener el valor de una función en un punto cuando no se conoce dicha función o, en algunos casos, cuando la función es difícil de evaluar. En otras palabras, si tenemos un conjunto de datos previamente validados, el problema de la interpolación consiste en hallar una función válida para esos datos y fácil de evaluar para el punto que origina el problema. En ingeniería y en algunas ciencias de base empírica es común obtener, por experimentación, un número de puntos que dan respuesta al comportamiento del fenómeno estudiado, a partir de estos datos se procede a construir una función llamada interpolante que los ajuste.

Se trata, entonces, a partir de un conjunto de puntos o nodos (x_i,y_i), obtener una función f que valide f(x_i) = y_i, con i = 1, ..., n. Algunas funciones interpolantes se obtienen por interpolación lineal, por interpolación polinómica de Newton o Lagrange (como el caso de la imagen), la interpolación por medio de splines (a trozos) o la interpolación polinómica de Hermite.

128

9. Polinomios interpolantes

La estimación de valores intermedios entre datos definidos por puntos, se puede lograr con una interpolación polinomial, teniendo en cuenta que la fórmula general para un polinomio de n-ésimo grado es

f(x) = a₀ + a₁x + a₂x² + ··· + a_nxⁿ

Si tenemos n + 1 puntos, podemos hallar un polinomio de grado n que contiene esos puntos. Para dos puntos, por ejemplo, la interpolación será lineal (polinomio de primer grado); para tres puntos será una interpolación cuadrática, pues una parábola los une. En este apartado veremos algunos métodos de interpolación polinómica.

9.1 Método de Newton - Diferencias divididas

El polinomio de interpolación de Newton en diferencias divididas, en general, se puede representar por la siguiente expresión (Véase Chapra y Canale, pág 506):

f_n(x) = b₀ + b₁(x - x₀) + b₂(x - x₀)(x - x₁)+ ··· + b_n(x - x₀)(x - x₁)···(x - x_n)

Para el caso de una interpolación cuadrática f₂(x) = b₀ + b₁(x - x₀) + b₂(x - x₀)(x - x₁), que al multiplicar los términos y agrupar términos semejantes, podemos hallar los coeficientes del polinomio, así:

a₀ = b₀ – b₁x₀ + b₂x₀x₁
a₁ = b₁ – b₂x₀ – b₂x₁
a₂ = b₂

130

La primera diferencia dividida finita, en forma general, se puede representar así:

La segunda diferencia dividida finita representa la diferencia de las dos primeras diferencias divididas:

La n-ésima diferencia dividida finita es:

Con estas diferencias divididas obtenemos el polinomio de interpolación de Newton en diferencias divididas:

En la página siguiente, hemos diseñado un algoritmo para la determinación de los polinomios de interpolación de grado uno hasta el cuatro en wxMaxima, utilizando datos de la función logaritmo natural. Las diferencias divididas las hemos representado como f012···n, por ejemplo:

f012 = f[x₂, x₁, x₀]. Luego, en Descartes, usamos el mismo algoritmo, permitiendo, además, ingresar un dato a interpolar.

132

Polinomio interpolante de Newton en Descartes

En la siguiente escena puedes observar polinomios desde el grado uno hasta el siete. La función es el logaritmo natural, practica con un valor a interpolar (dos, por ejemplo) y observa que la mejor aproximación es en el polinomio de grado cinco.

134

Para el polinomio de orden 2, la matriz obtenida es:

De la fila 1, concluimos: b₀ = 0, b₁ = 0.4621 y b₂ = -0.0519 por lo tanto:

p(x) = .4621(x-1) - 0.0519(x-1)(x-4).

En la escena de la siguiente página hemos usado el algoritmo para calcular los coeficientes del polinomio interpolante de Newton (b₀, b₁, ···, b_n) y, en consecuencia, muestra el polinomio por diferencias divididas. Por la extensión del polinomio, sólo mostramos hasta el grado tres.

Se incluyen cuatro ejemplos adicionales. El primero corresponde a los datos de la función f(x) = ln x, que hemos venido trabajando; en el último, con una serie de datos experimentales, podrás notar que los polinomios de orden 1 y 2 son lineales.

¿Qué concluyes con el polinomio de orden 3?

136

Polinomios de interpolación con GeoGebra

GeoGebra ofrece el comando Polinomio[ Lista de Puntos ], el cual crea la interpolación polinomial de grado n-1 a través de los n puntos listados. La documentación ofrecida por GeoGebra no explica el método utilizado para calcular el polinomio.

Hemos diseñado una escena utilizando este comando (observa la página siguiente), con los datos de la función f(x) = ln x, los cuales puedes cambiar en la hoja de cálculo de la derecha. Para interactuar con la escena, se han habilitado la barra de entrada y el clic derecho, de tal forma que puedas:

1. Incluir un punto a interpolar, para ello, escribe el valor en la celda A9 de la hoja de cálculo, y en la celda B9: +g(A9), en el caso del polinomio de orden 1, o +h(A9) para el polinomio de orden 2, y así sucesivamente.

2. Selecciona las dos celdas y haz clic derecho, luego escoge la opción Crea Lista de puntos, para que se muestre en la gráfica.

3. En la casilla de entrada, puedes crear más polinomios con el comando polinomio[{A,B,C,D,E,...}].

4. Puedes ampliar el ancho de la hoja de cálculo y con clic derecho abrir la ventana de propiedades para modificar colores, tamaños, etcétera, hazlo, preferiblemente con la escena ampliada. Ahora, si tienes instalado el software de GeoGebra, puedes abrir esta escena desde la carpeta interactivos y el archivo geogebra13.html.

138

9.2 Polinomios de interpolación de Lagrange

El polinomio de interpolación de Lagrange es una variante del polinomio de Newton, cuya representación es (Véase Chapra y Canale, pág 516):

en la cual

Para el polinomio de primer orden:

Para el polinomio de segundo orden:

140

La salida de la instrucción pol1(x);,es el polinomio interpolante de Lagrange::

que en forma expandida, sería:

Para los datos dados en el algoritmo, la gráfica del polinomio que muestra wxMaxima es:

142

9.3 Interpolación iterada de Neville

El método de Neville se suele usar cuando nos interesa obtener únicamente el valor numérico en un punto dado. Este método interpola un valor dado con polinomios que crecen hasta que los valores estén suficientemente cercanos, lo que nos permite, por inspección, estimar el valor que presente menor variación.

Tomando la notación de Mora (2013, pág. 36Mora, W. (2013). Introducción a los métodos numéricos. Obtenido de https://tecdigital.tec.ac.cr/revistamatematica/Libros /WMora_MetodosNumericos /WMora-ITCR-MetodosNumericos.pdf) con P_0,1 como polinomio que pasa por (x₀, y₀), (x₁, y₁); P_0,1,2 como polinomio que pasa por (x₀, y₀), (x₁, y₁), (x₂, y₂), P_0,1,2,3 como polinomio que pasa por (x₀, y₀), (x₁, y₁), (x₂, y₂)(x₃, y₃),..., el algoritmo sería el siguiente:

La matriz P es triangular superior, en la cual la diagonal principal contiene los valores de y_i, la columna 0 el polinomio P₀ (P_0,1, P_0,1,2, ...), la columna 1 los polinomios P₁, ···, la columna 4 los polinomios P₄.

Hemos diseñado una escena interactiva de Descartes que muestra la matriz P, en color azul los términos y_i y en rojo la mejor aproximación que, en la mayoría de los casos, corresponde a P_0,4 o, en la notación de Mora, P_0,1,2,3,4. Cuando n es igual a cuatro, mostrará el valor interpolado P(x) calculado por método de Neville y por Lagrange (L(x)). El último ejemplo corresponde a una ecuación cuadrática que se desarrolla con GeoGebra en la siguiente escena interactiva.

144

Método de Neville con GeoGebra

Haz clic para ampliar la escena

146

Los polinomios S_i−1 y S_i interpolan el mismo valor en el punto x_i, es decir, se cumple:
S_i-1(x_i) = y_i = S_i(x_i), para (1 ≤ i ≤ n-1), que garantiza la continuidad de f en todo el intervalo. Suponiendo, además, que S' y S" son continuas, es posible deducir una expresión como:

En la cual h_i = x_i+1 - x_i y z₀, z₁, ···, z_n es un vector de incognitas que se debe resolver con alguno de los métodos de resolución de sistemas de ecuaciones lineales. No obstante, Mora presenta el siguiente algoritmo, que permite calcular los coeficientes a_i, b_i, c_i y d_i de la función splin antes definida.

148

Imagen del primer ejemplo. Observa las oscilaciones de Lagrange (color verde) en los puntos extremos (haz clic sobre la imagen).

Imagen del cuarto ejemplo. Grandes oscilaciones de Lagrange (color verde) en los últimos nodos (haz clic sobre la imagen).

150

Splines cúbicos en GeoGebra

GeoGebra ofrece un comando para calcular una función splin, dada una nube de puntos; sin embargo, no se sabe qué tipo de trazador o qué método es el que se usa para construirlo, sólo se tiene la siguiente información: "Crea una spline que pase por todos los puntos". Hemos creado una escena con los 10 puntos del primer ejemplo de la escena anterior y usado el comando spline[{A,B,C,D,E,F,G,H,I,L}]. Observa el resultado y concluye (puedes mover los puntos).

152

Tres de los splines calculados son los siguientes:

La gráfica muestra los splines, la función original y, para comparar, el polinomio interpolante de Lagrange. Observa las oscilaciones presentadas en este último (haz clic en la imagen para agrandarla).

154

Un resumen del vídeo anterior lo puedes observar en esta presentación:

En la página siguiente, puedes interactuar con una escena diseñada por José Luis Abreu.

156

La introducción anterior la hicimos desde la estadística. Ahora, presentamos una explicación desde los métodos numéricos, realizada por Elena Álvarez:

158

Proyección de un vector sobre un subespacio

Dado un vector v y un subespacio S, el vector v se puede descomponer de manera única de la forma;

v = v₁ + v₂ con v₁ ∈ S y v₂ ortogonal a S

La componente v₁ se llama proyección ortogonal de v sobre S, la cual se denota por proy_S v

Se calcula de la forma

proy_Sv = proy_u1v + ... + proy_unv

siendo u₁, ..., u_n una base ortogonal de S.

Teniendo en cuenta como se obtiene la proyección de un vector sobre otro, la proyección se puede obtener de la forma

160

162

164

Ejemplo 2. Según la primera ley de Kepler, un cometa debe tener una órbita elíptica o hiperbólica. En coordenadas polares, la posición de un cometa, (r, θ), satisface la ecuación r = β + e(r cos φ), donde β es una constante y e la excentricidad de la cónica (si 0 ≤ e ≤ 1 es una elipse, si e=1 es una parábola y si e>1 es una hipérbola).

Suponiendo que las observaciones de un cometa proporcionan los datos que se indican en la escena de la página siguiente,

¿Qué órbita describe el planeta?

166

Mínimos cuadrados en wxMaxima

En wxMaxima resolveremos el primer ejemplo relacionado con el muelle o resorte.

Los datos los podemos ingresar como una matriz, sin embargo, se recomienda hacerlo en dos vectores, así:

F : [0.2,0.4,0.5,0.7,0.9,1];
L : [0.06,0.12,0.15,0.21,0.26,0.29];

que corresponden a los del ejemplo. Luego construimos la matriz de esta forma: A : transpose(matrix(F, L));. Finalmente, con la instrucción simple_linear_regression(A);, obtenemos el modelo de regresión, cuyo resultado sería:

168

En general, wxMaxima estima por mínimos cuadrados modelos no lineales, a través de la expresión: lsquares_estimates(A, [X, Y], Y = a*X^b, [a, b]);, que aplicada al ejemplo anterior, nos daría como resultados a = 0.2895 y b = 0.9515, es decir, un modelo posible para el ejemplo sería: L = 0.2895F^0.9515. Si dibujamos esta ecuación, con la expresión: wxdraw2d( points(A), color=red, explicit(0.2895*X^0.9515, X, 0, 1.2),title="Ejemplo 1",xlabel="Fuerza",ylabel="Alargamiento" )$, obtendríamos:

que es muy cercana al modelo de regresión lineal.

170

Introducción

Observa el siguiente vídeo:

Haz clic para ampliar el vídeo

174

En muchas situaciones nos encontramos con una serie de datos experimentales, que configuran curvas de carga que no tienen una función asociada, como en el ejemplo anterior. En muchos casos se idealiza la curva, con el fin de poder estimar, en forma algo aproximada, el comportamiento de un fenómeno dado. No obstante, en otros casos se requiere de un mayor rigor en los cálculos, pues de ello depende que no fallen los diseños en la ingeniería, los cálculos financieros, cálculo de probabilidades, entre muchos otras situaciones. Nos interesa, entonces, estimar los valores de la derivada o la integral de la función que responde a estos fenómenos, así la función sea difícil de determinar analíticamente y sea conveniente conocer su valor lo más aproximado posible. Es decir, estamos obligados a estimar y a aproximar pero, con un mayor rigor. En los estudios del comportamiento de materiales para estructuras y de las estructuras mismas, por ejemplo, se suelen usar factores de seguridad, pues las aproximaciones generan incertidumbre.

En esta parte del libro conoceremos y aplicaremos fórmulas de integración y derivación numérica, que dan respuesta a las aproximaciones requeridas. Contrario a lo que ocurre en los cursos regulares de Cálculo, iniciaremos con la integración numérica con las fórmulas de Newton-Cotes, las cuales usan como estrategia el reemplazo de una función no determinada o una nube de puntos por un polinomio interpolante. Algunas de estas fórmulas son la regla del trapecio (polinomio de primer grado) y la regla de Simpson (polinomios de segundo y tercer grado) y, para una mejor aproximación, el método de Romberg.

Por otra parte, las fórmulas de diferenciación numérica se utilizan en el desarrollo de algoritmos, para resolver numéricamente problemas de contorno de ecuaciones diferenciales ordinarias y ecuaciones en derivadas parciales. Para este caso, mostraremos diferentes fórmulas de derivación numérica (diferencias progresivas, regresivas y centradas), justificando su validez a partir de desarrollos en series de Taylor.

176

10.2 Fórmula del punto medio

En un curso de cálculo integral utilizamos la suma de Riemann para estimar el área de la región limitada por gráfica de la función f, el eje x y las rectas cuyas ecuaciones son
x = a y x = b.

Para encontrar esta área, dividimos el intervalo [a, b] en n subintervalos y, luego, calculamos los n rectángulos de área f(c_k)•Δx_k. Tras un análisis simple, deducimos la siguiente expresión:

En las siguiente escena interactiva de Descartes, observa cómo se calcula la integral utilizando la fórmula del punto medio; obviamente, con f conocida.

178

10.3 Fórmula del trapecio

La fórmula de punto medio lo que hace es corregir, en parte, los errores de la suma de Riemann cuando se configuran rectángulos por exceso o por defecto respecto al área delimitada por la curva de la curva. Observa, en la escena de GeoGebra, la suma de rectángulos con extremos izquierdos o derechos

Se podría pensar que bastaría con sumar los valores y promediar, para obtener el mismo resultado de la fórmula del punto medio, pero lo que realmente se obtendría es la fórmula del trapecio.

180

Fórmula de cuadratura

Ahora, se podría pensar que podríamos sumar los valores obtenidos con las fórmulas del punto medio y del trapecio y luego promediar, para obtener una mejor aproximación. Este razonamiento nos acerca a la fórmula de Simpson.

182

En la siguiente página presentamos el algoritmo en wxMaxima, para este mismo ejemplo. Observa los resultados para las fórmulas del punto intermedio, la del trapecio y la de Simpson.

184

Vamos a practicar con las fórmulas del trapecio y la de Simpson. En las siguiente escena puedes ingresar datos, para realizar cálculos y comparaciones, para ello:

1. Define el conjunto de puntos equidistantes x₀, x₁, ..., x_n, de la forma: x_k = x₀ + kh, siendo k un número natural desde 0 hasta n.

2. Introduce los valores de las ordenadas correspondientes a los puntos x₀, x₁, ..., x_n. Puedes generar estos datos, a partir de una función, pulsando el botón Generar datos

186

Observa los resultados con las diferentes fórmulas estudiadas en la escena de GeoGebra.

188

10.6 La derivada y la secante

Dada una función de una variable, y = f(x), se define la derivada primera en un punto a de su dominio como:

Se puede obtener una aproximación de la derivada en el punto a mediante el cociente incremental, tomando un valor de h pequeño:

En la escena interactiva de la siguiente página, puedes practicar aproximando el cociente incremental (de color azul) al reducir el tamaño de h. Puedes, también, pulsar el botón Play para ver cómo al hacer h pequeño, se aproxima mejor al valor de la derivada.

Prueba con otras funciones y sus derivadas, ingresándolas en los cuadros de texto correspondientes.

La escena presenta la posibilidad de verse en tamaño ampliado, con el fin de manipular mejor los controles en las opciones gráficas.

190

10.7 Aproximación de la primera derivada

Suponemos conocido el valor de la función en un número de nodos equiespaciados una distancia h>0.

El objetivo es desarrollar fórmulas que permitan aproximar la derivada k-ésima de f(x) en un punto x_i, mediante una expresión de la forma:

f^(k(x_i) ≅ A₁f(x₁) + ... + A_nf(x_n)

donde los puntos x₁, ..., x_n, que se llaman nodos, debe cumplir que están equiespaciados una distancia h>0.

Obtención: Como para la obtención de las distintas aproximaciones de f^(k(x_i), se utilizará la fórmula de Taylor en el punto (x_i), se supondrá que la función cumple las condiciones de regularidad que se precisen.

En la escena interactiva de la siguiente página, comenzamos con la obtención de las fórmulas de aproximación de la primera derivada.

Debes hacer clic sobre uno de los tres números, para ver la información sobre la fórmula.

192

Haz clic para ampliar la escena

194

de acuerdo a las expresiones anteriores, el resultado numérico puede escribirse, así:

en rojo está marcado el error de redondeo y en azul el de truncamiento.

Para valores pequeños de h se tiene que:

Una acotación aproximada del error total es entonces:

El análisis realizado es teórico, ya que para datos obtenidos, a partir de una tabla, la acotación realizada no es de utilidad (si no se conoce la derivada primera, no se conocerá, tampoco, la tercera).

A continuación, con fines de comparación con la escena interactiva anterior, presentamos el algoritmo en wxMaxima para el cálculo de la derivada primera con las diferencias progresiva, regresiva y central. Posteriormente, presentamos una escena interactiva para el cálculo de la derivada primera con diferencias de orden h².

196

10.9 Derivada primera: otras diferencias de orden h²

Observa e interactúa con las fórmulas de diferencias de segundo orden, con la escena de la siguiente página.

198

10.10 Fórmulas de diferencias para la derivada segunda

Observa e interactúa con las fórmulas de diferencias para la derivada segunda, con la escena de la siguiente página.

200

Practica con las diferencias de orden h², en la escena de la siguiente página.

202

Derivada segunda con GeoGebra

En las siguientes escenas reproducimos la escena de Descartes para f(x) = sen x + 2

204

Introducción

En el área de la Tecnología, en especial la Ingeniería, las ecuaciones diferenciales son de gran importancia para modelar y explicar el comportamiento de un fenómeno dado. Igualmente, las leyes físicas, generalmente, se expresan como una razón de cambio, por ejemplo, la caída de un cuerpo, el cambio de temperatura, las variaciones del caudal de una corriente de agua, la deformación de un elemento de hormigón, etcétera.

Las ecuaciones diferenciales ordinarias (EDO) surgen como problemas de valor inicial o problemas con condiciones de frontera. En esta última parte del libro, nos dedicaremos al primer tipo de problemas, los cuales se pueden escribir con la siguiente notación:

y(a) = y₀ es una condición inicial del problema.

Se trata, entonces, de hallar una solución verificando la condición inicial dada. Esta solución es una función y(t) definida en un intervalo, que pase por el punto (a, y₀).

Los métodos que vamos a usar para resolver las EDO se suelen llamar métodos o técnicas de Runge-Kutta, entre los cuales se encuentran: el método de Euler, de Heun, métodos de Taylor y los de Runge-Kutta propiamente.

208

El método de Euler inicia con el cálculo de la primera pendiente o derivada:
φ = f(t_i, y_i+1), donde f(t_i, y_i+1) es la pendiente de la recta tangente en t_i . Esta estimación da origen a la fórmula de punto pendiente o de Euler-Cauchy, conocida como método de Euler:

En la escena interactiva de la página siguiente, presentamos el ejemplo de Mora (2013, pág. 192Mora, W. (2013). Introducción a los métodos numéricos. Obtenido de https://tecdigital.tec.ac.cr/revistamatematica/Libros/ WMora_MetodosNumericos/ WMora-ITCR-MetodosNumericos.pdf), el cual parte del problema de valor inicial:

La escena inicia con 11 puntos (n = 10) en el intervalo [1, 2], lo que genera un valor de h de 0.10. Estos valores se pueden variar para observar que a un mayor valor de n la aproximación a la solución (marcada en azul) es mejor. Esta solución es
y(t) = 1.42857t² + 4.08163t − 4.42583e^0.7t + 4.40233. Es de observar que la aproximación depende del tamaño del intervalo y del tipo de ecuación, por ejemplo:

210

Método de Euler con wxMaxima

Presentamos el algoritmo para el método de Euler, en el cual hemos incorporado la ecuación diferencial ordinaria presentada en la escena de Descartes. Puedes ampliar la imagen, y descargar y ejecutar el archivo de wxMaxima, para que observes la gráfica comparativa de y(t) aproximada con y(t) exacta.

212

11.2 Mejoras del método de Euler

El error de aproximación en el método de Euler se genera al suponer que la derivada al inicio de un intervalo es la misma a lo largo de todo el intervalo. Existen dos métodos llamados correctivos de Euler: método de Heun y el método del punto medio, en este apartado explicaremos el primero. El método de Euler que predice la pendiente al inicio del intervalo, la podemos expresar así (el superíndice cero significa que la pendiente es una predicción intermedia):

esta ecuación permite estimar la pendiente al final del intervalo:

Ahora, combinadas las dos pendientes obtenemos una pendiente promedio que al reemplazar en la fórmula de Euler (predictor) obtenemos el método propuesto por Heun (corrector):

En la siguiente escena de Descartes puedes ver los dos métodos comparados. Puedes cambiar la ecuación y la solución y(t) (si no la conoces, déjala en blanco), para ello, ingresa primero las ecuaciones y luego cambia los valores del intervalo. Usa clic izquierdo sostenido para mover la gráfica y el derecho para zoom.

214

En la escena hemos realizado el ejemplo de Mathews y Fink (2000, pág. 476Mathews, j., & Fink, K. (2000). Métodos numéricos con Matlab (Tercera ed.). Madrid: Prentice Hall.), con los siguientes datos: y' = (t-y)/2 en [0, 3] y y(0) = 1. La solución exacta para esta ecuación es y(t) = 3e^-t/2 - 2 + t. En la escena de GeoGebra puedes ingresar estos datos con
y' = (x - y)/2 y desplazando el punto A a las coordenadas (0, 1).

Para h = 0.5, estos fueron los resultados:

216

Método de Heun con wxMaxima

Este es el algoritmo para el método de Heun, en el cual hemos incorporado la ecuación diferencial ordinaria presentada en las escenas enteriores. Puedes ampliar la imagen, y descargar y ejecutar el archivo de wxMaxima, para que observes la gráfica comparativa de y(t) Heun con y(t) exacta.

218

Método de Taylor con wxMaxima

Con el método de Taylor vamos a solucionar la ecuación y′(t) = e^-2t – 4y, desde t = 0 hasta t = 1. La condición inicial es t = 0, y(0) = 1.

Inicialmente hacemos las siguientes declaraciones:

Observa que, como en los casos anteriores, hemos incluido la solución que para este problema es y(t) = 0.5e^-2t + 0.5e^-4t, esto lo hacemos con fines didácticos, para verificar el nivel de aproximación de cada método. Es obvio que recurrimos a una solución numérica cuando la solución análitica es compleja o difícil de evaluar.

A continuación calculamos las derivadas sucesivas y"(t) hasta el máximo orden que queremos evaluar que, para nuestro ejemplo, será cuarto y^iv(t)

El siguiente paso es crear un ciclo o bucle que calcule los valores de la solución, de acuerdo a la expansión de Taylor descrita anteriormente, estos valores los almacenamos en un vector y[k], donde k toma valores de 0 hasta N

220

La imagen de salida muestra que para orden cuatro (m = 4), se obtiene una muy buena solución. Haz clic sobre la imagen para observar el código completo o sobre el ícono de wxMaxima para descargarlo.

En las siguientes escenas de Descartes y GeoGebra, presentamos la solución para la misma ecuación. Para el caso de Descartes, sólo hemos incluido los métodos de Taylor de orden uno (m=1) (Euler) y dos (m=2), con el fin de comparar el nivel de aproximación. Puede observarse que para m=2, la aproximación es superior a la obtenida por el método de Heun.

222

Método de Taylor con GeoGebra

Tanto GeoGebra como wxMaxima, por ser programas de cálculo simbólicos (CAS), permiten calcular las derivadas requeridas para la expansión de Taylor, por ello, se facilita la solución a un orden superior. Es de observar que GeoGebra obliga a usar variables en notación f(x, y). Una ventaja adicional, en las últimas versiones, es la disposición de una hoja de cálculo con funcionamiento similar a una hoja de Excel, lo cual permite la generación sencilla de los puntos de la solución. Por ejemplo, la fórmula de la celda E4 es:

que corresponde a un orden cuatro del método de Taylor.

Otro comando útil es el campo de direcciones, que no hemos incluido en la escena para evitar retardamientos en la interacción, pues dicho comando consume bastantes recursos. El comando para y'(t) que en la escena hemos denotado como f0(x, y), sería:

CampoDirecciones[ f0(x, y), 10]

En la ventana gráfica de la escena sólo hemos habilitado las soluciones hasta el orden tres (m=3), pero en la hoja de cálculo lo hicimos hasta cuarto (m=4).

224

11.4 Método de Runge-Kutta

Como lo enunciamos en la introducción, los métodos anteriores se suelen llamar métodos del tipo Runge-Kutta¹ o de un sólo paso (los métodos de varios pasos no serán tratados en este libro), pues su origen surge de la utilización de un número finito de términos de la serie de Taylor, de la forma:

y_i+1 = y_i + hf(t_i, y_i) + (h²/2!)f'(t_i, y_i) + (h³/3!)f"(t_i, y_i) + ...

Tanto Euler (1745) como Runge (1885), Heun (1900) y Kutta (1901) formulan sus métodos con fundamento en las series de Taylor (1715). Sin entrar en demostraciones o mayores detalles, que sugerimos consultar en algunos de los libros reseñados en la bibliografía, es posible llegar a la siguiente fórmula:

la cual se puede expresar como:

Esta fórmula se conoce como método de Rung-Kutta de segundo orden, que coinciden con los tres primeros términos de la serie de Taylor.

226

Método de Runge-Kutta en Descartes

El propósito de este libro no fue explicar como se desarrollan las escenas de Descartes, no obstante, en esta última escena interactiva, presentamos el algoritmo que permite el cálculo del método de Runge-Kutta diseñados en el editor de Descartes.

Algoritmo Runge-Kutta

Se puede apreciar la simpleza del algoritmo, el cual incluye los vectores t[i] y yrk4[i], en los que se almacenan los puntos de la solución Runge-Kutta. Descartes permite el diseño de ciclos o bucles tipo do-while.

En la siguiente escena interactiva puedes practicar con el método.

228

Método de Runge-Kutta con GeoGebra

GeoGebra, como lo hemos dicho antes, pemite utilizar una hoja de cálculo con funciones similares a las de Excel. Esto posibilita programar, fácilmente, el algoritmo a través de las celdas de la hoja de cálculo.

Runge-Kutta y Heun en la hoja de cálculo

En la celda B3, por ejemplo, se incluye la fórmula simplificada de Runge-Kutta.

Puedes practicar en la siguiente escena interactiva de GeoGebra.

230

Método de Runge-Kutta con wxMaxima

Este es el algoritmo para el método de Runge-Kutta. No se muestran las instrucciones para la gráfica, sin embargo, puedes descargar el archivo y ejecutarlo, para observar los resultados.

232

234