1

Solución Numérica de Ecuaciones Diferenciales

Parciales
Manuela Bastidas* and Freddy H. Marín**

Abstract—Durante los últimos años los problemas de ingeniería finitas, las derivadas que son parte fundamental de un modelo
de todas las índoles se han ido complicando de manera exponen- en forma de EDP son remplazadas por aproximaciones en
cial y cada vez se dificulta más encontrar soluciones analíticas diferencias, convirtiendo el planteamiento en un problema
para éstos; los métodos numéricos son entonces la manera de
encontrar una solución teórica y práctica a problemas más algebraico.
complejos, con condiciones y geometrías tan variadas como se Son de nuestro interés avances y desarrollos en los plan-
quiera. En éste trabajo se muestran dos aplicaciones del método teamientos matemáticos de las diferencias finitas y sobreto-
de diferencias finítas explicitas para encontrar aproximaciones a do en la computación y la complejidad algorítmica de los
la solución, pero además se realiza el análisis de las características mismos, basados en trabajos que se centran en solucionar
del problema y del esquema numérico que garantizan una
correcta implementación. ecuaciones diferenciales particulares. Ejemplos de lo aterior
son la ecuación de Helmholtz que Erlangga,Y. solocionó
Index Terms—Diferencias finitas, Esquema numérico, Consis- proponiendo un esquema numérico de diferencias finitas de
tencia, Estabilidad, Convergencia, Ecuación KdV, Ecuación de
Schrödinger. orden superior con dominios externos y haciendo un aná-
lisis estricto de las condiciones de frontera necesarias para
lograr una discretización estable, y donde además se hace
I. I NTRODUCTION una transformación y simplificiación del modelo absorbiendo
NA cantidad significativa de problemas de orden cien- propiedades y considerando condiciones de radiación. Por otra
U tífico y tecnológico se nos presentan diariamente y la
modelización de fenómenos es precisamente el área que se
parte, uno más de los trabajos que es interesante resaltar, y en
este caso por razones de interés general en el problema, es la
encarga de encontrar las soluciones a dichos problemas, eso publicación de Atanassov, E.I y su desarrollo numérico para
hace que ésta sea rama parte vital de la ciencia actual y es solucionar el sistema de ecuaciones estocásticas planteadas
aquí donde son útiles las ecuaciones diferenciales parciales, por Heston para calcular el precio de una opción y por otra
modelan sistemas de evolución en los que se describe una parte, el trabajo de Gandarias, M.L que por medio de modelos
dinámica a lo largo del tiempo y además pueden representar matematicamente simples establece y generaliza algunas ca-
diversos objetos que van desde la posición de un satélite en rácteristicas y particularidades de la ecuación KdV (Korteweg-
el espacio hasta la dinámica de un átomo, pasando por los de Vries), ecuación que físicamente describe la propagación
índices bursátiles o el grado en que una enfermedad afecta a de ondas en medios dispersivos como la superficie del agua
la población. en canales poco profundos, este trabajo está ampliamente
La complejidad de los modelos que actualemente se han relacionado con los desarrollos hechos por Duruflé, M.a. que
propuesto para reproducir la concepción humana de los siste- utiliza e implementa el método de diferencias finitas para
mas dinámicos y que son descritos por una cantidad diversa solucionar la ecuacion KdV y propone un esquema numérico
de EDP’S (Ecuaciones diferenciales parciales) hace que la que considere cualidades físicas como la conservación de
solución de dichos problemas sea analíticamente más com- energía y además hace una valiosa comparación de dicho
plicada cada vez, además, la diversidad de los problemas esquema con un planteamiento numérico basado en el método
obliga a que los métodos numéricos que encuentran soluciones de Galerkin discontinuo que es usual en estudios relacionados
aproximadas a los modelos sean más utilizados, tanto por la a los elementos finítos y otros métodos más sofisticados con
facilidad computacional como teórica y además, la precisión fuertes bases del análisis numérico.
de dichos métodos ha ido aumentando significativamente de la Éste trabajo pretende estudiar el método de diferencias
mano de estudios que hacen que estos sean también confiables finítas explícitas y su posible aplicación en la solución de
dadas ciertas características matemáticas específicas. ecuaciones diferenciales parciales, realizar un par de im-
En particular, el método de diferencias finítas es un clásica plementaciones computacionales de dicho método numérico
aproximación para encontrar la solución numérica de las para resolver modelos típicos en la literatura y analizar las
ecuaciones que plantean un modelo matemático de un sistema carácteristicas teóricas de los esquemas propuestos.
continuo, es decir, se plantea un problema aproximado que Luego de hacer una revisión del estado del arte y estu-
tiene características discretas. En una solución por diferencias diar el método de diferencias finitas explícitas para resolver
ecuaciones diferenciales parciales, se proponen un par de
*Manuela Bastidas- Departamento de Ciencias Básicas, Universidad problemas que es posible resolver utilizando discretizaciones
EAFIT, Medellín, Colombia, e-mail: mbastida@eafit.edu.co.
**Freddy H. Marín- Departamento de Ciencias Básicas, Universidad simples, hacer un análisis de las caracteristicas geométicas y
EAFIT, Medellín, Colombia, e-mail: fmarinsa@eafit.edu.co. matemáticas de los modelos, crear esquemas numéricos para
 2

lograr resultados significativamente buenos, analizar todas las interés particular de la ingeniería y la modelación ma-
condiciones del esquema numérico: consistencia, estabilidad, temática. De manera que las ecuaciones que el modelo
convergencia, entre otras e implementar computacionalmente que se plantea se aproxime correctamente a la realidad.
el método de diferencias finítas para resolver los modelos Primero es importante introducir el concepto de aproxi-
propuestos y analizar los resultados encontrados apartir de las mación que implementaremos y el método de diferencias
simulaciones. finitas donde mediante un proceso de discretización, el
Por otra parte, el alcance de éste trabajo como se mencionó conjunto infinito de números que representan la fúncion
anteriormente se basa en la importancia que desde hace varios o funciones incógnitas en el continuo es reemplazado
años la modelización de sistemas dinámicos ha tomado en por un número finito de parámetros, para hacer una
las ciencias actuales, lo que implica la estructuración fuerte aproximación útil del método se emplean las diferencias
de ecuaciones diferenciales que abarquen en gran parte el finitas que son expresiones matemáticas de la forma
comportamiento y cambios de estado de los modelos y sus
sistemas asociados. f (x + b) − f (x + a)
La solución de dichas ecuaciones diferenciales al igual que
Ahora bien, si una diferencia finita se divide por (b − a)
su planteamiento, progresivamente se ha hecho más compleja
se obtiene una expresión similar al cociente diferencial,
y es en este punto donde las soluciones numéricas se hacen
que se diferencia en que se emplean cantidades finitas en
fundamentales y los métodos sencillos que implican bajos
lugar de infinitesimales, es decir, las diferencias finitas
costos computacionales toman fuerza, siendo el caso de las
se pueden orientar como un método numérico de aproxi-
diferencias finitas, que durante muchos años han sido la forma
mación a las derivadas. Normalmente se consideran solo
de aproximarse a la solución exacta de una manera flexible y
tres formas de diferencias finitas: la anterior, la posterior
cómoda.
y la central.
Desde la concepción de ingeniería matemática se entiende
la importancia de los problemas de modelización aplicados a a) Una diferencia progresiva, adelantada o posterior
los fenomenos naturales, financieros, físicos, biológicos entre es una expresión de la forma
otros que son de nuestro interés. Para el desarrollo de éste
trabajo son útiles e importantes las áreas de modelación y la ∆f (x) = f (x + h) − f (x)
fuerte componente matemática que implica conocimiento del
b) Una diferencia regresiva, atrasada o anterior es de
planteamiento del problema, significado y solución de EDP’s,
la forma
además de nociones acerca de la discretización de problemas
de ingeniería ya sea con fundamentos básicos del cálculo ∆f (x) = f (x) − f (x − h)
diferencial o más profundamente como con la simulación de
sistemas dinámicos discretos y por último un alto componente c) Una diferencia central es la media de las diferen-
computacional que permita la modelación y simulación de cias anteriores y posteriores y está dada por
sistemas y soluciones numéricas de ésta índole.
Luego, el interés de éste trabajo más que en encontrar f (x + h) − f (x − h)
∆f (x) =
soluciones numéricas está en las condiciones naturales de 2
las ecuaciones, es decir, en el análisis de las ecuaciones, el Diferencias finitas para resolver ecuaciones diferencia-
esquema numérico y sus caracteristicas particulares. les: Para resolver ecuaciones diferenciales con condi-
Se considerarán algunos problemas y se espera solucionar- ciones iniciales por el método de diferencias finitas se
los numéricamente, en particular de estudiará el modelo plan- diferencia la variable independiente x con dominio [0,L]
teado por la ecuación KdV que modela el comportamiento de , es decir, se construye un conjunto (grilla o malla)
solitones en medios no lineales y la ecuación de Schrodinger de L+1 puntos discretos igualmente espaciados xl (l =
en una caja de potencial de la que resulta la amplitud de 0, 1, ...L) , con x0 = 0, xL = L, y xl+1 − xl = ∆x.
probabilidad de presencia para una partícula, ecuaciones que Luego se reemplazan aquellos términos de la ecuación
son de nuestro interés porque la cantidad de aplicaciones de diferencial que involucren diferenciación, por términos
éstas cuando se hace una correcta simulación en ciencia y en que contengan operaciones algebraicas. Este proceso
ingeniería son inimaginables para hacer desarrollos en base a trae implícito una aproximación y puede efectuarse
las aproximaciones de la solución. mediante la utilización de aproximaciones en diferencias
finitas para las derivadas de la función.
Aproximaciones a la derivada con diferencias finitas:
II. M ETODOLOGÍA
Para hacer las primeras aproximaciones a derivadas con
Para lograr el alcance de los objetivos que se han propuesto diferencias finitas es posible utilizar la serie de Taylor,
en este trabajo, la metodología planteada y seguida es la donde:
siguiente:
f 00 (x)(h)2
1) Hacer una revisión del estado del arte, es decir, básados f (x + h) = f (x) + f 0 (x)(h) + + ... (1)
2!
en una cantidad confiable de trabajos anteriores enca- f 00 (x)(h)2
minar la investigación hacia áreas del conocimiento del f (x − h) = f (x) − f 0 (x)(h) + − ... (2)
2!
 3

Y usando unicamente los primeros tres términos de cada

~2 ∂ 2 Ψ ∂ 2 Ψ
 
expresión y restando (1) de (2) resulta: ∂Ψ
i~ =− + +VΨ
∂t 2m ∂x2 ∂y 2
f (x + h) − f (x − h)
f 0 (x) = Donde ~ (Constante de Planck), m es la masa, V es
2h
una función de energía potencial.
Con un error del orden: O(h2 )
O para la segunda derivada podría resultar:
3) Elegimos el esquema numérico adecuado que se
00 f (x + h) − 2f (x) + f (x − h) adapte al modelo, a las capacidades técnicas y a
f (x) =
2h2 los alcances del proyecto, además de estudiar el
Con un error del orden: O(h3 ) método de diferencias finitas y como ha sido utilizado
Análogamente se opera para ecuaciones de mayor orden para solucionar las ecuaciones diferenciales parciales
en una dimensión y para ecuaciones con derivadas par- anteriormente propuestas.
ciales en más de una dimensión. Al aproximar al mayor
orden posible una ecuación diferencial con una ecuación 4) Análisis de las características y condiciones del esquema
en diferencias, se pueden obtener mejores resultados, numérico: consistencia, estabilidad, convergencia (entre
pero no se garantiza que el resultado sea “aceptable”, en otras, si es necesario) que son básicas para establecer una
un sentido amplio. Para precisar conceptos, es necesario solución realmente aproximada y garantizar resultados
formular algunas definiciones que son la clave de el de calidad, estás condiciones además son caracteristicas
análisis matemático que se realiza para cada uno de los básicas de las diferencias finitas y de las estructuras que
problemas que proponemos. se basan en la discretización de los modelos.
Consistencia: Un esquema en diferencias se dice con- El tipo de problemas que se proponen para ser resuel-
sistente si la ecuación discretizada tiende a la ecuación tos con métodos numéricos basados en discretizaciones
diferencial cuando h tiende a cero. Esto es relativamente deben cumplir ciertas carácteristicas que garantizan el
fácil de verificar, dado que plantear el desarrollo en serie correcto análisis y el buen comportamiento de las solu-
de Taylor es siempre posible. ciones.
Estabilidad: El esquema se dice estable si la diferen- La notación empleada usualmente es la siguiente:
cia entre la solución exacta y la numérica permanece • El problema se considera en una region Ω cerrada

acotada, ésta condición garantiza que los errores no se en un espacio abierto de una, dos o más dimensiones
amplifican con el tiempo, es decir, controlar los errores con coordenadas cartesianas, polares, cilíndricas,
de redondeo. etc...
Convergencia: Consistencia y Estabilidad garantizan • La región Ω tiene frontera δΩ

convergencia en un problema en el cual la solución en • La solución a encontrar es una función u de varia-

todo punto del dominio depende en forma continua de bles de espacio y el tiempo definida en Ω × [0, T ]
las condiciones iniciales, lo que implica que pequeñas • Las condiciones de frontera son valores de la fun-

perturbaciones en estas, producen pequeñas discrepan- ción y las derivadas en toda o parte de la frontera
cias en la solución de Ω y la cantidad y orden de dichas condiciones
hace que el problema tenga o no tenga solución.
2) Las ecuaciones que se consideran apropiadas dados los
Típicamente las condiciones de frontera son de tres
objetivos y espectativas del trabajo y conociendo bien la
tipos:
teoría de las diferencias finítas son:
• Condición de frontera de Dirichlet : Se especifica
• Ecuación KdV (Korteweg-de Vries) valores de la función en segmentos o puntos de la
A grandes rasgos es la ecuación que describe la frontera.
dinámica de las ondas solitarias (Solitones) que se • Condición de frontera de Neumann : Se especifican
propagan sin deformarse en un medio no lineal. los valores de la derivada de una solución tomada
sobre la frontera, indican entonces un flujo o cambio
en dichos segmentos.
du(x, t) du(x, t) d3 u(x, t)
+ u(x, t) +µ =0 • Condición de frontera mixta : Combinaciones pon-
dt dx dx3 deradas o no ponderas, por segmentos o puntos de
la frontera, de las condiciones de Dirichlet (sobre
Donde  y µ son reales ponderadores de la no la función) y condiciones de Neumann (sobre la
linealidad y la dispersión respectivamente. derivada
• Ecuación de Schrödinger
Las condiciones iniciales indican un valor de la función
Ecuación que soluciona la función de onda Ψ, que
u en el tiempo t = 0 sobre todo el dominio Ω.
es la función que contiene toda la información que
puede conocerse sobre la partícula y de la que Consistencia, convergencia y estabilidad: Cuando se
resulta la aplitud de probabilidad de presencia de propone un problema para solucionar por medio de un
dicha partícula. esquema numérico con buenas y confiables cualidades
 4

se debe primero garantizar que el problema está bien amplifican con el tiempo, es decir, controlar los errores
definido o bien propuesto (en el sentido de Hadamard), de redondeo.
lo anterior es que una ecuación diferencial con valores Una interpretación de la estabilidad de un esquema
iniciales tenga propiedades analíticas adecuadas y sus de diferencias finitas es que: Para un esquema estable
soluciones tengan una estructura conveniente. variaciones pequeñas en los valores de las condiciones
En general, los problemas se dicen bien propuestos en iniciales originan errores pequeños en la solución.
el sentido de Hadarmad si: En este caso lo relevante es asegurar la estabilidad
a) Unicidad: las soluciones estrictas están determina- de la solución numérica respecto de la analítica, que
das unívocamente por las condiciones iniciales. resulta ser, entonces, el criterio más general. Hablar de
b) Conjunto denso: el conjunto de todas las condi- la estabilidad de la solución numérica respecto de la
ciones iniciales correspondientes a las soluciones analítica,
 significa acotar
 de alguna manera la diferencia
posibles es denso en el espacio de Banach en el U (xi , yj , tk ) − U k  < ∞cuando k → ∞.
i,j
que se plantea problema. Para garantizar la estabilidad, sobretodo, cuando la solu-
c) Acotación local: Para algún intervalo finito [0, t0 ] ción analítica no se conoce a ciencia cierta, se implemen-
existe una constante K tal que cada solución es- tan los analisis de estabilidad basados en funciones de
tricta satisface la desigualdad: kut k ≤ K ku0 k. lyapunov, el método de Von Neumann basádo en series
Tener un problema bien definido entonces significa que de Fourier (que se usará en éste trabajo) u otros métodos
la solución depende continuamente de las condiciones de acotamiento que aprovechan los espacios normados
impuestas sobre la frontera y es uniformemente cerrada a los que pertenecen en general las funciones que son
en un intervalo compacto lo que hace entonces que solución de los problemas bien definidos en el sentido
todos los análisis posteriores sean más cómodos y tengan de Hadamard.
además un sentido más matemático garantizando así
el correcto comportamiento de todo el esquema, para
Convergencia: Se dice que una solución numérica con-
encontrar soluciones más aproximadas a las que se
verge a la solución analítica cuando, para cada punto (en
esperan de un desarrollo analítico pero de una manera
el espacio de 1as variables independientes Ω), la primera
más flexible y computacional.
tiende a la segunda al refinar la malla en la que se están
Luego bien, las carácteristicas que se deben garantizar
haciendo los cálculos del esquema numérico propuesto.
para hacer confiables las soluciones obtenidas de un
Es importante resaltar que la consistencia del esquema
desarrollo numérico basado en diferencias, tratan del
numérico es necesaria para que el error de truncamiento
análisis asintótico del error e indican que un esquema
efectivamente disminuya al refinar la malla pero no es
numérico se comporta bien en todos los instantes y
suficiente y para complementar el análisis es necesario
no hay posibilidades de soluciones inesperadas o no
demostrar que el problema en diferencias permanece
aproximadas correctamente.
estable.
Consistencia: Un esquema en diferencias se dice con- Estos conceptos se expresan en el Teorema de Equiva-
sistente si la ecuación discretizada tiende a la ecuación lencia de Lax-Richtmyer que se enunciará más adelante
diferencial cuando ∆ tiende a cero. En particular ga- La impotancia entonces de éste resultado es que dado
rantizar este correcto comportamiento de las ecuaciones un problema de valores iniciales bien planteado y una
es relativamente fácil de verificar, dado que plantear el aproximación en diferencias finitas que satisface la con-
desarrollo en serie de Taylor es siempre posible. dición de consistencia, la estabilidad es una condición
Dada una ecuación diferencial parcial y un esquema de necesaria y suficiente para la convergencia que es por
diferencias finitas es posible asegurar que éste es un es- si sola una condición complicada dado el análisis de la
quema consistente con la ecuación diferencial parcial si malla que requiere.
para cualquier función suave Φ(x, t) la diferencia entre
la ecuación diferencial y la evaluación en la ecuación
5) Implementación computacional del método de
en diferencias converge puntualmente a 0 cuando ∆x,
diferencias finítas, en Matlab, software que la
∆t convergen a 0. De donde resulta entonces que la
universidad EAFIT dispone para los desarrollos
consistencia significa que si u es una solución (suave)
investígativos de ésta indole, generando algoritmos de
de una ecuación diferencial parcial, entonces u es una
baja complejidad y fácil divulgación.
aproximación de la solución del esquema de diferencias
finitas que aproxima la EDP.
6) Análisis los resultados encontrados apartir de las si-
En este punto es importante anotar que la consistencia
mulaciones de los modelos propuestos, comparar los
es una condición necesaria para la convergencia de un
resultados obtenidos con los esperados ya sea con cono-
esquema de diferencias finitas, sin embargo no todo
cimiento previo del comportamiento de los sistemas, o
esquema consistente será convergente.
con modelos de prueba como datos de validación, para
Estabilidad: El esquema se dice estable si la diferen- de ésta misma manera garantizar experimentalmente que
cia entre la solución exacta y la numérica permanece las aproximaciones son correctas.
acotada, ésta condición garantiza que los errores no se
 5

III. R ESULTADOS
du(x, t)
j
U j − Ui−1
Para presentar los resultados de éste trabajo es muy impor- = i+1
tante tener claro que la implementación y el análisis se realizó dx 2∆x
de manera disjunta para los dos problemas ya mencionados, j j j j
entonces se presentan los resultados separados. d3 u(x, t) Ui+2 − 2Ui+1 + 2Ui−1 − Ui−2
=
dx3 2∆x3
A. Resultados KdV De lo anterior se obtiene entonces un esquema numérico
para la ecuación KdV en el dominio numérico propuesto,
La ecuación de Korteweg-de Vries o KdV es una ecuación
descrito por la ecuación
en derivadas parciales que incluye efectos de no linealidad y
dispersión a la vez. Físicamente es un modelo que describe, en
una dimensión espacial, la propagación de ondas de longitud Uij+1 = Uij−1 − AUij (Ui+1
j j
− Ui−1 ) (4)
de onda larga en medios dispersivos. − j
B(Ui+2 − j
2Ui+1 + j
2Ui−1 − j
Ui−2 ) (5)
3
du(x, t) du(x, t) d u(x, t) donde
+ u(x, t) +µ =0 (3)
dt dx dx3
∆t µ∆t
Donde  y µ son reales ponderadores de la no linealidad y A= y B=
3∆x ∆x3
la dispersión respectivamente.
Con condiciones inciales y de frontera : Gráficamente es posible ver que cada solución temporal
y espacial de la ecuación KdV que resulta del esquema
u(x, 0) = f (x) u(xmin , t) = g(t) propuesto en la ecuación(5) depende de 6 putnos en el espacio
u(xmax , t) = h(t) (Ver. fig 1).
u(x, t) describe la elongación de la onda en el lugar y en el
tiempo. KDV es no lineal debido al producto que se muestra en
el segundo término y de tercer orden es a causa de la tercera
derivada en el tercero. El término no lineal, es similar a la
ecuación de onda término usual e implica que cuando u(x, t)
no varía demasiado, la onda se propaga con una velocidad
proporcional; éste término además introduce la posibilidad
de ondas de choque en la solución. Y por otra parte el
término de tercer orden produce dispersión de ampliación que
exactamente puede compensar el estrechamiento provocado Figure 1. Esquema numérico - KdV
por el término no lineal en condiciones adecuadas.
Siendo (5) un esquema numérico relativamente sencillo en
Para solucionar la ecuación KdV se propone entonces
el dominio discretizado, es computacionalemente fácil de im-
discretizar todo el dominio y encontrar por medio de un
plementar, pero dado que la ecuación KdV tiene restricciones
esquema numérico propio de las diferencias finítas una buena
por las condiciones de frontera y la geometría del problema,
aproximación a la solución del comportamiento de las ondas.
es recomendable analizar esquemas alternos para las fronteras,
Es posible escribir el dominio numérico como una dis-
lo anterior se clarifica en la figura 2 donde se muestra que
cretización del espacio entre xmin y algún xmax elegidos y
es posible crear esquemas con menos dependencias dada la
además discretizar el tiempo hasta T que será el instante donde
rigidez del problema en las fronteras, y éstos se establecen en
nos interese la solución.
las ecuaciones (6) y (7).
(x, t) ∈ [xmin , xmax ] × [0, T ] Es importante resaltar que dichas modificaciones del es-
quema en las fronteras se logran implementando de otra
manera las diferencias finítas, es decir, no haciendo siempre
xi = xmin + i∆x i = 0, 1, ..., S implementaciones centradas.
tj = j∆t j = 0, 1, ..., P
La solución entonces se denotará u(x, t) = U(xi , tj ) = Uji
.
Para lograr una aproximación en diferencias de la ecuación
KdV se propone aproximar cada una de las derivadas par-
ciales con aproximaciones en diferencias como se definieron
antes entonces se tienen las aproximaciones a las derivadas
centradas:

du(x, t) U j+1 − Uij−1

= i Figure 2. Esquema numérico (Fronteras) - KdV
dt 2∆t
 6

Proof: Para probar que el esquema numérico prop-

∆t f f uesto para solucionar la ecuación KdV es convergente, se
Uif = Uif −1 − U (U f
− Ui−1 ) utiliza el siguiente teorema.
6∆x i i+1
µ∆t Teorema 1. (Teorema de Equivalencia de Lax-
− (U f − 2Ui+1
f f
+ 2Ui−1 f
− Ui−2 ) (6)
2∆x3 i+2 Richtmyer) Un esquema de diferencias finitas Consistente
para una ecuación diferencial parcial tal que el problema
∆t j j de valor inicial es bien planteado es Convergente si y sólo
Ufj+1 = Ufj−1 − U (U − Ufj−1 ) si es Estable.
3∆x f f +1
2µ∆t j Es decir, dado que ya se han garantizado el planteamiento,
− (U − 2Ufj+1 + 2Ufj−1 − Ufj−2 ) (7)
∆x3 f +2 la consistencia y la estabilidad del esquema, empleando el
donde f = 1 ó s − 1 teorema de Lax, garantiza la convergencia de la solución
numérica
Dada la metodología propuesta para éste trabajo, se estipula
que luego de encontrar esquemas visiblemente buenos para
lograr aproximaciones a la solución de la ecuación, se propone, 1) Caso particular : Luego de haber garantizado, como
como parte vital del trabajo, el desarrollo teórico de las buena práctica teórica todas las condiciones del esquema, para
propiedades de consistencia, estabilidad y consistencia (en éste realizar la implementación se tomó el caso particular de la
caso) que fueron descritas anteriormente y que además son ecuación KdV
propias de los esquemas en diferencias.
du(x, t) du(x, t) d3 u(x, t)
• Consistencia − 6u(x, t) + =0 (10)
Un esquema en diferencias se dice consistente si la dt dx dx3
ecuación discretizada tiende a la ecuación diferencial con condiciones de frontera
cuando ∆ tiende a cero.
Proof: Trivial (Por la misma construcción del es-
quema apartir de la ecuación diferencial, ésta propiedad u(0, t) = u(S, t) = 0
se puede obviar) √ 2
0,2 0,2
u(x, 0) = sech x
• Estabilidad 2 2
Un esquema en diferencias se dice estable si la diferencia
Se obtuvo como resultado la superficie de la figura 3. que
entre la solución exacta y la numérica permanece acotada,
muestra el comportamiento espacial-temporal de los solitones,
ésta condición garantiza que los errores no se amplifican
con una discretización del espacio ∆x = 0,5 y del tiempo
con el tiempo.
∆t = 0,0625.
Proof: Criterio de estabilidad de Von Neumann
Se considera una expansión en serie de Fourier del error
en un punto concreto de la malla.
U n (i) = An exp I[k1 m(i∆x)]


(8)
donde: A es el factor de amplificación o variación tem-
poral, k1 es un modulo de cambio en la coordenada x.
Con lo anterior se puede construir una expresión recursiva
para la solución en (n+1) como sigue

U n+1 (i) = An U n (i) (9)

Al garantizar que la amplificación o variación temporal
permanezca acotada y luego de realizar la respectiva
expansión de (9) en términos del esquema numérico (5), Figure 3. Superficie - KdV
se tiene que:
Siempre que |A| ≤ 1 se obtiene una relación de estabili- Es importante además que los resultado obtenidos con el
dad para las disctretizaciones del tiempo y del espacio esquema numérico concuerden con los resultados obtenidos
3∆x con otros métodos numéricos y aproximaciones a la solución
∆t ≥ exacta que se encontraron durante la revisión de la literatura,
2
como la solución encontrada por Kolebaje, Olusola et.al. 2012.

• Convergencia B. Resultados Schrödinger

Se dice que una solución numérica converge a la solución La ecuación de Schrödinger fue desarrollada por el físico
analítica cuando la primera tiende a la segunda al refinar austríaco Erwin Schrödinger en 1925. Describe la evolu-
la malla. ción temporal de una partícula masiva no relativista. Es de
 7

importancia central en la teoría de la mecánica cuántica, de Runge-Kutta que al igual que las diferencias finitas resultan
donde representa para las partículas microscópicas un papel de expansiones en series de Taylor y con éstos se obtienen
análogo a la segunda ley de Newton en la mecánica clásica. muy buenos resultados en términos de costo computacional
Generalmente se describe como sigue y convergencia de la solución, en nuestro caso, la derivada
temporal se considera como una discretización simple (hacia
~2 ∂ 2 Ψ ∂ 2 Ψ
 
∂Ψ atrás) del tiempo dado que es de nuestro interés particular el
i~ =− + +VΨ
∂t 2m ∂x2 ∂y 2 trabajo de las diferencias finitas sin la intervención de otros
Donde ~ (Constante de Planck), m es la masa, V es una métodos, aunque cabe anotar que sí garantizamos que, dada
función de energía potencial. la revisión de la literatura, no nos alejamos mucho de las
precisiones que se obtienen con métodos más sofisticados (R.
En nuestro caso, implementamos el método de diferencias Becerril, et.al 2008).
finítas en dos casos, para cajas de pontencial cero en dos y
tres dimensiones. Gráficamente es posible ver que cada solución temporal y
La partícula en una caja (también conocida como pozo de espacial de la ecuación KdV resultanto del esquema propuesto
potencial infinito) es un problema muy simple que consiste de en la ecuación (11) depende de 3 puntos en el espacio cuando
una sola partícula que rebota dentro de una caja inmóvil de la se analiza en 1 dimensión (Ver. figura 4). Resaltando que en
cual no puede escapar, y donde no pierde energía al colisionar éste caso, se construyeron las ecuaciones en diferencias de
contra sus paredes. manera que no haya problemas en las fronteras empleando
Por lo anterior se tienen las siguientes condiciones de fron- diferencias hacia atrás y no centradas como en la imple-
tera iguales a cero expresando el hecho de que la probabilidad mentación anterior.
de encontrar la partícula fuera de una caja de la que la partícula Analogaménte funciona el esquema numérico de manera
no puede escapar es cero. gráfica para la ecuación en 2 dimensiones.

Ψ(0, y, t) = Ψ(Lx , y, t) = 0
Ψ(x, 0, t) = Ψ(x, Ly , t) = 0
Para discretizar el problema e implementar el método de
diferencias finítas explícitas, se propone (como antes) un
dominio numérico con el fin de discretizar el espacio y el
tiempo

(x, y, t) ∈ [0, Lx ] × [0, Ly ] × [0, T ] Figure 4. Esquema numérico 1D - Schrödinger

xj = j∆x j = 0, 1, ..., S Luego de tener el esquema numérico se realiza entonces

yk = k∆y k = 0, 1, ..., P (de manera análoga a la ecuación KdV), las pruebas de
tn = n∆t n = 0, 1, ..., N consistencia, de estabilidad y de convergencia del esquema
y de la solución propuesta que son vitales para éste trabajo.
En particular, se construye una prueba generalizada en 2
denotando entonces la solución como sigue dimensiones para la estabilidad
• Estabilidad
Ψ(x, y, t) = Ψ(xj , yk , tn ) = Ψnj,k
Proof: Criterio de estabilidad de Von Neumann
Se Considera una expansión en serie de Fourier del error en
Con el mismo desarrollo anterior de cada una de las aprox- un punto concreto de la malla como se había indicado antes
imaciones de las derivadas parciales, se obtiene el siguiente
Ψn (i) = An exp I[k1 m(i∆x) + k2 m(j∆y)]


(13)
esquema numérico
• En 1 Dimensión donde: A es el factor de amplificación o variación temporal,
∂Ψ i Ψnj+1,k − 2Ψnj,k + Ψnj−1,k
  ki es un modulo de cambio en la coordenada x o y respectiva-
= + O(h)(11) mente.
∂t 2 ∆x2
Con lo anterior se puede construir una expresión recursiva
• En 2 Dimensiones para la solución en (n+1) como sigue
i Ψnj+1,k − 2Ψnj,k + Ψnj−1,k

∂Ψ
=
∂t 2 ∆x2 Ψn+1 (i) = An Ψn (i) (14)
Ψnj,k+1 − 2Ψnj,k + Ψnj,k−1

+ + O(h)(12) Al garantizar que la amplificación o variación temporal
∆y 2 permanezca acotada y luego de realizar la respectiva expansión
Para solucionar el esquema numérico (12) en la literatura se de (14) en términos del esquema numérico (12), se tiene
propone resolver la derivada temporal por medio de métodos que naturalmente el esquema mismo garantiza la estabilidad
 8

en términos de ∆x y ∆y para valores significativamente

pequeños.
Resaltamos que como condición inicial de cada uno de los Los resultados obtenidos para la ecuación de Schrödinger en
casos particulares que se mostrarán, se tomó la solución exacta 1D en una caja de potencial coinciden además con la solución
dado que, no existe otra manera de simular el comportamiento exacta y con los resultados aproximado que se revisaron en la
de una particula en el tiempo t = 0, y cabe anotar que las literatura.
soluciones analíticas inciales fueron construidas por el autor.
2) Caso particular : Solución numérica para la ecuación
1) Caso particular : Solución numérica para la ecuación de Schrodinger 2D en una caja de potencial
de Schrodinger 1D en una caja de potencial.
Condicion inicial (t = 0) Condicion inicial (t = 0)
√
Ψnx (x) = 2 cos(Enx t) sin(nx πx) (15)
Ψnx,ny (x, y) = 2 cos(En t) sin(nx πx) sin(ny πx) (16)
Solución Numérica (Particula 1D−Box)
1.5

0.5
ΨR

−0.5

−1

−1.5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x

Figure 5. Solución numérica nx = 2, ∆x = 0.01, ∆t = 0.0001

En la figura 5. se muestra entonces en diferentes colores

cada una de las amplitudes de probabilidad o funciones de
onda que son solución en cada instante del tiempo de la parte
real de la ecuación de Schrödinger, resaltando además que
se ha tomado nx = 2 es decir, 2 modos de vibración para la
condición inicial del problema y de donde surgen ambos picos Figure 7. Solución numérica nx = ny = 3, ∆x = ∆y = 0.01, ∆t =
de la gráfica. 0.0001
La interpreatación física y matemática del resultado anterior
se dificulta bastante dado que su significado toma mucha más
fuerza al leer la densidad de probabilidad asociada al espacio,
que además es independiente del tiempo y se encuentra como
En la figura 7. se muestra en 3 tiempos (uno inicial, uno
el complejo conjugado de la función de onda.
medio, y uno final) la solución de cada una de las amplitudes
ρ = |Ψ? Ψ| de probabilidad o funciones de onda que son solución de la
ecuación de Schrödinger en 2 dimensiones, resaltando además
Dicha función de densidad de probabilidad se muestra en la
que se ha tomado nx = ny = 3 es decir, 3 modos de vibración
figura 6. y efectivamente resulta ser independiente del tiempo
en cada una de las dimenciones espaciales para la condición
como se esperaba.
inicial del problema y de donde surgen los picos de la gráfica.
Función de densidad de probabilidad
La interpreatción física del resultado anterior se dificulta
2 bastante dado que su significado toma mucha más fuerza al
1.8
leer la densidad de probabilidad asociada al espacio y que
1.6

1.4
además es independiente del tiempo y se encuentra como el
1.2 complejo conjugado de la función de onda.
ρ

0.8

0.6
ρ = |Ψ? Ψ|
0.4

0.2

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

x Dicha función de densidad de probabilidad se muestra en la

figura 8. y efectivamente resulta ser independiente del tiempo
Figure 6. ρ Independiente del tiempo como se esperaba.
 9

analíticas es en muchos casos no triviales tediosa y

extensa.
• La importancia de garantizar siempre que los esquemas
numéricos cumplan con las condiciones básicas de con-
strucción y comportamiento radica en que la comodidad
al momento de la implementación es mucho mayor, se
pueden asegurar condiciones sobre las mallas y sobre
los parámetros y además se puede defender la solución
teóricamente.
• La ecuación KdV ha sido resuelta por medio de muchos
métodos numéricos y dada la estricta revisión bibliográ-
Figure 8. ρ Independiente del tiempo fica que se hizo para realizar éste trabajo, es posible
asegurar que las diferencias finítas explícitas en términos
de simplicidad tanto matemática como computacional,
de esfuerzo y de resultados, son uno de los métodos
Los resultados obtenidos para la ecuación de Schrödinger en más eficientes, que (por lo menos en el caso particular
2D en una caja de potencial tambien coinciden con la solución aquí expuesto) logra interpretar la difusión y todas las
exacta y con los resultados aproximado que se revisaron en la no linealidades en pasos muy pequeños del espacio y del
literatura. tiempo.
• Por otro lado, la ecuación de Schrödinger que como se
mencionaba, es una de las ecuaciones más importantes
IV. C ONCLUSIONES
de la física moderna, también suele resolverse por medio
• Los métodos numéricos realmente son una manera apli- de muchos métodos numéricos, y además tiene solución
cable de encontrar aproximaciones a la solución y se exacta en los casos particulares que se trabajaron, y
convierten en una herramienta muy útil cuando se trata conociendo éstas condiciones, el método de diferencias
de comprender algunos comportamientos o fenómenos. finitas resultó ser eficiente y eficaz para solucionarla, ya
Como resultado de éste trabajo se tienen un conjunto de que los resultados obtenidos son altamente semejantes
desarrollos en términos de la formación matemática y a lo analítico y además el esquema propuesto cumple
de la aplicabilidad de la ingeniería a favor de resolver con las condiciones para garantizar el comportamiento
problemas de manera teórica y computacional. correcto de las soluciones.
• El método de diferencias finitas a pesar de ser un método • Se puede entonces asegurar que con el método de difer-
simple, basado en la teoría de las series de Taylor y encias finítas explícitas se pueden construir soluciones
la definición de derivada, resulta ser muy bueno para muy aproximadas a las analíticas, y se pueden simular
representar las variaciones temporales y espaciales de los comportamiento lineales y no lineales siempre y cuando
fenómenos físicos. se garantice previamente la consitencia, convergencia y
• Es muy importante, además de tener aproximaciones estabilidad de la solución y del esquema.
buenas a la solución analítica por medio de métodos
numéricos, analizar dichos esquemas, es decir, aprovechar R EFERENCES
la flexibilidad matemática de los métodos numéricos y
anticiparse a mejoras en las implementaciones computa- [1] J. Strikwerda, Finite difference schemes and partial differential equa-
tions. Society for Industrial Mathematics, 2004.
cionales, lo anterior es una buena práctica teórica y [2] A. Capella Kort and J. Sarrate Ramos, “Estudio numérico del control
además minimiza las probabilidades de errores al mo- exacto desde el contorno de la ecuación korteweg de vries lineal,”
mento de simular los fenómenos. Revista internacional de métodos numéricos, 2003.
[3] W. Hui-Ping, W. Yu-Shun, and H. Ying-Ying, “An explicit scheme for
• No todas las carácteristicas de los métodos numéricos the KdV equation,” Chinese Physics Letters, vol. 25, no. 7, 2008.
(positividad, monotonicidad...) se deben garantizar en [4] R. Becerril, F. S. Guzmán, A. Rendón-Romero, and S. Valdez-Alvarado,
cada problema a resolver, dado que la naturaleza del “Solving the time-dependent schrodinger equation using finite difference
methods,” Revista mexicana de física E, vol. 54, no. 2, 2008.
modelo indica que la solución cumple o no con determi- [5] M. G. Vázquez and Capetillo., “Simulaciones numéricas de la ecuación
nadas condiciones, por otra parte, es de vital importancia del calor: diferencias finitas.”
entonces el análisis del planteamiento del problema, el [6] M. Duruflé and S. Israwi, “A numerical study of variable depth KdV
conocimiento de la EDP que lo describe, las condiciones equations and generalizations of CamassaâĂŞHolm-like equations,”
Journal of Computational and Applied Mathematics, 2012.
y la geometría sobre las que se trabaja. [7] A. C. Vliegenthart, “On finite-difference methods for the korteweg de
• Ecuaciones que son muy populares en la física y vries equation,” Journal of Engineering Mathematics, vol. 5, no. 2, 1971.
que tienen una cantidad innumerable de aplicaciones [8] B. F. Feng and T. Mitsui, “A finite difference method for the korteweg-
de vries and the kadomtsev-petviashvili equations,” Journal of compu-
como las aquí resueltas (Ecuación KdV y Ecuación de tational and applied mathematics, vol. 90, no. 1, 1998.
Schrödinger) se interpretan en muchos casaos de una [9] P. Mercader, J. Torres, M. J. Núñez, J. M. Zamarro, E. Martín, and G. J.
forma más cómoda cuando se obtiene para ellas una Molina-Cuberos, “Resolución numérica de la ecuación de Schrodinger.
Introducción al estudio de la teoría de bandas.”
aproximación y sobretodo una forma de visualizarlas ya [10] J. R. Nagel, “The one-dimensional finite-difference time-domain
que la matemática que se suele utilizar en las soluciones (FDTD) algorithm applied to the schrodinger equation.”
 10

[11] Y. Erlangga and E. Turkel, “Iterative schemes for high order compact
discretizations to the exterior helmholtz equation,” ESAIM: Mathemati-
cal Modelling and Numerical Analysis, vol. 46, no. 03, 2012.
[12] C. E. Froberg and C. E. Frhoberg, Introduction to numerical analysis.
Addison-Wesley Publishing Company, 1969.
[13] A. Iserles, A first course in the numerical analysis of differential
equations. Cambridge University Press, 2008, vol. 44.
[14] J. C. Strikwerda, Finite Difference Schemes and Partial Differential
Equations. SIAM, Philadelphia, 1994.
[15] K. W. Morton and D. F. Mayers, Numerical Solution of Partial Differ-
ential Equations: An Introduction. Cambridge University Press, Apr.
2005.
[16] S. Jamrud Aminuddin, “Numerical solution of the korteweg de vries
equation.”
[17] O. T. Kolebaje and O. E. Oyewande, “Numerical solution of the ko-
rteweg de vries equation by finite difference and adomian decomposition
method,” International Journal of Basic and Applied Sciences, vol. 1,
no. 3, 2012.
[18] M. L. Gandarias and M. S. Bruzón, “Some conservation laws for a forced
KdV equation,” Nonlinear Analysis: Real World Applications, 2012.
[19] E. Atanassov and S. Ivanovska, “Sensitivity study of heston stochas-
tic volatility model using GPGPU,” Large-Scale Scientific Computing,
2012.