ecuaciones diferenciales

Vista previa del material en texto

Ecuación del calor: El modelo 
CONCEPTOS PREVIOS 
Masa y densidad: Si V es un volumen ocupado por un material (generalmente no homogéneo) y 
Ω es una parte de V, la propiedad extensiva, masa de Ω, se obtiene por integración sobre Ω de la 
magnitud intensiva ρ(x, y, z) llamada densidad de masa en el punto (x, y, z) de Ω , es decir 
 m(Ω) = ∭ ρ�x, y, z�dxΩ dy	dz			 
Cilindro unidimensional: Si V es un cilindro de longitud L, diremos que V puede modelarse en una 
dimensión o que V es unidimensional, si todas sus propiedades físicas son constantes en cada 
sección transversal al eje del cilindro. En tal caso tendremos por ejemplo que, ρ(x, y, z) = ρ(x) y 
entonces 	m(Ω) = A
 ρ�x�dx		�� cuando Ω es la porción de cilindro comprendida entre a < x < b y 
A es una constante positiva con la cual indicamos el área de cualquier sección transversal de V. 
Calor específico y densidad de energía térmica: El calor específico es la cantidad de energía 
térmica (e.t.) necesaria para aumentar en un unidad de temperatura, una unidad de masa del 
material que constituye el cuerpo. Si denotamos con c(x) al calor específico del material que 
ocupa V en cualquier punto de la sección en x, y con u(x,t) a la temperatura en cualquier punto 
de dicha sección en el instante t , entonces la densidad de energía térmica, e(x,t), es 
e(x,t) = c(x) ρ(x) u(x,t) 
La energía térmica total en Ω en el instante t, se calcula y la denotamos de la siguiente manera: 
 					E�Ω, t� = 	� 
 c�x��� ρ�x�	u�x, t�dx				 
Vector de flujo de calor: es un campo vectorial definido dentro de V y para cada instante t, que 
denotamos con �����, �, �; ��, que mide la energía térmica que atraviesa el plano normal a ��� por el 
punto (x,y,z), por unidad de área y de tiempo. 
Ley de Fourier (1807): El flujo es un múltiplo negativo del gradiente de temperatura y siendo K0 la 
conductividad térmica del material, se cumple que 
 �����, �, �; �� = 	−!"��, �, ��∇���$��, �, �; �� 
También conocida como Ley de Ohm si u potencial electrostático y como Ley de Fick cuando u es 
concentración química. 
Caso de flujo unidimensional: Si aislamos la cara curva del cilindro unidimensional, resulta que ��� 
tiene sólo componente horizontal y denotamos con ���, �� = 	−!	"(x) %&%' ��, ��	 a la cantidad 
de energía térmica que fluye de izquierda a derecha por cualquier punto de la sección en x, por 
unidad de área y de tiempo. 
DEDUCCIÓN DE LA ECUACIÓN UNIDIMENSIONAL DEL CALOR 
Balance de energía térmica: aquí se usa el Principio de Conservación de la Energía que establece 
( )*	+*�ó-	./	0*1234		�/154+*6	./	/. �. /-	Ω	/-	/6	3-8�*-�/	� 9 = 	 (
)*	/. �. :$/	/-�+*	54+	648	/��+/148	./	Ω		54+	$-3.*.	./	�3/154	/-	/6	3-8�*-�/	� 9 + (
)*	/. �. </-/+*.*	./-�+4	./	Ω		54+	$-3.*../	�3/154	/-	/6	3-8�*-�/	�9 
Reeplazando en términos de la notación introducida anteriormente, resulta 
																																						=>=? 	�Ω, �� 	= 			 @	�	��*, �� 	− 	�	��2, ��	A 				+ 					�	 
 B��, ��.��� 
donde la función Q(x,t) representa la fuente interna de calor por unidad 
 de área en el instante t, que puede deberse a reacciones químicas o calentamiento por medio de 
electricidad. 
Luego, reemplazando la e.t. total E y aplicando el teorema fundamental del cálculo para los 
términos en � (suponiendo que %C%' es continua para � ∈ @*, 2A	) tenemos, 
�	E 0���F�����
G$G� ��, ��.�	 = 	−�	 H		E G�G� ��, ��.�
�
� 	I 				+ 					�	 EB��, ��.�
�
� 
Simplificando A, reemplazando el flujo � y escribiendo una sola integral en el 2do miembro, 
queda 
	E 0���F�����
G$G� ��, ��.�	 = 	E J GG� J!	"���G$G� ��, ��K + 	B��, ��K .��� 					 
Como esta igualdad de integrales es válida toda elección que se haga de a y b con 0 < a < b < L , 
podemos concluir que necesariamente los integrandos deben ser iguales, es decir 
0���F��� G$G� ��, �� = 	 GG� L!	"���G$G� ��, ��M	+ 	B��, �� 
Ecuación que se conoce como ecuación unidimensional del calor. 
El caso más simple para esta ecuación ocurre en un medio homogéneo (propiedades físicas 
constantes en V) y sin fuentes internas, es decir 0��� = 0	, F��� = F	, !	"��� = !	"	 y B��, �� = 0. 
Así la ecuación se puede escribir, denotando con k = 
OPQR	 a la llamada constante de difusividad 
térmica, como 	%&%? ��, �� = S	 %T&%'T ��, �� (*) 
Nota: La ecuación del calor también se llama ecuación de difusión cuando u representa la 
concentración de alguna sustancia química o ecuación del telégrafo cuando u es potencial 
eléctrico. 
Observación: Notemos que la ecuación del calor (*) que también se escribe como $? = S. $'' 
tiene infinitas soluciones, por ejemplo: 
u(x,t) = 0 , u(x,t) = x , u(x,t) = 2x + 1 , u(x, t) = a.x + b , u(x, t) = t + x2/(2k) , etc . Por lo cual para 
obtener unicidad de la solución es necesario imponer condiciones en los bordes y una condición 
inicial (3 condiciones deben agregarse puesto que la suma de los órdenes de las derivadas que 
 intervienen en la ecuación diferencial parcial es 3) 
CONDICIONES DE BORDE: Consideraremos 4 tipos de condiciones de borde (frontera o contorno) 
1) Condición de Dirichlet, cuando se prescribe en los bordes la variable dependiente u. 
Ejemplo: Si u es temperatura y u(0,t) = A(t) y u(L, t) = B(t) , se dice que las temperaturas 
prescritas en los bordes izquierdo y derecho del cilindro V son A(t) y B(t) respectivamente. 
2) Condición de Neuman, cuando se prescribe en los bordes el flujo. 
En el caso unidimensional, es suficiente con dar el valor de la derivada parcial de u con respecto a 
x en x = 0 y x = L, por ejemplo 
%&%' �0, �� = 	�"	 , %&%' �), �� = 	0	. En este último caso se dice que 
el borde derecho está aislado ya que el flujo resulta ser 0. 
3) Condición de Robin, es una combinación de las anteriores. 
 Por ejemplo, si los extremos del cilindro V están en contacto con un fluido (o aire) corresponde 
modelar la situación con la Ley de enfriamiento de Newton que establece “el flujo de calor que 
sale por cada extremo es proporcional a la diferencia entre la temperatura en ese extremo y la 
temperatura del medio externo (fluido o aire)”. Así las ecuaciones en x=0 y x= L serán, 
−!	"�0� G$G� �0, �� =	−U@$�0, ��− $V���A	 y −	!	"�)� G$G� �), �� 	= 		U@$�), ��− $V���A . 
La constante positiva H se llama coeficiente de transferencia de calor o de convección. Notar que 
el signo menos en el extremo izquierdo, responde al hecho de que si el cilindro está más caliente 
que el medio, pierde calor por los extremos, por lo cual en x = 0 la dirección del flujo saliente es de 
derecha a izquierda (o sea � < 0�. 
4) Contacto térmico perfecto: Se dice que dos cilindros unidimensionales de materiales 
diferentes, unidos en x= x0 , están en contacto térmico perfecto cuando la temperatura es continua 
en x = x0 , es decir u(x0
-, t) = u(x0
+, t) , y no se pierde calor en x = x0 (o sea, el flujo de calor que sale 
de uno fluye dentro del otro). Esto último, se modela con ��	x"X, �� = 	��	x"Y, ��	. 
Problemas de valores en el borde 
Un problema a valores en el borde (PVB) para la temperatura en un cilindro unidimensional 
homogéneo, sin fuentes internas de calor, con temperaturas prescritas en los bordes y 
distribución inicial de temperatura dada por la función f(x), puede escribirse como 
�Z[V� ∶ 	 ] ^_Z:				$? = 	S. $''				,						5*+*			0	 < 	� < 	)aV:		$�0, �� = ����	; 	$�), �� = V���	; 	�	 ≥ 0ac:			$�	�, 0� = 	d���		,									5*+*		0	 ≤ �	 ≤ )f 
Antes de abordar un método analítico para resolver estos problemas con ecuaciones diferenciales 
parciales (EDP) y valor inicial (CI) analizaremos una cuestión física relacionada con ecuaciones 
diferenciales ordinarias (EDO). 
Distribución de temperatura en equilibrio 
Cuando transcurre cierto tiempo, teóricamente infinito, que en la práctica depende del tipo de 
material que empleamos, se establece un estado estacionario en el que la temperatura de cada 
punto de la barra cilíndrica no varíacon el tiempo, matemáticamente lim?→k $��, �� = $��� y $? = 0 pues u = u(x). 
En estos casos, si se conocen las temperaturas en los bordes, corresponde considerarlas 
constantes y el PVB será: 
�PVB� ∶ 	
opq
pr^_s:		 .t$.�t = −B���!" 		 , 0 < � < )			aV1:			$�0� = 	vw																															aV2:			$�	)� = 	vt																													
f 
siendo Q(x) el término fuente y K0 la conductividad térmica. 
Llamaremos solución en estado estacionario o estado de equilibrio a la función u(x) solución de 
cualquier PVB independiente del tiempo. 
Ejemplo 1: Cuando no hay fuentes internas y las temperaturas están prescritas en los extremos, 
y 	=T&='T = 0																																			$�0� = 	vw																															$�	)� = 	vt																													 fintegrando 2 veces la EDO obtenemos su solución general, 
u(x) = a x + b , donde los coeficientes a y b se obtienen aplicando las condiciones de borde, 
u(0)= b ⇒ b = vw y u(L) = a L + vw =	vt ⇒ a = ⇒ 		b	 = �vt − 	vw)/L . Por lo tanto , 
	$��� = vw +	|TX	|}~ � es la solución estacionaria (UNICA) cuando las temperaturas están 
prescritas en los extremos x = 0 y x = L. 
Ejemplo 2: Cuál será la solución cuando los bordes están aislados? 
a) Plantear el problema y resolver. 
b) Notar que hay infinitas soluciones de equilibrio, u(x) = b (b constante) 
Comentamos el resultado: Luego de un tiempo suficientemente largo, cuando los bordes también 
están aislados la temperatura en la barra será constante, sin embargo no es razonable que sea 
cualquier constante. En efecto, tenemos lim?→k $��, �� = $��� = b y como la energía térmica 
total es constante pues 
=>=? = ����0� − 	��)�� = 0 (ya que los flujos son cero en los bordes) 
podemos escribir E(0) = E (∞ ) = lim?→k ^���. Calculando primer y tercer miembro a partir de 
la fórmula para la energía térmica total,			E�	t� = 	� 
 c�x�~" ρ�x�	u�x, t�dx		,		 con calor específico y 
densidad constantes, resulta que en t =0 y t = ∞, tenemos la siguiente igualdad 
 			�0	ρ 
 u�x, 0�dx	~" 	=			 �0	ρ 
 lim?→ku�x, t�dx	~" . Cancelando factores comunes y reemplazando 
por el valor constante b de la solución estacionaria, obtenemos 
 u�x, 0�dx	 =~" 2	L . 
Por lo tanto, la solución de equilibrio para bordes aislados también es ÚNICA, u(x) = 
w~ 
 u�x, 0�dx	~" 
y es el promedio de la temperatura inicial. 
Comentario: Cuando la barra está totalmente aislada, la barra “no se olvida” por completo de la 
condición inicial. 
Ejercicios propuestos: Sección 1.4 pág 19 (Haberman): 1.4.1 (f) y (g) – 1.4.2 – 1.4.4 – 1.4.5 – 1.4.7 . 
------------------------------------------------- 
BIBLIOGRAFÍA: Haberman,R. “Applied Partial Differential Equations” Prentice-Hall (2004)-Capítulo 1 
Zill y Cullen “Ecuaciones diferenciales con aplicaciones de valores en la frontera” Thomson (2006) 
 
Material elaborado por Adriana Frausin (2014) 
Método de Separación de variables para EDP 
Para resolver un problema de valores en el borde con ecuaciones diferenciales parciales, muchas 
veces se puede aplicar directamente esta técnica (inventada por Daniel Bernoulli en los 1700s) 
conocida como método de separación de variables. 
Concretamente este método puede usarse cuando la ecuación diferencial parcial (EDP) y las 
condiciones de borde (CB) del problema son lineales (ecuaciones de 1er grado en la variable 
dependiente y sus derivadas) y homogéneas (admiten la solución trivial pues cada término 
contiene la variable dependiente o una derivada). 
La propiedad fundamental de una ecuación diferencial lineal y homogénea es que cualquier 
combinación lineal de soluciones conocidas también la satisface (Principio de Superposición) 
Existen problemas que NO CUMPLEN con estas condiciones de linealidad y homogeneidad pero 
pueden convertirse, a partir de cambios de variables adecuados, en (PVB) donde la EDP y sus CB 
resultan lineales y homogéneas, pudiendo aplicar el método a las “nuevas variables”, para luego 
expresar la solución del problema original en término de las variables originales. 
Ejemplos: 
1) La ecuación no lineal �� =	��� +	��� puede reducirse a la ecuación del calor, �� =	��� , 
mediante el cambio de variables � =		
 (VERIFICARLO). 
2) Aplicación en economía: El valor de un derivado financiero se basa en el precio de otro activo 
(subyacente) S y depende del tiempo t. Cualquier derivado sobre S, V(S, t), debe verificar la 
siguiente ecuación diferencial parcial no lineal conocida como Modelo de Black-Scholes: 
 
���� + 	
���� ���� +	�� 	����, ���� + ������ − 	
� = 0 donde r es la tasa de interés y � la desviación 
de los cambios de las tasas de cambio. Esta ecuación también se reduce a la ecuación del calor, �� =	��� , si se le aplican los siguientes cambios de variables : � = 		�����������	������������ , !�,"� = 2��
	,			� = $	� ,			� = % − 	2��!		, siendo K el precio de ejercicio de la opción de compra y 
T su vencimiento. 
Separación de variables para un (PVB) con ecuación del calor 
 Interesados en poder predecir cómo cambia la energía térmica inicial (representada en la 
condición inicial CI ) a medida que pasa el tiempo en una situación física relativamente sencilla, 
vamos a aplicar la técnica del método de separación de variables, a la ecuación del calor en un 
cilindro unidimensional con temperatura prescripta cero en los extremos y sin fuentes internas de 
calor. Esta descripción nos permite formular el (PVB) con las siguientes ecuaciones (EDP y CB) 
lineales y homogéneas: 
�&�'� ∶ 	 ) *+&:				�� = 	". ���				,						./
/			0	 < 	 < 	12':		��0, �� = 0	; 									��1, �� = 0		; 					�	 ≥ 025:			��	 , 0� = 	6� �		,									./
/		0	 ≤ 	 ≤ 18 
El método de separación de variables consiste en determinar soluciones del problema que tengan 
forma de producto de funciones de “variables separadas”, es decir: 
Paso 1) Proponer 			�� , �� = 	∅� �:���	 donde ∅ es sólo función de x y G sólo función de t. 
Ahora como la función u propuesta debe satisfacer el (PVB), comenzamos por reemplazarla en la 
EDP, así obtenemos ∅� �:´��� = ". ∅´´� �:��� y para “separar variables” dividimos ambos 
miembros por ". ∅� �:��� ≠ 0 (ya que suponemos que existe solución no trivial). Entonces, 
: ´���". :��� = 	∅´´� �∅� � = 	−= 
ya que la única posibilidad de que una función solamente de t sea igual a una función solamente 
de x, es que ambas expresiones sean constantes. Constante que denotamos con – = ( =	 se conoce 
como constante de separación y el signo menos se introduce por conveniencia). 
De esa doble igualdad resultan las siguientes EDO, una para ∅ y otra para : : ∅´´� � = 	−=. ∅� � 
y : ´��� = −=. ". :���. 
Paso 2) Resolver el problema a valores en el borde para ∅� �: 
Ahora reemplazamos la función propuesta en las condiciones de borde CB del problema original, 
��0, �� = 	∅�0�:��� = 0				 ⇒ 		∅�0� = 0 pues :��� 	≠ 0		 sino u sería la solución trivial, del 
mismo modo , 
��1, �� = 	∅�1�:��� = 0		 		⇒ 			∅�1� = 0 . Por lo tanto el problema para ∅� � queda 
�∗� )∅´´� � = −=. ∅� �∅�0� = 0																∅�1� = 0																 8. Veamos para qué valores de = este problema tiene soluciones no triviales: 
• Si A < 0			 ⇒ 		∅� � = / cosh�√−= . � + G sinh�√−=. � 			 		 
Notar que podemos proponer esta combinación lineal de funciones hiperbólicas en vez de 
las funciones exponenciales linealmente independientes 	±	√�K	�. Por qué? 
Ahora reemplazamos esta ∅	en las condiciones de borde de (*) para hallar /		 y G ∶ 
														∅�0� = / cosh� 0� + G sinh�0� = /. 1 + G. 0		 ⇒ 		/ = 0 
∅�1� = 0. cosh�√−= . 1� + G sinhM√−=. 1N = 0		 ⇒ 		G = 0	, pues sinh(x) =0 sólo si x = 0. 
Por lo tanto cuando 	= < 0			la solución de (*) es ∅� � = 0	�OPQRSTóV		WXTYTZQ�. 
 
 
• Si A = [		 ⇒ 		∅� � = /	 + G		 ⇒		 ∅�0� = G		 				⇒ 		G = 0∅�1� = /. 1		 ⇒ 		/ = 0						⇒ ∅�\� = [	�OPQ. WXTYTZQ� 
 
Comentario: Se dice que los 	= ≤ 0 no son autovaloresdel problema a valores en el borde (*) 
• Si A > 0			 ⇒ 		∅� � = / cos�√= . � + G sinM√=. N			 
 Ahora reemplazamos esta ∅	en las condiciones de borde de (*) para hallar /		 y G ∶ 
									∅�0� = / cos� 0� + G sin�0� = /. 1 + G. 0		 ⇒ 		/ = 0 
									∅�1� = 0+	G sinM√=. 1N = 0		 ^_`abc		√=. 1 =d. e		 ⇒ 		A = �V.fg �h						para n = 0, 1, 2, 3, … 
son los autovalores del problema, o sea los que producen las soluciones no triviales 
 ∅� � = G. sin iV.fg j , llamadas autofunciones del problema a valores en el borde (PVB) 
 
Paso 3) Resolver la EDO para :��� para los autovalores hallados: 
 Reemplazamos A en la ecuación : ´��� = −=. ". :���		y resolvemos, 
 : ´��� = −iV.fg jh . ". :��� 	⇒					 	:��� = k. 	�iV.fg jh.�.�			 	 para cualquier constante k ∈ m. 
Paso 4) Aplicar el principio de superposición (extendido) a las soluciones producto halladas: 
 :���. ∅� � = G. k. 	�iV.fg jh.�.�			. sin iV.fg j para n = 0, 1, 2, 3, … 
Notamos que tenemos infinitas soluciones producto, una para cada n, por lo tanto será una 
combinación lineal infinita, la expresión para la solución u(x, t) del problema original (PVB); 
en la cual, incorporando b.c a los escalares, denotados con An, resulta que 
�� , �� = 	nop	�iV.fg jh.�.�			. sin id. e1 j
q
pr� 
Paso 5) Determinar los coef. An a partir de la aplicación de la condición inicial CI del (PVB): 
 Evaluamos la expresión anterior para t = 0 y la igualamos al dato inicial f (x), entonces 
																�� , 0� =	 	∑ op	�iV.fg jh.�.`			. sin ip.tu j			⇒ 		6� � = 		∑ op. sin ip.tu jqpr�qpr� 
 y notamos que ésta última expresión para f(x) es la serie de Fourier de senos de f(x) . 
 Sabemos que la serie de Fourier de los senos de f es igual a la serie de Fourier de la 
 extensión impar de f en el intervalo ( - L , L), por lo cual sus coeficientes (coeficientes de 
 Fourier para el desarrollo de semi-intervalo) se calculan resolviendo la siguiente integral: 
op =	21		v6� �. sin id. e1 j . w 
u
`
 
Por lo tanto, la solución analítica general del problema (PVB) es 
 
		�� , �� = 	n 21	�	v 6� � sin	 de L dx	�.u`
q
pr� 	�i
V.fg jh.�.�	 sin id. e1 j	 
OBSERVACIONES: 
1- Los coeficientes de Fourier op		no siempre son “fáciles” de evaluar, en cuyo caso se necesita 
aplicar integración numérica en una computadora. 
2- Se obtienen buenas aproximaciones para u(x, t) truncando la serie para “unos pocos términos” 
debido al rápido decrecimiento exponencial, cuadrático en n (VER Ejemplo 2.3.7 pág 52-53). 
3- Después de mucho tiempo, la condición inicial f(x) poco importa y las condiciones de borde van 
a dominar, en efecto es natural esperar que el perfil de temperatura se aproxime a la temperatura 
de equilibrio (convergencia al estado de equilibrio). La distribución de temperatura en estado 
estacionario para este problema de temperaturas cero en los bordes es la solución trivial, �� � = lim�→q �� , �� = 0	(la solución evolutiva o transiente converge a la solución estacionaria). 
4- Se obtiene SOLUCIÓN EXACTA para el (PVB) cuando la condición inicial f(x) es una combinación 
lineal de las autofunciones del problema. Ejemplos: RESOLVER Ejercicio 2.3.3 (a) y (b) de pág 54. 
Problemas propuestos: libro “Applied Partial Differential Equations” R. Haberman (2004) 
Pág 54-55 (EXERCISES 2.3) : 2.3.3 y 2.3.4 
Pág 59 (Sección 2.4 Worked Examples with de Heat Equation) : 2.4.1 (Extremos aislados ) y 2.4.2 
(Conducción del calor en un anillo circular delgado) 
Pág 69 (EXERCISES 2.4) : 2.4.1 (a) y (b) , 2.4.2 y 2.4.6 
-------------------------------------------------- 
 
ECUACIÓN de LAPLACE: Soluciones y propiedades cualitativas 
La ecuación del calor en 2 o 3 dimensiones para u = u(x,y,t) toma la siguiente forma 
k� , ~��� , ~� �
�� � , ~, �� = 	∇. M$	`� , ~�∇u� , ~, ��N 	+ 	�� , ~, �� donde ∇ es el operador nabla, 
∇u	es	el	gradiente	de	la	temperatura		y		∇. �$	`∇u�		es	la	divergencia	del	campo	vectorial	�$	`∇u�.	 
En los casos en que no hay fuentes internas de energía térmica (Q = 0) y los coeficientes térmicos 
son constantes, la ecuación queda:			 	�
	�� 		= 	"	∇�u , donde 	∇�u = uxx+ uyy es el Laplaciano de u. 
Ahora si las condiciones de borde son estacionarias (independientes de t) es posible obtener una 
solución de equilibrio, es decir una u = u (x, y), solución de la ecuación diferencial parcial 
 	∇�u = 0 conocida como Ecuación de Laplace o ecuación potencial (pues en ausencia de fuentes 
la satisfacen el potencial electrostático y el gravitacional) 
Comentario: Cuando hay un término fuente Q, la ecuación 	∇�u = -Q se conoce como ecuación 
de Poisson (para el caso de potencial electrostático, Q es la densidad de carga). Esta ecuación 
también aparece en el cálculo de deflexiones luego de aplicar el cambio de variables � = ��� +	���	 en la siguiente ecuación diferencial parcial que modela las deformaciones (VERIFICAR) 
������ 	+ 2	 ������	��� +	 ������ 	= 	 �			→		�����	���	
p����	��	á����		→		�������	������p��																		 
 
Separación de variables para la ecuación de Laplace en un rectángulo 
Sea u (x, y) la temperatura de equilibrio en un rectángulo R ={ (x,y): 0 ≤ x ≤ L , 0 ≤ y ≤ H } cuyos 
bordes laterales están aislados, con temperatura u(x, 0) = 0 en el borde inferior y u(x, H) = f(x) en 
el borde superior. El (PVB) que modela esta situación es (Ejercicio 2.5.1 (a) pág 85): 
�&�'� ∶ 	 ) *+&:				��� +	���	 = 0		,						./
/								0	 < 	 < 	1,				0 < ~ < �							2'1,2:			���0, ~� = 0	,								���1, ~� = 0		,											0 ≤ ~	 ≤ �																		2'3,4:			�� , 0� = 0			,								�� , �� = 6� �		,						0 ≤ 	 ≤ 1																				8 
Podemos resolverlo aplicando el método de separación de variables (por qué?) En efecto: 
Paso 1) Proponemos 			�� , ~� = 	∅� �	ℎ�~�	 reemplazamos en la EDP y separamos variables: 
∅´´� �∅� � = 	−	ℎ´´�~�ℎ�~� 	= 	−=		 ⇒  ∅´´� � = −=	∅� �ℎ´´�~� = =	ℎ�~� 8 
 
Paso 2) Resolvemos el problema a valores en el borde con 2 condiciones de borde homogéneas: 
En este caso, corresponden a x =0 y x = L las condiciones de borde homogéneas CB1,2 (lados 
aislados), entonces resolvemos primero el problema para ∅� �: 
�∗� )∅´´� � = −=. ∅� �∅´�0� = 0																∅´�1� = 0																8.Veamos para qué valores de = este problema tiene soluciones no triviales: 
• Si A < 0			 ⇒ 		∅� � = / cosh�√−= . � + G sinh�√−=. � 			 		 
Ahora ∅´�0� = √−=M/ senh�√−= . 0N + G coshM√−=. 0N =√−=. G = 	0		 ⇒ 	G = 0 ∅´�1� = √−=M/ senh�√−= . 1N = 0 sólo si a = 0. Por lo tanto cuando 	= < 0			la solución de 
(*) es ∅� � = 0	�OPQRSTóV		WXTYTZQ�.	 
 
• Si A = [		 ⇒ 		∅� � = /	 + G		 ⇒		 ∅′�0� = /		 	⇒ 		/ = 0∅′�1� = /			 ⇒ 		/ = 0	 	⇒∅�\� = ¢	�SPVOWZVW£� , 
entonces A = [		es autovalor de esta problema (pues produce solución no trivial) 
 
• Si A > 0			 ⇒ 		∅� � = / cos�√= . � + G sinM√=. N			 
 									∅¤�0� = 	√=�−	/ sen� 0� + G cos�0�� = √=�	G� = 0		 ⇒ 		G = 0 
									∅′�1� = √=	�−	/ senM√=. 1N = 0		 �_`abc		√=. 1 =d. e		 ⇒ 		A = �V.fg �h			para n = 0, 1, 2, 3, … 
Entonces los A > 0	 también son autovalores del problema, cuyas autofunciones son 
∅�\� = Z. ¥¦§ iV.fg \j para n = 0, 1, 2, 3, … 
 
Paso 3) Para los autovalores hallados resolver la EDO para h�¨�	con la restante CB homogénea: 
 Para A = [		 ⇒		 ©ℎ¨�~� = 0ℎ�0� = 0 8 			⇒	…		⇒							 	¬�¨� = Z	¨	 	 
 Para A = iV.fg jh 		⇒		 ­ℎ¨�~� = 	 ip.tu j
� 	ℎ�~�ℎ�0� = 0																	 8 		⇒	 
 	)			ℎ�~� = / cosh pt�u + 	G	®	dℎ	 pt�u 																																ℎ�0� = / cosh pt.`u + 	G	®	dℎ	 pt.`u = 0	 ⇒ 	/ = 0 	8 				⇒		 ¬�¨� = ¢. O£V¬	
Vf¨g 		 
 
 
Paso 4) Aplicar el principio de superposición (extendido) a las soluciones producto halladas: 
�� , ~� = 	o`~	̄Ar[ +		nop		®	dℎ	�
de~1 �cos id. e1 j
q
pr�°±±±±±±±±±²±±±±±±±±±³A´`
 
Paso 5) Determinar loscoef. An a partir de la aplicación de la condición borde NO homogénea: 
�� , �� = o`�	̄µ¶· +		nop		®	dℎ	�
de�1°±±±±²±±±±³µ·̧ �cos i
d. e1 j = 6� �
q
pr� 
 con o·̀ = 	 �u 		¹ 6� �w ù y op· =	 �u 		¹ 6� �. cos ip.tu j . w ù de donde 
�� , ~� = º 1�. 1 		v6� �w 
u
`
»~°±±±±±²±±±±±³A	r	[
+		n 	º 21. ®	dℎ	�de�1 �		v6� � cos i
d. e1 j w 	
u
`
» 		®	dℎ	�de~1 �cos id. e1 j
q
pr�°±±±±±±±±±±±±±±±±±±±±±±±²±±±±±±±±±±±±±±±±±±±±±±±³A	´	`
 
Esta es la expresión analítica para la solución de la ecuación de Laplace en un rectángulo. 
Observación: Siempre es posible formular los problemas “de Laplace” en un rectángulo 
descomponiéndolos en problemas como el anterior, en el cual de las cuatro CB sólo una sea no 
homogénea y se “deje para el final” es decir, “se la trate” como CI en el último paso del 
procedimiento, para determinar los coeficientes de la serie. 
Ejercicio propuesto: 2.5.1(b) de pág 85. Ayuda: al resolver el problema de valor inicial para ∅� � 
es conveniente proponer como solución, una combinación lineal de funciones trasladadas, esto es 
 ∅� � = / cosh√= � − 1� + b	senh√= � − 1�		pues la condición homogénea está dada en x= L. 
Comentario: Para el caso de Ecuación de Laplace en un círculo de radio R, también se puede 
aplicar el método de separación de variables pero luego de usar “las coordenadas polares (r, ½�” 
para reescribir el problema, que de esta manera pasa a tener un dominio matemático rectangular, 
con 0 < r < R y −	e < 	½ < e donde deben imponerse condiciones de borde apropiadas, como ser 
condiciones de regularidad (contacto térmico perfecto) en ½ = ±	e. (Sección 2.5.2: OPTATIVO) 
Propiedades cualitativas 
A menudo, el método de separación de variables no resulta apropiado. Si se desea información 
cuantitativa de la solución, suelen ser necesarios y mucho más eficientes los métodos numéricos 
que veremos más adelante. Sin embargo, al momento de determinar propiedades cualitativas de 
las soluciones de gran interés práctico, es el método de separación de variables el que nos permite 
obtener estas respuestas: 
1-Propiedad del valor medio para la ecuación de Laplace: Si u es solución de la ecuación de 
Laplace en una región arbitraria Ω, siempre es posible considerar para cualquier punto P dentro de 
Ω, un círculo de radio r0 (con r0 suficientemente chico para que el círculo esté contenido en Ω) y 
se puede demostrar (usando el teorema del valor medio) que 
En estado de equilibrio (estacionario)la temperatura en cualquier punto P interior a Ω es igual al 
promedio de la temperatura alrededor de cualquier círculo con centro en P contenido en Ω 
2-Principio del máximo: A partir de la propiedad anterior, se puede concluir que 
En estado de equilibrio, la temperatura no puede tener un máximo en ningún punto interior al 
dominio, es decir, las temperaturas máximas y mínimas se dan en el borde del dominio 
3-Condición de compatibilidad para la existencia de soluciones: En el caso de ecuación de 
Laplace con flujo prescrpto en la frontera, el problema podría no tener solución. En efecto, si u es 
solución de la ecuación de Laplace en Ω, entonces 	∇�u = 0 en toda la región Ω e integrando 
obtenemos 
0 = ∬ 	∇�u	dx	dy =	Ω ∬ ∇. �∇u�	dx	dy =ÀÁÂÃ.ÄÅÂÂƹ ∇u.ÇΩ nÈΩ w®	,			donde ∂Ω es la curva cerrada que 
rodea a	Ω y nÈ	es el vector normal exterior a ∂Ω . Esto nos dice que si u es solución de la ecuación 
de Laplace el flujo de calor neto a través de la frontera debe ser cero. Esto también está claro 
desde el punto de vista físico, puesto que en caso contrario habría cambio en el tiempo de la 
energía térmica total dentro de la región, violando la condición de estado estacionario. 
Por lo tanto la condición ¹ ∇u.ÇΩ nÈds = 	0 se llama condición de compatibilidad para la Laplace. 
Ejercicio propuesto: 2.5.2 (a) y (b) de pág 85 
--------------------------------------------------------------------------------------------------------------------------- 
Material elaborado por Adriana Frausin (2014) 
BIBLIOGRAFÍA: Haberman,R. “Applied Partial Differential Equations” Prentice-Hall (2004)-Capítulo 1 
Zill y Cullen “Ecuaciones diferenciales con aplicaciones de valores en la frontera” Thomson (2006) 
 
 
 
 
Ecuación de onda: 
Un problema físico importante y
que modela las vibraciones de una cuerda 
Vamos a deducir la EDP que gobierna las pequeñas vibraciones transversales de una cuerda 
elástica, la cual se estira hasta la longitud 
la cuerda se deforma transversalmente y después de cierto insta
deja vibrar. El problema es determinar las vibraciones de la cuerda, es decir su deformación 
cualquier punto x y en cualquier instante
Las suposiciones simplificadoras
demasiado complicada son: 
1- La cuerda es homogénea (la masa de la cuerda por unidad de longitud
perfectamente flexible (no ofrece resistencia alguna a la flexión y las tensiones serán tangenciales 
y constantes en todo punto). 
2- La tensión causada al estirar la cuerda antes de fijarla en los ptos extremos es tan grande que la 
acción de la fuerza gravitacional sobre la cuerda puede despreciarse
actúa sobre la cuerda. 
3- La deformación y la pendiente de la cuerda son en todo punto pequeñas en valor absoluto, lo 
cual nos permite despreciar la componente horizontal del movimiento, 
cuerda se mueve estrictamente en sentido vertical 
Denotando con � a la magnitud de la tensión 
pequeña porción de cuerda de longitud 
 
 u(x,t) 
 
 
Ecuación de onda: el modelo 
importante y diferente a los anteriores, que involucra otra clase de EDP es el 
que modela las vibraciones de una cuerda o membrana elástica. 
Vamos a deducir la EDP que gobierna las pequeñas vibraciones transversales de una cuerda 
elástica, la cual se estira hasta la longitud L y, a continuación se fija en los puntos extremos. Luego 
la cuerda se deforma transversalmente y después de cierto instante, digamos t= 0, se suelta y se 
deja vibrar. El problema es determinar las vibraciones de la cuerda, es decir su deformación 
y en cualquier instante t > 0. 
s que consideraremos para que la ecuación resultante no quede 
la masa de la cuerda por unidad de longitud, ρ es constante
(no ofrece resistencia alguna a la flexión y las tensiones serán tangenciales 
a tensión causada al estirar la cuerda antes de fijarla en los ptos extremos es tan grande que la 
acción de la fuerza gravitacional sobre la cuerda puede despreciarse y, ninguna otra fuerza externa 
endiente de la cuerda son en todo punto pequeñas en valor absoluto, lo 
cual nos permite despreciar la componente horizontal del movimiento, así cada partícula de la 
cuerda se mueve estrictamente en sentido vertical u = u(x,t) (llamada función de onda
a la magnitud de la tensión en cada punto, ilustramos la situación en una 
pequeña porción de cuerda de longitud ∆s en el intervalo [x, x + ∆x] . 
 cuerda 
Cuerda en equilibrio
, que involucra otra clase de EDP es el 
Vamos a deducir la EDP que gobierna las pequeñas vibraciones transversales de una cuerda 
y, a continuación se fija en los puntos extremos. Luego 
nte, digamos t= 0, se suelta y se 
deja vibrar. El problema es determinar las vibraciones de la cuerda, es decir su deformación u en 
resultante no quede 
es constante ) y 
(no ofrece resistencia alguna a la flexión y las tensiones serán tangenciales 
a tensión causada al estirar la cuerda antes de fijarla en los ptos extremos es tan grande que la 
ninguna otra fuerza externa 
endiente de la cuerda son en todo punto pequeñas en valor absoluto, lo 
cada partícula de la 
función de onda). 
en cada punto, ilustramos la situación en una 
 perturbada 
Cuerda en equilibrio 
 
Aplicando la segunda ley de Newton (m.a = F ) para las componentes verticales de las fuerzas que 
actúan sobre ese pequeño segmento de la cuerda, adimensionalizando el área de la sección 
transversal (A = 1) y usando 	∆�	 ≈ ∆� (por suposición3-) resulta: 
�	∆�	 �	
��	 ��, �� = 	�	�����´	� − �	�����	� 
Además, para � y �´ pequeños esta fuerza neta vertical que actúa en el segmento de cuerda ∆� 
puede expresarse en términos de las correspondientes funciones tangentes, con lo cual 
�	∆�	 �	
��	 ��, �� = 	�	������´	� − �����	�� = 	�	�	
��� + ∆�, �� − 
���, ���	 
Ahora dividiendo ambos miembros por ∆� y tomando el límite para ∆� → 0	, obtenemos la 
ecuación unidimensional de onda, donde 			�	 =	 �� es el cuadrado de la velocidad de 
propagación de la onda (rapidez con la que se propaga la perturbación por la cuerda) 
	�	
��	 =	�	 	�	
��	 
Observaciones: 
1- Esta ecuación no sólo gobierna las vibraciones de una cuerda perturbada (onda elástica) sino 
también las ondas mecánicas (sonido) y las ondas electromagnéticas (radio, televisión, internet) 
2- Conociendo la función u, podemos calcular el desplazamiento, la velocidad ut con que se 
mueven las partículas cuando son perturbadas por la onda (no confundir con c), su aceleración utt 
en cualquier instante de tiempo y la forma de la cuerda ( su pendiente= ux y su curvatura = uxx) 
3- Dado que la ecuación de onda tiene derivada de segundo orden para x, se deben imponer dos 
condiciones de borde (para el caso de extremos fijos son u(0,t) = 0 = u (L, t)) y como también hay 
derivada de segundo orden para t, son dos las condiciones iniciales necesarias (posición inicial 
u(x,0) y velocidad inicial ut(x,0)) 
4- Dado que la EDP de onda el lineal y homogénea, cuando las condiciones de borde también lo 
sean será posible obtener la solución del (PVB) aplicando la técnica de separación de variables. 
Más aún, si las condiciones iniciales fueran combinación lineal de las autofunciones del problema, 
obtendremos solución exacta para el problema. Por ejemplo: 
 
 
 
Separación de variables para un (PVB) con ecuación de onda 
Si la cuerda es de longitud L = π y es liberada desde el reposo dado por f(x) = sen(x), las 
ecuaciones que modelan la situación son: 
����� ∶ 	 !"#�	 ∶ 				 
$$ = 	�	. 
��				; 								0	 < 	� < 		π		; 			t > 0															+�1,2:			
�0, �� = 0			; 									
�	π	, �� = 0		; 					�	 ≥ 0																			+01,2:			
�	�, 0� = ������			; 							
$�	�, 0� = 0		; 		0	 ≤ �	 ≤ 2 3 
Paso 1) Proponer 			
��, �� = 	∅���5���	 reemplazando en la EDP y separando variables, 
5 ´´����	. 5��� = 	∅´´���∅��� = 	−7 ⇒ 9∅´´��� = −7∅���					5 ´´��	� = −7�	. 5���3 
Paso 2) Resolver el problema para ∅�:� con las cond de borde homogéneas dadas por CB1,2 
!∅´´��� = −7. ∅���∅�0� = 0																∅�π� = 0																3 . Ya resolvimos este problema para la ecuación del calor y obtuvimos que, 
∅��� = ;. �<�	����		�=�	�	 = 	1,2, … son las autofunciones no triviales correspondientes a los 
autovalores positivos ? = @A (dado que L = π� 
Paso 3) Resolver la EDO con 5��� para los autovalores hallados (o mejor, el PVI con la cond inicial 
que proporciona la CI homogénea del problema original): 
(PVI) : B5 ´´��	� = −�	�	. 5���5´�0� = 0																					 3 
Entonces 		5��� = � cos��. � . �� + ; sin��. �. ��			 ⇒ 		5´(0) = b. cos(n.c.0) = b = 0. 
Por lo tanto 5��� = � cos��. � . ��		�=�		� = 1,2, … . 
Paso 4) Aplicar el principio de superposición (extendido) a las soluciones producto halladas: 
��, �� = 	 H IJ cos��. � . ��. sin��. ��KJLM 
Paso 5) Determinar los coef. An a partir de la aplicación de la CI no homogénea del (PVB): 
�	�, 0� = 	 H IJ. sen��. ��KJLM = ������ OPQP		JLMRSSSSSST	IM = 1		U		IJ = 	0		V�W�	� = 2,3, …		 
Por lo tanto 		
��, �� = ������ cos��. �� es la solución exacta para este (PVB) 
Material elaborado por Adriana Frausin (2014) 
 
 
Método de las Diferencias Finitas 
Introducción: Hemos estudiado los tres tipos más simples de ecuaciones diferenciales parciales 
(EDP) como representantes de cada una de las tres clases en las que se clasifican las EDP de 
segundo orden: parabólica (calor), elíptica (Laplace) e hiperbólica (onda). Esta terminología 
proviene de “comparaciones” con las ecuaciones de las secciones cónicas (Weinberger, 1965). 
Hasta ahora, hemos aplicado el método de separación de variables (MSV) para obtener soluciones 
explícitas de algunos problemas de valores en el borde (PVB) de interés físico. Contar con una 
expresión analítica para la solución de un problema tiene la ventaja de permitir un estudio 
cualitativo de su solución, como la estabilidad o decaimiento de la misma cuando el tiempo crece 
(ecuación evolutiva), o el estudio de la regularidad de las soluciones para distintas condiciones de 
borde y/o iniciales. 
Sin embargo, si bien en muchos casos el MSV permite encontrar “soluciones exactas”, vimos que 
en general, lo que se obtiene es una expresión de la solución que involucra una serie infinita y a la 
hora de calcular efectivamente el valor de la solución en un punto del dominio es necesario 
truncarla y aproximar el valor de la suma infinita por una suma finita. Otra desventaja del MSV es 
que los dominios están restringidos a intervalos, rectángulos, círculos y regiones simples. 
Con el propósito de obtener soluciones numéricas para los casos en que el MSV resulta 
insatisfactorio, surgieron los métodos numéricos. Los más usados son, el Método de las 
Diferencias Finitas y el Método de los Elementos Finitos, ambos proporcionan soluciones 
numéricas aproximadas de problemas con EDP o EDO. 
Comenzaremos con el estudio de la técnica del método de las diferencias finitas (MDF) para 
discretizar un problema continuo (PVB) y obtener valores numéricos de la variable de estado en 
algunos puntos del dominio: 
Idea del MDF → Aproximar las derivadas por cocientes incrementales o 
diferencias finitas 
Para aplicar esta idea, en problemas con EDP o EDO y condiciones de borde que involucran 
derivadas, se usan las Series de Taylor que permiten aproximar estas derivadas de varias formas. 
���� = ����� + 	�´����	ℎ + 	�´´���� ��
! + 	�´´´���� ���! + 	 … +	 �������� ���! + 	 ��	 , donde el error 
de truncamiento �� = ��������� ���������! con �� < � < � = 	 �� + ℎ. 
 
 
Aproximaciones de la primera derivada
 Para aproximar derivadas de primer orden tenemos tres opciones
pequeño, aproximar la pendiente de la recta tangente a la curva y=f(x) en el punto (x, f(x))
pendiente de la secante que pasa por dicho punto y por (x+h, 
la secante que pasa por éstos dos últimos. Cada una de 
diferencia adelantada, diferencia atrasada 
diferentes algoritmos. 
Diferencia adelantada 	�´����
f(x), que el error de truncamiento en esta aproximación lineal es del orden de h y lo denotamos
Error = O(h). 
Diferencia atrasada: 	�´���� 	
 
Aproximaciones de la primera derivada: 
Para aproximar derivadas de primer orden tenemos tres opciones según se escoja, para h 
pequeño, aproximar la pendiente de la recta tangente a la curva y=f(x) en el punto (x, f(x))
pendiente de la secante que pasa por dicho punto y por (x+h, f(x+h)) o por (x-h, f(x
la secante que pasa por éstos dos últimos. Cada una de estas aproximaciones recibe 
diferencia adelantada, diferencia atrasada y diferencia centrada respectivamente, y conducen a 
 
� � 	 � �������	�	������ es fácil ver, a partir de la serie de Taylor
f(x), que el error de truncamiento en esta aproximación lineal es del orden de h y lo denotamos
 
� � � �����	�	�������� con Error = O(h) 
según se escoja, para h 
pequeño, aproximar la pendiente de la recta tangente a la curva y=f(x) en el punto (x, f(x)) por la 
-h)) o, por la de 
recibe el nombre de 
respectivamente, y conducen a 
serie de Taylor para 
f(x), que el error de truncamiento en esta aproximación lineal es del orden de h y lo denotamos, 
 
Diferencia centrada: Una mejor aproximación se obtiene promediando las anteriores 
 	�´���� 	 � �������	�	�������
� con Error = O(h2) . 
Ejercicio propuesto: Demostrar que el errorque se comete al aproximar la primer derivada �´����,		 con una diferencia centrada es O(h2) . Sugerencia: Restar las series de Taylor para 
f(�� + ℎ�	 y f(�� − ℎ� 
------------------------------ 
Aproximación de la segunda derivada: 
Como �´´���� 	 = "	�´����#´		 � �´����	–	�´������� 	 � %�&��'�	(	%�&��' 	�		%�&��	–	%�&�('�' 		� 	 , por lo tanto 
�´´���� � 	 ���� + ℎ� 	 − 2����� 	 + 	���� − ℎ�ℎ
 
Esta se llama diferencia centrada de segundo orden y se puede demostrar que el Error = O(h
2
). 
Sugerencia para la demostración: Sumar las series de Taylor para f(�� + ℎ�	 y f(�� − ℎ� 
-------------------------------- 
Malla de puntos en el plano y la aproximación de las derivadas parciales 
Todas las aproximaciones anteriores involucran puntos separados una distancia ∆� = ℎ		para la 
variable espacial �. Esto da lugar a una partición del intervalo [0, L] que denotaremos 0 = x0 , x1, 
x2, … , xN = L , donde xj = j. ∆x y es en estos puntos donde obtendremos el valor numérico de la 
solución de un problema con dominio [0, L] . 
Cuando la variable de estado u sea función de dos variables, por ejemplo u = u(x,t), al aproximar 
las derivadas con respecto a t , se tendrá una partición para el intervalo 0 < t < T : 0 = t0 , t1, t2, … , 
tM = T , donde tm = m. ∆, . De esta manera resulta una malla o reticulado de puntos (xj , tm) en el 
plano xt definida por los “pasos” ∆x, ∆t para espacio y tiempo respectivamente. 
Gráficamente, si por ejemplo dividimos los intervalos [0, L] y [0, T] en 10 subintervalos de 
longitud ∆x y ∆t , o sea L = 10∆x y T = 10∆t con lo cual N = 10 en el eje x (horizontal) y M = 
10 en el eje t (vertical) , tenemos la siguiente malla de 11 x 11 = 121 puntos: 
 
Donde el primer punto rojo es (x2, t2) = (2∆x, 2 ∆t) y el último (x9, t10) = (9∆x, 10∆t) 
Observación: En general dados dos enteros positivos N y M (no necesariamente iguales) el 
dominio continuo - = 	 .��, ,� ∶ 		0	 ≤ �	 ≤ 2	,			0	 ≤ ,	 ≤ 	34 se subdivide en NxM rectángulos ∆�	x ∆, y la malla puntos (dominio discreto) tendrá (N+1)x(M+1) puntos. 
Ahora para formular el problema discreto (p) correspondiente al siguiente problema continuo (P): 
 �5�	6 7-5:				9: = 	;. 9��				,						<=>=			0	 < 	� < 	2?@:		9�0, ,� = A�,�	; 	9�2, ,� = @�,�	; 	,	 ≥ 0?D:			9�	�, 0� = 	����		,									<=>=		0	 ≤ �	 ≤ 2E convenimos en denotar con ujm al valor 
numérico de la temperatura en cualquier punto de la malla (xj,tm) . Es decir, 
9"	�F	, ,G# ↔ 	 9FG	, <=>=				I	 = 0, 1, 2, … , K			L			M	 = 0, 1, 2, … , N 	
Si implementamos la notación anterior al aproximar la derivada parcial respecto de t, 9:�	�F	, ,G�, 
por una diferencia adelantada y a 9��	"	�F	, ,G#	 por la diferencia centrada de orden dos, tenemos 
9:	"	�F	, ,G# 	 � 	 9"�F		, ,G 	 + ∆,		# 	 − 	9"�F 	, ,G		#	∆, = 	 9FG�� 	 − 		 9FG			∆, 
 
9��	"	�F	, ,G# 	 � 	 9"�F + ∆�, ,G# − 29"�F	, ,G		# + 9��F − ∆�	, ,G��∆��
 	 = 9F��	G 	 − 29FG + 9F��G�∆��
 
Reemplazando en la 7-5 de �5� obtenemos una ecuación algebraica que llamamos edf 
(ecuación en diferencias finitas) para el problema discreto (p) 
OP�	 ∶ 	 9FG�� 	 − 		 9FG			∆, 	 = ;	 9F��	G 	 − 29FG + 9F��G�∆��
 
 
Despejando 	9FG�� y denotando con Q = ; ∆:�∆��� resulta que el problema discreto (p) queda así: 
�<�	ROP�:			9FG�� = Q. 9F��	G +	 �1 − 2Q�. 9FG + 	Q. 9F��G 		; 			I	 = 1, 2, … , K − 1			L			M	 = 0, 1, 2, … , N − 1ST:		9�G = A�m. ∆,		�			; 				9VG = @�m. ∆,		�	; 			<=>=			M	 = 1, 2, … , N																																																				SW:			9F� = 	��j. ∆��		; 										<=>=		I	 = 0, 1, 2, … , K																																																																																							 E 
 
Observación: Notemos que para calcular 		9FG��, valor numérico aproximado de la temperatura en 
el punto "	�F	, ,G��#, se utilizan los tres valores de u del instante anterior ,G a izquierda y derecha 
de 	�F	 y en 	�F	. Gráficamente, el esquema o átomo computacional para la ecuación del calor es 
										9FG��				 
 s 1- 2s s 
 													9F��G 																9FG 				9F��G 
 
¿Cómo saber qué tan buenas son las aproximaciones obtenidas para este tipo de (PVB)? 
Estabilidad y Convergencia del método: 
Denotemos con YFG = Z	9FG − 	9"	�F	, ,G#Z al error local de discretización en cada nodo de la 
malla (	9FG es el valor aproximado y 9"	�F	, ,G#		el valor exacto de la temperatura. 
Se dice que el método de diferencias finitas usado para calcular las aproximaciones 	9FG 
CONVERGE, si el max	]YFG^ → 0			cuando ∆�	 → 0		�L		∆, = Q. �∆��
 	 → 0	�; donde el máximo se 
toma sobre todos los puntos (j , m) de la grilla (malla). 
 
Puede demostrarse que para asegurar la convergencia y estabilidad (es decir, que el error 
disminuya conforme avanza el cálculo) es SUFICIENTE con que ` = a ∆b�∆c�d ≤ ed 		.		 
 
Advertencia: Esto significa que, cuando s > ½ el método puede (aunque NO necesariamente) 
volverse inestable (es decir oscilante en espacio y tiempo). VER en Ejemplo 1 (c) de página sigte. 
 
Observación: Dado un valor para ∆� , la condición s ≤ ed		 permite escoger ∆, de manera que “todo 
ande bien”. VER en Ejemplo 1 (b) de página siguiente. 
 
 
 
Teorema de estimación del error (ecuación del calor con temperatura prescripta en los bordes) 
Hipótesis: 
Sea (P) el problema continuo que modela la evolución de la temperatura u(x; t) en un “cilindro 
unidimensional” de longitud L; homogéneo, con difusividad k; temperaturas prescritas en los 
bordes dadas por A(t) y B(t) y temperatura inicial f(x). Con A, B y f “suficientemente” suaves. 
Asumimos que se cumplen las condiciones laterales, A(0) = f(0) y B(0) = f(L); de manera que (P) 
tiene una solución u(x; t). 
 
Sea		9FG	 la solución numérica aproximada de este problema por el método explícito de diferencias 
finitas en el punto (xj , tm) (o sea la solución de (p) en el nodo j,m de la malla) y Q = ; ∆:�∆��� ≤ �
 . 
 
Tesis: 
 
El error local de discretización en cada punto satisface para j = 0, 1, 2, … ,N y m = 0, 1, 2, … , M 
 
fgh 	 ≤ 	 i. jd 	k∆b + 	 �∆c�dl m 					nopnq				i = 	 rstu	v	c	v	wu	v	b	v	j]|y´´�b�|, |z´´�b�|, Z{|}�c�Z^ 
 
 
Es decir, si s = k.T.N2/(L2M) ≤ ½ , el error local de discretización tiende a cero cuando M , N → ∞; 
donde k es el coeficiente de difusividad de la ecuación del calor del problema original (P). 
 
Advertencia: no confundir el coef de difusividad k (que está en s) con el K del máximo para fgh. 
 
Observaciones: 
-La aproximación numérica converge en forma cuadrática para ∆�	 → 0		y lineal para ∆,	 → 0	. 
-La convergencia empeora en el tiempo, ya que la cota para el error depende de T (no de L). 
 
 
 
Ejemplo 1: Consideremos el sigte problema de fuga de calor (se pierde energía en los bordes) 
 
�5�	67-5:				9: = 	 9��				,						<=>=			0	 < 	� < 	5													?@:		9�0, ,� = 2,	, 9�5, ,� = 25 + 2,	; 	,	 ≥ 0?D:			9�	�, 0� = 	 �
		,									<=>=		0	 ≤ �	 ≤ 5										 E 
 
(a) Verificar que u(x,t) = x
2
 +2t es solución la solución exacta del problema (P) 
 
(b) Para una malla con ∆� = 1 elegir un valor “conveniente” de ∆, para calcular, usando el 
método explícito de diferencias finitas, u(1, 0.1) y u(4, 0.2). Comparar con los valores exactos. 
 
(c) Usar ∆� = 1 y ∆, = 2 para calcular u(1, 2) , u(2,2), u(3, 2) y u(4,2). Comparar con los valores 
exactos. Qué observas? Comentar. 
Solución (a) 
 
u(x,t) = x2 +2t es la solución exacta ya que satisface TODAS las ecuaciones de (P) : 
 7-5:				9: = 	2 = 	 9��	 
 ?@:		9�0, ,� = 0
 + 2,	, 9�5, ,� = 5
 + 2, 
 ?D:			9�	�, 0� = 	 �
 + 2�0� 
 
Solución (b) 
 
 
 
 0.2 O 
 
 O OO O 
 				∆, = 0.1 
 O O O O O O 
 ∆� = 1 2 3 4 5 = L 
 
Como 		∆� = 1 y 		; = 1	(difusividad) para asegurar estabilidad: Q = 1 ∆:���� 	 ≤ 	 �
 		 ⇒ 		 ∆,	 ≤ 	0.5	 
“alcanza para que todo ande bien”. Entonces teniendo en cuenta que se pide temperatura cuando , = 0.1 y , = 0.2 , escogemos ∆, = 0.1	 y luego calculamos todos los 		9FG	necesarios para 
obtener 		9�� (aproximación de u(1, 0.1)) y			9�
 (aproximación de u(4, 0.2) ) con Q = 1 �.����� = 0.1	. 
 
El esquema para este problema muestra que necesitamos calcular las 6 aproximaciones “para 
todos los nodos del eje x”, para esto usamos la ?D:			9�	�, 0� = 	 �
 y calculamos 		9F� para j = 0, 
1, 2, …, 5 (Notar que N coincide con L=5 pues ∆c = e	y por lo tanto g. ∆c = g		en este problema): 
 		9�� = 0 ; 		9�� = 1
		 = 1		; 		 		9
� = 2
		 = 4 ; 		9�� = 3
		 = 9 ; 		9�� = 4	
		 = 16		; 		 		9�� = 5
		 = 25 
 
Ahora		9�� 	 = 0.1		9�� 	 + "	1 − 2	�0.1�#	9�� 	 + 	0.1		9
� = 0.1�0� + 	0.8�1� + 	0.1	�4� = e. d = 		�ee 
 
Análogamente, 
 		9�� = 0.1	9
� + 0.8 	9�� +0.1 		9�� 	 = 0.1�4� + 		0.8	�9� + 	0.1	�16� = 		9.2 
 		9�� = 0.1	9�� + 0.8 	9�� +0.1 		9�� 	 = 0.1�9� + 		0.8	�16� + 	0.1	�25� = 		16.2 
 		9�� = 25 + 2	�0.1� = 25.2		 (pues es condición de borde derecho) 
 
Finalmente, con los tres valores anteriores, calculamos: 
 		9�
 = 0.1	9�� + 0.8 	9�� +0.1 		9�� 	 = 0.1�9.2� + 		0.8	�16.2� + 	0.1	�25.2� = 		 el. � = 		��d 
 
Calculamos los valores exactos usando la solución dada u(x,t) = x2 +2t : 9�1, 0.1� = 1
 + 2	�0.1� = 1.2 = 		�ee	! 9�4, 0.2� = 4
 + 2�0.1� = 16.2 = 	��d	! 
 
Notamos que los valores aproximados coinciden con los exactos dado que la cota que proporciona 
el teorema de estimación para los errores locales es cero pues 
 
 i = rstu	v	c	v	�u	v	b	v	u.d ]|y´´�b�|, |z´´�b�|, Z{|}�c�Z^ = 0 ya que |y´´�b�| = u	, |z´´�b�| = u	, Z{|}�c�Z = u. 
 
 
Solución (c) Acá nos piden calcular, con ∆� =1 y ∆, = 2, u(1, 2) , u(2,2) , u(3, 2) y u(4,2). 
 
Entonces tenemos 	Q = 1 
���� = 2 > �
 	 ! y sabemos que el método puede volverse inestable. 
 
Los cuatro valores pedidos corresponden a los cuatro nodos interiores del primer renglón de la 
malla definida por el “nuevo paso” ∆, = 2	, gráficamente 
 
 
 
 ∆, =	2 		9�� 		9
� 		9��														9��	 
 
 			 
 O O O O O O 
 ∆� = 1 2 3 4 5 = L 
 		9�� 	 = 2	9�� 	 + "	1 − 2	�2�#	9�� 	 + 	 	2	9
� = 2�0� + �−3��1� + 	2	�4� = � = 		�ee 
 		9
� 	 = 2	9�� 	 + "	1 − 2	�2�#	9
� 	 + 	 	2	9�� = 2�1� + �−3��4� + 	2	�9� = � = 		�de 
 		9�� 	 = 2	9
� 	 + "	1 − 2	�2�#	9�� 	 + 	 	2	9�� = 2�4� + �−3��9� + 	2	�16� = e� = ��e 
 		9�� 	 = 2	9�� 	 + "	1 − 2	�2�#	9�� 	 + 	 	2	9�� = 2�9� + �−3��16� + 	2	�25� = du = ��e 
 
Los valores exactos usando la solución dada u(x,t) = x2 +2t : 
 9�1, 2� = 1
 + 2	�2� = 5 = 		�ee	! 9�2, 2� = 2
 + 2	�2� = 8 = 	�de	! 9�3, 2� = 3
 + 2	�2� = 13 = ��e	! 9�4, 2� = 4
 + 2	�2� = 20 = 	��e	! 
 
Comentario: los valores exactos también coinciden con los “aproximados”. Este ejemplo muestra 
que la condición 	` ≤ 	 �
		es suficiente pero NO NECESARIA (esto ocurrirá siempre que las 
condiciones de borde sean funciones a lo sumo lineales y la condición inicial a lo sumo cúbica). 
 
 
 
Ejemplo 2: ¿Cómo modificarías el método para resolver el siguiente problema? 
 
�5�	67-5:				9: = 	;. 9��				,						<=>=			0	 ≤ 	� ≤ 	2	,				0 ≤ ,	 ≤ 3	?@:		9�0, ,� = A�,�	; 	 		9��2, ,� = 0	; 		,	 ≥ 0																												?D:			9�	�, 0� = 	����		,									<=>=		0	 ≤ �	 ≤ 2																									 E 
 
Solución: Notamos que ahora tenemos una condición de borde aislado en el extremo derecho x = 
L, las restantes ecuaciones son las mismas. Por lo cual sólo hay que modificar la ecuación para ese 
borde. Como en x = L hay una derivada, 9��2, ,� = 0	, debemos aproximarla por una diferencia 
finita (adelantada, atrasada o centrada). Aunque la “mejor” desde el punto de vista del error es la 
centrada, vamos a escoger usar diferencia atrasada para no “salirnos” del dominio [0,L] x [0, T], 
entonces 9��2, ,� = 9���V, ,G� = 	 ���	�		��(�� 			∆� = 	0			 ⇒ 			 9VG = 		 9V��G para m= 1, 2, …, M. 
 
Por lo tanto el problema discreto en este caso es 
 
�<�	R OP�:			9FG�� = Q. 9F��	G + 	 �1 − 2Q�. 9FG + 	Q. 9F��G 		; 			I	 = 1, 2, … , K − 1			L			M	 = 0, 1, 2, … , N − 1ST:		9�G = A�m. ∆,		�			; 				9VG = 		 9V��G 													; 			<=>=			M	 = 1, 2, … , N																																																				SW:			9F� = 	��j. ∆��		; 										<=>=		I	 = 0, 1, 2, … , K																																																																																							 E 
 
Propuesta: Aplicando lo recién deducido, resolver el Problema 1 (c), (d), (e) y Problema 4 (c) , (d) 
de la Guía de Problemas complementarios 
 
------------------------------------------------ 
 
Ejemplo 3: Considerar el siguiente problema para una cuerda vibrante de longitud π, fija en sus 
 extremos y liberada desde el reposo dado por f(x) = sen(x), 
 
�5� ∶ 	 67-5	 ∶ 				 9:: = 	 9��				; 								0	 < 	� < 		π		; 			t > 0																														?@:			9�0, ,� = 0			; 																			9�	π	, ,� = 0		; 					,	 ≥ 0																						?D:			9�	�, 0� = QO����			; 							9:�	�, 0� = 0		; 		0	 ≤ �	 ≤ 2															 E 
(a) Aplicar la técnica del método de diferencias finitas para formular el problema discreto para �5� 
 
(b) Implementar el esquema anterior con ∆� = �/4 y ∆, = 0.5 para aproximar u(�/4 , 1). 
 
(c) Comparar con el valor exacto proporcionado por la solución u(x, t) = (sen x)(cos t). Comentar. 
 
 
Solución (a) 
 
Notamos que la 7-5 de este problema de ondas tiene 2 derivadas parciales de 2do orden, por lo 
cual para aproximarlas usamos las correspondientes diferencias centradas de 2do orden. Así la 
edf (ecuación en diferencias finitas) para el problema discreto (p) queda 
OP�	 ∶ 	 9FG�� 	 − 29FG + 9FG���∆,�
 	 = 	 9F��	G 	 − 29FG + 9F��G�∆��
 
Despejando 	9FG�� y denotando con Q = �∆:���∆��� (en general, es Q = S
 �∆:���∆��� ) resulta que, 
OP�	 ∶	 9FG�� 	 =	 Q. 9F��	G + 2�1 − Q�	9FG + Q. 9F��G −		9FG�� 
Observación: Notemos que para calcular 		9FG��, valor numérico aproximado de la posición en el 
punto "	�F	, ,G��# se utilizan 4 “valores previos” (los tres valores “centrados” en el instante 
anterior ,G 	y el valor en "	�F	, ,G��# . Gráficamente, el esquema o átomo computacional para la 
ecuación de onda es 
										9FG��				 
 s 2(1-s) s 
 													9F��G 																	9FG 				9F��G 
 -1 
 
 		9FG�� 
 
Volviendo a las ecuaciones del problema continuo (P), tenemos otra derivada por aproximar, la 
velocidad inicial, lo cual hacemos con un diferencia adelantada (cuando t = 0 es m = 0) y 
tenemos, 9:�	�, 0� = 	 ���	�		���			∆: = 	0		 ⇒ 		 9F� = 		 9F� entonces 
 
�<�	R OP�:			9FG�� = Q. 9F��	G + 2�1 − Q�	9FG + Q. 9F��G 		−		9FG��		; 			I = 1, … , K − 1		L			M	 = 1, … , N − 1																									ST:			9�G= 0			; 				9VG = 0	; 			<=>=			M	 = 0, 1, 2	, … , N							"S9O>P=	�WI=	O�	��Q	O�,>OM�Q"																																								SW:			9F� = 	QO��j. ∆��		; 							9F� = 		 9F�			<=>=		I	 = 	1, 2, … , K − 1					"�O��SWP=P	W�WSW=�	SO>�"																																			
 
Finalmente este problema discreto (p) con S
 = 1,	 para L = �, cuando ∆� = �/4 y ∆, = 0.5 (o 
sea con Q = S
 �∆:���∆��� = ��� ) queda así: 
�<�	
���
��OP�:			9FG�� = 4�
 9F��	G + 2 �1 − 4�
 	9FG + 4�
 9F��G 		−		9FG��		; 			I	 = 1, 2, 3		L			M	 = 1, 2																													ST:		9�G = 0			; 				9�G = 0	; 			<=>=			M	 = 1, 2																																																																																																																SW:			9F� = 	QO��j. ∆��		; 								9F� = 		 9F�			<=>=		I	 = 	1, 2, 3																																																																																										
E 
 
Solución (b) 
 
Implementamos este algoritmo para calcular u(�/4 , 1) o sea 9�
 . Haciendo el esquema con 
todos los puntos necesarios notamos 
 
 
 
 
 9�
 
 
 				9�� 9�� 9
� 
 						9��				9��					9
� 
 
�� �
 ��� 	� 
 9�� 	 = 	 9�� = 0 (por ser borde izquierdo) 
 9�� = 	QO� ¡1. ��¢ = 	 √

 				y 9
� = 	QO� ¡2. ��¢ = 1 (por ser posiciones iniciales) 
 9�� = 9�� = √

 			y 9
� = 9
� = 1 (velocidades iniciales aproximadas con diferencia hacia adelante) 
 
		9�
 = ��� 9
	� + 2 ¡1 − ���¢	9�� + ��� 9��		−		9�� = 0.539 = 	9�
 
 
 
Solución (c) 
 
La solución exacta u(x,t) =(sen x)(cos t) puede hallarse aplicando separación de variables. 
 
Calculamos u(�/4 , 1) = sen(�/4�cos	�1� = 	 √

 0.5403 = 0.382 = 	u��/4	, 1� 
 
Comentario: Aunque el método es estable pues se cumple S	 ≤ ∆�∆: (es decir 1 < �/2� el error 
es grande pues usamos diferencia adelantada para aproximar ut con 	∆, = 0.5		(se ha introducido 
un error = O(∆,��. Para obtener una aproximación aceptable es necesario refinar la malla, por 
ejemplo ∆, = 1/100 y ∆� = �/10 (VERIFICAR, resolviendo en computadora) o tomar ∆, = 1/10 y 
usar diferencia centrada para ut (ver en Haberman pág 256) 
 
Nota: La condición de estabilidad establece que, la velocidad de propagación del método 
numérico, 
∆�∆: , debe ser mayor o igual que la velocidad de propagación de la onda, S. 
 
Ejercitación propuesta: en “Problemas Complementarios” 
 
 
 Material elaborado por A. Frausin - Agosto 2014 
 
BIBLIOGRAFIA: 
• Bleecker D., Csordas G. “Basic Partial Differential Equations” - Ed. Internacional - 1996 
• Haberman R. “Applied Partial Differential Equations” - Prentice Hall - 2004 
 
El M étodo de los Elementos Finitos (MEF)
1 Introducci ón
El Método de los Elementos Finitos (MEF) es un “nuevo” método numérico que se suma al tradicional Método
de Diferencias Finitas (MDF) empleado por ingenieros en an´alisis estructural, resistencia de materiales, mecánica
de fluı́dos, electromagnetismo, etc.
Las ideas básicas del método surgieron a partir de los avances en el área del análisis estructural relacionado
con la aviación (Hrenikoff en 1941 para problemas de elasticidad y Courant en 1943 para problemas de torsión).
Después, a mediados de los años 50 aparece Turner desarrollando matrices de rigidez para la solución de
problemas de elasticidad en barras y vigas, entre otros elementos.
Con grandes logros y siguiendo los pasos de Turner, las Corporaciones MacNeal-Schwendler and Computer
Sciences elaboran en la NASA el primer código de importancia para el análisis de elementos finitos, llamado
NASTRAN (usado en la industria aeroespacial y en áreas de laingenierı́a civil).
Pero no fue hasta 1960 cuando Clough utilizó por primera vezel término de elemento finito y en 1967 fue
publicado el primer libro de Elemento Finito por Zienkiewicz y Chung poniendo de manifiesto la versatilidad
del método, sus fuertes bases matemáticas y diversas aplicaciones.
Con los grandes avances tecnológicos que se han logrado en el área de la computación y sobre todo en los
sitemas de diseño asistido por computadoras, ahora es relativamente más fácil la modelización de prototipos,
en los cuales se pueden tener geometrı́as y superficies complicadas e irregulares, con aplicaciones de carga en
forma especı́fica para el estudio preciso de esfuerzos internos y lograr una modelización ajustada a los perfiles
y estructuras que se emplean teniendo en consideración ciertas caracterı́sticas como el cambio de secciones,
estructuras huecas con pared delgada y con caracterı́sticas en secciones transversales muy especı́ficas. Por ejem-
plo, la simulación de deformaciones de vehı́culos por impacto y el análisis de mecanismos de transportación
ósea para reducción de fracturas se pueden citar como dos aplicaciones concretas en áreas bien diferentes.
Actualmente se cuenta con una gran cantidad de paquetes computacionales que permiten realizar cálculos
con elementos finitos. Sin embargo, el manejo correcto de este tipo de herramientas exige un profundo
conocimiento no sólo del material con el que se trabaja, sino también de los principios en los que se basa
el método. Sólo en este caso se estará en condiciones de garantizar que los resultados obtenidos en los análisis
se ajustan a la realidad.
2 Fundamentos de los elementos finitos
Uno de los fundamentos del MEF está basado en la discretización de los cuerpos en estudio. De esta manera,
al igual que en el MDF, el método de los elementos finitos require de un espacio geométrico (o dominio) para
ser dividido en subregiones formando unared omalla.
En el MDF, la malla consiste de filas y columnas de ĺıneas ortogonales. En cambio en el MEF cada división
es única y no necesariamente ortogonal, lo cual es una ventaja ante sistemas de geometrı́a irregular.
Cada división o subregión (subintervalos enIR, triángulos o cuadriláteros enIR2, tetraedros o he xaedros
enIR3 ) se llamaelemento, y en cada uno de éstos se resuelve la ecuación diferencial,para luego generar la
solución totalensamblando las soluciones individuales.
La habilidad del MEF para representar con una colección de elementos finitos a dominios de
geometrı́a irregular, es lo que hace de este método una herramienta de mucho valor para resolver problemas a
valores en el borde (PVB) en diversos campos de la ingenierı́a.
Por otra parte, ası́ como el MDF consiste en discretizaciones convencionales, donde las derivadas se aprox-
iman por cocientes en diferencias (cocientes incrementales), el MEF basa su idea en la aproximación de la
función solución. El factor esencial es que la integral deuna función se puede escribir como suma de integrales
en dominios disjuntos cuya unión es el dominio original, entonces se puede hacer un análisis local del problema
y haciendo divisiones del dominio en dominios suficientemente pequeños se pueden elegir funciones sencillas
que sean adecuadas para una buena representación del comportamiento de la solución.
Veamos con un ejemplo unidimensional sencillo la implementación de las ideas expuestas.
3 El MEF para problemas unidimensionales
Consideremos el siguiente problema a valores en el borde (P)que modela los desplazamientos axiales de una
barra de longitud unitaria fija en los extremos y sometida a una carga tangencial f, o la temperatura en estado
estacionario en un cilindro unidimensional con fuente interna de calor f y temperatura cero en los extremos,
o potencial eléctrico en cualquier punto entre dos placas paralelas (ubicadas enx = 0 y x = 1) en un medio
homogéneo, con carga f y potencial nulo en ambas placas:
(P )



u′′(x) = −f(x) , 0 < x < 1
u(0) = 0
u(1) = 0
Idea : A partir de una partición del intervalo(0, 1), aproximar la soluciónu = u(x) por una funciónU(x)
que sea una combinación lineal de funciones “adecuadas”,Tj , llamadas funciones tests (o funciones de base).
La aproximación más simple esconsiderar funciones lineales a trozos definidas en[0, 1]. De esta maneraU(x)
es una poligonal que aproxima au(x)
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
1.2
x
y
exacta
aprox.
¿Cómo elegirTj ?
(1o) Se define una partición0 = x0 < x1 < x2 < ... < xn+1 = 1 del intervalo(0, 1) enn+ 1 subintervalos
(o elementos finitos)Ij = (xj−1, xj) de longitudhj = xj − xj−1. El númeroh = max {hj} será la
medida de cuan fina es la partición.
(2o) Se construyen las funcionesTj que permitirán expresar de manera única a cualquier función lineal a
trozos que aproxime a la solución del problema (P) y satisfaga sus condiciones de borde. A partir de
estas consideraciones se definen las siguientesn funciones lineales a trozos (una por cada nodo interior
de la malla) y tales que
Tj(xi) =
{
1 si i = j
0 si i 6= j
Ejercicio 1 : Representar graficamente las funcionesTj para una partición del(0, 1) en4 subintervalos de
igual longitud (partición equiespaciada).
Observacíon: Notemos que lasTj asi definidas son linealmente independientes y por lo tanto constituyen
una base del espacio vectorial que generan,Vh = gen {T1, T2, ..., Tn} . De esta manera cualquier función
U ∈ Vh tiene una representación única de la forma
U(x) =
∑n
j=1 ujTj(x) donde uj = U(xj)
¿Cómo calcular los coeficientesuj ?
Estosn coeficientesuj son precisamente los valores de la función “aproximante” evaluada en losn nodos
interiores.
(1o) Se reformula (P) obteniendo la llamadaforma débil o variacional (Galerkin):
Se multiplica la ecuación diferencial por una función de base arbitrariaTi,
Ti(x)u
′′(x) = −Ti(x)f(x) con Ti(0) = Ti(1) = 0 entonces,
∫
1
0
Ti(x)u
′′(x)dx = −
∫
1
0
Ti(x)f(x)dx , integrando por partes resulta
Ti(x)u
′(x)|10 −
∫
1
0
T ′i (x)u
′(x)dx = −
∫
1
0
Ti(x)f(x)dx , y comoTi satisface las condiciones de borde
del problema (P), el primer término se anula y ası́ se obtiene la
Forma débil o variacional : (D)
∫
1
0
T ′i (x)u
′(x)dx =
∫
1
0
Ti(x)f(x)dx
De esta manera se ha demostrado que siu es solución del problema diferencial (P) entonces,u satisface
la forma débil (D) para todaTi ∈ Vh.
Ejercicio 2: Demostrar que siu′′ y f son continuas también vale el recı́proco, es decir, siu es solución de
(D) para todaTi ∈ Vh entoncesu es solución de (P).
(2o) En (D) se reemplazau por su aproximaciónU obteniendo lamatriz de rigidez y el vector de cargas:
En efecto,
∫
1
0
T ′i (x)


n
∑
j=1
ujTj(x)


′
dx =
∫
1
0
Ti(x)f(x)dx , que se puede escribir en la forma
n
∑
j=1
uj
∫
1
0
T ′i (x)T
′
j(x)dx =
∫
1
0
Ti(x)f(x)dx para i = 1, 2, ...n,
que no es más que un sistema lineal den ecuaciones conn incógnitasu1, u2, ..., un.
En forma matricial ese sistema puede escribirse comoKU = F dondeK = (Kij) es una matriz
n×n con elementosKij =
∫
1
0
T ′i (x)T
′
j(x)dx , F = (Fi) es un vector columnan× 1 con elementos
Fi =
∫
1
0
Ti(x)f(x)dx y U = (u1, u2, ..., un)t es el vector de las incógnitas, que son precisamente
los valores aproximados de la variable de estadou en cada uno de los nodos interiores del dominio(0, 1).
La matrizK es conocida comomatriz de rigidez y F comovector de cargaterminologı́a que proviene
de las primeras aplicaciones del MEF en mecánica estructural.
Observaciones:
• Al ser las funcionesTi lineales a trozos, sus derivadasT ′i no son más que las pendientes de los corre-
spondientes segmentos de recta que las conforman.
• K es una matriz simétrica puesKij = Kji (son los productos de las correspondientes pendientes).
• En virtud de la elección hecha de las funciones de baseTj la matrizK es tridiagonal, es decir,Kij = 0
si |i− j| > 1 (puesTi y Tj tienen soportes disjuntos).
• Puede demostrarse queK es invertible con lo cual el sistema lineal tieneúnica solucíon U = K−1F
Estimación del error:
Si u′′ es continua en[0, 1] una estimación del error viene dada por
|u(x)− U(x)| ≤ h max
0≤x≤1
∣
∣u′′(x)
∣
∣ para todox ∈ [0, 1] ,
que demuestra la convergencia del método porque el error tiende a cero cuando el diámetroh de la partición
tiende a cero (pues al seru′′ continua en(0, 1) resulta acotada).
Ejercicio 3: Ejemplo de PVB con carga o término fuente ”f(x)” CONSTANTE
(a) Aplicar el MEF con una partición equiespaciada de4 elementos y funciones de base lineales, para obtener
la solución aproximada del siguiente problema a valores enel borde.
(P )



u′′(x) = −1 , 0 < x < 1
u(0) = 0
u(1) = 0
(b) Hallar la solución exacta de (P) y comparar con los valores aproximados.
(a) Para calcular los elementos de la matriz de rigidez, es útil representar graficamente las correspondientes
funciones tests lineales,T1, T2, T3 parah = 0.25 (HACERLO). LUEGO se calculan las correspondientes
integrales (NOTAR que el integrando de la fórmula paraKij contiene las derivadas de las funciones tests
T , que por ser lineales no son más que sus correspondientes PENDIENTES):
K11 =
∫
1
0
[
T ′1(x)
]2
dx =
∫
0.25
0
dx
(0.25)2
+
∫
0.5
0.25
dx
(−0.25)2
=
0.5
(0.25)2
= 8
Análogamente,K22 = K33 = 8 ( pues la partición es equiespaciada) .
K12 = K21 =
∫
1
0
T ′1(x)T
′
2(x)dx =
∫
0.5
0.25
[
−1
0.25
] [
1
0.25
]
=
−0.25
(0.25)2
= −4
Análogamente,K23 = K32 = −4 ( pues la partición es equiespaciada) .
K13 = K31 = 0 puesT1 y T3 tienen soportes disjuntos (la matriz es tridiagonal).
Las componentes del vector de carga serán iguales pues el t´ermino fuente es constante,f(x) = 1, y la
partición es equiespaciada, en efecto
parai = 1, 2, 3 Fi =
∫
1
0
(1)Ti dx = 1/4 ( área bajoTi para cualquieri).
Por lo tanto los valores aproximados deu en los nodos interiores,x1 = 0.25, x2 = 0.5, y x3 = 0.75,
son
U =

4


2 −1 0
−1 2 −1
0 −1 2




−1
1
4


1
1
1

 =


3/32
1/8
3/32

 =


u1
u2
u3


(b) Hallamos la solución exacta integrando dos veces
u′′(x) = −1 ⇒ u′(x) = −x+ c ⇒ u(x) = −x2/2 + cx+ d.
Luego por las condiciones de borde se tiened = 0 y c = 1/2 ⇒ u(x) = −x2/2 + x/2.
Notamos que los valores exactosu(1/4) = 3/32, u(1/2) = 1/8, u(3/4) = 3/32, coinciden con los
hallados aplicando el método de los elementos finitos.
A fines comparativos se ilustra en la siguiente figura, la solución anaĺıtica exacta del problema diferencial
(P),u = u(x), conjuntamente con la solución aproximada obtenida con MEF, la poligonalU = U(x),
0 0.2 0.4 0.6 0.8 1
0
5 · 10−2
0.1
x
y
exacta
aprox.
VER más ejercitación enGuı́a de Problemas complementarios(problemas 11 y 12).
4 El MEF para problemas bidimensionales
Se puede ver en la sección 6.7 (pág 267 a 274) de la bibliografı́a básica “Applied Partial Differential Equations”
de Richard Haberman (Cuarta edición2004) el procedimiento, totalmente análogo al caso unidimensional, que
se puede seguir para la construcción de las funciones de base lineales y la deducción de la forma débil, matriz
de rigidez y vector de carga correspondiente a un problema del tipo
(P )
{
∇2u = f(x, y) para(x, y) ∈ R
u(x, y) = 0 para(x, y) ∈ ∂R
donde conR se designa a cualquier región poligonal enIR2 y ∂R denota su frontera (curva cerrada)
BIBLIOGRAFIA:
• Reddy J.N. “An introduction of the Finite Element Method” - Mc Graw Hill - 1984
• Haberman R. “Applied Partial Differential Equations” - Prentice Hall -2004
Elaborado por A. Frausin - 2012 ( modificado 2014)