Definición de método numérico

Un método numérico es un procedimiento mediante el cual se obtiene de manera aproximada,

en la mayoría de los casos, la solución de problemas realizando cálculos lógicos consecutivos
(operaciones aritméticas elementales, cálculo de funciones, consulta de una tabla de valores,
cálculo preposicional, etc.). Consiste en una lista finita de instrucciones precisas que especifican
una secuencia de operaciones (algoritmo), que producen una aproximación de la solución del
problema (solución numérica). La eficiencia en el cálculo de dicha aproximación depende, de la
implementación del algoritmo y de las características y limitaciones de los instrumentos de cálculo
(los computadores).
Aunque existen muchos tipos de métodos numéricos, éstos comparten una característica común:
invariablemente requieren de un elevado número de cálculos. Por ello con el desarrollo de
computadoras digitales eficientes y rápidas, el papel de los métodos numéricos en la solución de
problemas en ingeniería ha aumentado de forma considerable en los últimos años.
Además de proporcionar un aumento en la potencia de cálculo, la disponibilidad creciente de las
computadoras y su asociación con los métodos numéricos han influido de manera muy
significativa en el proceso de la solución actual de los problemas en ingeniería.
Errores por truncamiento
Son los que se producen al pasar del problema matemático al numérico, por ejemplo, cuando se
desprecia un término complementario, suplantando una suma infinita por una finita.

Errores por redondeo

El redondeo está definido como la aproximación de un número a una determinada cantidad de
cifras significativas, donde la última cifra significativa tendrá un valor que podrá variar
dependiendo del valor de la siguiente cifra, es decir, si la siguiente cifra tiene un valor definido
por el intervalo [0 , 4] la última cifra significativa conservará su valor, en el caso que la siguiente
cifra esté en el intervalo [5 , 9] la última cifra significativa aumentará su valor en una unidad de su
respectivo decimal.
Ejemplo de redondeo y truncamiento:

√3 ≈ 1.7320508 Raíz cuadrada de 3 truncada a 9 cifras significativas

√3 ≈ 1.7320 Raíz cuadrada de 3 truncada a 5 cifras significativas

√3 ≈ 1.7321 Raíz cuadrada de 3 redondeada a 5 cifras significativas (la quinta cifra

significativa, es decir nuestra última cifra significativa tiene a su derecha, el siguiente dígito, el
dígito 5 antes del truncamiento, por lo que la última cifra significativa toma el valor de 1 en lugar
de 0 para seguir la regla del redondeo.
Principales errores en computación.

Normalmente, desde que nos enseñan los números utilizamos base 10, también llamado
sistema decimal, en el cual se debe realizar sumas consecutivas de cada una de sus cifras,
siguiendo la regla de las unidades, decenas, centenas, etc, que matemáticamente es una
sumatoria de estos valores, como se muestra:

327,302 = 3 ∙ 102 + 2 ∙ 101 + 7 ∙ 100 + 3 ∙ 10−1 + 0 ∙ 10−2 + 2 ∙ 10−3

Da la misma forma pasa para los números decimales en base 10:

d1 d2 d3 d𝑚
x = + 2
+ 3+∙∙∙+ 𝑚
10 10 10 10

• Donde d es el valor entero de la posición del decimal

Los computadores utilizan base 2 (digital; binaria; 0/1), la cual para los números enteros trabaja
sin problemas:
(8)10 = (1000)2 = 23 + 02 + 01 + 0
Sin embargo, para expresar un número decimal en base binaria
c1 c2 c3 c𝑡
x = + 2+ 3+∙∙∙+ β
2 2 2 2 𝑡
Donde c es 1 o 0

(0,1)10 = (0,00011001100110011. . . )2
Error absoluto
𝒆𝒂 = |𝒙𝒊 − 𝒙|
𝑒𝑎 ∶ 𝑒𝑟𝑟𝑜𝑟 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑜
𝑥𝑖 ∶ 𝑉𝑎𝑙𝑜𝑟 𝑖𝑑𝑒𝑎𝑙 𝑜 𝑟𝑒𝑎𝑙
𝑥 ∶ 𝑉𝑎𝑙𝑜𝑟 𝑚𝑒𝑑𝑖𝑑𝑜 𝑜 𝑐𝑎𝑙𝑐𝑢𝑙𝑎𝑑𝑜
Aunque el valor real del error sea un valor que puede ser calculado usualmente no se utiliza en
ingeniería ya que el error se toma como un intervalo, es decir:
• Después de tomar un número determinado de medidas de un elemento que idealmente
debería medir 1 m se determina que el elemento mide en promedio 1.010 m con un error
de ± 0.006
Ya que una medida no se puede tomar como exacta, así mismo la medida del error no es exacta,
se define su existencia en un intervalo definido por:
𝑥 − 𝑒𝑎 < 𝑥𝑖 < 𝑥 + 𝑒𝑎
Error Relativo (absoluto)
𝒆𝒂 |𝒙𝒊 − 𝒙|
𝒆𝒓 = =
𝒙𝒊 𝒙𝒊
Error Relativo Porcentual
𝒆𝒑 = 𝒆𝒓 ∙ 𝟏𝟎𝟎%

Propagación del error

Para una variable 𝑧 dependiente de 𝑥 y 𝑦 el valor estará definido por la función que defina su
dependencia de estas variables y el error será calculado como se muestra
𝑧 = 𝑓(𝑥, 𝑦)
𝜕𝑓(𝑥, 𝑦) 𝜕𝑓(𝑥, 𝑦)
𝑒𝑎𝑧 = | | 𝑒𝑎𝑥 + | | 𝑒𝑎𝑦
𝜕𝑥 𝜕𝑦
• Ejemplo
Calcule el valor y el error de 𝑧 si ésta está definida por 𝑧 = 𝑥𝑦 2 donde 𝑥 = 2.0 ± 0.1 y 𝑦 =
3.0 ± 0.2

𝑧 = (2.0)(3.0)2 = 18
𝑒𝑎𝑧 = |𝑦 2 |𝑒𝑎𝑥 + |2𝑥𝑦|𝑒𝑎𝑦 = |3.02 |(0.1) + |2(2.0)(3.0)|(0.2) = 3.3
𝑧 = 18 ± 3
Número de cifras significativas validas según el error
Si ∝𝒏 es una cifra significativa valida entonces el error absoluto y relativo está definido por

𝒆𝒂 ≤ 𝟎. 𝟓 ∙ 𝟏𝟎𝒎−𝒏+𝟏
𝟎. 𝟓
𝒆𝒓 ≤
∝𝟏 ∙ 𝟏𝟎𝒏−𝟏
∝1 ∶ 𝐷í𝑔𝑖𝑡𝑜 𝑑𝑒 𝑙𝑎 𝑝𝑟𝑖𝑚𝑒𝑟 𝑐𝑖𝑓𝑟𝑎 𝑠𝑖𝑔𝑛𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑣𝑎
𝑛 ∶ 𝑁ú𝑚𝑒𝑟𝑜 𝑣𝑎𝑙𝑖𝑑𝑜 𝑑𝑒 𝑐𝑖𝑓𝑟𝑎𝑠 𝑠𝑖𝑔𝑛𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑣𝑎𝑠
𝑚 ∶ 𝐸𝑥𝑝𝑜𝑛𝑒𝑛𝑡𝑒 𝑏𝑎𝑠𝑒 10 𝑑𝑒 𝑙𝑎 𝑝𝑟𝑖𝑚𝑒𝑟 𝑐𝑖𝑓𝑟𝑎 𝑠𝑖𝑔𝑛𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑣𝑎
• Ejemplo
Definir el número de cifras significativas necesarias para no sobrepasar un error relativo de 𝒆𝒓𝒃
En este caso la desigualdad se ve invertida ya que es el error relativo definido por la ecuación
debe ser menor que el error 𝒆𝒓𝒂 razón por la cual
𝟎. 𝟓
≤ 𝒆𝒓𝒃
∝𝟏 ∙ 𝟏𝟎𝒏−𝟏
𝟎. 𝟓
𝒏 ≥ 𝐥𝐨𝐠 ( )+𝟏
∝𝟏 𝒆𝒓𝒃
Inestabilidad numérica
Sería deseable que los algoritmos numéricos proporcionaran soluciones exactas a los problemas
numéricos. Sin embargo, debido a que las computadoras solo pueden representar un sistema
numérico discreto, esto no es posible en general. El error de redondeo juega un papel importante
en este caso. La noción de estabilidad es la forma común de caracterizar lo que es posible, es
decir, de obtener la respuesta correcta, aunque no sea la respuesta exacta. En un método estable,
los errores debidos a las aproximaciones se atenúan a medida que la computación procede. En
un método inestable, cualquier error en el procesamiento se magnifica conforme el cálculo
continúa.
Serie de Taylor
Solución de ecuaciones no lineales
Ecuaciones no lineales

En matemáticas, los sistemas no lineales representan sistemas cuyo comportamiento no es

expresable como la suma de los comportamientos de sus descriptores. Más formalmente, un
sistema físico, matemático o de otro tipo es no lineal cuando las ecuaciones de movimiento,
evolución o comportamiento que regulan su comportamiento son no lineales. En particular, el
comportamiento de sistemas no lineales no está sujeto al principio de superposición, como lo es
un sistema lineal.

Idea de los métodos:

Se realizan los cálculos del método y cada nuevo intento de encontrar un intervalo que contenga
la raíz de la función, es denominado “iteración”, para la búsqueda de aproximaciones cada vez
mejores para la raíz de la función.

A cada iteración verificamos si ya se tiene una buena aproximación a través de criterios de parada.

Cuando llevamos a cabo iteraciones a través de algún método, en general es necesario

proporcionar una estimación para el error cometido.

Es lógico que al hablar de métodos de resolución de ecuaciones no lineales se trate de

hallar la solución de la ecuación, no las raíces de una función f, sin embargo, esta relación
se halla al igualar la ecuación no lineal a 0 dejando todos los términos a un lado de la
ecuación y llamando estos términos como la función f.

Criterios de parada: sea 𝑥 una aproximación para la raíz de la función y 𝑥=𝜓 su valor exacto. Se
pueden establecer dos tipos de criterio de convergencia para una precisión definida ε:

El primer criterio está relacionado con el valor de la raíz de la función, se refiere a la precisión de
la abscisa definida en el método numérico, o convergencia en x.
El segundo criterio se establece para el valor de la función, que debe ser muy pequeño cuando
se esté cerca de 𝑓(𝜓), o convergencia en y.
Ambos criterios deben ser establecidos para que se pueda registrar una buena aproximación para
el valor de la raíz de la función.

En la práctica no es conocida la raíz exacta, razón por la cual ψ es reemplazada dependiendo del
método utilizado.
Teorema de Bolzano

Sea f una función real continua en un intervalo cerrado [a,b] con f(a) y f(b) de signos contrarios.
Entonces existe al menos un punto c del intervalo abierto (a, b) con f(c) = 0.

Método de bisección
1. Dividir un intervalo inicialmente establecido exactamente en su punto medio.
2. En seguida, se verifica en cuál lado de esa división se localiza la raíz de la función. Este
pasará a ser el nuevo intervalo de búsqueda.
3. El punto medio representa una nueva aproximación para la raíz de la función.
4. Así, se determina a cada iteración un intervalo con tamaño igual a la mitad de aquel de la
iteración anterior.
5. Sea un intervalo inicial [a,b] que contenga la raíz de la función. El valor estimado por la
bisección, a cada iteración es dado por:

• Criterio de parada en x:
|𝒃 − 𝒂| < 𝜺
Método de falsa posición
Puede ser considerado una variación del Método de la Bisección.
Es utilizada la información del valor de la función calculada en los extremos del intervalo [a,b]
para estimar de manera más elaborada la nueva posición de división del intervalo.
Se asume que la función pueda ser aproximada por una recta en el intervalo [a,b].

• Criterio de parada en x: Al menos una de las siguientes debe cumplirse.

𝑓(𝑎) ∙ 𝑓(𝑎 + 𝜀) < 0

𝑓(𝑏) ∙ 𝑓(𝑏 + 𝜀) < 0
Desde la primera iteración uno de los puntos de la recta permanece constante (el punto más
alejado de la raíz.
Método del punto fijo
Un punto fijo de una función 𝑔 , es un número 𝑝 en [𝑎, 𝑏] tal que 𝑔(𝑝) = 𝑝 . El problema de
encontrar las soluciones de una ecuación 𝑓(𝑥) = 0 y el de encontrar los puntos fijos de una
función ℎ(𝑥) son equivalentes en el siguiente sentido: dado el problema de encontar las
soluciones de una ecuación 𝑓(𝑥) = 0, podemos definir una función 𝑔 con un punto fijo 𝑝 de
muchas formas; por ejemplo, 𝑓(𝑥) = 𝑥 − 𝑔(𝑥). En forma inversa, si la función 𝑔 tiene un punto
fijo en 𝑝, entonces la función definida por 𝑓(𝑥) = 𝑥 − 𝑔(𝑥) posee un cero en 𝑝.
El método de punto fijo inicia con una aproximación inicial 𝑥0 y 𝑥𝑖+1 = 𝑔(𝑥𝑖 ) genera una sucesión
de aproximaciones la cual converge a la solución de la ecuación 𝑓(𝑥) = 0. A la función 𝑔 se le
conoce como función iteradora. Se puede demostrar que dicha sucesión 〈𝑥𝑛 〉 converge siempre
y cuando |𝑔′(𝑥)| < 1 para todo 𝑥 en [𝑎, 𝑏].
• Ejemplo
𝑥3
𝑓(𝑥) = − 𝑥2 − 𝑥
10
𝑥3
𝑓(𝑥) = 𝑥 − (− + 𝑥 2 + 2𝑥)
10
3𝑥 2
𝑔′(𝑥) = − + 2𝑥 + 2
10
En el intervalo que se define [10 , 12] |𝑔′(𝑥)| > 1 razón por la cual esa función g no nos sirve y
es necesario calcular una nueva.
𝑥2
𝑥 ( − 𝑥 − 1) = 0
10

𝑥 = ±√10(𝑥 + 1) = 𝑔(𝑥)
10
𝑔′(𝑥) =
2√10(𝑥 + 1)
Esta nueva función si cumple la condición |𝑔′(𝑥)| < 1 por ello, para el cálculo de la raíz se
utiliza esta función.

𝑥𝑖+1 = √10(𝑥𝑖 + 1)
Método Newton-Raphson
Este método utiliza una aproximación a la raíz de la función por medio del corte en x de la recta
tangente al punto que se está evaluando.

𝑓(𝑥𝑖 )
𝑥𝑖+1 = 𝑥𝑖 −
𝑓′(𝑥𝑖 )

• Para obtener una convergencia es necesario que el valor inicial sea cercano a la raíz
buscada.
• En la práctica el valor inicial para aplicar el método de Newton-Raphson puede
obtenerse tras ejecutar algunas iteraciones del método de bisección.
• Si no se aporta un valor lo suficientemente cercano o si la función tiene cambios de
concavidad cuando se acerca a la raíz, puede correr el riesgo de divergir.
 𝑚𝑔 (−𝑏𝑡) 𝑚𝑔
𝑣(𝑡) = (𝑣0 − )𝑒 𝑚 +
𝑏 𝑏
2. Una partícula de masa 𝑚 = 5 𝑘𝑔 es lanzada con una inclinación 𝜃 = 45° respecto al suelo
con una velocidad inicial |𝑣0 | este vector velocidad tiene una componente en 𝑥 y 𝑦 de las
𝑚
cuales se conoce 𝑣0𝑥 = 40 𝑠 (Velocidad inicial en dirección x). Durante su vuelo, la
partícula presenta una fuerza de fricción con el aire de valor ⃗⃗⃗
𝐹𝑟 = −𝑏𝑣 (Vectorialmente
esto indica que la fuerza de fricción siempre es opuesta a la velocidad por lo que este
vector no es constante ya que el vector velocidad varía de forma tangente a la trayectoria),
donde 𝑏 representa la constante de fricción de la partícula con el aire. (Suponga la
dirección en y positiva hacia arriba)
a. En el instante 𝑡 = 3 𝑠 la partícula se encuentra en una posición 𝑥 = 100 𝑚 .
Determine la magnitud y dirección del vector velocidad en ese instante.
(Se dejan las siguientes ecuaciones a las cuales se debe llegar en uno de los pasos durante la
resolución del ejercicio.)
𝑏𝑡
𝑣𝑥 (𝑡) = 𝑣0𝑥 𝑒 (− 𝑚 )
𝑚𝑔 (−𝑏𝑡) 𝑚𝑔
𝑣𝑦 (𝑡) = (𝑣0𝑦 + )𝑒 𝑚 −
𝑏 𝑏
Capítulo III. Sistemas de Ecuaciones Lineales
En este módulo trataremos el tema “Sistema de Ecuaciones Lineales”
Existen muchos problemas en ingeniería que culminan en sistemas lineales. Por ejemplo:

• Circuitos eléctricos
• Estática de estructuras (MEF)
• Conducción de calor (MDF)
• Problemas no lineales de modo general, son resueltos con linealizaciones sucesivas (por
ejemplo, a través del Método de Newton-Raphson)

Veremos dos enfoques para resolver sistemas lineales

• Forma directa
• Forma iterativa

Un sistema de ecuaciones es clasificado como lineal, si tiene la forma:

Un ejemplo de sistema mecánico muy simple es:

En el cual f1 y f2 son fuerzas constantes. Si se determina el sistema de ecuaciones que describe

el sistema mecánico, se tiene:

Se puede reorganizar el sistema separando x1 y x2:

Y si lo convertimos a su forma matricial se obtiene:

• La matriz que multiplica el vector de las variables de ese sistema lineal es denominada
“Matriz de Rigidez”.
• Esa matriz tiene una importancia muy grande principalmente para análisis estáticos.
• Toda la información de rigidez del sistema está contenida en esa matriz, incluyendo la
organización y acoplamiento entre elementos estructurales y su forma de conexión.
• Una característica importante de esta matriz es que en todo caso es simétrica.
Estos sistemas pueden ser resueltos de forma exacta de múltiples formas como por ejemplo por
Regla de Cramer o inversión de la matriz de rigidez. En este curso estudiaremos las siguientes
formas.

• Forma Directa: A partir del sistema en la forma Ax = b, reescribirlo a través de un proceso

de escalonamiento de la matriz A. Este proceso desea dejar la matriz A escrita en la forma
triangular, o sea, con un triángulo de elementos nulos. Principal método: Método de la
Eliminación Gaussiana.

• Forma Iterativa: Se despejan las incógnitas en todas las ecuaciones del sistema lineal,
una a una, despejando una incógnita diferente para cada ecuación. Es dada una
estimación inicial para la solución del sistema, y se calculan las incógnitas con base en
esa primera aproximación. Después, las incógnitas encontradas forman una segunda
aproximación para la solución del sistema. El proceso es repetido hasta que haya
convergencia.

Antes de iniciar la resolución de sistemas lineales es importante recordar la clasificación de los

mismos en relación al número de soluciones:
• Sistema compatible determinado (una única solución)

• Sistema compatible indeterminado (infinitas soluciones)

• Sistema incompatible (no tiene solución)

En este curso nos enfocaremos en los sistemas compatibles determinados.

Sistemas Compatibles Indeterminados pueden surgir en algunos tipos de problemas de

autovalores, como por ejemplo análisis modal y análisis de pandeo lineal.

Cuando el sistema sea compatible indeterminado, su solución no será única, y sí una familia de
posibles soluciones.

En los sistemas compatibles determinados la solución será única.

En algunas situaciones podemos resolver sistemas indeterminados en equilibrio utilizando un

método iterativo.

Métodos Exactos
Método de eliminación Gaussiana

Consiste en realizar el escalonamiento de una matriz (triangulación) para, posteriormente, realizar

la solución de las incógnitas del sistema linear (sustitución hacia atrás).

Idea central del método: sistemas lineales son equivalentes si presentan la misma solución. Por
tanto, existen algunas manipulaciones posibles que pueden ser realizadas entre las ecuaciones
del sistema, que no afectan su solución. A saber:
• Intercambiar dos ecuaciones de posición, es decir, cambiar el orden de las ecuaciones.
• Multiplicar una ecuación por un escalar no nulo.
• Adicionar un múltiplo de una ecuación, a otra.

Descripción del método

1. Se escribe la matriz cómo su forma extendida
 2. Se definirá en este método una iteración completa todos los cálculos necesarios para
hacer una columna, por debajo de la diagonal principal, 0. Para contar esta iteración se
define un contador 𝑘 el cual a su vez indica la columna que se está haciendo 0.

𝑘 = 1, … , 𝑛 − 1
3. Para cada iteración se calculan los factores 𝑚 que multiplicarán la fila 𝑘 que
posteriormente se sustraerán a cada fila.

𝑎𝑖𝑘 (𝑘)
𝑚𝑖𝑘 = 𝑖 = 𝑘 + 1, 𝑘 + 2, … , 𝑛
𝑎𝑘𝑘 (𝑘)
𝑎𝑖𝑗 (𝑘+1) = 𝑎𝑖𝑗 (𝑘) − 𝑚𝑖𝑘 𝑎𝑘𝑗 (𝑘) 𝑗 = 1, … , 𝑛
(𝑘)
𝑓𝑖 (𝑘+1) = 𝑓𝑖 (𝑘) − 𝑚𝑖𝑘 𝑓𝑘
4. Al completar la triangulación se realiza la sustitución hacia atrás para determinar la
solución o soluciones del sistema.

• Ejemplo
6 −2 2 4 𝑥1 12 6 −2 2 4 𝑥1 12
12 −8 6 10 𝑥2 34 0 −4 2 2 𝑥2 10
( ) (𝑥 ) = ( ) → ( ) (𝑥 ) = ( )
3 −13 9 3 3 27 0 −12 8 1 3 21
−6 4 1 −18 𝑥4 −38 0 2 3 −14 𝑥4 −26
6 −2 2 4 𝑥1 12 6 −2 2 4 𝑥1 12
0 −4 2 2 𝑥 2 10 0 −4 2 2 𝑥 2 10
( ) (𝑥 ) = ( ) → ( ) (𝑥 ) = ( )
0 0 2 −5 3 −9 0 0 2 −5 3 −9
0 0 4 −13 𝑥4 −21 0 0 0 −3 𝑥4 −3
La anterior implementación del método no utiliza la permutación de filas, o cambio del orden de
ecuaciones. Sin embargo, es necesario aclarar que su uso se lleva a cabo usualmente cuando
uno de los valores de la diagonal se hace 0 en algún momento del método, por lo que se hace
necesaria la modificación del sistema, intercambiando filas de la matriz principal así como las filas
correspondientes del vector columna.
Factorización PA=LU
Esta factorización es una técnica matemática para expresar una matriz como facturización de un
conjunto de matrices. Algunas implementación de códigos la utilizan para determinar la inversa
de la matriz 𝐴 valiéndose de la inversa de sus matrices factores los cuales sus inversas no
requieren un número de cálculos tan alto en comparación al número de cálculos necesario para
definir la matriz 𝐴.

𝑷𝑨 = 𝑳𝑼
Donde 𝐴 es la matriz inicial, 𝑈 es la matriz triangular superior obtenida por medio de la eliminación
Gaussiana, 𝑃 es una matriz de permutación que tiene en cuenta la reorganización de las filas que
se haya realizado durante el proceso de triangulación y 𝐿 es la matriz triangular inferior que se
forma con los factores 𝑚 obtenidos en el mismo procedimiento, esta matriz también tiene en
cuenta las permutaciones realizadas, se mostrará un ejemplo de ello posteriormente. Una matriz
𝐿 obtenida de una eliminación Gaussiana sin permutaciones se observa a continuación.
1 0 ⋯ 0 0
𝑚21 1 ⋯ 0 0
𝑳= ⋮ ⋮ ⋱ ⋮ ⋮
𝑚(𝑛−1)1 𝑚(𝑛−1)2 ⋯ 1 0
( 𝑚𝑛1 𝑚𝑛2 ⋯ 𝑚𝑛(𝑛−1) 1)
Posteriormente se aplica la función inversa a ambos lados de la ecuación matricial de arriba.

(𝑷𝑨)−𝟏 = (𝑳𝑼)−𝟏
𝑨−𝟏 𝑷−𝟏 = 𝑼−𝟏 𝑳−𝟏
Y se postmultiplica por 𝑃 a ambos lados

𝑨−𝟏 = 𝑼−𝟏 𝑳−𝟏 𝑷

• Ejemplo 1
Tomando como ejemplo la matriz triangulada por eliminación Gaussiana
1 0 0 0
2 1 0 0 6 −2 2 4
1 0 −4 2 2
𝐿= 3 1 0 𝑈=( )
2 0 0 2 −5
1 0 0 0 −3
(−1 − 2 2 1)

1 0 0 0
1 0 0 0 6 −2 2 4 2 1 0 0 6 −2 2 4
0 1 0 0 12 −8 6 10 1 0 −4 2 2
( )( )= 3 1 0 ( )
0 0 1 0 3 −13 9 3 2 0 0 2 −5
0 0 0 1 −6 4 1 −18 1 0 0 0 −3
(−1 −
2
2 1)

En este caso la matriz P es igual a la matriz identidad ya que no se realizó ningún cambio de filas
 • Ejemplo 2
Realizar la factorización 𝑃𝐴 = 𝐿𝑈 para la matriz que se tiene a continuación.
0 1 1
𝐴 = (−1 2 −4)
2 −5 1

En este caso se puede que ver que el elemento 𝑎11 = 0 lo cual para el método de eliminación
Gaussiana supone un problema a la hora de halla los factores 𝑚 por lo cual se realizará una
permutación de la matriz realizando un cambio de filas. Adicionalmente realizaremos otros
cambios en las filas para observar el cambio en la matriz L y la matriz de permutación P.
0 1 1 2 −5 1
(−1 2 −4) 𝐢𝐧𝐭𝐞𝐫𝐜𝐚𝐦𝐛𝐢𝐚𝐧𝐝𝐨 𝐟𝐢𝐥𝐚 𝟏 𝐲 𝟑 (−1 2 −4)
2 −5 1 0 1 1
1
Se obtiene m21 = − y se realiza la eliminación del nuevo término a21
2
2 −5 1 1 0 0
(0 −1/2 −7/2) El estado de la matriz L, tendríamos L = (−1/2 1 0)
0 1 1 0 0 1
2 −5 1 2 −5 1
(0 −1/2 −7/2) 𝐢𝐧𝐭𝐞𝐫𝐜𝐚𝐦𝐛𝐢𝐚𝐧𝐝𝐨 𝐟𝐢𝐥𝐚 𝟐 𝐲 𝟑 (0 1 1 )
0 1 1 0 −1/2 −7/2
1 0 0
Ahora L = ( 0 1 0) a causa del cambio entre las filas 2 y 3
−1/2 0 1
1
Se obtiene m32 = − y se realiza la eliminación del nuevo término a32
2
2 −5 1 1 0 0
U = (0 1 1) y L=( 0 1 0)
0 0 −3 −1/2 −1/2 1
En este caso 𝑃 ≠ 𝐼, para determinarlo se realizan las mismas permutaciones realizadas sobre
la matriz identidad.

De esta forma 𝑃𝐴 = 𝐿𝑈 queda de la siguiente forma.

Método de Gauss-Jordan

Este método complementa la eliminación Gaussiana con el objetivo de evitar la necesidad de

utilizar la sustitución hacia atrás, por medio de la diagonalización de la matriz. Después de realizar
la eliminación Gaussiana se continúa hasta dejar la matriz cómo se muestra.

𝑎11 0 ⋯ 0 0 | 𝑓1
0 𝑎22 ⋯ 0 0 | 𝑓2
⋮ ⋮ ⋱ ⋮ ⋮ | ⋮
0 0 ⋯ 𝑎(𝑛−1)(𝑛−1) 0 | 𝑓𝑛−1
( 0 0 ⋯ 0 𝑎𝑛𝑛 | 𝑓𝑛 )

Esto se logra utilizando como “pivote” los elementos de la diagonal comenzando con
𝒂𝒏𝒏 hasta 𝒂𝟐𝟐 .

𝑘 = 𝑛, 𝑛 − 1, … ,2
𝑎𝑖𝑘 (𝑘)
𝑚𝑖𝑘 = 𝑖 = 1,2, … , 𝑘 − 1
𝑎𝑘𝑘 (𝑘)
Ya que todos los elementos por debajo de la diagonal son 0 los únicos elementos de la matriz
extendida que van a ver involucrados son los términos 𝑓. Ya que la resta entre ecuaciones solo
está involucrando la variable de la ecuación a hacer 0 ya que todos los demás elementos de la
fila de abajo son cero, a excepción de la diagonal por lo que la diagonal mantendrá el valor de
sus términos y todos los términos por encima de la diagonal se harán 0 uno por uno.
(𝑘−1)
𝑓𝑖 (𝑘) = 𝑓𝑖 (𝑘−1) − 𝑚𝑖𝑘 𝑓𝑘 𝑖 = 1,2, … , 𝑘 − 1
• Ejemplo
Métodos aproximados
Método de Jacobi

𝑛
1
𝑥𝑖 (𝑘+1) = 𝑏 − ∑ 𝑎𝑖𝑗 𝑥𝑗 (𝑘)
𝑎𝑖𝑖 𝑖
𝑗=1
𝑗≠𝑖
( )

1
𝑥1 (𝑘+1) = (𝑏1 − 𝑎12 𝑥2 𝑘 − 𝑎13 𝑥3 𝑘 − ∙∙∙ −𝑎1𝑛 𝑥𝑛 𝑘 )
𝑎11
1
𝑥2 (𝑘+1) = (𝑏 − 𝑎21 𝑥1 𝑘 − 𝑎23 𝑥3 𝑘 − ∙∙∙ −𝑎2𝑛 𝑥𝑛 𝑘 )
𝑎22 2
⋮
1
𝑥𝑛 (𝑘+1) = (𝑏𝑛 − 𝑎𝑛1 𝑥1 𝑘 − 𝑎𝑛2 𝑥2 𝑘 − ∙∙∙ −𝑎𝑛𝑛−1 𝑥𝑛−1 𝑘 )
𝑎𝑛𝑛
Método Gauss-Seidel
A diferencia del método de Jacobi, no es necesario calcular todas las 𝑥𝑖 para calcular las
siguientes 𝑥𝑖 . En este caso se pueden calcular una por una lo que disminuye la memoria de
máquina necesaria para calcular 𝑥𝑖 .
𝑖−1 𝑛
1
𝑥𝑖 (𝑘+1) = (𝑏𝑖 − ∑ 𝑎𝑖𝑗 𝑥𝑗 (𝑘+1) − ∑ 𝑎𝑖𝑗 𝑥𝑗 (𝑘) )
𝑎𝑖𝑖
𝑗=1 𝑗=𝑖+1

1
𝑥1 (𝑘+1) = (𝑏 − 𝑎12 𝑥2 (𝑘) − 𝑎13 𝑥3 (𝑘) − ⋯ − 𝑎1𝑛 𝑥𝑛 (𝑘) )
𝑎11 1
1
𝑥2 (𝑘+1) = (𝑏 − 𝑎21 𝑥1 (𝑘+1) − 𝑎23 𝑥3 (𝑘) − ⋯ − 𝑎2𝑛 𝑥𝑛 (𝑘) )
𝑎22 2
1
𝑥3 (𝑘+1) = (𝑏3 − 𝑎31 𝑥1 (𝑘+1) − 𝑎32 𝑥2 (𝑘+1) − 𝑎34 𝑥4 (𝑘) − ⋯ − 𝑎3𝑛 𝑥𝑛 (𝑘) )
𝑎33
⋮
1
𝑥𝑛−1 (𝑘+1) = (𝑏 − 𝑎(𝑛−1)1 𝑥1 (𝑘+1) − ⋯ − 𝑎(𝑛−1)(𝑛−2) 𝑥𝑛−2 (𝑘+1)
𝑎(𝑛−1)(𝑛−1) 𝑛−1
− 𝑎(𝑛−1)𝑛 𝑥𝑛 (𝑘) )
1
𝑥𝑛 (𝑘+1) = (𝑏 − 𝑎𝑛1 𝑥1 (𝑘+1) − 𝑎𝑛2 𝑥2 (𝑘+1) − ⋯ − 𝑎𝑛𝑛−1 𝑥𝑛−1 (𝑘+1) )
𝑎𝑛𝑛 𝑛

Convergencia del método

La convergencia de este método depende de la dominancia de la diagonal de la matriz, es decir,
se asegura la convergencia si la matriz se considera diagonalmente dominante, de lo contrario,
el método no converge o tomará un número mayor de iteraciones para alcanzar la convergencia.
Una matriz es diagonalmente dominante, si en cada una de las filas, el valor absoluto del
elemento de la diagonal principal es mayor que la suma de los valores absolutos de los elementos
restantes de la fila. En ocasiones la matriz de un sistema de ecuaciones no es diagonalmente
dominante pero cuando se modifica el orden de las ecuaciones el nuevo sistema puede tener
matriz de coeficientes diagonalmente dominante.
𝑛

|𝑎𝑘𝑘 | > ∑|𝑎𝑘𝑖 |

𝑖=1
𝑖≠𝑘
Ejercicio Gauss Seidel

∑𝐹1 = 𝑓1 − 𝑘1 𝑥1 + 𝑘2 (𝑥2 − 𝑥1 ) = 0
∑𝐹2 = 𝑓2 − 𝑘2 (𝑥2 − 𝑥1 ) + 𝑘3 (𝑥3 − 𝑥2 ) = 0
∑𝐹3 = 𝑓3 − 𝑘3 (𝑥3 − 𝑥2 ) + 𝑘4 (𝑥4 − 𝑥3 ) = 0
∑𝐹4 = 𝑓3 − 𝑘4 (𝑥4 − 𝑥3 ) + 𝑘5 (𝑥5 − 𝑥4 ) = 0
∑𝐹5 = 𝑓5 − 𝑘5 (𝑥5 − 𝑥4 ) = 0

(𝑘1 + 𝑘2 )𝑥1 − 𝑘2 𝑥2 = 𝑓1
−𝑘2 𝑥1 + (𝑘2 + 𝑘3 )𝑥2 − 𝑘3 𝑥3 = 𝑓2
−𝑘3 𝑥2 + (𝑘3 + 𝑘4 )𝑥3 − 𝑘4 𝑥4 = 𝑓3
−𝑘4 𝑥3 + (𝑘4 + 𝑘5 )𝑥4 − 𝑘5 𝑥5 = 𝑓4
−𝑘5 𝑥4 + 𝑘5 𝑥5 = 𝑓5

(𝑘1 + 𝑘2 ) −𝑘2 0 0 0 𝑥1 𝑓1
−𝑘2 (𝑘2 + 𝑘3 ) −𝑘3 0 0 𝑥2 𝑓2
0 −𝑘3 (𝑘3 + 𝑘4 ) −𝑘4 0 𝑥3 = 𝑓3
0 0 −𝑘4 (𝑘4 + 𝑘5 ) −𝑘5 𝑥4 𝑓4
( 0 0 0 −𝑘5 𝑥
𝑘5 ) ( 5 ) (𝑓5 )

[𝐾 ][𝑋] = [𝐹 ]
Matriz de rigidez multiplicada por el vector de posición es igual al vector de fuerzas.
1. Darle desarrollo al sistema por medio del método Gauss-Seidel para unos
𝑁
valores: 𝑘1 = 𝑘2 = 𝑘3 = 𝑘4 = 𝑘5 = 100 𝑚 𝑓1 = 10𝑁 𝑓2 = 𝑓3 = 20𝑁 𝑓4 =
25𝑁 𝑓5 = 30𝑁
Método de Newton para sistemas de Ecuaciones no lineales
El método de Newton para la resolución de ecuaciones no lineales puede ser extrapolado para
sistemas de ecuaciones no lineales. Dado un sistema de dos ecuaciones, se desea hallar el punto
que hace que dichas ecuaciones sean iguales a cero.

Aplicando a este caso el método de Newton, se puede generar un sistema matricial de la siguiente
forma.

Para comenzar, hay que definir los valores para los puntos iniciales y calcular 𝐽−1 que
corresponde a la matriz inversa del jacobiano del sistema de ecuaciones, el cual se calcula de la
siguiente forma.

𝜕𝑓1 𝜕𝑓1
𝜕𝑥 𝜕𝑦
𝑱=
𝜕𝑓2 𝜕𝑓2
( 𝜕𝑥 𝜕𝑦 )

A partir de entonces, se calcula el valor siguiente de la aproximación y se repite iterativamente

este procedimiento. Hasta encontrar el vector 𝑿𝑲 tal que el vector 𝑭(𝑿𝑲 ) ≈ 𝟎 .

𝑿𝑲+𝟏 = 𝑿𝑲 − 𝑱−𝟏 |𝑿𝑲 𝑭(𝑿𝑲 )

Descripción del método
1. Se determina el jacobiano de F determinando las derivadas parciales necesarias
2. Se define un vector inicial 𝑋0
3. Se reemplaza en el jacobiano a partir de 𝑋0
4. Se determina la inversa del jacobiano (esta inversa puede ser calculada cuando el
jacobiano aún está en su forma de ecuaciones, sin embargo, la matriz quedaría escrita en
términos de los productos y divisiones entre las derivadas parciales, razón por la cual hace
los términos muy extensos).
5. Se determina 𝐹(𝑋0 ) y se realizan las operaciones descritas en la ecuación de iteración del
método y se continúa con las iteraciones.

• Ejemplo
Dado el siguiente sistema de ecuaciones, hallar por el método de Newton para sistemas de
ecuaciones un punto solución:

𝑥 2 − 10𝑥 = −8 − 𝑦 2
𝑥𝑦 2 + 𝑥 = −8 + 10𝑦
 𝑓1 (𝑥, 𝑦) = 𝑥 2 − 10𝑥 + 𝑦 2 + 8 = 0
𝑓2 (𝑥, 𝑦) = 𝑥𝑦 2 + 𝑥 − 10𝑦 + 8 = 0
𝜕𝑓1 𝜕𝑓1
𝜕𝑥 𝜕𝑦 2𝑥 − 10 2𝑦
𝑱= =( 2 )
𝜕𝑓2 𝜕𝑓2 𝑦 + 1 2𝑥𝑦 − 10
( 𝜕𝑥 𝜕𝑦 )

Renombrando como 𝑥1 𝑦 𝑥2
2𝑥1 − 10 2𝑥2
𝐽=( )
𝑥2 2 + 1 2𝑥1 𝑥2 − 10
−1
(𝑘+1) (𝑘) (𝑘) (𝑘) (𝑘) 2 (𝑘) (𝑘) 2
𝑥 𝑥 2𝑥1 − 10 2𝑥2 𝑥1 − 10𝑥1 + 𝑥2
+8
( 1(𝑘+1) ) = ( 1(𝑘) ) − ( (𝑘) 2 (𝑘) (𝑘)
) ( 2 )
𝑥2 𝑥2 (𝑘) (𝑘) (𝑘) (𝑘)
𝑥2 + 1 2𝑥1 𝑥2 − 10 𝑥1 𝑥2 + 𝑥1 − 10𝑥2 + 8
Por medio de una modificación a la ecuación se puede reescribir para evitar la necesidad de
calcular la inversa de la matriz Jacobiana.

𝑌 = 𝐽−1 𝐹
𝐽𝑌 = 𝐹
Esta nueva ecuación toma la forma de un sistema lineal, ya que al reemplazar los valores k-
esimos de 𝑥 los resultados serán números reales dando como resultado un sistema como se
observa a continuación.
(𝑘) (𝑘) (𝑘) (𝑘) 2 (𝑘) (𝑘) 2
2𝑥1 − 10 2𝑥2 𝑦1 𝑥1 − 10𝑥1 + 𝑥2
+8
( (𝑘) 2 (𝑘) (𝑘)
)( (𝑘)
) =( 2 )
𝑦2 (𝑘) (𝑘) (𝑘) (𝑘)
𝑥2 +1 2𝑥1 𝑥2 − 10 𝑥1 𝑥2 + 𝑥1 − 10𝑥2 + 8
Este sistema puede solucionarse por medio de otro método numérico y la solución se
reemplazará en el método de Newton.

𝑋 (𝑘+1) = 𝑋 (𝑘) − 𝑌 (𝑘)

Ejercicio Newton para sistemas de ecuaciones

3
Donde 𝐹 = 𝑘 √𝑢 siendo u la diferencia entre la longitud original del resorte y la longitud del
resorte.
 ∑𝐹1 = −𝑘1 3√𝑥1 + 𝑘2 3√𝑥2 − 𝑥1 + 𝑓1 = 0

∑𝐹2 = −𝑘2 3√𝑥2 − 𝑥1 + 𝑓2 = 0

Se multiplican ambas ecuaciones por -1 y se definen cómo las funciones 𝑔 que harán parte del
método de Newton

𝑔1 (𝑥1 , 𝑥2 ) = 𝑘1 3√𝑥1 − 𝑘2 3√𝑥2 − 𝑥1 − 𝑓1

𝑔2 (𝑥1 , 𝑥2 ) = 𝑘2 3√𝑥2 − 𝑥1 − 𝑓2
𝑘1 𝑘2 −𝑘2
3
+ 3 3
3 ∙ √(𝑥1 )2 3 ∙ √(𝑥2 − 𝑥1 )2 3 ∙ √(𝑥2 − 𝑥1 )2
𝐽= −𝑘2 𝑘2
3 3
( 3 ∙ √(𝑥2 − 𝑥1 )2 3 ∙ √(𝑥2 − 𝑥1 )2 )
3 3
𝑘1 √𝑥1 − 𝑘2 √𝑥2 − 𝑥1 − 𝑓1
𝐹=( 3
)
𝑘2 √𝑥2 − 𝑥1 − 𝑓2

1. Dar desarrollo al sistema no lineal por medio del método de Newton, utilizando el
método Gauss-Seidel para determinar 𝑌 (𝑘) para el método de Newton sin
necesidad de determinar la inversa de la matriz Jacobiana de cada iteración del
método de Newton.

𝑋 (𝑘+1) = 𝑋 (𝑘) − 𝑌 (𝑘)

𝑌 = 𝐽−1 𝐹
𝐽𝑌 = 𝐹
𝑁
Para unos datos 𝑘1 = 𝑘2 = 10 𝑚 𝑓1 = 20𝑁 𝑓2 = 30𝑁
Interpolación
“En el subcampo matemático del análisis numérico, se denomina interpolación a la obtención de
nuevos puntos partiendo del conocimiento de un conjunto discreto de puntos.”

Con frecuencia se encontrará con que tiene que estimar valores intermedios entre datos definidos
por puntos. El método más común que se usa para este propósito es la interpolación polinomial.
Recuerde que la fórmula general para un polinomio de n-ésimo grado es:

𝑓(𝑥) = 𝑎0 + 𝑎1 𝑥 + 𝑎2 𝑥 2 + · · · + 𝑎𝑛 𝑥 𝑛
Dados n + 1 puntos, hay uno y sólo un polinomio de grado n que pasa a través de todos los
puntos. Por ejemplo, hay sólo una línea recta (es decir, un polinomio de primer grado) que une
dos puntos. De manera similar, únicamente una parábola une un conjunto de tres puntos. La
interpolación polinomial consiste en determinar el polinomio único de n-ésimo grado que se ajuste
a n + 1 puntos. Este polinomio, entonces, proporciona una fórmula para calcular valores
intermedios. Aunque hay uno y sólo un polinomio de n-ésimo grado que se ajusta a n + 1 puntos,
existe una gran variedad de formas matemáticas en las cuales puede expresarse este polinomio.

El método de los elementos finitos, entre otros métodos, usan polinomios para realizar la
interpolación de funciones.
En el método de elementos finitos, los polinomios se utilizan como base para la aproximación de
soluciones de ecuaciones diferenciales parciales.
Por lo tanto, es muy importante realizar un estudio de la interpolación con polinomios, porque
esta teoría es necesaria, por ejemplo, para la formulación de elementos de celosía, viga o sólido,
utilizados en aplicaciones prácticas de elementos finitos.

Interpolación polinomial de Newton

El polinomio de interpolación de Newton en diferencias divididas es una de las formas más
populares y útiles. Antes de presentar la ecuación general, estudiaremos las versiones de primero
y segundo grados por su sencilla interpretación visual.
Interpolación lineal
La forma más simple de interpolación consiste en unir dos puntos con una línea recta. Dicha
técnica, llamada interpolación lineal, se ilustra de manera gráfica en la figura utilizando triángulos
semejantes.
Cuanto menor sea el intervalo entre los datos, mejor será la aproximación.
Interpolación cuadrática
Es común encontrar altos errores resultantes de aproximar una curva mediante una línea recta.
En consecuencia, una estrategia para mejorar la estimación consiste en introducir alguna
curvatura a la línea que une los puntos. Si se tienen tres puntos como datos, éstos pueden
ajustarse en un polinomio de segundo grado (también conocido como polinomio cuadrático o
parábola).

Forma general de los polinomios de interpolación de Newton

𝑓𝑛 : 𝑃𝑜𝑙𝑖𝑛𝑜𝑚𝑖𝑜 𝑑𝑒 𝑖𝑛𝑡𝑒𝑟𝑝𝑜𝑙𝑎𝑐𝑖ó𝑛 𝑑𝑒 𝑔𝑟𝑎𝑑𝑜 𝑛

Donde las evaluaciones de la función colocadas entre corchetes son diferencias divididas finitas.
Por ejemplo, la primera diferencia dividida finita en forma general se representa como:

Errores de la interpolación polinomial de Newton

Obsérvese que se requiere calcular la diferencia dividida n+1 por lo que para un polinomio de
grado n el cual se ha calculado con n+1 puntos necesita de n+2 puntos conocidos para calcular
el error aproximado para cualquier punto x.
Interpolación polinomial de Lagrange
El polinomio de interpolación de Lagrange es simplemente una reformulación del polinomio de
Newton que evita el cálculo de las diferencias divididas, y se representa de manera concisa como:
Interpolación Polinómica
Es necesario aclarar que el polinomio de grado 𝑛 que pasa por 𝑛 + 1 puntos El polinomio
interpolador de grado 𝑛 determinado a partir de 𝑛 + 1 puntos (𝑥𝑖 , 𝑓(𝑥𝑖 )) pasa por todos esos
puntos. Lo anterior se ve de esta forma.
𝑓(𝑥0 ) = 𝑎0 𝑥0 0 + 𝑎1 𝑥01 + 𝑎2 𝑥0 2 + · · · + 𝑎𝑛 𝑥0 𝑛
𝑓(𝑥1 ) = 𝑎0 𝑥1 0 + 𝑎1 𝑥11 + 𝑎2 𝑥1 2 + · · · + 𝑎𝑛 𝑥1 𝑛
𝑓(𝑥2 ) = 𝑎0 𝑥2 0 + 𝑎1 𝑥21 + 𝑎2 𝑥2 2 + · · · + 𝑎𝑛 𝑥2 𝑛
⋮
𝑓(𝑥𝑛 ) = 𝑎0 𝑥𝑛 0 + 𝑎1 𝑥𝑛1 + 𝑎2 𝑥𝑛 2 + · · · + 𝑎𝑛 𝑥𝑛 𝑛
Si lo anterior se cumple se tiene el siguiente sistema de ecuaciones lineales, teniendo como
incógnitas a los coeficientes 𝑎𝑖 del polinomio.
1 𝑥0 𝑥0 2 ⋯ 𝑥0 𝑛−1 𝑥0 𝑛 𝑎0 𝑓(𝑥0 )
1 𝑥1 𝑥1 2 ⋯ 𝑥1 𝑛−1 𝑥1 𝑛 𝑎 1 𝑓(𝑥1 )
⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ = ⋮
1 𝑥𝑛−1 𝑥𝑛−1 2 ⋯ 𝑥𝑛−1 𝑛−1 𝑥𝑛−1 𝑛 𝑎 𝑛−1 𝑓(𝑥𝑛−1 )
(1 𝑥𝑛 𝑥𝑛 2 ⋯ 𝑥𝑛 𝑛−1 𝑥𝑛 𝑛 ) ( 𝑎𝑛 ) ( 𝑓(𝑥𝑛 ) )
𝑋𝐴 = 𝐹
𝐴 = 𝑋 −1 𝐹
Interpolación trigonométrica
Las anteriores interpolaciones tienen la particularidad de que son únicamente validas en un
intervalo determinado o incluso en las cercanías externas a los límites del intervalo, pero cuando
se aleja del intervalo de la interpolación, naturalmente los valores se alejan de la función original
o de la tendencia de los datos, ya que no se posee información externa al intervalo para el cálculo
del polinomio.
Hay muchos casos en ingeniería en los que se trabaja con funciones que son periódicas, o con
datos que tienen una tendencia periódica, como es común en máquinas rotatorias, esta
característica periódica hace que la aproximación de su comportamiento por medio de un
polinomio no sea tan apreciada. Por ello se han creado métodos para determinar una
aproximación a estas funciones periódicas, el más básico de ellos es la interpolación
trigonométrica, la cual se vale de series trigonométricas para expresar un “polinomio”
trigonométrico interpolador que pueda expresar la tendencia periódica de una función.
𝑛

𝑝𝑛 (𝑥) = 𝑎0 + ∑(𝑎𝑘 𝑐𝑜𝑠(𝑘𝑥) + 𝑏𝑘 𝑠𝑖𝑛(𝑘𝑥))

𝑘=1
Esta ecuación tiene una particularidad ya que el periodo de la función es 2𝜋, sin embargo es claro
que un conjunto de datos periódicos seguramente no tendrán este mismo periodo, por lo cual,
hay una forma de generalizarla para una función de periodo 𝑇 de la siguiente forma.
𝑛
2𝜋 2𝜋
𝑝𝑛 (𝑥) = 𝑎0 + ∑ (𝑎𝑘 𝑐𝑜𝑠 ( 𝑘𝑥) + 𝑏𝑘 𝑠𝑖𝑛 ( 𝑘𝑥))
𝑇 𝑇
𝑘=1
Ejercicios
1. Un sensor de determinó la posición de una partícula p en diferentes tiempos
como se observa en la tabla.
t x
0,21 2,5
0,38 6,7
0,79 9,1
1,01 14,5
1,23 17,2
Determine la posición de la partícula en el instante 𝑡 = 0,5 𝑠.
a. Por medio del polinomio interpolador de mayor grado posible a partir de los
datos.
b. Por medio de la serie trigonométrica de mayor grado posible a partir de los
datos.
2 1
2. Una partícula se mueve a una velocidad 𝑣(𝑡) = 𝑒 −𝑡 + 3
𝑡 +1
a. Determine la distancia recorrida por la partícula entre el instante 𝑡 = 0 y 𝑡 =
1 𝑠.
Consejo: Utilizar un polinomio de interpolación de grado n, a placer, para
expresar la función de la velocidad.
Diferenciación Numérica
Diferencias finitas
El método consiste en una aproximación de las derivadas parciales por expresiones algebraicas
con los valores de la variable dependiente en un número finito de puntos seleccionados en el
dominio.
- Diferencias finitas hacia adelante

𝑑𝐹 𝐹(𝑥 + ℎ) − 𝐹(𝑥)
≈
𝑑𝑥 ℎ
𝜕𝐹 𝐹(𝑥 + ℎ, 𝑦) − 𝐹(𝑥, 𝑦)
≈
𝜕𝑥 ℎ
𝜕𝐹 𝐹(𝑥, 𝑦 + ℎ) − 𝐹(𝑥, 𝑦)
≈
𝜕𝑦 ℎ
- Diferencias finitas hacia atrás

𝑑𝐹 𝐹(𝑥) − 𝐹(𝑥 − ℎ)
≈
𝑑𝑥 ℎ
𝜕𝐹 𝐹(𝑥, 𝑦) − 𝐹(𝑥 − ℎ, 𝑦)
≈
𝜕𝑥 ℎ
𝜕𝐹 𝐹(𝑥, 𝑦) − 𝐹(𝑥, 𝑦 − ℎ)
≈
𝜕𝑦 ℎ
- Diferencias finitas centrales

𝑑𝐹 𝐹(𝑥 + ℎ) − 𝐹(𝑥 − ℎ)
≈
𝑑𝑥 2ℎ
𝜕𝐹 𝐹(𝑥 + ℎ, 𝑦) − 𝐹(𝑥 − ℎ, 𝑦)
≈
𝜕𝑥 2ℎ
𝜕𝐹 𝐹(𝑥, 𝑦 + ℎ) − 𝐹(𝑥, 𝑦 − ℎ)
≈
𝜕𝑦 2ℎ
- Diferencias finitas, segundo orden

𝑑 2 𝐹 𝐹(𝑥 + ℎ) − 2𝐹(𝑥) + 𝐹(𝑥 − ℎ)

≈
𝑑𝑥 2 ℎ2
Integración numérica
Regla del trapecio

La regla del trapecio es un método de integración numérica, es decir, un método para calcular
aproximadamente el valor de la integral definida. El método consiste en determina el área
delimitada por la función 𝑓(𝑥) por medio del área del trapecio delimitado por la recta que pasa
por los puntos (𝑎, 𝑓(𝑎)) 𝑦 (𝑏, 𝑓(𝑏)) .
𝑏
𝑓(𝑎) + 𝑓(𝑏)
∫ 𝑓(𝑥) ≈ (𝑏 − 𝑎)
𝑎 2
Regla del trapecio compuesta

Es una modificación a la regla del trapecio que consiste en dividir el intervalo (𝑎, 𝑏) en 𝑛
subintervalos iguales a los cuales se les aplicará la regla del trapecio, de esta forma se aproximará
de manera más acertada la integral de 𝑓(𝑥).
 𝑏
𝑓 (𝑎 ) 𝑓 (𝑏 )
∫ 𝑓 (𝑥)𝑑𝑥 ≈ ℎ ( + 𝑓 (𝑎 + ℎ) + 𝑓 (𝑎 + 2ℎ) + ⋯ + 𝑓 (𝑎 + (𝑛 − 1)ℎ) + )
𝑎 2 2
Donde
𝑏−𝑎
ℎ=
𝑛
𝑏 𝑛−1
𝑓(𝑎) + 𝑓(𝑏)
∫ 𝑓(𝑥) ≈ ℎ ( + ∑ 𝑓(𝑎 + 𝑖ℎ))
𝑎 2
𝑖=1

Regla de Simpson (1/3)

El método consiste en determina el área delimitada por la función 𝑓(𝑥) por medio del área
delimitada por la parábola que pasa por los puntos (𝑎, 𝑓(𝑎)), (𝑚, 𝑓(𝑚)), (𝑏, 𝑓(𝑏)) al dividir el
intervalo (𝑎, 𝑏) en dos sub-intervalos iguales.

𝑏−𝑎
𝑚=
2
𝑏
𝑏−𝑎
∫ 𝑓(𝑥)𝑑𝑥≈ (𝑓(𝑎) + 4𝑓(𝑚) + 𝑓(𝑏))
𝑎 6

Regla de Simpson (1/3) compuesta

Así como la regla del trapecio, la regla de Simpson puede ser utilizada para aproximar de manera
más acertada la integral por medio de la partición del intervalo (𝑎, 𝑏), en 𝑛 subintervalos de esta
forma se tienen (𝑛 + 1) funciones cuadráticas para interpolar la función y definir con estas una
mejor aproximación de la integral. Por esto 𝑛 debe ser un número par.
𝑏
ℎ ℎ ℎ
∫ 𝑓(𝑥)𝑑𝑥≈ (𝑦0 + 4𝑦1 + 𝑦2 ) + (𝑦2 + 4𝑦3 + 𝑦4 )+. . . . + (𝑦𝑛−2 + 4𝑦𝑛−1 + 𝑦𝑛 )
𝑎 3 3 3

𝑦𝑖 = 𝑓(𝑥𝑖 )
𝑥𝑖 = 𝑎 + 𝑖ℎ ; 0 < 𝑖 < 𝑛
 𝑏−𝑎
ℎ =
𝑛
𝑛 𝑛
−1
𝑏 2 2
ℎ
∫ 𝑓(𝑥)𝑑𝑥 ≈ (𝑦 + 4 ∑(𝑦2𝑖−1 ) + 2 ∑ (𝑦2𝑖 ) +𝑦𝑛 )
𝑎 3 0
𝑖=1 𝑖=1
Cuadratura de Gauss
- Cuadratura:

La cuadratura en matemáticas indica el problema de la determinación del área de formas

geométricas por medio de la determinación del cuadrado de área igual a la figura
geométrica en cuestión.

- Cuadratura numérica

Se trata de la determinación de la integral de una función 𝑓(𝑥) por medio de una sumatoria
en la que es necesaria la determinación de los 𝑥𝑖 , o puntos del intervalo, y los 𝑤𝑖 , también
llamados pesos, más adecuados, para calcular la integral. Este problema se vuelve
entonces la determinación de los 𝑤𝑖 y los 𝑥𝑖 .
Se ha determinado que para un intervalo determinado estos 𝑤𝑖 y 𝑥𝑖 son constantes sin
importar las funciones concernientes si se ha de tener en cuenta que estas funciones
pueden ser expresadas con un polinomio interpolador.
Por ello es válido afirmar que si 𝑓(𝑥) es un polinomio hasta de grado 2𝑛 − 1 la integral
puede ser determinada de forma exacta, siendo el error igual a 0.
1 𝑛

∫ 𝑓(𝑥)𝑑𝑥 = ∑ 𝑤𝑖 𝑓(𝑥𝑖 )
−1 𝑖=1
De la anterior ecuación hay un total de 2𝑛 incognitas. Que a continuación se ejemplificará
la determinación de ellas por medio de funciones polinómicas, para un 𝑛 = 2.
1
∫ 𝑓(𝑥)𝑑𝑥 = 𝑤1 𝑓(𝑥1 ) + 𝑤2 𝑓(𝑥2 )
−1

1
𝑓(𝑥) = 1 ; ∫ 𝑑𝑥 = 2 = 𝑤1 + 𝑤2
−1
1
𝑓(𝑥) = 𝑥 ; ∫ 𝑥𝑑𝑥 = 0 = 𝑤1 𝑥1 + 𝑤2 𝑥2
−1
1
2
𝑓(𝑥) = 𝑥 2 ; ∫ 𝑥 2 𝑑𝑥 = = 𝑤1 𝑥1 2 + 𝑤2 𝑥2 2
−1 3
 1
3
𝑓(𝑥) = 𝑥 ; ∫ 𝑥 3 𝑑𝑥 = 0 = 𝑤1 𝑥1 3 + 𝑤2 𝑥2 3
−1

𝑤1 + 𝑤2 = 2
𝑤1 𝑥1 + 𝑤2 𝑥2 = 0
2
𝑤1 𝑥1 2 + 𝑤2 𝑥2 2 =
3
{𝑤1 𝑥1 3 + 𝑤2 𝑥2 3 = 0

Por lo que resulta un sistema de no lineal de dimensiones (2𝑛)𝑥(2𝑛), razón por la cual se
hace compleja la determinación de la integral si es necesario primeramente la
determinación de estas incognitas. Sin embargo, su utilidad se ve reflejada en el hecho
de que estos valores son constantes sin importar la función, a pesar de que solo es exacto
para funciones polinómicas las integrales de otras funciones pueden ser aproximadas por
medio de este método, utilizando los valores como constantes.
1 1
Para este caso 𝑤1 = 𝑤2 = 1 y 𝑥1 = − y 𝑥2 =
√3 √3

En la tabla anterior se observan los valores de estas constantes hasta el 𝑛 = 5.

La dificultad que supone la determinación de estas constantes al ser necesaria la
resolución del sistema no lineal dio paso a la relación existente en estas ecuaciones con
otras en matemáticas, de la cual nació la relación con los polinomios de Legendre, donde:
1
𝑤𝑖 =
(1 − 𝑥 2 )[𝑃𝑛′ (𝑥𝑖 )]2
Donde 𝑃𝑛′ es la primera derivada el 𝑛 − 𝑒𝑠𝑖𝑚𝑜 polinomio de Legendre. Por esta relación,
al método se le ha dado el nombre de Cuadratura de Gauss-Legendre.
- Polinomios de Legendre

En el análisis de ecuaciones diferenciales ordinarias, las funciones de Legendre son las

soluciones de las ecuaciones diferenciales de Legendre:

𝑑 𝑑
[(1 − 𝑥 2 ) 𝑃𝑛 (𝑥)] + 𝑛(𝑛 + 1) 𝑃𝑛 (𝑥) = 0
𝑑𝑥 𝑑𝑥
Donde los polinomios de Legendre son definidos por:

1 𝑑𝑛
𝑃𝑛 (𝑥) = 𝑛 𝑛
[(𝑥 2 − 1)𝑛 ]
2 𝑛! 𝑑𝑥

- Traslación de la cuadratura de Gauss

Se puede observar que en la ecuación para la determinación de la integral y las

constantes determinadas solo serían validas en el caso de que se tuviera el mismo
intervalo de integración (−1,1). Para generalizar el problema Hay una forma de definir una
ecuación a partir de la anterior, conservando los valores de las constantes. Evaluando la
𝑏−𝑎 𝑏+𝑎
función en 2 𝑥𝑖 + 2 .
𝑏 1
𝑏−𝑎 𝑏+𝑎 𝑏−𝑎
∫ 𝑓(𝑥)𝑑𝑥 = ∫ 𝑓 ( 𝑡+ ) 𝑑𝑡
𝑎 −1 2 2 2

𝒃 𝒏
𝒃−𝒂 𝒃−𝒂 𝒃+𝒂
∫ 𝒇(𝒙)𝒅𝒙 = ∑ 𝒘𝒊 𝒇 ( 𝒙𝒊 + )
𝒂 𝟐 𝟐 𝟐
𝒊=𝟏
Integrales múltiples
𝑏 𝑑(𝑥)

∫ ∫ 𝑓(𝑥, 𝑦)𝑑𝑦𝑑𝑥
𝑎 𝑐(𝑥)

• Por Método de Simpson 1/3

Utilizando la regla de Simpson desde la integral externa hacia la interna se puede determinar una
aproximación a esta integral
𝑏 𝑑(𝑥) 𝑑(𝑥0 ) 𝑑(𝑥1 ) 𝑑(𝑥2 )
𝑥2 − 𝑥0
∫ ∫ 𝑓(𝑥, 𝑦)𝑑𝑦𝑑𝑥≈ [ ∫ 𝑓(𝑥0 , 𝑦)𝑑𝑦 + 4 ( ∫ 𝑓(𝑥1 , 𝑦)𝑑𝑦) + ∫ 𝑓(𝑥2 , 𝑦)𝑑𝑦]
6
𝑎 𝑐(𝑥) 𝑐(𝑥0 ) 𝑐(𝑥1 ) 𝑐(𝑥2 )

𝑏−𝑎
Donde 𝑥0 = 𝑎 , 𝑥1 = y 𝑥2 = 𝑏
2
𝑑(𝑥0 )
𝑦20 − 𝑦00
∫ 𝑓(𝑥0 , 𝑦)𝑑𝑦 ≈ [𝑓(𝑥0 , 𝑦00 ) + 4𝑓(𝑥0 , 𝑦10 ) + 𝑓(𝑥0 , 𝑦20 )]
6
𝑐(𝑥0 )

𝑑(𝑥0 )−𝑐(𝑥0 )
Donde 𝑦00 = 𝑐(𝑥0 ) , 𝑦10 = y 𝑦20 = 𝑑(𝑥0 )
2
𝑑(𝑥1 )
𝑦21 − 𝑦01
∫ 𝑓(𝑥0 , 𝑦)𝑑𝑦 ≈ [𝑓(𝑥1 , 𝑦01 ) + 4𝑓(𝑥0 , 𝑦11 ) + 𝑓(𝑥0 , 𝑦21 )]
6
𝑐(𝑥1 )

𝑑(𝑥1 )−𝑐(𝑥1 )
Donde 𝑦01 = 𝑐(𝑥1 ) , 𝑦11 = y 𝑦21 = 𝑑(𝑥1 )
2
𝑑(𝑥2 )
𝑦22 − 𝑦02
∫ 𝑓(𝑥2 , 𝑦)𝑑𝑦 ≈ [𝑓(𝑥0 , 𝑦02 ) + 4𝑓(𝑥0 , 𝑦12 ) + 𝑓(𝑥0 , 𝑦22 )]
6
𝑐(𝑥2 )

𝑑(𝑥2 )−𝑐(𝑥2 )
Donde 𝑦02 = 𝑐(𝑥2 ) , 𝑦12 = y 𝑦22 = 𝑑(𝑥2 )
2

• Por cuadratura de Gauss

𝑏 𝑑(𝑥) 𝑛 𝑚
𝑏−𝑎
∫ ∫ 𝑓(𝑥, 𝑦)𝑑𝑦𝑑𝑥 ≈ ∑ ∑ 𝑤𝑖 𝑤𝑗 [𝑑(𝑃𝑥𝑖 ) − 𝑐(𝑃𝑥𝑖 )]𝑓(𝑃𝑥𝑖 , 𝑃𝑦𝑗 )
4
𝑎 𝑐(𝑥) 𝑖=1 𝑗=1

𝑏−𝑎 𝑏+𝑎
𝑃𝑥𝑖 = 𝑥𝑖 +
2 2
𝑑(𝑃𝑥𝑖 ) − 𝑐(𝑃𝑥𝑖 ) 𝑑(𝑥𝑖 ) + 𝑐(𝑥𝑖 )
𝑃𝑦𝑗 = 𝑥𝑗 +
2 2
Ejercicios
1. Si el coeficiente de fricción cinética entre el embalaje de 200 𝑘𝑔 y el suelo es 𝜇𝑐 =
0.25, determine la rapidez del embalaje cuando 𝑡 = 4 𝑠. El embalaje comienza a
𝑡2
(− )
moverse desde el punto de reposo y lo remolca la fuerza de 𝐹𝑒 16 , con 𝐹 =
10𝑘𝑁. 𝑡2
(− )
𝐹𝑒 16

Sumatoria de fuerzas en y para determinar la normal:

𝑡2
(− )
𝛴𝐹𝑦 = 0 = 𝑁 + sin(30°) 𝐹𝑒 16 − 𝑚𝑔
𝑡2
(− )
𝑁 = 𝑚𝑔 − sin(30°) 𝐹𝑒 16

Sumatoria de fuerzas en x para determinar las fuerzas involucradas

𝑡2
(− )
𝛴𝐹𝑥 = cos(30°) 𝐹𝑒 16 𝑑𝑡 − 𝑁𝜇𝑐
𝑡2 𝑡2
(− ) (− )
𝛴𝐹𝑥 = cos(30°) 𝐹𝑒 16 𝑑𝑡 − (𝑚𝑔 − sin(30°) 𝐹𝑒 16 ) 𝜇
𝑐

Principio de conservación de cantidad de movimiento

4 2
4 2
√3 𝑡 1 (− 𝑡 )
(− )
𝐹 ∫ 𝑒 16 𝑑𝑡 − 𝜇𝑐 ∫ (𝑚𝑔 − 𝐹𝑒 16 ) 𝑑𝑡 + 𝑚𝑣1 = 𝑚𝑣2
2 2
0 0
4 2
𝐹 𝑡
(− )
(√3 + 𝜇𝑐 ) ∫ 𝑒 16 𝑑𝑡 − 4𝜇𝑐 𝑚𝑔 + 𝑚𝑣1 = 𝑚𝑣2
2
0

2. Cuál es la temperatura promedio de una placa cuadrada de acero de lado 2𝜋 𝑚

de la cual se conoce la distribución de temperatura bidimensional en un instante
2
sin(𝑥) −𝑦
de tiempo. 𝑇(𝑥, 𝑦) = 300 ( 𝑒 10 + 1) °𝐶
𝑥+1

2𝜋 2𝜋
1 1 sin(𝑥) −𝑦 2
𝑇𝑃𝑟𝑜𝑚 = ∬ 𝑇(𝑥, 𝑦)𝑑𝐴 = ∫ ∫ 300 ( 𝑒 10 + 1) 𝑑𝑦𝑑𝑥
𝐴 (2𝜋)2 𝑥+1
0 0
3. Integre la función siguiente tanto en forma analítica como numérica. Emplee las
reglas del trapecio y de Simpson 1/3 con 𝑛 = 4 y cuadratura de Gauss-Legendre
con 𝑛 = 3 para integrar numéricamente la función. Calcule los errores relativos
porcentuales para los resultados numéricos.
𝜋/2

∫ cos(𝑥) 𝑒 𝑥 𝑑𝑥
0

ln(𝑠+1)
4. Una grúa levanta una viga de 2.50 𝑀𝑔 con una fuerza 𝐹 = (30 + 3 ) 𝑘𝑁.
𝑠3 +2
Determine la velocidad de la viga cuando alcanza 𝑠 = 3 𝑚. También, ¿cuánto
tiempo se requiere para que alcance esta altura a partir del punto de reposo?
Nota : La energía cinética inicial del sistema, además del trabajo realizado por
todas las fuerzas externas e internas que actúan en el sistema, es igual a la
energía cinética final del sistema.
𝐸𝑐1 + 𝑊1→2 = 𝐸𝑐2
𝑠
ln(𝑠 + 3) 1
0 + 1000 ∫ (30 + 3 3
) 𝑑𝑠 − 𝑚𝑔𝑠 = 𝑚𝑣22
𝑠 +2 2
0
Ecuaciones diferenciales ordinarias
(Ordinal Differential Equations) ODE

En todas las áreas de ingeniería las ecuaciones diferenciales son una herramienta matemática
que permite modelar el comportamiento de diversos fenómenos bien sea físicos, químicos,
económicos entre otros.
Una categoría particular de estas ecuaciones diferenciales son las ecuaciones diferenciales
ordinarias, que relacionan una función desconocida de una variable independiente con sus
derivadas. Es necesario recalcar el hecho de que solo hay variable independiente a diferencia de
las ecuaciones diferenciales parciales, que involucran derivadas parciales de varias variables y
una o más de sus derivadas respecto de tal variable.
De forma general una ecuación diferencial ordinaria de grado n se escribe de la forma

𝑑𝑛 𝑦 𝑑𝑦 𝑑 2 𝑦 𝑑 𝑛−1 𝑦
= 𝑓 (𝑥, 𝑦, , 2 , ⋯ , 𝑛−1 )
𝑑𝑥 𝑛 𝑑𝑥 𝑑𝑥 𝑑𝑥
Donde es claro que la variable independiente es x.
En general, las ecuaciones diferenciales de bajo orden, con coeficientes constantes y/o donde
las funciones involucradas respecto a 𝑥 pueden ser fácilmente identificadas, separadas entre si
y la variable dependiente pueden tener una solución analítica, sin embargo son casos particulares,
cuando el caso es diferente la solución de estas ecuaciones se deja exclusivamente a los
métodos numéricos.
Inicialmente analizaremos las ecuaciones diferenciales ordinarias de primer orden
𝑑𝑦
= 𝑓(𝑥, 𝑦)
𝑑𝑥
Método de Euler Explícito
Este método utiliza la aproximación de la función por medio de su derivada 𝑓(𝑥𝑘 , 𝑦𝑘 ) y un “paso”
ℎ utilizando la siguiente ecuación.

𝑦𝑘+1 = 𝑦𝑘 + 𝑓(𝑥𝑘 , 𝑦𝑘 )ℎ
Es necesario el conocimiento del valor inicial o de frontera 𝑦0 = 𝑓(𝑥0 , 𝑦0 ) y la definición del paso
ℎ el cual a medida que se hace más pequeño la aproximación de la función 𝑦 se hace más precisa.

Ejercicios
1. Dar solución a la ecuación diferencial teniendo como valor de frontera 𝑦(0) = 4.
Determinar el error relativo porcentual a partir de la solución exacta de la ecuación.
𝑑𝑦
= 4𝑥 3 − 9𝑥 2 − 6𝑥 − 20
𝑑𝑥
𝑦(𝑥) = 𝑥 4 − 3𝑥 3 − 3𝑥 2 − 20𝑥 + 4
 𝑦𝑘+1 = 𝑦𝑘 + (4𝑥𝑘3 − 9𝑥𝑘2 − 6𝑥𝑘 − 20)ℎ
2. Dar solución a la ecuación diferencial teniendo como valor de frontera 𝑦(0) = 1.
Determinar el error relativo porcentual a partir de la solución exacta de la ecuación.
𝑑𝑦
+ 𝑦 = 2𝑥 − 10
𝑑𝑥
𝑦(𝑥) = 13𝑒 −𝑥 + 2𝑥 − 12
𝑦𝑘+1 = 𝑦𝑘 + (2𝑥𝑘 − 𝑦𝑘 − 10)ℎ
Método de Euler Implícito
A diferencia del método explícito, el método implícito utiliza la información del siguiente punto
para determinarlo en sí mismo, es decir, no se vale únicamente de la información previa o “externa”
sino que utiliza información aún desconocida o “interna” para dar solución, razón por la cual el
método toma una forma implícita ya que la aproximación no es determinada de forma directa.

𝑦𝑘+1 = 𝑦𝑘 + 𝑓(𝑥𝑘+1 , 𝑦𝑘+1 )ℎ

Ya que se desconoce 𝑦𝑘+1 la determinación de este k-ésimo término debe ser determinado por
medio de la solución de la ecuación (seguramente no lineal) que resulta del método, es decir:

𝑦𝑘+1 − 𝑦𝑘 − 𝑓(𝑥𝑘+1 , 𝑦𝑘+1 )ℎ = 0

En el caso en el cual la función 𝑓 no dependa de 𝑦 no será necesaria la resolución de la ecuación,
ya que no será muy diferente del método explícito.

Ejercicios
1. Dar solución a la ecuación diferencial teniendo como valor de frontera 𝑦(0) = 4.
Determinar el error relativo porcentual a partir de la solución exacta de la ecuación
utilizando tanto el método implícito como el explícito.
𝑑𝑦
= 4𝑥 3 − 9𝑥 2 − 6𝑥 − 20
𝑑𝑥
𝑦(𝑥) = 𝑥 4 − 3𝑥 3 − 3𝑥 2 − 20𝑥 + 4
Explicito:
𝑦𝑘+1 = 𝑦𝑘 + (4𝑥𝑘3 − 9𝑥𝑘2 − 6𝑥𝑘 − 20)ℎ
Implícito
3 2
𝑦𝑘+1 = 𝑦𝑘 + (4𝑥𝑘+1 − 9𝑥𝑘+1 − 6𝑥𝑘+1 − 20)ℎ
2. Dar solución a la ecuación diferencial teniendo como valor de frontera 𝑦(0) = 1.
Determinar el error relativo porcentual a partir de la solución exacta de la ecuación.
𝑑𝑦
+ 𝑦 = 2𝑥 − 10
𝑑𝑥
 𝑦(𝑥) = 13𝑒 −𝑥 + 2𝑥 − 12
Explícito
𝑦𝑘+1 = 𝑦𝑘 + (2𝑥𝑘 − 𝑦𝑘 − 10)ℎ
Implícito
𝑦𝑘 + (2𝑥𝑘+1 − 10)ℎ
𝑦𝑘+1 =
1+ℎ
Tarea
Una barra de metal, cuya temperatura inicial es de 25°𝐶, se mete en un horno eléctrico que se
1
encuentra a 300°𝐶. 𝑘 = −0.05
𝑠

Calcular la temperatura de la barra en 𝑡 = 12 𝑠

𝑑𝑇
= 𝑘(𝑇 – 𝑇𝛼)
𝑑𝑡
Sabiendo que
𝑇𝛼 = 300°
𝑇(0) = 25°
Métodos de Runge-Kutta
Para un mismo valor de paso ℎ, los métodos de Runge-Kutta (RK) logran una mayor exactitud en
comparación a los métodos anteriores. Existen muchas variantes, pero todas tienen la forma
generalizada de la ecuación.

𝑦𝑖+1 = 𝑦𝑖 + 𝜙(𝑥𝑖 , 𝑦𝑖 , ℎ)ℎ

Donde 𝜙(𝑥𝑖 , 𝑦𝑖 , ℎ) se conoce como función incremento, la cual puede interpretarse como una
pendiente representativa en el intervalo. La función incremento se escribe en forma general
como
𝜙 = 𝑎1 𝑘1 + 𝑎2 𝑘2 + · · · + 𝑎𝑛 𝑘𝑛
Donde las a son constantes y las k son
𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 )
𝑘2 = 𝑓(𝑥𝑖 + 𝑝1 ℎ , 𝑦𝑖 + 𝑞11 𝑘1 ℎ)
𝑘3 = 𝑓(𝑥𝑖 + 𝑝2 ℎ , 𝑦𝑖 + 𝑞21 𝑘1 ℎ + 𝑞22 𝑘2 ℎ)
⋮
𝑘𝑛 = 𝑓(𝑥𝑖 + 𝑝𝑛–1 ℎ , 𝑦𝑖 + 𝑞(𝑛–1)1 𝑘1 ℎ + 𝑞(𝑛–1)2 𝑘2 ℎ + ··· + 𝑞(𝑛–1)(𝑛−1) 𝑘𝑛–1 ℎ)
Donde las 𝑝 y las 𝑞 son constantes. Observe que las k son relaciones de recurrencia. Es decir,
𝑘1 aparece en la ecuación 𝑘2 , la cual aparece en la ecuación 𝑘3 , etcétera. Como cada k es una
evaluación funcional, esta recurrencia vuelve eficientes a los métodos RK para cálculos en
computadora.

Las constantes 𝑎, 𝑝 y 𝑞 dependen del grado 𝑛 del método de Runge-Kutta que se desea realizar
y estas están estrechamente relacionadas con los coeficientes que se encuentran en la Serie
de Taylor en dos variables de orden 𝑛 de la función 𝑓. Estas constantes deben cumplir que:
𝑛 𝑛

𝑝𝑘 = ∑ 𝑞𝑘𝑗 𝑦 ∑ 𝑎𝑗 = 1
𝑗=1 𝑗=1

Los métodos de Runge-Kutta de alto orden no serán ejemplificados a causa de la necesidad de

calcular las constantes 𝑎, 𝑝, 𝑞. Las constantes están determinadas para los métodos de bajo orden
y en ocasiones son modificadas al determinar los matemáticos que ciertos valores mejoran la
aproximación de las ecuaciones diferenciales. Los grados más utilizados computacionalmente
son del grado 1 al 5.

Runge-Kutta grado 1 (Método de Euler)

Observe que el método de Runge-Kutta (RK) de primer orden con n = 1 es, de hecho, el método
de Euler.
𝑦𝑖+1 = 𝑦𝑖 + 𝑘1 ℎ
𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 )

Runge-Kutta grado 2 (Heun con un corrector)

1 1
𝑦𝑖+1 = 𝑦𝑖 + ( 𝑘1 + 𝑘2 ) ℎ
2 2
𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 )
 𝑘2 = 𝑓(𝑥𝑖 + ℎ , 𝑦𝑖 + 𝑘1 ℎ)

Runge-Kutta grado 2 (Ralston)

1 2
𝑦𝑖+1 = 𝑦𝑖 + ( 𝑘1 + 𝑘2 ) ℎ
3 3
𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 )
3 3
𝑘2 = 𝑓 (𝑥𝑖 + ℎ , 𝑦𝑖 + 𝑘1 ℎ)
4 4

Runge-Kutta grado 3
1
𝑦𝑖+1 = 𝑦𝑖 + (𝑘1 + 4𝑘2 + 𝑘3 )ℎ
6
𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 )
1 1
𝑘2 = 𝑓 (𝑥𝑖 + ℎ , 𝑦𝑖 + 𝑘1 ℎ)
2 2
𝑘3 = 𝑓(𝑥𝑖 + ℎ , 𝑦𝑖 − 𝑘1 ℎ + 2𝑘2 ℎ)

Runge-Kutta grado 4
1
𝑦𝑖+1 = 𝑦𝑖 + (𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4 )ℎ
6
𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 )
1 1
𝑘2 = 𝑓 (𝑥𝑖 + ℎ , 𝑦𝑖 + 𝑘1 ℎ)
2 2
1 1
𝑘3 = 𝑓 (𝑥𝑖 + ℎ , 𝑦𝑖 + 𝑘2 ℎ)
2 2
𝑘4 = 𝑓(𝑥𝑖 + ℎ , 𝑦𝑖 + 𝑘3 ℎ)

Runge-Kutta grado 5 (Método de Butcher)

1
𝑦𝑖+1 = 𝑦𝑖 + (7𝑘1 + 32𝑘3 + 12𝑘4 + 32𝑘5 + 7𝑘6 )ℎ
90
𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 )
1 1
𝑘2 = 𝑓 (𝑥𝑖 + ℎ , 𝑦𝑖 + 𝑘1 ℎ)
4 4
1 1 1
𝑘3 = 𝑓 (𝑥𝑖 + ℎ , 𝑦𝑖 + 𝑘1 ℎ + 𝑘2 ℎ)
4 8 8
1 1
𝑘4 = 𝑓 (𝑥𝑖 + ℎ , 𝑦𝑖 − 𝑘2 ℎ + 𝑘3 ℎ)
2 2
3 3 9
𝑘5 = 𝑓 (𝑥𝑖 + ℎ , 𝑦𝑖 + 𝑘1 ℎ + 𝑘4 ℎ)
4 16 16
3 2 12 12 8
𝑘6 = 𝑓 (𝑥𝑖 + ℎ , 𝑦𝑖 − 𝑘1 ℎ + 𝑘2 ℎ + 𝑘3 ℎ − 𝑘4 ℎ + 𝑘5 ℎ)
7 7 7 7 7
Sistemas de ecuaciones diferenciales ordinarias
Los sistemas de ecuaciones diferenciales ordinarias son sistemas de ecuaciones en los que todas
las variables del sistema dependen de una misma variable independiente.

𝑑𝑦1
= 𝑓1 (𝑥, 𝑦1 , 𝑦2 , ⋯ , 𝑦𝑛 )
𝑑𝑥
𝑑𝑦2
= 𝑓2 (𝑥, 𝑦1 , 𝑦2 , ⋯ , 𝑦𝑛 )
𝑑𝑥
⋮
𝑑𝑦𝑛−1
= 𝑓𝑛−1 (𝑥, 𝑦1 , 𝑦2 , ⋯ , 𝑦𝑛 )
𝑑𝑥
𝑑𝑦𝑛
{ 𝑑𝑥
= 𝑓𝑛 (𝑥, 𝑦1 , 𝑦2 , ⋯ , 𝑦𝑛 )

𝑦(1𝑘+1) = 𝑦(1𝑘) + 𝑓1 (𝑥, 𝑦(1𝑘) , 𝑦(2𝑘) , ⋯ , 𝑦𝑛(𝑘) ) ℎ

𝑦(2𝑘+1) = 𝑦(2𝑘) + 𝑓2 (𝑥, 𝑦(1𝑘) , 𝑦(2𝑘) , ⋯ , 𝑦𝑛(𝑘) ) ℎ
⋮
𝑦𝑛−1 = 𝑦𝑛−1 + 𝑓𝑛−1 (𝑥, 𝑦(1𝑘) , 𝑦(2𝑘) , ⋯ , 𝑦(𝑛𝑘) ) ℎ
(𝑘+1) (𝑘)

{ 𝑦(𝑛𝑘+1) = 𝑦𝑛(𝑘) + 𝑓𝑛 (𝑥, 𝑦1(𝑘) , 𝑦(2𝑘) , ⋯ , 𝑦(𝑛𝑘) ) ℎ

Para generalizar se puede escribir de forma vectorial.

𝒀𝑘+1 = 𝒀𝑘 + 𝑭(𝑥, 𝒀𝑘 )ℎ
Donde 𝒀𝑘 es el vector columna que define el valor de cada variable en el paso 𝑘 y 𝑭 es la
función vectorial, como vector columna que contiene todas las funciones 𝑓𝑖 .
Ecuaciones diferenciales ordinarias de orden superior
Una ecuación diferencial de orden 𝑛 está definida por la siguiente ecuación.

𝑑𝑛 𝑦 𝑑𝑦 𝑑 2 𝑦 𝑑 𝑛−1 𝑦
= 𝑓 (𝑥, 𝑦, , 2 , ⋯ , 𝑛−1 )
𝑑𝑥 𝑛 𝑑𝑥 𝑑𝑥 𝑑𝑥
Con el objetivo de utilizar los métodos ya estudiados, es necesario reescribir la ecuación
utilizando múltiples cambios de variables, comenzando con un cambio de variable trivial.

𝑦 = 𝑦1
Posteriormente este cambio de variable se deriva y se nombra una nueva variable a partir de esta
derivada y el proceso se repite hasta llegar a la derivada de grado 𝑛. El último cambio de variable,
llamado 𝑦𝑛 , no tendrá un cambio de variable para su derivada, por el contrario llevará la función
𝑑𝑦 𝑑 2 𝑦 𝑑 𝑛−1 𝑦
𝑓 (𝑥, 𝑦, 𝑑𝑥 , 𝑑𝑥 2 , ⋯ , 𝑑𝑥 𝑛−1 ) ahora evaluada en los cambios de variable, reemplazando todos los
términos de derivadas contenidos en la función.
𝑑𝑦1
= 𝑦2
𝑑𝑥
𝑑 2 𝑦1 𝑑𝑦2
= = 𝑦3
𝑑𝑥 2 𝑑𝑥
⋮
𝑑 𝑛 𝑦1 𝑑 𝑛−1 𝑦2 𝑑 2 𝑦𝑛−1 𝑑𝑦𝑛
= =⋯= = = 𝑓(𝑥, 𝑦1 , 𝑦2 , 𝑦3 , ⋯ , 𝑦𝑛 )
𝑑𝑥 𝑛 𝑑𝑥 𝑛−1 𝑑𝑥 2 𝑑𝑥
Después de ello las igualdades se organizan en forma de sistema de ecuaciones ascendente.
𝑑𝑦1
= 𝑦2
𝑑𝑥
𝑑𝑦2
= 𝑦3
𝑑𝑥
⋮
𝑑𝑦𝑛−1
= 𝑦𝑛
𝑑𝑥
𝑑𝑦𝑛
{ 𝑑𝑥
= 𝑓(𝑥, 𝑦1 , 𝑦2 , 𝑦3 , ⋯ , 𝑦𝑛 )
Lo anterior genera un sistema de ecuaciones diferenciales de primer orden al cual se le da
solución por medio de los métodos ya explicados.

𝑦(𝑘+1)
1
= 𝑦(𝑘)
1
+ 𝑦(𝑘)
2
ℎ
𝑦(𝑘+1)
2
= 𝑦(𝑘)
2
+ 𝑦(𝑘)
3
ℎ
⋮
(𝑘+1)
𝑦𝑛−1 = 𝑦(𝑘)
𝑛−1
+ 𝑦𝑛(𝑘) ℎ
(𝑘+1) = 𝑦(𝑘) + 𝑓 (𝑥, 𝑦(𝑘) , 𝑦(𝑘) , 𝑦(𝑘) , ⋯ , 𝑦(𝑘) ) ℎ
{𝑦𝑛 𝑛 1 2 3 𝑛
Tarea

Se construye un modelo equivalente de un

automóvil conduciendo en una vía con múltiples
resaltos. Este sistema se compone de un sistema
masa resorte amortiguador como el de la figura,
suponga que el sistema es excitado desde el
reposo, y la fuerza que genera los resaltos en la
amortiguación del auto se trata de la fuerza 𝐹(𝑡).
a. Determine en el instante 𝑡 = 1 𝑠 la posición de
la masa 𝑚 equivalente.

𝑘𝑔 𝑘𝑁
𝑚 = 1000 𝑘𝑔 𝑏 = 1000 𝑘 = 10
𝑠 𝑚
𝐹(𝑡) = 5000𝑒−2𝑡 sin 𝑡

𝑑2𝑥 𝑑𝑥
𝑚 2 +𝑏 + 𝑘𝑥 = 𝐹(𝑡)
𝑑𝑡 𝑑𝑡
𝑑 2 𝑥 𝐹(𝑡) 𝑏 𝑑𝑥 𝑘
= − − 𝑥
𝑑𝑡 2 𝑚 𝑚 𝑑𝑡 𝑚
𝑑𝑥
=𝑣
{ 𝑑𝑡
𝑑𝑣 𝐹(𝑡) 𝑏 𝑘
= − 𝑣− 𝑥
𝑑𝑡 𝑚 𝑚 𝑚
Donde 𝑥1 es la posición de la masa y 𝑥2 su velocidad
Aplicando el método de Euler explícito

𝑥𝑘+1 = 𝑥𝑘 + 𝑣𝑘 ℎ
{ 5000𝑒−2𝑡𝑘 sin 𝑡𝑘 𝑏 𝑘
𝑣𝑘+1 = 𝑣𝑘 + ( − 𝑣𝑘 − 𝑥𝑘 ) ℎ
𝑚 𝑚 𝑚
Aplicando Runge-Kutta Grado 3

1
𝑥𝑘+1 = 𝑥𝑘 + (𝑝1 + 4𝑝2 + 𝑝3 )ℎ
{ 6
1
𝑣𝑘+1 = 𝑣𝑘 + (𝑞1 + 4𝑞2 + 𝑞3 )ℎ
6
Teniendo en cuenta la forma en que se calcula las constantes que representan la pendiente
ponderada de una función.
𝑘1 = 𝑓(𝑡𝑘 , 𝑥𝑘 , 𝑣𝑘 )
 1 1 1
𝑘2 = 𝑓 (𝑡𝑘 + ℎ , 𝑥𝑘 + 𝑘1 ℎ , 𝑣𝑘 + 𝑘1 ℎ)
2 2 2
𝑘3 = 𝑓(𝑡𝑘 + ℎ , 𝑥𝑘 − 𝑘1 ℎ + 2𝑘2 ℎ , 𝑣𝑘 − 𝑘1 ℎ + 2𝑘2 ℎ)
Entonces
𝑝1 = 𝑣𝑘
1
𝑝2 = 𝑣𝑘 + 𝑝1 ℎ
2
𝑝3 = 𝑣𝑘 − 𝑝1 ℎ + 2𝑝2 ℎ

5000𝑒 −2𝑡𝑘 sin 𝑡𝑘 𝑏 𝑘

𝑞1 = − 𝑣𝑘 − 𝑥𝑘
𝑚 𝑚 𝑚
1 1
5000𝑒 −2(𝑡𝑘+2ℎ ) sin (𝑡𝑘 + 2 ℎ ) 𝑏 1 𝑘 1
𝑞2 = − (𝑣𝑘 + 𝑞1 ℎ) − (𝑥𝑘 + 𝑞1 ℎ)
𝑚 𝑚 2 𝑚 2
5000𝑒 −2(𝑡𝑘+ℎ ) sin(𝑡𝑘 + ℎ ) 𝑏 𝑘
𝑞3 = − (𝑣𝑘 − 𝑞1 ℎ + 2𝑞2 ℎ) − (𝑥𝑘 − 𝑞1 ℎ + 2𝑞2 ℎ)
𝑚 𝑚 𝑚
Elemento Tipo Cercha o Barra (Truss)
El problema de una barra de sección transversal 𝐴 constante, módulo de elasticidad 𝐸, densidad
𝜌 y longitud 𝐿, sometida a cargas (de cuerpo, distribuidas y/o concentradas) de tracción y/o
compresión. Se puede considerar que los desplazamientos en las direcciones transversales al
eje de la barra son despreciables cuando se comparan con los desplazamientos en la dirección
axial. En otras palabras, se supone que todos los puntos de una sección transversal de la barra
experimentan un mismo desplazamiento en la dirección axial.

Para un elemento de tipo resorte se sabe que la constante de rigidez tiene la forma
𝐹
𝑘=
𝑢
Así mismo se sabe que para una barra de área 𝐴 de longitud 𝐿 y módulo de elasticidad 𝐸 a la que
se le aplica una carga axial 𝐹, su deformación se calcula como se muestra a continuación:
𝐹𝐿
𝑢=
𝐴𝐸
Por consiguiente, si se relaciona esta deformación con la constante de rigidez de un elemento
tipo resorte, se puede determinar la constante de rigidez para un elemento de tipo barra.
𝐴𝐸
𝑘=
𝐿
Entonces la matriz de rigidez de este elemento puede ser escrita como.

𝐴𝐸 1 −1
𝑲= ( )
𝐿 −1 1
Adicionalmente la finalidad del uso de este tipo de elemento es el cálculo del esfuerzo que soporta
el elemento y para su cálculo se vale de las deformaciones calculadas y el conocimiento de:

𝜎 = 𝐸𝜖
Donde épsilon (𝜖) es la deformación unitaria del elemento, entonces:
𝑢2 − 𝑢1
𝜎 = 𝐸( )
𝐿
Ejercicio
Se tiene la barra empotrada de sección transversal como se observa en la figura.
Determinar los esfuerzos en los elementos.

𝐴2 = 350 𝑚𝑚2
𝐸 = 2 ∗ 104 𝑁/𝑚𝑚2 𝑃 = 100𝑘𝑁 𝐿1 = 𝐿2 = 150 𝑚𝑚 𝐿3 = 𝐿4 = 200 𝑚𝑚
𝐴1 = 300 𝑚𝑚2 𝐴2 = 200 𝑚𝑚2 𝐴3 = 250 𝑚𝑚2 𝐴2 = 200 𝑚𝑚2
Esta se puede discretizar como se muestra.

1 2 3 4 5

𝐴1 𝐸 1 −1 𝑢1 𝑓1,1 𝐴2 𝐸 1 −1 𝑢2 𝑓1,2
𝐸𝑙𝑒𝑚𝑒𝑛𝑡𝑜 1: ( ) (𝑢 ) = ( ) 𝐸𝑙𝑒𝑚𝑒𝑛𝑡𝑜 2: ( ) (𝑢 ) = ( )
𝐿1 −1 1 2 𝑓2,1 𝐿2 −1 1 3 𝑓2,2
𝐴3 𝐸 1 −1 𝑢3 𝑓1,3 𝐴4 𝐸 1 −1 𝑢4 𝑓1,4
𝐸𝑙𝑒𝑚𝑒𝑛𝑡𝑜 3: ( ) (𝑢 ) = ( ) 𝐸𝑙𝑒𝑚𝑒𝑛𝑡𝑜 4: ( ) (𝑢 ) = ( )
𝐿3 −1 1 4 𝑓2,3 𝐿4 −1 1 5 𝑓2,4

𝐴1 𝐴1
− 0 0 0
𝐿1 𝐿1
𝐴1 𝐴1 𝐴2 𝐴2
− + − 0 0 0 𝑅𝐴
𝐿1 𝐿1 𝐿2 𝐿2
𝐴2 𝐴2 𝐴3 𝐴3 𝑢2 0
𝐸 0 − + 0 𝑢3 = 𝑄
𝐿2 𝐿2 𝐿3 𝐿3 𝑢4 0
𝐴3 𝐴3 𝐴4 𝐴4 𝑢
0 0 + − ( 5) ( 𝑃 )
𝐿3 𝐿3 𝐿4 𝐿4
𝐴4 𝐴4
0 0 0 −
( 𝐿4 𝐿4 )
Elementos uniaxiales en 2 dimensiones

Nota: En el gráfico no es ilustrada la fuerza que se aplica en el nodo 2 ni la deformación del nodo
1.
Ya que todas las deformaciones y fuerzas a través del elemento ocurren en su dirección axial, la
matriz de rigidez, vector de deformación y fuerzas en el elemento se escribe de la siguiente forma.
′
1 −1 𝑢1 𝑓1′
𝑘( ) ( ′ ) = ( ′)
−1 1 𝑢2 𝑓2
Sin embargo se desea determinar una expresión en términos no de 𝑢′ sino en términos de 𝑢 y 𝑣
ya que 𝑢 y 𝑣 están alineados con los ejes 𝑥, 𝑦 globales del sistema. Para ello debemos tener en
cuenta que el sistema coordenado local del elemento está alineado con su eje 𝑥 local y las
deformaciones en el eje 𝑦 local se les dará el nombre de 𝑣 ′ . De esta forma se puede escribir la
deformación 𝑢′ y 𝑣 ′ en términos de 𝑢 y 𝑣

𝑢1′ = 𝑢1 cos (𝜃) + 𝑣1 sin (θ)

{ ′
𝑣1 = −𝑢1 sin (θ) + 𝑣1 cos (𝜃)
𝑢1′ cos(𝜃 ) sin(𝜃 ) 𝑢1
( ′) = ( ) (𝑣 ) → 𝑈 ′ = 𝑇𝑈
𝑣1 − sin(𝜃 ) cos(𝜃 ) 1
A la matriz de senos y cosenos se le llama matriz de transformación o matriz de rotación del
elemento y la denotaremos como 𝑇.
𝑇 ∶ 𝑀𝑎𝑡𝑟𝑖𝑧 𝑜𝑟𝑡𝑜𝑔𝑜𝑛𝑎𝑙 𝑑𝑒 𝑡𝑟𝑎𝑛𝑠𝑓𝑜𝑟𝑚𝑎𝑐𝑖ó𝑛
Por convención escribiremos 𝐶 = cos (𝜃) y 𝑆 = 𝑠𝑖𝑛 (𝜃).
Esta matriz de transformación también se puede utilizar para escribir la proyección de la fuerza
en cada nodo en términos de las coordenadas globales.
′
𝑓1𝑥 𝐶 S 𝑓1𝑥
( ′ )=( ) ( ) → 𝐹 ′ = 𝑇𝐹
𝑓1𝑦 −S C 𝑓1𝑦
′
Donde claramente 𝑓1𝑥 = 𝑓1′ y 𝑓1𝑦
′
=0
Y se puede escribir la matriz que expresa la transformación de ambos nodos del elemento así:

𝑢1′ 𝐶 𝑆 0 0 𝑢1
′ 𝑣1
𝑣1 −𝑆 𝐶 0 0
′ =( ) (𝑢 )
𝑢2 0 0 𝐶 𝑆 2
′
(𝑣2 ) 0 0 −𝑆 𝐶 𝑣2
Así mismo, si se realiza la extensión de la matriz de rigidez, vector de desplazamientos y fuerzas
en el sistema bidimensional local del elemento, escribiéndola de la siguiente forma.

1 0 −1 0 𝑢1′ 𝑓1′ 𝑥1 𝑙𝑜𝑐𝑎𝑙

′ 𝑦1 𝑙𝑜𝑐𝑎𝑙
0 0 0 0 𝑣1 0
𝐾𝑈 ′ = 𝐹 ′ → 𝑘 ( ) = ( ′ ) (𝑥 )
−1 0 1 0 𝑢2′ 𝑓2 2 𝑙𝑜𝑐𝑎𝑙
0 0 0 0 (𝑣2′ ) 0 𝑦2 𝑙𝑜𝑐𝑎𝑙
Donde claramente no existe fuerza, rigidez en la dirección 𝑦 local del elemento y por consiguiente
el desplazamiento así mismo será 𝑣1′ = 𝑣2′ = 0
Si se reemplaza 𝑈′ con 𝑇𝑈 se obtiene.

𝐾𝑇𝑈 = 𝐹 ′
Por las propiedades de las matrices ortogonales la transpuesta de 𝑇 es igual a su inversa 𝑇 𝑡 , y
es claro que 𝑇 −1 𝐹 ′ = 𝑇 𝑡 𝐹 ′ = 𝐹 entonces si se premultiplica la transpuesta de 𝑇 a ambos lados de
la ecuación:

𝑇 𝑡 𝐾𝑇𝑈 = 𝑇 𝑡 𝐹 ′
𝑇 𝑡 𝐾𝑇𝑈 = 𝐹
Y 𝑇 𝑡 𝐾𝑇 es la matriz de rigidez del elemento en el sistema global de coordenadas. Multiplicando
estas matrices se obtiene la siguiente expresión

𝐶2 𝑆𝐶 −𝐶 2 −𝑆𝐶 𝑢1 𝑓1𝑥

𝑘 ( 𝑆𝐶2 𝑆2 −𝑆𝐶 −𝑆 2 ) ( 𝑣1 ) = 𝑓1𝑦

−𝐶 −𝑆𝐶 𝐶2 𝑆𝐶 𝑢2 𝑓2𝑥
−𝑆𝐶 −𝑆 2 𝑆𝐶 𝑆 2 𝑣2 (𝑓2𝑦 )
𝐴𝐸
Siendo 𝑘 = si se utiliza un elemento tipo cercha.
𝐿
Entonces los esfuerzos en el elemento se calcularán como:
𝑢1
𝐸 𝑣1
𝜎 = (−𝐶 −𝑆 𝐶 )
𝑆 𝑢 )
(
𝐿 2
𝑣2
Ejercicio

𝑃𝑦

𝑃𝑥
3

𝑃𝑥 = 𝑃𝑦 = 100𝑘𝑁 𝐴 = 6 ∗ 10−4 𝑚2 𝐸 = 210𝐺𝑃𝑎 𝐿 = 1𝑚

Tabla de conectividad de elementos

Elementos Nodo i Nodo j
1 1 3
2 2 3
Elemento 1
El ángulo desde la horizontal hasta el elemento con respecto al nodo i es
𝜃 = 45°
1 1 1 1
− −
2 2 2 2
1 1 1 1 𝑢1 1 1 −1 −1 𝑢1
𝐴𝐸 2 − − 𝐴𝐸
2 2 2 (𝑣1 ) → (
1 1 −1 −1 𝑣1
) (𝑢 )
𝐿 1 1 1 1 𝑢3 2𝐿 −1 −1 1 1 3
− − 𝑣3 𝑣3
2 2 2 2 −1 −1 1 1
1 1 1 1
−
( 2 −
2 2 2 )
Elemento 2
El ángulo desde la horizontal hasta el elemento con respecto al nodo i es
𝜃 = −45°
1 −1 −1 1 𝑢2
𝐴𝐸 −1 1 1 −1 𝑣2
( ) (𝑢 )
2𝐿 −1 1 1 −1 3
1 −1 −1 1 𝑣3

Ensamble de ecuaciones

1 1 0 0 −1 −1 𝑢1 𝐹1𝑥
1 1 0 0 −1 −1 𝑣1 𝐹1𝑦
𝐴𝐸 0 0 2 0 −1 1 𝑢2 𝐹2𝑥
𝑣2 = 𝐹2𝑦
2𝐿 0 0 0 2 1 −1
−1 −1 −1 1 2 0 𝑢3 𝐹3𝑥
(−1 −1 1 −1 0 2 ) (𝑣3 ) (𝐹3𝑦 )

1 1 0 0 −1 −1 0 𝑅1𝑥
1 1 0 0 −1 −1 0 𝑅1𝑦
𝐴𝐸 0 0 2 0 −1 1 0 𝑅2𝑥
= 𝑅
2𝐿 0 0 0 2 1 −1 0 2𝑦
−1 −1 −1 1 2 0 𝑢3 𝑃𝑥
(−1 −1 1 −1 0 2 ) (𝑣3 ) ( 𝑃𝑦 )
Esfuerzos
0
𝐸 0
𝜎1 = (− √2 − √2 √2 √2) ( )
𝐿 2 2 2 2 𝑢3
𝑣3
0
𝐸 0
𝜎2 = (− √2 √2 √2 − √2) ( )
𝐿 2 2 2 2 𝑢3
𝑣3
Ejercicio

𝑃𝑥 = 𝑃𝑦 = 100𝑘𝑁 𝐴 = 6 ∗ 10−4 𝑚2 𝐸 = 210𝐺𝑃𝑎 𝐿 = 1𝑚

Elemento Tipo Viga
Los elementos viga ofrecen resistencia tanto a las fuerzas como a los momentos; la diferencia
entre las vigas y las barras es que las vigas soportan esfuerzos de torsión y flexión.
• Para aplicar correctamente un elemento viga se debe cumplir con lo siguiente:
• La longitud del elemento es mucho mayor que su ancho.
• Es constante en sus secciones y en sus propiedades.
• Puede trasferir momentos.
• Es capaz de manejar cargas distribuidas en su longitud.

La deformación vertical de la viga entre dos nodos, con fuerzas verticales y momentos
únicamente aplicados en sus nodos estará definida por un polinomio de grado 3 como se muestra.

𝑦(𝑥) = 𝑎1 + 𝑎2 𝑥 + 𝑎3 𝑥 2 + 𝑎4 𝑥 3
Y se sabe que el ángulo de rotación de la viga estará definido por la derivada de la deformación
vertical, entonces:

𝜃(𝑥) = 𝑎2 + 2𝑎3 𝑥 + 3𝑎4 𝑥 2

Teniendo en cuenta las condiciones de frontera de la viga (deformación y rotación en los nodos)

𝑦(0) = 𝑎1 = 𝑦𝑖
𝜃(0) = 𝑎2 = 𝜃𝑖
𝑦(𝐿) = 𝑎1 + 𝑎2 𝐿 + 𝑎3 𝐿2 + 𝑎4 𝐿3 = 𝑦𝑗
2
{ 𝜃(𝐿) = 𝑎2 + 2𝑎3 𝐿 + 3𝑎4 𝐿 = 𝜃𝑗
1 0 0 0 𝑎1 𝑦𝑖
0 1 0 0 𝑎2 𝜃𝑖
( ) (𝑎3 ) = (𝑦𝑗 )
1 𝐿 𝐿2 𝐿3
0 1 2𝐿 3𝐿2 𝑎4 𝜃𝑗

𝑎1 𝐿3 0 0 0 𝑦𝑖
𝑎2 1
(𝑎 ) = 3 ( 0 𝐿3 0 0 ) (𝜃𝑖 )
3 𝐿 −3𝐿 −2𝐿2 3𝐿 −𝐿2 𝑦𝑗
𝑎4 2 𝐿 −2 𝐿 𝜃𝑗

Adicionalmente se sabe que:

 𝑑2𝑦 𝑑3𝑦
𝐸𝐼 2 = 𝑀(𝑥) y 𝐸𝐼 3 = 𝑉(𝑥)
𝑑𝑥 𝑑𝑥
𝑀(𝑥) = 𝐸𝐼(2𝑎3 + 6𝑎4 𝑥) 𝑉(𝑥) = 𝐸𝐼(6𝑎4 )
𝑉(0) = 𝐸𝐼(6𝑎4 ) = 𝐹𝑖
𝑀(0) = 𝐸𝐼(2𝑎3 ) = 𝑀𝑖
𝑉(𝐿) = 𝐸𝐼(6𝑎4 ) = 𝐹𝑗
{𝑀(𝐿) = 𝐸𝐼(2𝑎3 + 6𝐿𝑎4 ) = 𝑀𝑗
0 0 0 6 𝑎1 𝐹𝑖
0 0 2 0 𝑎2 𝑀𝑖
𝐸𝐼 ( ) (𝑎 ) = ( 𝐹 )
0 0 0 6 3 𝑗
0 0 2 6𝐿 𝑎4 𝑀𝑗
Reemplazando el vector de coeficientes A:

0 0 0 6 𝐿3 0 0 0 𝑦𝑖 𝐹𝑖
𝐸𝐼 0
(
0 2 0
)( 0 𝐿3 0 0 ) (𝜃𝑖 ) = (𝑀𝑖 )
𝐿3 0 0 0 6 −3𝐿 −2𝐿2 3𝐿 −𝐿2 𝑦𝑗 𝐹𝑗
0 0 2 6𝐿 2 𝐿 −2 𝐿 𝜃𝑗 𝑀𝑗
El sistema que resulta de este producto es:

12 6𝐿 −12 6𝐿 𝑦𝑖 𝐹𝑖
𝐸𝐼 6𝐿 4𝐿2 −6𝐿 2𝐿2 𝜃𝑖 𝑀𝑖
( ) (𝑦𝑗 ) = ( 𝐹 )
𝐿3 −12 −6𝐿 12 −6𝐿 𝑗
6𝐿 2𝐿2 −6𝐿 4𝐿2 𝜃𝑗 𝑀 𝑗
Ejercicio 1

𝑙𝑏𝑓
𝑃1 = 1000 𝑙𝑏𝑓 𝐸 = 30 ∗ 106 𝐼 = 100 𝑖𝑛4 𝐿 = 20 𝑖𝑛
𝑖𝑛2
Elemento 1
12 6𝐿 −12 6𝐿 𝑦1
𝐸𝐼 6𝐿 4𝐿2 −6𝐿 2𝐿2 𝜃1
( ) (𝑦 )
𝐿3 −12 −6𝐿 12 −6𝐿 2
6𝐿 2𝐿2 −6𝐿 4𝐿2 𝜃2
Elemento 2
12 6𝐿 −12 6𝐿 𝑦2
𝐸𝐼 6𝐿 4𝐿2 −6𝐿 2𝐿2 𝜃2
( ) (𝑦 )
𝐿3 −12 −6𝐿 12 −6𝐿 3
6𝐿 2𝐿2 −6𝐿 4𝐿2 𝜃3
Ensamble de ecuaciones
12 6𝐿 −12 6𝐿 0 0 𝑦1 𝐹1
6𝐿 4𝐿2 −6𝐿 2𝐿2 0 0 𝜃1 𝑀1
𝐸𝐼 −12 −6𝐿 24 0 −12 6𝐿 𝑦2 𝐹2
𝐿3 6𝐿 2𝐿2 0 8𝐿2 −6𝐿 2𝐿2 𝜃2 = 𝑀2
0 0 −12 −6𝐿 12 −6𝐿 𝑦3 𝐹3
( 0 0 6𝐿 2𝐿2 −6𝐿 4𝐿 2 ) (𝜃 )
3 (𝑀 3)

𝐹1 = −𝑃 = −1000 𝑙𝑏𝑓 𝐹2 = 𝑅2 𝐹3 = 𝑅3
𝑦2 = 𝑦3 = 𝜃3 = 0 𝑀1 = 𝑀2 = 0
12 6𝐿 −12 6𝐿 0 0 𝑦1 −𝑃
6𝐿 4𝐿2 −6𝐿 2𝐿2 0 0 𝜃1 0
𝐸𝐼 −12 −6𝐿 24 0 −12 6𝐿 0 𝑅𝐴
=
𝐿3 6𝐿 2𝐿2 0 8𝐿2 −6𝐿 2𝐿 2 𝜃2 0
0 0 −12 −6𝐿 12 −6𝐿 0 𝑅𝐵
( 0 0 6𝐿 2𝐿2 −6𝐿 4𝐿2 ) ( 0 ) (𝑀𝐵 )
𝑦1 −𝑃
𝐸𝐼 12 6𝐿2 6𝐿
2 ) (𝜃 ) = (
(6𝐿 4𝐿 2𝐿 1 0 )
𝐿3 𝜃2
6𝐿 2𝐿2 8𝐿2
0
Ejercicio 2
Se tiene la viga doblemente empotrada de la figura.

𝑦1 = 𝑦3 = 𝜃1 = 𝜃3 = 0
Elemento 1

𝑤𝐿
−
2
12 6𝐿 −12 6𝐿 𝑦1 𝑤𝐿2 𝑅1
𝐸𝐼 6𝐿 𝜃1 −
4𝐿2 −6𝐿 2𝐿2
12 + ( 𝑀1 )
( ) (𝑦 ) =
𝐿3 −12 −6𝐿 12 −6𝐿 2 𝑤𝐿 𝑅21
𝜃2 −
6𝐿 2𝐿2 −6𝐿 4𝐿2
2 𝑀2
𝑤𝐿2
( 12 )
Elemento 2

𝑃
−
2
12 6𝐿 −12 6𝐿 𝑦2 𝑃𝐿 −𝑅21
𝐸𝐼 6𝐿 𝜃2 −
4𝐿2 −6𝐿 2𝐿2
8 + ( −𝑀2 )
( ) (𝑦 ) =
𝐿3 −12 −6𝐿 12 −6𝐿 3 𝑃 𝑅3
𝜃3 −
6𝐿 2𝐿2 −6𝐿 4𝐿2
2 𝑀3
𝑃𝐿
( 8 )
Ensamble de ecuaciones
𝑤𝐿 0
−
2 0
12 6𝐿 −12 6𝐿 0 0 𝑦1 𝑤𝐿 2 𝑃 𝑅1
− −
6𝐿 4𝐿2 −6𝐿 2𝐿2 0 0 𝜃1 12 2 𝑀1
𝐸𝐼 −12 −6𝐿 24 0 −12 6𝐿 𝑦2 𝑤𝐿 𝑃𝐿 0
𝜃 = − + − +
𝐿3 6𝐿 2𝐿2 0 8𝐿2 −6𝐿 2𝐿 2
2 2 8 0
0 0 −12 −6𝐿 12 −6𝐿 𝑦3 𝑤𝐿 2 𝑃 𝑅3
( 0 −
0 6𝐿 2𝐿2 −6𝐿 4𝐿2 ) (𝜃3 ) 12 2 (𝑀3 )
0 𝑃𝐿
( 0 ) ( 8 )
𝑤𝐿
−+ 𝑅1
2
𝑤𝐿2
− + 𝑀1
12 6𝐿 −12 6𝐿 0 0 0 12
6𝐿 4𝐿2 −6𝐿 2𝐿2 0 0 0 𝑤𝐿 𝑃
𝐸𝐼 −12 −6𝐿 − −
24 0 −12 6𝐿 𝑦2 2 2
=
𝐿3 6𝐿 2𝐿2 0 8𝐿2 −6𝐿 2𝐿 2
𝜃2 𝑤𝐿 2
𝑃𝐿
0 0 −12 −6𝐿 12 −6𝐿 0 −
12 8
( 0 0 6𝐿 2𝐿2 −6𝐿 4𝐿 2 )( )
0 𝑃
− + 𝑅3
2
𝑃𝐿
( 8 + 𝑀3 )
 𝑤𝐿 𝑃
𝐸𝐼 24 0 𝑦2 − −
2 2
( 2 ) ( 𝜃2 ) = ( 2 )
𝐿3 0 8𝐿 𝑤𝐿 𝑃𝐿
−
12 8
𝑤𝐿4 + 𝑃𝐿3
𝑦2 = −
48𝐸𝐼
𝑤𝐿3 𝑃𝐿2
𝜃2 = −
96𝐸𝐼 64𝐸𝐼
Reacciones
𝑤𝐿
−
2
12 6𝐿 −12 6𝐿 𝑦1 𝑤𝐿 2 𝑅1
−
𝐸𝐼 6𝐿
(
4𝐿2 −6𝐿 2𝐿2 𝜃1
) (𝑦 ) − 12 = 𝑀1
𝐿3 −12 −6𝐿 12 −6𝐿 2 𝑤𝐿 𝑅21
𝜃2 −
6𝐿 2𝐿2 −6𝐿 4𝐿2 2 ( 𝑀2 )
𝑤𝐿2
( 12 )
𝑤𝐿
0 2
12 6𝐿 −12 6𝐿 0 𝑤𝐿2 𝑅1
4 3
𝐸𝐼 6𝐿 4𝐿2 −6𝐿 2𝐿2 𝑤𝐿 + 𝑃𝐿 12 𝑀
( ) − + = ( 1)
𝐿3 −12 −6𝐿 12 −6𝐿 48𝐸𝐼 𝑤𝐿 𝑅21
3 2
6𝐿 2𝐿2 −6𝐿 4𝐿2 𝑤𝐿 𝑃𝐿 2 𝑀2
−
(96𝐸𝐼 64𝐸𝐼 ) 𝑤𝐿2
−
( 12 )
13𝑤𝐿 5𝑃
𝑅1 = +
16 32
11𝑤𝐿2 3𝑃𝐿
𝑀1 = +
48 32
3𝑤𝐿 5𝑃
𝑅21 = −
16 32
𝑤𝐿2 𝑃𝐿
𝑀2 = +
12 16