T-ESPE-048495

Vista previa del material en texto

VICERRECTORADO DE INVESTIGACIÓN, INNOVACIÓN Y
TRANSFERENCIA DE TECNOLOGÍA
CENTRO DE POSGRADOS
Programa: Maestría en la Enseñanza de la Matemática
Promoción Primera
TESIS DE GRADO:
Previa a la Obtención del Título de:
Magíster en la Enseñanza de la Matemática
TEMA
Desarrollo de las Ecuaciones Diferenciales Parciales Parabólicas
mediante Diferencias Finitas, Elementos Finitos y Meshless
Ing. Manuel Patricio Pugarín Díaz
Director: Ph.D. Nelson Subía Cepeda
Quito - Ecuador
2015
I
UNIVERSIDAD DE LAS FUERZAS ARMADAS - ESPE
VICERRECTORADO DE INVESTIGACIÓN, INNOVACIÓN Y
TRANSFERENCIA DE TECNOLOGÍA
CERTIFICACIÓN DEL DIRECTOR
Ph.D. Nelson Edmundo Subía Cepeda
CERTIFICA:
Que el trabajo titulado “Desarrollo de las Ecuaciones Diferenciales Parciales
Parabólicas mediante Diferencias Finitas, Elementos Finitos y Meshless”, rea-
lizado por el maestrante Pugarín Díaz Manuel Patricio, ha sido guiado y re-
visado periódicamente, cumpliendo con las normas establecidas por el De-
partamento de Ciencias Exactas de la Universidad de las Fuerzas Armadas –
ESPE, por tanto, se autoriza su presentación para los fines legales pertinen-
tes.
Sangolquí, marzo 2015
PhD. Nelson Subía C.
DIRECTOR DE TESIS
II
UNIVERSIDAD DE LAS FUERZAS ARMADAS - ESPE
VICERRECTORADO DE INVESTIGACIÓN, INNOVACIÓN Y
TRANSFERENCIA DE TECNOLOGÍA
DECLARACIÓN DE RESPONSABILIDAD
Yo: Pugarín Díaz Manuel Patricio
DECLARO QUE:
La tesis de grado titulada “Desarrollo de las Ecuaciones Diferenciales Parcia-
les Parabólicas mediante Diferencias Finitas, Elementos Finitos y Meshless”,
ha sido desarrollada con base a una profunda investigación, respetando dere-
chos intelectuales de terceros, conforme a las citas correspondientes, cuyas
fuentes constan en la bibliografía. Consecuentemente este trabajo es de mi
autoría.
En virtud de esta declaración, me responsabilizo del contenido, veracidad
y alcance científico de la tesis de grado en mención.
Sangolquí, marzo 2015
Patricio Pugarín D.
C.C: 1708038318
III
UNIVERSIDAD DE LAS FUERZAS ARMADAS - ESPE
VICERRECTORADO DE INVESTIGACIÓN, INNOVACIÓN Y
TRANSFERENCIA DE TECNOLOGÍA
AUTORIZACIÓN
Yo: Pugarín Díaz Manuel Patricio
Autorizo a la Universidad de las Fuerzas Armadas – ESPE, la publicación en
la biblioteca virtual de la institución, de mi trabajo denominado: “Desarrollo de
las Ecuaciones Diferenciales Parciales Parabólicas mediante Diferencias Fini-
tas, Elementos Finitos y Meshless”, cuyo contenido, ideas y criterios son de
mi exclusiva responsabilidad y autoría.
Sangolquí, marzo 2015
Patricio Pugarín D.
C.C: 1708038318
IV
AGRADECIMIENTO
Al Ph.D. Nelson Subía Cepeda y al Ph.D. Paul Medina Vásquez quienes
con sus vastos conocimientos y experiencia me guiaron, confiaron y enseña-
ron de la mejor manera en el desarrollo de este trabajo.
A mis compañeros del Departamento de Ciencias Exactas, de la Univer-
sidad de las Fuerzas Armadas – ESPE, por la paciencia y aliento brindado
durante el desarrollo de la presente investigación.
V
DEDICATORIA
A mi DIOS por estar siempre a mi lado y proveerme de la fuerza, salud y
sabiduría necesaria, a mi madre por su ejemplo de superación y trabajo y a
mi esposa e hijos por su apoyo incondicional; quienes me permitieron ahora
alcanzar este logro tan especial.
Índice general
Certificación del Director I
Declaración de Responsabilidad II
Autorización III
Agradecimiento IV
Dedicatoria V
Índice de Figuras XI
Índice de Cuadros XII
Notación utilizada XIII
Resumen XIV
Abstract XV
Introducción 1
1. Ecuaciones Parabólicas 3
1.1. Ecuaciones Diferenciales Parciales . . . . . . . . . . . . . . . . 3
1.2. Importancia y Selección . . . . . . . . . . . . . . . . . . . . . . . 3
1.3. Ecuaciones del tipo parabólico . . . . . . . . . . . . . . . . . . . 6
1.4. Métodos de solución . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4.1. Diferencias finitas . . . . . . . . . . . . . . . . . . . . . . 7
1.4.2. Método de Elementos Finitos . . . . . . . . . . . . . . . . 8
VI
VII
1.4.3. Método de Meshless . . . . . . . . . . . . . . . . . . . . . 9
1.4.4. Formulación y Objetivo General . . . . . . . . . . . . . . 9
2. Diferencias Finitas 11
2.1. Ecuación Diferencial Parcial Unidimensional . . . . . . . . . . . . 11
2.1.1. Ecuaciones Homogéneas . . . . . . . . . . . . . . . . . . 11
2.1.2. Ecuaciones No Homogeneas . . . . . . . . . . . . . . . . 19
2.2. Otros esquemas numéricos . . . . . . . . . . . . . . . . . . . . . 20
2.2.1. Esquema de Richardson . . . . . . . . . . . . . . . . . . 20
2.2.2. Esquema de Crank - Nicolson . . . . . . . . . . . . . . . 20
2.3. Ecuación Bidimensional . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.1. Análisis de la estabilidad del método . . . . . . . . . . . . 24
2.4. Ejercicios de aplicación . . . . . . . . . . . . . . . . . . . . . . . 26
2.4.1. Problema tipo . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4.2. Problema tipo. Esquema explícito . . . . . . . . . . . . . 27
2.4.3. Problema tipo. Esquema de Crank - Nicolson . . . . . . . 34
2.4.4. Problema 2 (No Homogéneo) . . . . . . . . . . . . . . . . 40
3. Elementos Finitos 50
3.1. El Método de Elementos Finitos. . . . . . . . . . . . . . . . . . . 50
3.1.1. Métodos existentes . . . . . . . . . . . . . . . . . . . . . 50
3.1.2. Utilidad del MEF . . . . . . . . . . . . . . . . . . . . . . . 51
3.1.3. Un poco de historia del MEF . . . . . . . . . . . . . . . . 51
3.1.4. Descripción general del MEF . . . . . . . . . . . . . . . . 52
3.1.5. Descripción matemática del MEF . . . . . . . . . . . . . . 53
3.2. Desarrollo de las EDP’s del tipo parabólico, unidimensional, por
el MEF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.2.1. Formulación débil . . . . . . . . . . . . . . . . . . . . . . 54
3.2.2. Construcción de la formulación débil . . . . . . . . . . . . 55
3.2.3. Formulación de Galerkin . . . . . . . . . . . . . . . . . . 58
3.2.4. Funciones lineales por partes . . . . . . . . . . . . . . . . 61
3.2.5. Ensamblaje de la matriz de elementos finitos . . . . . . . 64
3.3. Ejercicio de aplicación . . . . . . . . . . . . . . . . . . . . . . . . 65
VIII
3.3.1. Problema tipo. Método de elementos finitos . . . . . . . . 65
4. Meshless 79
4.1. El método Meshless . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.1.1. Un poco de historia del MM . . . . . . . . . . . . . . . . . 79
4.1.2. La Función de Base Radial (FBR) . . . . . . . . . . . . . 80
4.1.3. El Método Meshless y las Funciones de Base Radial . . 80
4.2. Desarrollo de una EDP unidimensional . . . . . . . . . . . . . . . 82
4.2.1. Descripción Matemática del Método Meshless . . . . . . 82
4.3. Ejercicio de aplicación . . . . . . . . . . . . . . . . . . . . . . . . 85
4.3.1. Problema tipo. Método de Meshless . . . . . . . . . . . . 85
5. Error de la aproximación numérica y tiempos de cómputo 91
5.1. Cálculo del error . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
5.1.1. Error absoluto (Err) . . . . . . . . . . . . . . . . . . . . . . 91
5.1.2. Error cuadrático medio (Ecm) . . . . . . . . . . . . . . . . 92
5.1.3. Error medio cuadrático relativo (Emcr) . . . . . . . . . . . 92
5.2. Cálculo del error cometido . . . . . . . . . . . . . . . . . . . . . . 93
5.2.1. Método de las diferencias finitas . . . . . . . . . . . . . . 93
5.2.2. Método de los elementos finitos . . . . . . . . . . . . . . 96
5.2.3. Método de Meshless . . . . . . . . . . . . . . . . . . . . . 99
5.3. Análisis del error . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.3.1. Error absoluto . . . . . . . . . . . . . . . . . . . . . . . . 103
5.3.2. Error cuadrático medio . . . . . . . . . . . . . . . . . . . . 104
5.3.3. Error medio cuadrático relativo . . . . . . . . . . . . . . . 105
5.4. Tiempos de cómputo . . . . . . . . . . . . . . . . . . . . . . . . . 106
6. Conclusiones y Recomendaciones 108
6.1. Respuesta a las Interrogantes Planteadas para este trabajo . . . 108
6.2. Conclusiones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
6.3. Recomendaciones . . . . . . . . . . . . . . . . . . . . . . . . . . 111
IX
Tabla de Algoritmos 112
Algoritmo 1, problema tipo, solución exacta . . . . . . . . . . . . . . . 112
Algoritmo 2, problema tipo, esquema explícito . . . . . . . . . . . . . . 112
Algoritmo 3, problema tipo, esquema de Crank-Nicolson . . . . . . . . 114
Algoritmo 4, problema 2 no homogeneo, Crank-Nicolson . . . . . . . . 117
Algoritmo 5, problema tipo, elementos finitos, n = 6 . . . . . . . . . . . 119
Algoritmo 6, problema tipo, elementos finitos, todo n . . . . . . . . . . 121
Algoritmo 7, problema tipo, método Meshless . . . . . . . . . . . . . . 124
Algoritmo 8, error absoluto . . . . . . . . . . . . . . . . . . . . . . . . 127
Algoritmo 9, error cuadrático medio . . . . . . . . . . . . . . . . . . . 129
Algoritmo 10, error medio cuadrático relativo . . . . . . . . . . . . . . 130
Algoritmo 11, polinomio interpolador Lagrange . . . . . . . . . . . . . 130
Algoritmo 12, sistema matricial tridiagonal de ecuaciones . . . . . . . 131
Tabla de Matrices 133
Solución exacta, problema tipo . . . . . . . . . . . . . . . . . . . . . . 133
Solución aproximada, problema tipo, esquema explícito . . . . . . . . 137
Solución aproximada, problema tipo, esquema explícito estable . . . . 142
Solución aproximada, problema tipo, literal a). Crank-Nicolson . . . . 143
Solución aproximada, problema tipo, literal b). Crank-Nicolson . . . . . 148
Solución exacta, problema tipo, literal b) . . . . . . . . . . . . . . . . . 150
Solución aproximada, problema tipo comprobación . . . . . . . . . . . 151
Solución aproximada, problema 2 no homogéneo . . . . . . . . . . . . 156
Solución aproximada, problema tipo, elementos finitos, n = 6 . . . . . 160
Solución exacta, problema tipo, n = 6 . . . . . . . . . . . . . . . . . . 161
Solución aproximada, problema tipo, elementos finitos, n = 19 . . . . 162
Solución aproximada, problema tipo, método Meshless . . . . . . . . 167
Anexo A 172
Problema 3 bidimensional, método de diferencias finitas . . . . . . . . 172
Bibliografía 189
Índice de figuras
2.1. Mallado espacio - tiempo. . . . . . . . . . . . . . . . . . . . . . . 13
2.2. Esquema de diferencias progresivas. . . . . . . . . . . . . . . . . 14
2.3. Esquema del punto (i, j + 1
2
). . . . . . . . . . . . . . . . . . . . . 21
2.4. Esquema Implícito. . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.5. Malla discretizada bidimensional. . . . . . . . . . . . . . . . . . . 24
2.6. Esquema Explícito bidimensional. . . . . . . . . . . . . . . . . . 24
2.7. Solución exacta. Problema tipo . . . . . . . . . . . . . . . . . . . 27
2.8. Mallado del problema tipo. . . . . . . . . . . . . . . . . . . . . . . 28
2.9. Solución aproximada. Problema tipo, explícito . . . . . . . . . . . 30
2.10.Soluciones aproximada y exacta. Problema tipo, comparación. . 31
2.11.Solución aproximada. Problema tipo explícito, estable. . . . . . . 33
2.12.Solución aproximada. Problema tipo, m = 5776. . . . . . . . . . . 34
2.13.Solución aproximada. Problema tipo, literal a), Crank-Nicolson . 39
2.14.Solución aproximada. Problema tipo, literal b).Crank −Nicolson 40
2.15.Solución exacta. Problema tipo, literal b). . . . . . . . . . . . . . . 40
2.16.Solución aproximada. Problema tipo (algoritmo 4), comprobación. 48
2.17.Solución aproximada. Problema 2 no homogéneo (algoritmo 4). . 48
3.1. Elemento lineal unidimensional. . . . . . . . . . . . . . . . . . . . 62
3.2. Soluciones aproximada y exacta. Problema tipo, n = 6. . . . . . 73
3.3. Solución aproximada. MEF. Problema tipo, n = 19. . . . . . . . . 78
4.1. Solución aproximada. Problema tipo, método Meshless. . . . . . 89
5.1. El Err máximo en el MDF . . . . . . . . . . . . . . . . . . . . . . 94
5.2. Comparación entre U y Uexat en el MDF . . . . . . . . . . . . . 94
X
XI
5.3. Comportamiento del error, nivel t = 0,0263, en el MDF . . . . . . 94
5.4. El Ecm máximo en el MDF . . . . . . . . . . . . . . . . . . . . . 95
5.5. El Emcr máximo en el MDF . . . . . . . . . . . . . . . . . . . . . 96
5.6. El Err máximo en el MEF . . . . . . . . . . . . . . . . . . . . . . 97
5.7. Comparación entre U y Uexat en el MEF . . . . . . . . . . . . . 97
5.8. Comportamiento del error, nivel t = 0,0263, en el MEF . . . . . . 97
5.9. El Ecm máximo en el MEF . . . . . . . . . . . . . . . . . . . . . 98
5.10.El Emcr máximo en el MEF . . . . . . . . . . . . . . . . . . . . . 99
5.11.El Err máximo en el MM . . . . . . . . . . . . . . . . . . . . . . . 100
5.12.Comparación entre U y Uexat en el MM . . . . . . . . . . . . . . 100
5.13.Comportamiento del error, nivel t = 0,026316, en el MM . . . . . 100
5.14.El Ecm máximo en el MM . . . . . . . . . . . . . . . . . . . . . . 101
5.15.El Emcr máximo en el MM . . . . . . . . . . . . . . . . . . . . . . 102
5.16.El Emcr máximo, nivel t = 0,013158, en el MDF . . . . . . . . . . 102
5.17.El Emcr máximo, nivel t = 0,013158, en el MEF . . . . . . . . . . 102
5.18.El Emcr máximo, nivel t = 0,48684, en el MM . . . . . . . . . . . 103
6.1. Solución aproximada. Problema 3, hipermatriz k = 1. . . . . . . . 180
6.2. Solución aproximada. Problema 3, hipermatriz k = 10. . . . . . . 181
6.3. Solución aproximada. Problema 3, hipermatriz k = 20 . . . . . . 183
6.4. Solución aproximada. Problema 3, hipermatriz k = 50. . . . . . . 184
6.5. Solución aproximada. Problema 3, hipermatriz k = 100. . . . . . 185
6.6. Solución aproximada. Problema 3, hipermatriz k = 250. . . . . . 187
6.7. Solución aproximada. Problema 3, hipermatriz k = 401. . . . . . 188
Índice de cuadros
2.1. Sistema matricial, problema tipo, Crank - Nicolson . . . . . . . . 38
2.2. Sistema matricial, problema 2 no homogéneo. . . . . . . . . . . 46
3.1. Coeficientes de las matrices AKi y MKi. . . . . . . . . . . . . . . 77
3.2. Coeficientes generalizados para n. . . . . . . . . . . . . . . . . . 77
4.1. Funciones de Base Radial . . . . . . . . . . . . . . . . . . . . . . 82
5.1. Error absoluto en tres métodos de aproximación . . . . . . . . . 104
5.2. Error cuadrático medio en tres métodos de aproximación . . . . 105
5.3. Error cuadrático relativo en tres métodos de aproximación . . . . 106
5.4. Tiempos de cómputo, problema tipo, diferencias finitas. . . . . . 106
5.5. Tiempos de cómputo en los tres métodos aproximados . . . . . 107
XII
XIII
Notación utilizada
EDP’s ecuaciones diferencial parcial (1.1)
PVF problemas de valor de frontera (1.3)
EDP ecuación diferencial parcial (1.4)
MDF método de diferencias finitas (2.1.1)
Uexat solución exacta (2.4.1)
MEF método de elementos finitos (3.1.1)
MM método Meshless (4.1)
FBR función de base radial (4.1.2)
MQ multicuádrica (4.1.3)
Err error absoluto (5.1.1)
Ecm error cuadrático medio (5.1.2)
Emcr error medio cuadrático relativo (5.1.3)
XIV
Resumen
Muchos problemas presentados en las ramas de la ingeniería, a través de la
modelación matemática, presentan ecuaciones diferenciales parciales del ti-
po parabólico. La solución analítica de una ecuación parcial proporciona una
mayor comprensión del comportamiento del fenómeno estudiado, ya que per-
mite ser determinada en todo instante de tiempo; pero, no siempre es posi-
ble encontrar esta solución. Es aquí en donde los métodos aproximados son
una herramienta de solución a este problema presentado. El presente traba-
jo pretende ser considerado como un material didáctico, guía para cualquier
docente, en donde paso a paso se mostrarán y desarrollarán, en forma pe-
dagógica, los métodos aproximados clásicos diferencias finitas y elementos
finitos y de vanguardia como el método Meshless, en la solución de estas
ecuaciones; utilizando para ello aplicaciones orientadas a la ingeniería y para
la enseñanza de la matemática a nivel superior. Se propone en las Diferencias
Finitas, utilizando fórmulas de aproximación para las derivadas, encontrar la
ecuación discretizada parcial de un problema en estudio,construyendo el ma-
llado correspondiente y encontrando su solución. Para los Elementos Finitos,
partiendo de una reformulación del problema en estudio, construir la formu-
lación débil; que, combinada con el método de Galerkin, permite plantear el
problema en un espacio de funciones finito que garantiza la existencia de la
solución. Para el Método Meshless se propone el desarrollo del mismo, ba-
sado en la colocación de puntos del dominio con Funciones de Base Radial
Multicuádrica. Adicionalmente, a más de desarrollar estos métodos se realiza
un análisis del error cometido en la aproximación, comparándole con la solu-
ción analítica conocida de un problema en estudio. Esto tiene el propósito de
demostrar la validez y confianza en el método. Además, se construyeron pro-
gramas computacionales en, Matlab, para cada uno de los métodos tratados.
PALABRAS CLAVES:
APROXIMACIÓN DE LA DERIVADA
FORMULACIÓN DÉBIL
FUNCIONES DE BASE RADIAL
DERIVADA TEMPORAL
SOLUCIÓN NUMÉRICA
PROBLEMAS DE FRONTERA
XV
Abstract
Many engineering problems are expressed and modeled using parabolic Par-
tial Differential Equations (PDE’s). Although its analytical solution provides grea-
ter understanding for the behavior of certain phenomena at any instant of ti-
me, it is not always possible to find these solutions. In accordance of this, the
approximate methods are a useful tool to numerically solve a problem with
PDE’s. By using classical approximate methods of finite difference and finite
edge elements such as the Meshless Method, this work aims to be considered
as a teaching aid guide for any instructor. Indeed, these methods are explai-
ned through a very detailed and pedagogical way. Moreover, it includes many
examples of applications oriented to engineering and mathematics teaching
at a higher level. As far as Finite Differences, I used approximation formulae
of derivatives to find the partial equation discretized, and thereby building the
corresponding meshing on which someone would find a numerical solution. In
the other hand, Finite Elements are based on a reformulation of the problem
using the Weak Derivative Formulation; with which combined the Galerkin Met-
hod allows setting the problem in a finite space that guarantees existence and
uniqueness. The Meshless method is based on a dot placement domain with
Radial Basis Function Multi-quadrics. Furthermore, I provide error analysis for
all the methods mentioned, comparing them with the known analytical solution.
All in all, this work is intended to demonstrate the validity and confidence of all
of these methods. Finally, I show several computer programs using Matlab for
each of the methods discussed.
KEY WORDS:
APPROXIMATION OF THE DERIVATIVE
WEAK DERIVATIVE FORMULATION
RADIAL BASIS FUNCTION
TEMPORARY DERIVATIVE
NUMERICAL SOLUTION
PROBLEMS OF BORDER
1
Introducción
Desarrollo de las Ecuaciones Diferenciales Parciales Parabólicas mediante
Diferencias Finitas, Elementos Finitos y Meshless
Con frecuencia no es posible encontrar la solución analítica de una ecuación
diferencial parcial (EDP) debido, como por ejemplo, a la no linealidad de la
ecuación, por el alto grado de la EDP o por tener coeficientes variables, que
constituye el modelo matemático; o debido al dominio donde se estudia el
problema. Es aquí en donde los métodos numéricos, o aproximados, son una
herramienta que permite convertir modelos matemáticos en procedimientos
computacionales, cuyos resultados pueden ser contrastados con las solucio-
nes analíticas, en los casos en los que éstas existan.
De la experiencia vivida en la profesión, material didáctico en donde se
desarrolle, paso a paso, la metodología que se aplica para resolver una ecua-
ción diferencial, es muy escasa. Como ejemplo extremo, para el método Mesh-
less, los textos de matemática avanzada para Ingeniería no consideran aún,
en sus contenidos, este método de solución.
Por tanto, este trabajo pretende incentivar el estudio de estas ecuaciones,
presentado como un material didáctico, en donde paso a paso se mostrarán
y desarrollarán en forma pedagógica métodos aproximados de solución clási-
cos como diferencias finitas y elementos finitos y de vanguardia como el de
Meshless; encaminados todos a la programación computacional en Matlab.
La aproximación por Diferencias Finitas, basada en fórmulas de aproxi-
mación para las derivadas, se presenta con fundamento matemático para la
discretización de las ecuaciones requeridas, para resolver un problema de
frontera determinado y considerando criterios de estabilidad. Se desarrollan
problemas con EDP’s del tipo homogénea y no homogénea, unidimensional,
mediante los métodos Explícito y Crank-Nicolson. Adicionalmente, se desa-
rrolla un problema bidimensional por el método Explícito. Las soluciones se
presentan en forma numérica y gráfica.
2
Los Elementos Finitos, como aproximación numérica para resolver proble-
mas de frontera, parten de una reformulación de este problema; para lo cual
se debe construir la formulación débil, en la que se realiza la discretización del
tiempo continuo, para obtener la formulación discreta. Esta formulación última,
combinada con el método de Galerkin, permite plantear la ecuación diferencial
en un espacio de funciones adecuado que permite encontrar y verificar la exis-
tencia de la solución de esta formulación; correspondiendo ésta a la solución
del problema inicial.
El Método Meshless, a diferencia de los métodos diferencias finitas y ele-
mentos finitos, no requiere de un mallado para obtener la solución numérica
de una ecuación diferencial parcial (EDP). Se basa en la colocación de puntos
del dominio del problema, con Funciones de Base Radial del tipo Multicuádrico
[23], realizando previamente una discretización de la variable temporal. Todo
se encamina a resolver, finalmente, un sistema lineal de ecuaciones del cual
se obtiene la solución numérica del problema.
Este trabajo se compone de 6 Capítulos. El Capítulo 1 resalta la importan-
cia de trabajar con este tipo de ecuaciones y hace una pequeña descripción
de cómo operan estos métodos aproximados de solución. En el Capítulo 2 se
trata, en forma detallada, el método de las Diferencias Finitas, con fundamen-
to matemático, aplicado a resolver problemas de frontera con condición inicial
en 1 y 2 dimensiones. El Capítulo 3 se enfoca en el método de los Elemen-
tos Finitos para resolver un problema homogéneo unidimensional, partiendo
de una base conceptual con fundamento matemático. El Capítulo 4 trata el
método sin mallado o de Meshless, para un problema homogéneo unidimen-
sional, presentando el fundamento matemático correspondiente hasta llegar a
obtener la solución numérica del mismo. En el Capítulo 5 se hace el cálculo
y análisis del error cometido en la aproximación, para los tres métodos; esta-
bleciendo finalmente las Conclusiones y Recomendaciones pertinentes en el
Capítulo 6.
Capítulo 1
Ecuaciones Parabólicas
1.1. Ecuaciones Diferenciales Parciales
La teoría clásica de las Ecuaciones Diferenciales en Derivadas Parciales
(EDP’s) las clasifica en tres grandes grupos: elípticas, parabólicas e hiperbóli-
cas. En el presente trabajo se trabaja sobre las del tipo parabólico, ya que sus
características matemáticas son muy distintas a las otras. La ecuación del ca-
lor es el prototipo de ecuación de evolución de tipo parabólico cuyas variantes
están presentes de manera sistemática en todos los modelos matemáticos de
la Difusión y de la Mecánica de Fluidos [4]. Se trata de un modelo fuertemente
irreversible en tiempo en que la información se propaga rápidamente.
1.2. Importancia y Selección
En problemas de Ingeniería es común encontrar modelos matemáticos
que incluyen ecuaciones diferenciales en derivadas parciales del tipo parabóli-
co [1]. Este trabajo pretende exponer los conocimientos, en los métodos numé-
ricos, que se consideran básicos para que un estudiante de ingeniería pueda
entender, plantear y desarrollar un modelo matemáticobasado en EDP’s; el
cual servirá de base para formulaciones más complejas y reales en proble-
mas que se presentan al estudiar los procesos de conductibilidad térmica y
difusión.
3
4
Como ejemplo se tiene el modelo unidimensional del flujo de calor en un
alambre aislado de longitud a [1]. La ecuación del calor, que nos da la tempe-
ratura U(x, t) en la posición x del alambre y en el instante t, es
cUxx(x, t) = σ · δ · Ut(x, t) para 0 ≤ x ≤ a y 0 ≤ t <∞,
cuya condición inicial de temperatura es
U(x, 0) = f(x) para t = 0 y 0 < x < a,
y sus condiciones de contorno a considerar son
U(0, t) = c1 para x = 0 y 0 < t <∞,
U(a, t) = c2 para x = a y 0 < t <∞,
siendo c, el coeficiente de conducción térmica, σ, el calor específico, y, δ, la
densidad del material, constantes a considerar.
La selección de este tipo de EDP’s se efectuó tomando en consideración
las múltiples aplicaciones como por ejemplo en aerodinámica, acústica, elas-
ticidad, transferencia de calor, metereología, física del plasma, etc.; presen-
tando, a continuación, algunos ejemplos de gran relevancia en la Física y la
Ingeniería.
Las ecuaciones de Maxwell, que describen la propagación de las ondas
electromagnéticas en un determinado medio, siendo una de ellas
∇x ~H = 4π
c
~J,
donde ~H es el campo magnético, ~J la densidad de flujo de carga [17] y c
la velocidad de la luz.
La ecuación de Schrödinger de la mecánica cuántica, que describe la
evolución de los sistemas microscópicos como son los átomos o las mo-
léculas [18], siendo esta
− h
2
2m
∇2 · U + V (~r) · U = i · ~ · Ut,
donde U es la función de onda, V la energía potencial, ~r es la energía
de posición, ~ la constante de Planck e i la unidad imaginaria.
5
Las EDP’s del tipo parabólico son muy utilizadas en el estudio del creci-
miento de tumores y la forma en que se difunden sobre los tejidos que
los rodean. En la utilización de estas ecuaciones se supone un sistema
descrito bajo un comportamiento mecánico, donde el sistema puede ser
un fluido o una mezcla entre líquido y sólido (los fluidos son normalmen-
te los nutrientes). Se presentan fenómenos de difusión y transporte de
nutrientes teniendo en cuenta efectos de concentración, tamaño y velo-
cidad de células. En este modelo se aplican condiciones iniciales y de
frontera que están relacionadas con el tamaño, la permeabilidad del me-
dio, la geometría y las dimensiones del sistema. Para migraciones no
estacionarias debe satisfacer la EDP
c = ∇2n+ V (c, n),
siendo D el coeficiente de difusión, V (c, n) la función específica de pro-
ducción de la señal química [19], n la densidad y c la concentración.
En el caso particular de la Ingeniería de Alimentos, es necesario me-
jorar los procesos tradicionales de conservación, y de optimización en
el uso de los recursos energéticos; siendo necesario estudiar aplicacio-
nes matemáticas y computacionales para el análisis del calentamiento
convectivo de alimentos envasados. Este es un ejemplo de aplicaciones
ingenieriles a temas vinculados con los procesos biotecnológicos, parti-
cularmente en la conservación de alimentos. Así, para la determinación
de coeficientes de difusión, de agua y de sólidos, se considera la segun-
da ley de Fick para la difusión unidireccional
Ct = Def · C2xx,
siendo C la concentración y Def el coeficiente de difusión efectivo[20].
Los fluidos perfectos o ideales, si bien constituyen modelos matemáticos
interesantes, no dejan de ser un tanto irrealistas en la medida en que
todo fluido posee un cierto grado de viscosidad. La ecuación de Burgers
viscosa, se presenta en un sentido más estricto y riguroso, adoptando la
6
forma:
Ut − v · Uxx + U · Ux = 0,
donde U es la propagación de la onda y v > 0 la constante de viscosidad
[16].
Aunque la evolución de un sistema cancerígeno está limitada y condicio-
nada a parámetros químicos y biológicos, su diseminación a través de
los tejidos y órganos contaminados obedece a procesos estocásticos,
es decir, que describen la evolución temporal de una variable continua;
las cuales, en la modelación desarrollan EDP’s de este tipo. Ánderson y
Chaplain plantearon una EDP para describir la dinámica de la densidad
de células endoteliales que migran a través de un tumor
n = D · ∇2n −∇(X(c) · n · ∇c) + g(n, c),
siendo D el coeficiente de difusión, X(c) el parámetro químico táctico,
c(x, t) la concentración de TAF químico específico [19] y g(n, c)=la fun-
ción de proliferación.
Estos pocos ejemplos muestran el interés de su estudio; ahora, se pondrá
énfasis en las EDP’s de tipo parabólico.
1.3. Ecuaciones del tipo parabólico
Requiere la solución de este tipo de EDP’s de un dominio espacio-temporal
abierto. La solución está sujeta a un conjunto de condiciones iniciales y de
frontera, correspondiente a Problemas de Valor de Frontera (PVF) con condi-
ción inicial. La solución se calculará yendo hacia adelante a partir de la condi-
ción inicial y satisfaciendo siempre las condiciones de contorno. Como ejem-
plo se tiene la ecuación de la conducción del calor, con coeficiente térmico c
constante [4].
Ut(x, t) = cUxx(x, t) +Q(x, t) para 0 ≤ x ≤ a y 0 ≤ t <∞,
7
en donde Q(x, t) representa las fuentes de calor existentes, sujeta a la condi-
ción inicial
U(x, 0) = f(x), para 0 < x < a,
y a dos condiciones de contorno
U(0, t) = g1(t) para t > 0,
U(a, t) = g2(t) para t > 0.
1.4. Métodos de solución
En este trabajo se considerarán tres métodos aproximados para encontrar
la solución de una EDP; siendo estos: Diferencias Finitas, Elementos Finitos
y el método Meshless. Es importante recalcar que muchos PVF no se pueden
resolver satisfactoriamente mediante técnicas analíticas, por lo que es nece-
sario obtener aproximaciones numéricas de las soluciones.
1.4.1. Diferencias finitas
Se basa en la utilización de fórmulas para aproximar las derivadas de una
función. Estas fórmulas de aproximación de las derivadas de una función U
pueden ser: centradas, progresivas o regresivas, con un orden de la aproxi-
mación O(hn) [1], n = 1, 2, ... Como ejemplo se tienen las siguientes fórmulas:
Para la primera derivada de la función U(x, t), con respecto a x, utilizando
un orden de aproximación O(h2) y O(h)
Diferencias centradas
Ux(x, t) =
U(x+ h, t)− U(x− h, t)
2h
+O(h2).
Diferencias progresivas
Ux(x, t) =
U(x+ h, t)− U(x, t)
h
+O(h).
8
Diferencias regresivas
Ux(x, t) =
U(x, t)− U(x− h, t)
h
+O(h).
Otro ejemplo es la fórmula de la segunda derivada de la función U(x, t), cen-
trada y con respecto a x, con orden de aproximación O(h2)
Uxx(x, t) =
U(x+ h, t)− 2U(x, t) + U(x− h, t)
h2
+O(h2).
Fácilmente se puede aproximar la derivada n-ésima de U(x, t) con un orden
de aproximación O(hn) [1][6].
1.4.2. Método de Elementos Finitos
Es una técnica numérica muy importante que se ha destacado especial-
mente por modelar en dominios irregulares, condiciones de contorno, sistemas
de cargas complejos, etc.; así como también por la facilidad en la selección del
mecanismo de aproximación de las variables involucradas en un problema es-
pecífico.
El análisis de un problema mediante el método de elementos finitos par-
te de la definición de su dominio, el cual se debe discretizar en regiones de
tamaño finito. Se debe a continuación identificar la variable de estado, esta-
bleciendo sistemas de referencia: locales para los elementos y globales para
el sistema completo; los cuales permiten construir las funciones de aproxima-
ción de cada elemento, y, a partir de estas determinar las ecuaciones a nivel
de cada elemento mediante los métodos: directo, variacional o de residuos
ponderados[10]. Luego introduciendo las condiciones de contorno, se debe
resolver el sistema de ecuaciones resultante, obteniendo valores aproximados
de las variables de estado en los nodos del dominio. Estos valores deben ser
interpretados y pueden usarse para el cálculo de otras cantidades físicas tales
como flujos de calor, esfuerzos, etc. En la práctica,se obtienen resultados con-
fiables comparando diferentes análisis, basados en diferentes discretizaciones
del dominio, para el mismo problema. Se profundizará sobre este método al
tratarlo con detalle en el Capítulo III, de este trabajo.
9
1.4.3. Método de Meshless
Comúnmente, encontrar una aproximación a la solución de una EDP se
realiza empleando los métodos de Diferencias Finitas y Elementos Finitos. Es-
tos métodos tienen en común, que requieren de un mallado que da el soporte
a la solución numérica con un orden de convergencia algebraico.
Para reducir los tiempos de cómputo se han desarrollado los métodos li-
bres de malla (Meshless), y, como su nombre lo indica no requieren de un
mallado para obtener la solución numérica de la EDP. Los métodos libres de
mallas están basados en un conjunto independiente de puntos o nodos y el
costo de la generación de la malla es eliminado[22]. Se presenta la solución
numérica de las Ecuaciones Diferenciales Parciales empleando funciones de
base radial, mediante el método de colocación asimétrico; el cual está basado
en un conjunto de nodos aleatorios y no requiere de un proceso de mallado
como ocurre en los métodos de Diferencias Finitas y Elementos Finitos, obje-
tos de estudio de este trabajo. Se profundizará sobre este método al tratarlo
con detalle en el Capítulo IV.
1.4.4. Formulación y Objetivo General
Considerando los criterios de claridad, síntesis y operatividad, el problema
o preguntas que se planteó contestar con este trabajo son: ¿Cómo resolver
una EDP que no tiene solución analítica?, ¿Los métodos aproximados son
una alternativa de solución para las EDP’s que no tienen solución analítica?,
¿Los métodos aproximados son realmente confiables?, ¿Qué método aproxi-
mado es el más recomendable para resolver EDP’s del tipo parabólico? Con
estas interrogantes planteadas el problema específico, objeto de este trabajo,
se formuló en los términos: “¿Cómo desarrollar las Ecuaciones Diferenciales
Parciales Parabólicas mediante Diferencias Finitas, Elementos Finitos y Mesh-
less?”. A partir de aquí se plantean los objetivos de este trabajo.
10
Objetivo General
Resolver las ecuaciones diferenciales parciales parabólicas aplicando los
métodos: Diferencias Finitas, Elementos Finitos y Meshless, determinando su
grado de confiabilidad y, presentándolos en forma didáctica como un material
científico, útil para el aprendizaje a nivel de pregrado en los cursos de ingenie-
ría.
Objetivos Específicos
Revisar las metodologías de solución aproximada de las EDP’s parabó-
licas .
Determinar la validez de las soluciones numéricas de las EDP’s, contras-
tando con las soluciones analíticassobre problemas con solución analíti-
ca conocida.
Comparar los resultados obtenidos con los tres métodos sobre un pro-
blema con solución analítica conocida.
Presentar un material didáctico para la enseñanza de este tipo de EDP’s
a nivel superior.
Capítulo 2
Diferencias Finitas
2.1. Ecuación Diferencial Parcial Unidimensional
Se trabajará en este Capítulo con las ecuaciones diferencialles parcia-
les homogéneas y no homogéneas. Se define como una EDP homogénea a
aquella de la forma f(Ut, Uxx, U) = 0; mientras que, una ecuación de la forma
f(Ut, Uxx, U, x, t) = 0 será una EDP no homogénea.
2.1.1. Ecuaciones Homogéneas
Se inicia el análisis de este método encontrando la solución aproximada de
la ecuación parabólica típica, la del Calor y sin fuentes, en el intervalo acotado
0 ≤ x ≤ a, aplicando el Método de Diferencias Finitas (MDF).
Sea
Ut(x, t) = cUxx(x, t), (2.1)
con la condición inicial
U(x, 0) = f(x), 0 < x < a, t = 0, (2.2)
y las condiciones de frontera
U(0, t) = 0, para x = 0 y 0 < t ≤ b, (2.3)
U(a, t) = 0, para x = a y 0 < t ≤ b; (2.4)
11
12
en donde c es una constante térmica que depende del problema particular.
Se trata de sustituír la EDP (2.1), dada en un punto (x, t) cualquiera, por
una aproximación basada en fórmulas en diferencias finitas para las derivadas.
En forma arbitraria, tomando una diferencia progresiva en el tiempo para
Ut(x, t), discretizando en intervalos iguales de magnitud k, se tiene[4]
Ut(x, t) =
U(x, t+ k)− U(x, t)
k
− 1
2
Utt(x, η1)k, (2.5)
en donde t < η1 < t+ k.
Para la derivada espacial, la fórmula en diferencias espaciales, centrada y
discretizada en intervalos iguales de magnitud h es[4]
Uxx(x, t) =
U(x+ h, t)− 2U(x, t) + U(x− h, t)
h2
− 1
12
Utttt(ζ1, t)h
2, (2.6)
siendo x < ζ1 < x+ h
Remplazando las ecuaciones (2.5) y (2.6) en (2.1), la ecuación del calor
en un punto cualquiera (x,t) es
U(x, t+ k)− U(x, t)
k
= c
U(x+ h, t)− 2U(x, t) + U(x− h, t)
h2
+ E, (2.7)
siendo el error de truncamiento
E =
1
2
Utt(x, η1)k,−
1
12
Utttt(ζ1, t)h
2. (2.8)
Ahora, no se puede resolver la ecuación (2.7) por ser el error E desconocido.
Se procede, entonces, a ignorar el error de truncamiento E, obteniendo de
esta forma una aproximación de la ecuación (2.7), que es
U(x, t+ k)− U(x, t)
k
= c
U(x+ h, t)− 2U(x, t) + U(x− h, t)
h2
. (2.9)
Es importante hacer notar que el error de truncamiento local, dado por la ecua-
ción (2.8), es el máximo entre O(k) y O(h2). Como el error E → 0, cuando
h→ 0 y k → 0, se puede decir que la aproximación dada por la ecuación (2.9)
es consistente con la ecuación en derivadas parciales dada por la ecuación
(2.1).
Se espera que la solución U(x, t) aproximada, encontrada al resolver la
ecuación (2.9), se aproxime de forma precisa a la solución exacta que se en-
contraría al resolver (2.7).
13
Por otra parte, la ecuación (2.9) involucra puntos separados que permiten
introducir una malla espacio - tiempo, representado en la Figura 2.1, mediante
una discretización constante en el tiempo (k).
Figura 2.1: Mallado espacio - tiempo.
En este mallado, se divide la magnitud espacial x, de longitud a, en n
intervalos iguales; cada uno de longitud h = a/n, obteniendo (n + 1) puntos,
que en forma general se pueden representar
xi = (i− 1)h i = 1, 2, ...(n+ 1).
En forma análoga para el tiempo b, k = b
m
, siendo m el número de intervalos
iguales en que se divide el tiempo b
tj = (j − 1)k, j = 1, 2, ..., (m+ 1)
Por tanto, la temperatura exacta en un punto de la malla se aproxima por
U(xi, tj), que satisface la ecuación (2.9).
Ahora, se debe pensar ya en llevar todo este proceso a la programación;
para lo cual se introduce la notación
U(xi, tj) = U(i, j) o U(xi−h, tj+2k) = U(i− 1, j + 2);
bajo esta notación, la ecuación (2.9) toma la forma
U(i, j + 1)− U(i, j)
k
= c
U(i+ 1, j)− 2U(i, j) + U(i− i, j)
h2
, (2.10)
para i = 2, 3, ..., n y j > 1.
14
Adicionalmente, se puede verificar que la ecuación (2.10) satisface la con-
dición inicial en los puntos de la malla
U(x, 0) = U(i, 1) = f(xi), i = 2, 3, ..., n; (2.11)
y las condiciones de frontera
U(0, t) = U(1, j) = 0, j > 1 (2.12)
U(a, t) = U(n+ 1, j) = 0, j > 1. (2.13)
Ahora, si se define
r =
ck
h2
, (2.14)
remplazando en la ecuación (2.10) y simplificando la misma se tiene
U(i, j + 1) = U(i, j) + r[U(i+ 1, j)− 2U(i, j) + U(i− i, j)], (2.15)
siendo r un parámetro adimensional.
El valor U(i, j + 1) es una combinación lineal de tres valores en el nivel
temporal j, debidamente especificados (Ver Figura 2.2).
Figura 2.2: Esquema de diferencias progresivas.
El cálculo de los valores de la solución aproximada, para la primera colum-
na de la malla, se encuentra aplicando la condición inicial
U(i, 1) = f(xi), i = 2, 3, ..., n.
Luego, por la ecuación (2.15), se obtienen los valores de la solución U(i, 2),
i = 2, 3, ..., n; es decir, se calculan los valores de la segunda columna de la
15
malla; permitiendo continuar con el cálculo en forma análoga para las demás
columnas.
Los valores de la solución en las fronteras de la malla; es decir, los de las
filas 1 y (n + 1), se obtienen aplicando las condiciones de frontera dadas por
las ecuaciones (2.12) y (2.13).
De esta forma se puede resolver esta EDP, discreta, numéricamente. El
método propuesto se programafácilmente en un computador personal, el mis-
mo que se desarrollará y explicará con detalle, más adelante al resolver los
problemas propuestos para este trabajo.
Como se puede observar, la ecuación en diferencias parciales (2.15) de-
pende principalmente de r = ck
h2
. Se debe, por tanto, entender como escoger
r para poder obtener soluciones numéricas razonables, asociadas con la EDP
en estudio y, que sean una buena aproximación a la solución exacta del pro-
blema.
Análisis de la estabilidad del método
El método de separación de variables, utilizado para resolver analíticamen-
te una EDP, se puede aplicar también a una ecuación en diferencias parciales;
es decir, suponiendo que la ecuación (2.15) tiene soluciones producto de la
forma U(i, j) = θiϕj. Se pueden obtener soluciones homogéneas sustituyen-
do en esta θi = Qi[4].
No es de interés, en el presente trabajo, profundizar el cálculo analítico;
pero analizando en forma rápida, las condiciones de contorno φ0 = φn = 0
sugieren que la solución pueda oscilar. Esto ocurre cuando Q es complejo1,
con |Q| = 1, por lo que una sustitución equivalente es [4]
φi = (|Q|eîθ)i = eîθi = eîθ(
xi
h
) = eîαxi
Sobre la base de lo último anotado y para analizar la estabilidad del método ex-
plícito, se toman soluciones producto especiales, con número de onda α, con-
siderando tj = (j−1)k, j = 1, 2, ..., (m+1) y xi = (i−1)h, i = 1, 2, ..., (n+1);
1La notación utilizada en este trabajo para un número complejo es z = a+ b̂i, para diferen-
ciarle del contador i.
16
según el mallado propuesto (ver Figura 2.1)
U(i, j) = eîαxiQ
tj
k = eîα(i−1)hQ(j−1). (2.16)
Sustituyendo (2.16) en (2.15)
eîα(i−1)hQ(j) = eîα(i−1)hQ(j−1) + r[eîα(i)hQ(j−1) − 2eîα(i−1)hQ(j−1) + eîα(i−2)hQ(j−1)]
Sacando el factor común eîα(i−1)hQ(j−1) y simplificando en la última expresión
Q = 1 + r[eiαh − 2 + e−iαh]
Q = 1− 2r[1− ( eiαh+e−iαh
2
)];
pero
( e
iαh+e−iαh
2
) = cos(αh),
entonces, el valor de Q que se obtiene es
Q = 1− 2r[1− cos(αh)]. (2.17)
Ahora, se puede utilizar una combinación lineal de las funciones e±îαx para la
ecuación (2.16), que se ajusten a las condiciones del PVF.
La condición de frontera U(0, t) = 0 se cumple con la función sen(αx).
La condición U(a, t) = 0 implica que α = n̂π
a
, n̂ = 1, 2, 3, ....
Por tanto, hay solución para el PVF, dado por las ecuaciones (2.15), (2.12) y
(2.13), de la forma
U(i, j) = sen
(
n̂πxi
a
)
Q
tj
k ,
en donde Q se determina a partir de la ecuación (2.17), obteniendo una solu-
ción producto U(x, t) de la ecuación en derivadas parciales
U(i, j) = sen
(
n̂πxi
a
)[
1− 2r
(
1− cos
(
n̂πh
a
))] tj
k
, n̂ = 1, 2, ..., (n− 1),
(2.18)
siendo r = ck
h2
. El número máximo de auto funciones independientes para esta
ecuación, en diferencias parciales, es[4] n̂ = (n− 1).
En la ecuación (2.18), la dependencia en el tiempo está dada por
Q
tj
k =
[
1− 2r
(
1− cos
(
n̂π
n
))] tj
k
. (2.19)
17
Pero, el valor que toma Q es fundamental para lograr oscilaciones convergen-
tes, en el tiempo, para la solución[4].
Si Q > 1 o Q < −1, existe un crecimiento exponencial en el tiempo de la
solución buscada.
Si 0 < Q < 1 o −1 < Q < 0, se presenta un decaimiento exponencial.
Si Q = ±1, la solución es constante en el tiempo.
Ahora, es posible que haya una oscilación convergente en el tiempo si |Q| ≤ 1.
El valor de Q, por tanto, determina la estabilidad. En resumen, si |Q| ≤ 1
para todas las soluciones, el esquema numérico es estable; caso contrario es
inestable.
Por tanto, si |Q| ≤ 1, la solución no será nunca una exponencial perma-
nente creciente en el tiempo; pero, a medida que decae exponencialmente, la
solución puede ser una oscilación convergente o divergente. Lo que se pre-
tende es que el esquema numérico tenga solo oscilaciones convergentes en
el tiempo. Aclarando esto último, las oscilaciones convergentes no reproducen
el comportamiento de la ecuación en derivadas parciales; pero garantizan que
decaen.
Por otra parte, el valor de r juega un papel muy importante para Q. Si r
es demasiado grande Q puede llegar a ser demasiado negativa; pero, debe
cumplirse para que la solución sea estable |Q| ≤ 1; es decir −1 ≤ Q ≤ 1, de la
cual Q ≥ −1.
Remplazando Q en la ecuación (2.19)[
1− 2r
(
1− cos
(
n̂π
n
))]
≥ −1, n̂ = 1, 2, ..., (n− 1);
despejando r
r ≤ 1
1− cos
(
n̂π
n
) ;
pero, r debe ser menor o igual al valor más pequeño que se obtiene cuando
n̂ = (n − 1), ya que en el denominador de expresión anterior cos( n̂π
n
) → −1;
por tanto
0 < 1− cos
(
n̂π
n
)
< 2.
18
Interesa, por simplificar el problema, tomar el límite superior.
r ≤ 1
2
<
1
1− cos
(
(n−1)π
n
) .
Por tanto, existe la garantía de que la solución numérica será estable si
r ≤ 1
2
. (2.20)
Se realiza, a continuación, algunas interpretaciones al valor de r:
Si r > 1
2
, la solución numérica contendrá una oscilación divergente. Se
llama a este efecto Inestabilidad Numérica.
Si r > 1
2
, la solución que crece más rápidamente corresponde a una
oscilación rápida en el espacio n̂ = n− 1.
Como r = ck
h2
y r ≤ 1
2
, se puede escribir
k ≤
1
2
h2
c
(2.21)
Como h debe ser pequeño, entonces k debe ser demasiado pequeño.
Si por el contrario, los pasos k en el tiempo son demasiado grandes, el
esquema se convierte en inestable.
En resumen, para obtener soluciones razonables en el método explícito se
debe cumplir r ≤ 1/2; para lo cual:
1. Para un número de puntos n de discretización de la variable espacial
(0 < x < a), calculamos el espaciamiento h con la fórmula h = a/(n− 1).
2. Si es dato del problema el espaciamiento k, de la discretización de la
variable temporal (0 < t < b), calculamos el valor de r con la fórmula
r = ck/h2, siendo c una constante de la EDP en estudio.
3. Si r ≤ 1/2 el método es estable y la solución es razonable. Pero, si
r ≥ 1/2 se debe calcular un nuevo valor para k aplicando la fórmula
k ≤ h2/(2c).
4. Finalmente, con lo valores de h, k y c se calcula el nuevo valor de r que
hace que el método explícito sea estable.
El esquema planteado, como se manifestó al inicio, se conoce como explícito.
19
2.1.2. Ecuaciones No Homogeneas
Las EDP’s parabólicas como la del calor, con fuentes, es un caso típico de
estas ecuaciones. Podemos resolver numéricamente esta ecuación en forma
análoga a las EDP’s homogéneas.
Sea
Ut(x, t) = cUxx(x, t) + q(x, t), (2.22)
con la condición inicial
U(x, 0) = f(x), 0 < x < a, t = 0,
y las condiciones de frontera
U(0, t) = g1(t), para x = 0 y 0 < t ≤ b,
U(a, t) = g2(t), para x = a y 0 < t ≤ b;
en donde c es una constante térmica en la EDP.
Utilizando diferencias progresivas en el tiempo y diferencias centradas en
el espacio, ecuaciones (2.5) y (2.6) respectivamente, reemplazando en (2.22)
se obtiene la siguiente aproximación numérica.
U(i, j + 1)− U(i, j)
k
= c
U(i+ 1, j)− 2U(i, j) + U(i− 1, j)
h2
+ q(xi, tj), (2.23)
para i = 2, 3, ..., n y j > 1; sujeta a las condiciones
U(i, 1) = f(xi), j = 0, i = 2, 3, ..., n
U(1, j) = g1(tj), i = 1, j > 1
U(n+ 1, j) = g2(tj), i = (n+ 1), j > 1
La solución se obtiene fácilmente calculando U(i, j+1), a partir de la ecua-
ción (2.23).
El análisis de estabilidad, realizado anteriormente, es válido también para
problemas no homogéneos (se demuestra en forma análoga); es decir, para
estos cálculos se tomará
r =
ck
h2
≤ 1
2
20
2.2. Otros esquemas numéricos
El esquema explícito, analizado anteriormente, utiliza diferencias progresi-
vas en el tiempo con espaciamiento k muy pequeño y diferencias centradas
en el espacio. Además, es estable si ck
h2
≤ 1
2
.
Este esquema numérico tiene un costo computacional alto y, el error de
truncamiento corresponde a la suma de O(k) y O(h2). Ahora, solo si r = 1
2
se
consigue ambos errores con un orden de aproximación O(h2), ya que en este
caso k = rh
2
c
.
Entonces, es necesario contar con un esquema numérico de menor costo
computacional; así por ejemplo: están el esquema de Richardson y Crank -
Nicolson.
2.2.1. Esquema de RichardsonLa propuesta es intentar obtener unas diferencias en el tiempo más pre-
cisas. Es así que Richardson, en 1927, propuso utilizar diferencias centradas
para el tiempo y para el espacio; para lo cual planteó la ecuación
U(i, j + 1)− U(i, j − 1)
2k
= c
U(i+ 1, j)− 2U(i, j) + U(i− i, j)
h2
,
para i = 2, 3, ..., n y j > 1.
Simplificando esta última expresión y llamando r = ck
h2
, se tiene la ecuación
U(i, j + 1) = U(i, j − 1) + 2r[U(i+ 1, j)− 2U(i, j) + U(i− 1, j)]
En este caso, el error de truncamiento es la suma de los errores O(k2) y O(h2).
Este esquema no se analizará en el presente trabajo; pero, aunque es más
preciso que el esquema explícito, por lo general, es un método numérico muy
inestable [4] [5].
2.2.2. Esquema de Crank - Nicolson
Crank y Nicolson (en 1947) sugirieron una nueva alternativa de utilizar las
diferencias centradas. El artificio es interpretar una diferencia progresiva en el
21
tiempo
Ut(i, j) =
U(i, j + 1)− U(i, j)
k
,
como una diferencia centrada alrededor del punto (i, j + 1
2
) (Figura 2.3).
Expliquemos esto con el PVF dado por las ecuaciones (2.1), (2.2), (2.3),
y (2.4), despreciando el error de truncamiento. Para el nivel j, la diferencia
progresiva en el tiempo es
Figura 2.3: Esquema del punto (i, j + 1
2
).
Ut(i, j) =
U(i, j + 1)− U(i, j)
k
; (2.24)
y la diferencia centrada en el espacio es
Uxx(i, j) =
U(i− 1, j)− 2U(i, j) + U(i+ 1, j)
h2
. (2.25)
Remplazando las ecuaciones (2.24) y (2.25) en la ecuación (2.1) nos queda
U(i, j + 1)− U(i, j)
k
= c
U(i− 1, j)− 2U(i, j) + U(i+ 1, j)
h2
. (2.26)
Ahora, aplicando el esquema regresivo en el nivel j + 1, considerando la dife-
rencia regresiva en el tiempo
Ut(i, j + 1) =
U(i, j + 1)− U(i, j)
k
, (2.27)
y la diferencia centrada espacial
Uxx(i, j + 1) =
U(i− 1, j + 1)− 2U(i, j + 1) + U(i+ 1, j + 1)
h2
, (2.28)
22
al reemplazar (2.27) y (2.28) en (2.11), se tiene
U(i, j + 1)− U(i, j)
k
= c
U(i− 1, j + 1)− 2U(i, j + 1) + U(i+ 1, j + 1)
h2
. (2.29)
Para cumplir con el esquema de Crank-Nicholson, se sumaron las ecuaciones
(2.26) y (2.29); obteniendo
2
U(i, j + 1)− U(i, j)
k
= c
[
U(i− 1, j)− 2U(i, j) + U(i+ 1, j)
h2
(2.30)
+
U(i− 1, j + 1)− 2U(i, j + 1) + U(i+ 1, j + 1)
h2
]
.
El error al aproximar Ut(i, j+ 12) es del orden O(k
2). La discretización realizada
para Uxx(i, j + 12) con esquemas de diferencias centradas es del orden O(h
2).
La ecuación (2.30) corresponde al valor de la función discretizada en este
punto intermedio del tiempo, correspondiendo a la media entre t y t + k. Esto
representa el esquema de Crank-Nicolson[4] [1].
Entre las ventajas más importantes del método se pueden resaltar:
El error de truncamiento corresponde a la suma de O(h2) y O(k2).
Este método tiene la ventaja de ser estable[4] para todo r = ck
h2
.
El valor de k puede ser tan grande como queramos ya que el error es del
orden O(h2), permitiendo realizar un menor costo computacional; que
hera lo que se quería alcanzar.
Se conoce al esquema de Crank-Nicolson como Implícito.
Como desventaja del método, como se puede ver en la Figura 2.4, utiliza
seis puntos; tres de los cuales corresponden al tiempo posterior. Esto impi-
de avanzar directamente hacia adelante en el tiempo aplicando la ecuación
(2.30), sino que para cada nivel de tiempo calculado se requiere resolver un
sistema de (n-2) ecuaciones, con una matriz de coeficientes tridiagonal.
Este método implícito, con detalle y con la programación correspondiente,
se desarrollará más adelante al resolver los problemas seleccionados para
este trabajo.
23
Figura 2.4: Esquema Implícito.
2.3. Ecuación Bidimensional
Aplicamos criterios análogos a las EDP’s unidimensionales para encontrar
soluciones numéricas a los PVF bidimensionales.
Sea la EDP
Ut(x, t) = α(Uxx(x, t) + Uyy(x, t)). (2.31)
Se requiere la construcción de una malla bidimensional, discretizada, según
se indica en la Figura 2.5.
Utilizando diferencias progresivas para Ut y diferencias centradas en x y y
para el Laplaciano, y despreciando el error de truncamiento, se tiene
Ukt (i, j) =
Uk+1(i, j)− Uk(i, j)
p
(2.32)
Ukxx(i, j) =
Uk(i− 1, j)− 2Uk(i, j) + Uk(i+ 1, j)
h21
(2.33)
Ukyy(i, j) =
Uk(i, j − 1)− 2Uk(i, j) + Uk(i, j + 1)
h22
(2.34)
Reemplazando (2.32), (2.33) y (2.34) en (2.31) se tiene
Uk+1(i, j)− Uk(i, j)
p
= α
[
Uk(i− 1, j)− 2Uk(i, j) + Uk(i+ 1, j)
h21
(2.35)
+
Uk(i, j − 1)− 2Uk(i, j) + Uk(i, j + 1)
h22
]
24
Figura 2.5: Malla discretizada bidimensional.
La ecuación (2.35) nos permite avanzar directamente hacia adelante con
el tiempo; es decir encontramos el valor de la solución Uk+1(i, j), utilizando 5
puntos conocidos del nivel temporal anterior k (Ver Figura 2.6). Pero, por ser
un método explícito es importante considrar criterios de estabilidad.
Figura 2.6: Esquema Explícito bidimensional.
2.3.1. Análisis de la estabilidad del método
Como se analizó anteriormente, este esquema explícito puede ser ines-
table. Ignorando las condiciones de contorno para este análisis, nos interesa
investigar el posible crecimiento de las ondas periódicas en la variable espa-
cial. Tomando soluciones producto especiales, con números de ondas γ1 y γ2,
25
para x y y respectivamente
Uk(i, j) = eî(γ1xi)+γ2yjQ
tk
p
Considerando que:
tk = (k − 1)P, k = 1, 2, ..., (l + 1)
xi = (i− 1)h1, i = 1, 2, ..., (n+ 1)
yj = (j − 1)h2, j = 1, 2, ..., (m+ 1);
según el mallado propuesto (ver Figura 2.5), se tiene
Uk(i, j) = eî(γ1(i−1)h1)+γ2(j−1)h2Qk−1 (2.36)
Por otra parte, por facilidad para este análisis de estabilidad, se considera
h1 = h2 = h; por lo que la ecuación (2.35) toma la forma
Uk+1(i, j)− Uk(i, j)
p
=
α
h2
[Uk(i− 1, j) + Uk(i+ 1, j) (2.37)
+Uk(i, j − 1) + Uk(i, j + 1)− 4Uk(i, j)]
Sustituyendo la ecuación (2.36) en (2.37)
eî(γ1(i−1)h)+γ2(j−1)hQk − eî(γ1(i−1)h)+γ2(j−1)hQk−1 = r[eî(γ1(i)h)+γ2(j−1)hQk−1
+eî(γ1(i−2)h)+γ2(j−1)hQk−1
+eî(γ1(i−1)h)+γ2(j)hQk−1
+eî(γ1(i−1)h)+γ2(j−2)hQk−1
+eî(γ1(i−1)h)+γ2(j−1)hQk−1]
Sacando factor común eî(γ1(i−1)h)+γ2(j−1)hQk−1 y simplificando la expresión an-
terior, se tiene:
Q = 1 + 2r[ e
iγ1h+e−iγ1h
2
+ e
jγ2h+e−jγ2h
2
− 2],
es decir
Q = 1 + 2r[cos(γ1h) + cos(γ2h)− 2]. (2.38)
Para garantizar la estabilidad se sabe que |Q| ≤ 1; es decir −1 ≤ Q ≤ 1.
Vamos a considerar Q ≥ −1. Entonces, la ecuación (2.38) toma la forma
1 + 2r[cos(γ1h) + cos(γ2h)− 2] ≥ −1,
26
y
r < − 1
cos(γ1h) + cos(γ2h)− 2
. (2.39)
Además, se puede afirmar con exactitud
−4 < cos(γ1h) + cos(γ2h)− 2 < 0
Interesa, por simplificar el problema, tomar el de mayor valor; es decir
cos(γ1h) + cos(γ2h)− 2 = −4
Considerando la expresión anterior, la ecuación (2.39) se puede escribir
r ≤ 1
4
;
es decir
r =
αp
h2
≤ 1
4
(2.40)
2.4. Ejercicios de aplicación
Se pretende estudiar la ecuación del calor, con temperatura cero en los
extremos, por ser un problema físico relevante [4][1] correspondiente a una
varilla unidimensional (0 < x < a), sin fuentes (homogéneo) y con fuente (no
homogéneo), con ambos extremos sumergidos en un depósito a 0◦ de tempe-
ratura; y, de esta forma descubrir cómo cambia la energía térmica inicial dada
por la condición inicial. Además, resolver el problema homogéneo permitirá
poder resolver el problema no homogéneo.
2.4.1. Problema tipo
Sea la ecuación diferencial del calor
Ut(x, t) = Uxx(x, t) para 0 ≤ x ≤ 1 y 0 ≤ t <∞,
con la condición inicial de temperatura
U(x, 0) = sen(πx) · (1 + 2cos(πx)) para t = 0 y 0 < x < 1,
27
y sus condiciones de frontera (o de contorno) a considerar
U(0, t) = 0 para x = 0 y 0 < t <∞,
U(1, t) = 0 para x = 1 y 0 < t <∞.
Se tiene como solución exacta de este problema la función U(x, t) = sen(πx)e−π2t+
sen(2πx)e−4π
2t. Para un mallado de (nxm) puntos de discretización de las va-
riables espacial y temporal, las líneas de programación, en Matlab, que per-
miten representar gráficamente la solución exacta U(x, t) se encuentran en la
tabla de algoritmos, en el programa denominado Algoritmo 1, problema tipo,
solución exacta.
Ejecutandoesta función de Matlab, para los datos del problema, se tiene la
solución exacta U(x, t) que se encuentra en la tabla de matrices con el nombre
Solución exacta, problema tipo.
La gráfica de la solución exacta, del problema en estudio, se representa
en la Figura 2.7.
Figura 2.7: Solución exacta. Problema tipo
2.4.2. Problema tipo. Esquema explícito
Use el esquema explícito para resolver el problema tipo, usando mallados
correspondientes. Compare los resultados obtenidos de la aproximación con
la solución exacta de este problema en los niveles de tiempo indicados:
28
a). n=19, m=39 y en el nivel de tiempo t = 0,5.
b). n=39, m=99 y en el nivel de tiempo t = 2.
Para resolver este problema, es necesario recurrir a la programación; es
decir, se va a desarrollar un programa computacional (en Matlab) que permita
hacer diferentes corridas para poder visualizar de mejor manera los resultados
de la solución aproximada U(x, t) obtenida.
Representemos el PV F dado, de la siguiente manera:
Ut(x, t) = c
2Uxx, para 0 < x < a y 0 < t ≤ b, (2.41)
U(x, 0) = f(x) para 0 < x < a y t = 0, (2.42)
U(0, t) = g1(t) para x = 0 y 0 ≤ t ≤ b, (2.43)
U(a, t) = g2(t) para x = a y 0 ≤ t ≤ b. (2.44)
La discretización de este problema se realiza mediante el mallado de la Figura
2.8, en donde n es el número de puntos en que se divide la variable espacial a
ym el número de puntos en que se divide b, espaciados h y k respectivamente.
Figura 2.8: Mallado del problema tipo.
29
El cálculo de los valores de la solución aproximada U(x, t) se realiza de la
siguiente manera:
Cálculo de la primera columna
Para t1 = 0; es decir si j = 1, se aplica la condición inicial dada por la
ecuación (2.42)
U(i, 1) = f(xi), i = 2, 3, ..., (n− 1). (2.45)
Cálculo de la segunda columna y las restantes
Las fórmulas en diferencias finitas para las derivadas primera y segunda
en t y x respectivamente son:
Ut(i, j) =
U(i, j + 1)− U(i, j)
k
+O(k) (2.46)
Uxx(i, j) =
U(i− 1, j)− 2U(i, j) + U(i+ 1, j)
h2
+O(h2) (2.47)
Remplazando las ecuaciones (2.46) y (2.47) en la ecuación (2.41); y,
despreciando el error de truncamiento se tiene:
U(i, j + 1)− U(i, j)
k
= c2
U(i− 1, j)− 2U(i, j) + U(i+ 1, j)
h2
Ahora, definiendo r = kc
2
h2
, en la expresión anterior, se tiene
U(i, j + 1)− U(i, j) = rU(i− 1, j)− 2rU(i, j) + rU(i+ 1, j).
Simplificando y ordenando esta última expresión
U(i, j+1) = rU(i−1, j)+(1−2r)U(i, j)+rU(i+1, j), i = 2, 3, ..., (n−1).
(2.48)
La ecuación (2.48) es la ecuación en diferencias explícita y se emplea
para calcular las aproximaciones en la columna (j + 1)_ésima de la ma-
lla, a partir de las aproximaciones de la columna j_ésima, permitiéndo
avanzar en el tiempo.
Programando las ecuaciones (2.45) y (2.48), en Matlab, se obtienen el pro-
grama que se encuentra en la tabla de algoritmos con el nombre Algoritmo 2,
problema tipo, método explícito (ver tabla de algoritmos).
30
Es importante indicar la forma de ingresar los datos para ejecutar el pro-
grama en mención. La función, en Matlab, que ejecuta el programa tiene la
forma:
>> U=Problema_1_Explicito(f,g1,g2,a,b,c,n,m);
en donde: f es la función f(x) dada de la condición inicial, g1 y g2 son los
valores de las funciones de las condiciones de frontera, a es la magnitud de
la variable espacial x, b es el tiempo a considerar, c2 es la constante térmica
de la ecuación diferencial, n es el número de puntos en que se divide x y m el
número de puntos en que se divide el tiempo b.
Ahora, ejecutando el programa para los datos de este problema: fx =
sen(πx) · (1 + 2cos(πx)), g1(t) = g2(t) = 0, a = 1 y c = 1.
Para n = 19, m = 39 y t = 0,5
>> U=Problema_1_Explicito(’sin(pi*x)*(1+2*cos(pi*x))’,’0’,’0’,1,
0.5,1,19,39)
La matriz U resultante, evaluada en cada uno de los puntos de la malla, se en-
cuentra en la tabla de matrices con el nombre Solución aproximada, problema
tipo, método explícito (ver tabla de matrices).
La gráfica de la solución aproximada U(x, t), de este problema, se repre-
senta en la Figura 2.9.
Figura 2.9: Solución aproximada. Problema tipo, explícito
31
Observando, por simple inspección (Figura 2.10), las soluciones aproxima-
da y exacta son completamente diferentes. Pero ¿Que está pasando?
Figura 2.10: Soluciones aproximada y exacta. Problema tipo, comparación.
Recordemos que en el punto 2.1.1, de este trabajo, se estudió la estabi-
lidad del método explícito y se llegó a concluir que este método es estable
(ecuación (2.20)) si
r ≤ 1
2
.
El programa ejecutado nos proporciona el valor de r para este problema y,
según se puede ver en la corrida de este programa
r = 4,2632 >
1
2
Por tanto, la solución aproximada contendrá una oscilación divergente; es decir
esta solución crece rápidamente cuando t → m, produciéndose un efecto de
inestabilidad numérica. Esto se puede apreciar en el gráfico de esta solución
(Figura 2.9) y en los valores muy grandes de la matriz solución aproximada,
columnas 37, 38 y 39
También se dijo que para que se produzca esta inestabilidad numérica,
el paso k de la variable temporal no es demasiado pequeño. El conseguir la
estabilidad del método producirá un costo computacional muy alto.
Ahora, se va a conseguir la estabilidad del método explícito. Para esto, se
debe cumplir
r =
c2k
h2
≤ 1
2
,
32
de la cual,
k ≤
1
2
h2
c2
,
es decir
b
m− 1
≤
1
2
h2
c2
.
Despejando m de la última expresión
m ≥ bc
2
1
2
h2
+ 1.
Por tanto, con los datos del problema se tiene
m ≥ 325.
Es por el valor de m, muy grande, que se concluye que el costo computacional
es muy alto. Nuevamente se realiza la corrida de este ejercicio, tomando el
valor de m = 339.
La solución U(x, t), para m = 339, se encuentra en la tabla de matrices
con el nombre Solución aproximada, problema tipo, método explícito estable.
El tiempo requerido por este programa, para encontrar la solución aproxi-
mada U(x, t), en segundos, es
tiempo =
0.3552
>>
Esto permite concluir que el costo computacional es alto.
La gráfica de la solución aproximada, bajo estas nuevas consideraciones,
se representa en la figura 2.11.
Como se puede observar, la solución aproximada encontrada se aproxima
mucho a la solución exacta representada matricialmente y graficada, ya que
se ha conseguido un valor de
r = 0,4793 <
1
2
.
33
Figura 2.11: Solución aproximada. Problema tipo explícito, estable.
Ahora, para responder el literal b), es decir para n = 39 m = 99 y t = 2,
se tiene un valor de r = 29,17 > 1
2
; por lo que el método explícito es inesta-
ble. Con la finalidad de conseguir la estabilidad de este método se procede a
calcular el valor de m adecuado, de manera que se cumpla con el criterio de
estabilidad del método explícito r ≤ 1
2
.
Del cálculo realizado en el literal a), de este problema
m ≥ bc
2
1
2
h2
+ 1.
Con los datos del literal b), se tiene
m ≥ 5776.
Tomando m = 5776 se procede a ejecutar el programa del algoritmo 1.
La gráfica de la solución aproximada U , del problema en estudio, bajo
estas nuevas consideraciones se representa en la Figura 2.12.
Como se puede observar, el valor de m es demasiado grande; es decir, el
espaciamiento de la variable temporal k es demasiado pequeño y por tanto el
costo computacional muy alto.
Es necesario por tanto, para reducir este costo, encontrar un nuevo método
aproximado de solución. El método que se propone a continuación desarrollar
es el de Crank-Nicolson, explicado en la sección 2.2 de este trabajo.
34
Figura 2.12: Solución aproximada. Problema tipo, m = 5776.
2.4.3. Problema tipo. Esquema de Crank - Nicolson
Use el esquema de Crank - Nicolson para resolver el problema tipo, usan-
do mallados correspondientes. Compare los resultados obtenidos de la apro-
ximación con la solución exacta de este problema en los niveles de tiempo
indicados:
a). n=19, m=39 y en el nivel de tiempo t = 0,5.
b). n=39, m=99 y en el nivel de tiempo t = 2.
Para resolver este PVF, es necesario desarrollar el método de Crank-Nicolson
y llevarlo a la programación.Como se explicó anteriormente, en el punto 2.2.2 de este trabajo y, toman-
do en consideración el mallado de la Figura 2.8, la diferencia progresiva en el
tiempo, en el nivel j,es
Ut(i, j) =
U(i, j + 1)− U(i, j)
k
+O(k); (2.49)
y la diferencia centrada en el espacio es
Uxx(i, j) =
U(i− 1, j)− 2U(i, j) + U(i+ 1, j)
h2
+O(h2). (2.50)
35
Remplazando las ecuaciones (2.49) y(2.50) en la ecuación (2.41), consideran-
do c = 1 y despreciando el error de truncamiento, se tiene
U(i, j + 1)− U(i, j)
k
=
U(i− 1, j)− 2U(i, j) + U(i+ 1, j)
h2
. (2.51)
En forma análoga, en el nivel j + 1, para las diferencias regresiva y centrada
se tiene
Ut(i, j + 1) =
U(i, j + 1)− U(i, j)
k
+O(k), (2.52)
Uxx(i, j + 1) =
U(i− 1, j + 1)− 2U(i, j + 1) + U(i+ 1, j + 1)
h2
+O(h2), (2.53)
respectivamente. Remplazando las ecuaciones (2.52) y(2.53) en (2.41), con la
misma consideración para c y para el error de truncamiento
U(i, j + 1)− U(i, j)
k
=
U(i− 1, j + 1)− 2U(i, j + 1) + U(i+ 1, j + 1)
h2
(2.54)
Ahora, para cumplir con el esquema de Crank-Nicholson se deben sumar las
ecuaciones (2.51) y (2.54)
2
U(i, j + 1)− U(i, j)
k
=
U(i− 1, j)− 2U(i, j) + U(i+ 1, j)
h2
+
U(i− 1, j + 1)− 2U(i, j + 1) + U(i+ 1, j + 1)
h2
.
Haciendo r = k
h2
en la última expresión, dividiendo para r y simplificando se
tiene:
U(i− 1, j + 1) + 2 + 2r
r
U(i, j + 1)− U(i+ 1, j + 1) = U(i− 1, j) (2.55)
+
2− 2r
r
U(i, j)
+U(i+ 1, j),
para i = 2, 3, ..., (n− 1).
Por tanto, la ecuación (2.55) permite, a partir de los valores de la solución
U del nivel j, conocer las incógnitas de la solución U en el nivel j+1, formando,
para cada columna de cálculo, un sistema de ecuaciones.
Cálculo de la primera columna de la representación matricial del
sistema.
36
Para t = 0, es decir si j = 1 y para 0 < xi < 1, i = 2, 3, ...(n−1), aplicamos
la condición inicial, dada por la ecuación (2.42)
U(i, 1) = sen(πxi)(1 + 2cos(πxi)) = f(xi).
Cálculo de la primera y última filas de la representación matricial
del sistema.
Se aplican las condiciones de frontera para j = 1, 2, ...,m
U(1, j) = 0
U(n, j) = 0
Cálculo de las filas segunda en adelante de la representación ma-
tricial del sistema.
En la ecuación (2.55), para valores de i = 2, 3, ...(n− 1), haciendo
S1 =
2+2r
r
, S2 = 2−2rr ; y, considerando que
U(1, j + 1) = 0, U(1, j) = 0, U(n, j + 1) = 0 y U(n, j) = 0, se tiene:
Para i = 2
−U(1, j + 1) + S1U(2, j + 1)− U(3, j + 1) = U(1, j) + S2U(2, j) + U(3, j)
Para i = 3
−U(2, j + 1) + S1U(3, j + 1)− U(4, j + 1) = U(2, j) + S2U(3, j) + U(4, j)
Para i = 4
−U(3, j + 1) + S1U(4, j + 1)− U(5, j + 1) = U(3, j) + S2U(4, j) + U(5, j)
Para i = 5
−U(4, j + 1) + S1U(5, j + 1) − U(6, j + 1) = U(4, j) + S2U(5, j) + U(6, j)
En forma análoga, continuando con este proceso, para i = n− 1
−U(n− 2, j + 1) + S1U(n− 1, j + 1)− U(n, j + 1) = U(n− 2, j)
+S2U(n− 1, j)
+U(n, j)
37
Las ecuaciones anteriormente obtenidas forman el sistema matricial de ecua-
ciones, con las incógnitas U(2, j+1), U(3, j+1), U(4, j+1), U(5, j+1), ..., U(n−
1, j + 1), representado en el cuadro 2.1.
Para la programación de este método, se tomaron las siguientes conside-
raciones:
Las incógnitas U(i, j + 1), i = 2, 3, ..., (n − 1) de la columna j + 1 se
encuentran resolviendo el sistema matricial de ecuaciones, en el que los
términos independientes están formados por los valores conocidos de la
función solución U de la columna anterior j. Las incógnitas de la colum-
na j + 2 se obtendrán resolviendo el sistema que tiene como términos
independientes a los valores de la solución U de la columna j + 1. Es
decir, se debe resolver el sistema para cada columna de incógnitas que
se quiera encontrar.
El sistema que se forma siempre es un sistema tridiagonal
Sobre la base de lo anterior indicado, el programa desarrollado para el proble-
ma en estudio, se encuentra en la tabla de algoritmos con el nombre Algoritmo
3, problema tipo, Crank - Nicolson.
En primer lugar se va a responder el literal a), del Problema tipo, aplicando
el método de Crank-Nicolson; es decir:
Para t = 0,50, n = 19 y m = 39
Ejecutando el programa, para los datos del problema, se tiene la solución
U(x, t) que se encuentra en la tabla de matrices con el nombre Solución apro-
ximada, problema tipo, literal a).CrankNicolson.
El tiempo requerido por este programa, para encontrar la solución aproxi-
mada U(x, t), en segundos, es
tiempo =
0.3045
>>
La gráfica de la solución aproximada U, del problema en estudio, con las con-
sideraciones señaladas se encuentra en la Figura 2.13.
38
Cuadro 2.1: Sistema matricial, problema tipo, Crank - Nicolson
                         1
0
0
0
0
0
0
0
0
0
0
S
1
−
1
0
0
0
0
0
0
0
0
−
1
S
1
−
1
0
0
0
0
0
0
0
0
−
1
S
1
−
1
0
0
0
0
0
0
0
0
−
1
S
1
−
1
0
0
0
0
0
0
0
0
−
1
S
1
−
1
0
0
0
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
0
0
0
0
0
0
−
1
S
1
−
1
0
0
0
0
0
0
0
0
−
1
S
1
0
0
0
0
0
0
0
0
0
0
1
                                                  
U
(1
,j
+
1)
U
(2
,j
+
1)
U
(3
,j
+
1)
U
(4
,j
+
1)
U
(5
,j
+
1)
U
(6
,j
+
1)
. . . U
(n
−
2,
j
+
1)
U
(n
−
1,
j
+
1)
U
(n
,j
+
1)
                         =
                         
g 1
(t
1
)
S
2
U
(2
,j
)
+
U
(3
,j
)
U
(2
,j
)
+
S
2
U
(3
,j
)
+
U
(4
,j
)
U
(3
,j
)
+
S
2
U
(4
,j
)
+
U
(5
,j
)
U
(4
,j
)
+
S
2
U
(5
,j
)
+
U
(6
,j
)
U
(5
,j
)
+
S
2
U
(6
,j
)
+
U
(7
,j
) . . .
U
(n
−
4,
j)
+
S
2
U
(n
−
3,
j)
+
U
(n
−
2,
j)
U
(n
−
3,
j)
+
S
2
U
(n
−
2,
j)
+
U
(n
−
1,
j)
U
(n
−
2,
j)
+
S
2
U
(n
−
1,
j)
                         
39
Figura 2.13: Solución aproximada. Problema tipo, literal a), Crank-Nicolson
Se puede observar que la solución aproximada U(x, t), obtenida con el
método de Crank-Nicolson, es bastante aproximada a la solución exacta en-
contrada en el Problema tipo (Figura 2.7), de este trabajo. En el Capítulo 5 se
realizará el cálculo del error cometido en la aproximación.
Además, podemos comprobar lo anotado anteriormente, que este método
es estable para cualquier valor de r (r = 4,2632).
Para n = 39, m = 99 y t = 2 Aplicando el esquema de Crank-Nicolson y
ejecutando el programa para los datos de este problema, se tiene la solución
U(x, t) que se encuentra en la tabla de matrices con el nombre Solución apro-
ximada, problema tipo, literal b), Crank −Nicolson.
La gráfica de la solución aproximada U del problema en estudio, bajo estas
nuevas consideraciones se representa en la Figura 2.14.
Ahora, se debe conocer la solución exacta y su representación gráfica
correspondiente; para lo cual, se debe ejecutar el programa Algoritmo 2, para
los datos de este problema.
La solución exacta U(x, t) se encuentra en la tabla de matrices con el
nombre Solución exacta, problema tipo, literal b).
40
Figura 2.14: Solución aproximada. Problema tipo, literal b).Crank −Nicolson
La representación gráfica de la solución exacta, para los datos del proble-
ma en estudio, se representan en la Figura 2.15.
Figura 2.15: Solución exacta. Problema tipo, literal b).
Comparando los valores de las soluciones aproximada y exacta, se puede
concluir que con el método de Crank-Nicolson se tiene una buena aproxima-
ción a la solución exacta.
2.4.4. Problema 2 (No Homogéneo)
Use el esquema de Crank - Nicolson para resolver la siguiente EDP, no
homogénea, usando mallados correspondientes, si:
Ut(x, t) = Uxx(x, t) + q(x, t), 0 < x < 1 y 0 < t <∞ y q(x, t) = 10x+ t.
41
con la condición inicial de temperatura
U(x, 0) = sen(πx) · (1 + 2cos(πx)) para t = 0 y 0 < x < 1,
y con las condiciones de frontera (o de contorno) a considerar
U(0, t) = 0 para x = 0 y 0 < t <∞,
U(1, t) = −t para x = 1 y 0 < t <∞,
Desarrollo de un programa computacional para PVF no Homogéneos.
Los diferentes problemas unidimensionalesvarían por lo general por las con-
diciones inicial, de frontera y por la EDP misma; dando lugar a un PVF No
Homogéneo.
Es necesario por tanto desarrollar un programa computacional que resuel-
va algunas formas de la EDP, con sus condiciones de frontera e inicial. Se
aplicará, por tanto, el esquema de Crank-Nicolson, para resolver un PVF no
homogéneos.
Sea
Ut(x, t) = c
2Uxx(x, t) + q(x, t), 0 < x < a, 0 < t ≤ b. (2.56)
Con la condición inicial
U(x, 0) = f(x), t = 0, 0 < x < a, (2.57)
y, con las condiciones de frontera
α1U(0, t) + β1Ut(0, t) = g1(t), x = 0, 0 < t ≤ b, (2.58)
α2U(a, t) + β2Ut(a, t) = g2(t), x = a, 0 < t ≤ b. (2.59)
Siendo: c, α1, α2, β1, β2 constantes numéricas, q(x, t) una fuente externa, f(x)
una función dada por la condición inicial y g1(t), g2(t) funciones de t para las
condiciones de frontera.
Considerando el mallado de la Figura 2.8 y aplicando la diferencia progre-
siva, en el tiempo, para el nivel j
Ut(i, j) =
U(i, j + 1)− U(i, j)
k
+O(k), (2.60)
42
y la diferencia centrada para la variable espacial, en el mismo nivel de tiempo
Uxx(i, j) =
U(i− 1, j)− 2U(i, j) + U(i+ 1, j)
h2
+O(h2) (2.61)
Se tiene remplazando las ecuaciones (2.60) y (2.61) en (2.56), despreciando
el error de truncamiento
U(i, j + 1)− U(i, j)
k
= c2
U(i− 1, j)− 2U(i, j) + U(i+ 1, j)
h2
+ q(xi, tj) (2.62)
Ahora, aplicando la diferencia regresiva en el tiempo y centrada en la variable
espacial, para el nivel j + 1, se tiene
Ut(i, j + 1) =
U(i, j + 1)− U(i, j)
k
+O(k) (2.63)
Uxx(i, j + 1) =
U(i− 1, j + 1)− 2U(i, j + 1) + U(i+ 1, j + 1)
h2
+O(h2) (2.64)
Remplazando las ecuaciones (2.63) y(2.64) en (2.56) y realizando la misma
consideración para el error de truncamiento, se tiene
U(i, j + 1)− U(i, j)
k
=
U(i− 1, j + 1)− 2U(i, j + 1) + U(i+ 1, j + 1)
h2/c2
+q(xi, tj+1) (2.65)
Para cumplir con el esquema de Crank-Nicolson sumamos las ecuaciones
(2.62) y (2.65)
2
U(i, j + 1)− U(i, j)
k
= c2
U(i− 1, j)− 2U(i, j) + U(i+ 1, j)
h2
+c2
U(i− 1, j + 1)− 2U(i, j + 1) + U(i+ 1, j + 1)
h2
+q(xi, tj) + q(xi, tj+1)
En la expresión anterior haciendo r = c
2k
h2
, se tiene
−rU(i− 1, j + 1) + (2 + 2r)U(i, j + 1)− rU(i+ 1, j + 1) = rU(i− 1, j)
+(2− 2r)U(i, j)
+rU(i+ 1, j)
+k[q(xi, tj)
+q(xi, tj+1)]
43
Dividiendo para r en la última expresión
−U(i− 1, j + 1) + (2
r
+ 2)U(i, j + 1)− U(i+ 1, j + 1) = U(i− 1, j)
+(
2
r
− 2)U(i, j)
+U(i+ 1, j)
+
k
r
[q(xi, tj)
+q(xi, tj+1)]
Ahora, llamemos S1 = 2r + 2, S2 =
2
r
− 2 y S3 = kr y remplazando en la última
expresión se tiene
−U(i− 1, j + 1) + S1U(i, j + 1)− U(i+ 1, j + 1) = U(i− 1, j) (2.66)
+S2U(i, j) + U(i+ 1, j)
+S3[q(xi, tj)
+q(xi, tj+1)],
para i = 2, 3, ...(n− 1)
Cálculo de la primera columna de la representación matricial del
sistema
Para t = 0, es decir si j = 1 y para 0 < xi < a, i = 2, 3, ...(n − 1),
aplicamos la condición inicial dada por la ecuación (2.57). Se tiene
U(i, 1) = f(xi)
Cálculo de la primera y última filas de la representación matricial
del sistema
Para este problema, las condiciones de frontera izquierda y derecha, en
forma general se pueden escribir, respectivamente
α1U(i, j) + β1Ut(i, j) = g1(tj) (2.67)
α2U(i, j) + β2Ut(i, j) = g2(tj) (2.68)
Ahora, se desarrollará la discretización de estas condiciones de frontera:
44
Para la frontera izquierda (ecuación (2.67)), aplicando las diferencias
progresivas a Ut se tiene
α1U(1, j) + β1
(
U(1,j+1)−U(1,j)
k
)
= g1(tj);
es decir
kα1 + β1U(1, j + 1)− β1U(1, j) = kg(tj).
De esta última expresión se tiene
U(1, j + 1) = −(kα1 − β1)
β1
U(1, j) +
k
β1
g1(tj). (2.69)
Ahora, en forma análoga para la frontera derecha (ecuación (2.68)), aplica-
cando diferencias regresivas a Ut, se tiene
α2U(n, j) + β2
(
U(n,j)−U(n,j−1)
k
)
= g2(tj)
De esta última expresión
U(n, j + 1) =
(
β2
kα2 + β2
)
U(n, j) +
(
k
kα2 + β2
)
g2(tj+1) (2.70)
Cálculo de las filas segunda en adelante de la representación ma-
tricial del sistema
De la ecuación (2.66), para valores de i = 2, 3, ...(n− 1), se tiene:
Para i = 2
−U(1, j + 1) + S1U(2, j + 1)− U(3, j + 1) = U(1, j) + S2U(2, j)
+U(3, j) + +S3[q(x2, tj)
+q(x2, tj+1)]
Remplazando la ecuación (2.69) en la última expresión y simplificando
se tiene
S1U(2, j + 1)− U(3, j + 1) =
(
1− kα1 − β1
β1
)
U(1, j) + S2U(2, j)
+U(3, j) + S3[q(x2, tj) + q(x2, tj+1)]
+
k
β1
g1(tj)
Para i = 3
45
−U(2, j + 1) + S1U(3, j + 1)− U(4, j + 1) = U(2, j) + S2U(3, j)
+U(4, j) + S3[q(x3, tj)
+q(x3, tj+1)]
Para i = 4
−U(3, j + 1) + S1U(4, j + 1)− U(5, j + 1) = U(3, j) + S2U(4, j)
+U(5, j) + S3[q(x4, tj)
+q(x4, tj+1)]
Para i = 5
−U(4, j + 1) + S1U(5, j + 1)− U(6, j + 1) = U(4, j) + S2U(5, j)
+U(6, j) + S3[q(x5, tj)
+q(x5, tj+1)]
En forma análoga, continuando con este proceso hasta i = n − 1, se
tiene
−U(n− 2, j + 1) + S1U(n− 1, j + 1)− U(n, j + 1) = U(n− 2, j)
+S2U(n− 1, j)
+U(n, j)
+S3[q(xn−1, tj)
+q(xn−1, tj+1)]
Remplazando la ecuación (2.70) en la última expresión nos queda
−U(n− 2, j + 1) + S1U(n− 1, j + 1) = U(n− 2, j) + S2U(n− 1, j)
+
(
1 +
β2
kα2 + β
)
U(n, j)
+S3[q(xn−1, tj) + q(xn−1, tj+1)]
+
(
k
kα2 + β2
)
g2(tj+1)
Las ecuaciones anteriormente obtenidas, para i = 2, 3, ...(n − 1), forman el
sitema matricial de ecuaciones, representado en el cuadro 2.2, con las incóg-
nitas:
46
Cuadro 2.2: Sistema matricial, problema 2 no homogéneo.
                         1
0
0
0
0
0
0
0
0
0
0
S
1
−
1
0
0
0
0
0
0
0
0
−
1
S
1
−
1
0
0
0
0
0
0
0
0
−
1
S
1
−
1
0
0
0
0
0
0
0
0
−
1
S
1
−
1
0
0
0
0
0
0
0
0
−
1
S
1
−
1
0
0
0
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
0
0
0
0
0
0
−
1
S
1
−
1
0
0
0
0
0
0
0
0
−
1
S
1
0
0
0
0
0
0
0
0
0
0
1
                                                  
U
(1
,j
+
1)
U
(2
,j
+
1)
U
(3
,j
+
1)
U
(4
,j
+
1)
U
(5
,j
+
1)
U
(6
,j
+
1)
. . . U
(n
−
2,
j
+
1)
U
(n
−
1,
j
+
1)
U
(n
,j
+
1)
                         =
                         
−
k
α
1
−
β
1
β
1
U
(1
,j
)
+
k β
1
g 1
(t
j
)
S
4
U
(1
,j
)
+
S
2
U
(2
,j
)
+
U
(3
,j
)
+
S
3
[q
(x
2
,t
j
)
+
q(
x
2
,t
j+
1
)]
+
k β
1
g 1
(t
j
)
U
(2
,j
)
+
S
2
U
(3
,j
)
+
U
(4
,j
)
+
S
3
[q
(x
3
,t
j
)
+
q(
x
3
,t
j+
1
)]
U
(3
,j
)
+
S
2
U
(4
,j
)
+
U
(5
,j
)
+
S
3
[q
(x
4
,t
j
)
+
q(
x
4
,t
j+
1
)]
U
(4
,j
)
+
S
2
U
(5
,j
)
+
U
(6
,j
)
+
S
3
[q
(x
5
,t
j
)
+
q(
x
5
,t
j+
1
)]
U
(5
,j
)
+
S
2
U
(6
,j
)
+
U
(7
,j
)
+
S
3
[q
(x
6
,t
j
)
+
q(
x
6
,t
j+
1
)] . . .
U
(n
−
3,
j)
+
S
2
U
(n
−
2,
j)
+
U
(n
−
1,
j)
+
S
3
[q
(x
n
−
2
,t
j
)
+
q(
x
n
−
2
,t
j+
1
)]
U
(n
−
2,
j)
+
S
2
U
(n
−
1,
j)
+
S
5
U
(n
,j
)
+
S
3
[q
(x
n
−
1
,t
j
)
+
q(
x
n
−
1
,t
j+
1
)]
+
S
6
g 2
(t
j+
1
)
β
2
k
α
2
+
β
2
U
(n
,j
)
+
k
k
α
2
+
β
2
g 2
(t
j+
1
)
                         
47
U(1, j + 1), U(2, j + 1), U(3, j + 1), U(4, j + 1), U(5, j + 1), ..., U(n − 1, j +
1), U(n, j + 1),
Para la programación de este método aproximado, se tomaron las siguien-
tes consideraciones:
Las incógnitas U(i, j + 1), i = 2, 3, ..., (n − 1) de la columna j + 1 se
encuentran resolviendo el sistema matricial de ecuaciones, en el que los
términos independientes están formados por los valores conocidos de la
función solución U de la columna anterior j.
Las incógnitas de la columna j + 2 se obtendrán resolviendo el sistema
que tiene como términos independientes a los valores de la solución U
de la columna j+1. Es decir, se debe resolver un sistema de ecuaciones
para cada columna de incógnitas que se quiera encontrar.
Los valores de U(1, j+1) y U(n, j+1) corresponden a las fórmulas (2.69)
y (2.70), respectivamente.
La programación se realizará aplicando comandos de Matlab.
El programa que se desarrolló, para el PVF planteado, se encuentra en la tabla
de algoritmos con el nombre Algoritmo 4, problema 2 no homogéneo, Crank