1. No category

Práctica 8: Ecuaciones diferenciales. Problemas de
contorno
1
Introducción. El método de disparo
En esta práctica consideramos ecuaciones diferenciales de la forma
y 00 =
d2 y
= f (x, y) ,
dx2
(1)
con la condición de contorno (o condición en la frontera)
y(a) = y1 ; y(b) = y2
(2)
Para resolver este problema vamos a emplear el método de disparo. Esquemáticamente el
procedimiento comienza escribiendo la ecuación diferencial (1) como el sistema de ecuaciones diferenciales de primer orden:
dy
= z
dx
dz
= f (x, y)
dx
(3)
Este sistema debe resolverse con la condición inicial y(a) = y1 , pero se desconoce la
condición inicial para z. En el método de disparo se lleva a cabo un proceso iterativo
en el que se va variando la condición inicial z(a) hasta que se cumple la condición de
contorno y(b) = y2 .
Consideremos el siguiente ejemplo: una barra metálica de longitud L está colocada
entre dos cuerpos a temperaturas T1 y T2 . Esta barra no está aislada y puede intercambiar
calor con el entorno, que se encuentra a una temperatura Tα . En estado estacionario la
temperatura de la barra verifica la ecuación:
d2 T
+ h(Tα − T ) = 0
dx2
(4)
donde h es un coeficiente de transferencia de calor. Las condiciones de contorno son:
T (0) = T1 ; T (L) = T2
(5)
Escribimos la ecuación (4) en la forma:
dT
= z
dx
dz
= h(T − Tα )
dx
1
(6)
Resolveremos este sistema con la condición inicial T (0) = T1 ; z(0) = α e iremos variando
α hasta que se cumpla T (L) = T2 .
Consideremos el caso concreto: L=10m, h=0.01m−2 , T (0) = 400 C, T (10) = 2000 C,
Tα = 200 C. Para ilustrar el procedimiento, se dispone del programa barra en el que
hemos programado las ecuaciones (6). Como primer ”disparo”, tomemos un valor inicial
α = 10, y resolvamos el sistema utilizando el comando ODE45 empleado en la práctica
anterior.
y0=[40 10];
ode45 (’barra’, [0,10], y0)
En la figura resultante observamos que hemos obtenido un valor de la temperatura en
el extremo de la barra inferior al valor de contorno (200). El valor numérico se determina
escribiendo:
[x,T]= ode45 (’barra’, [0,10], y0)
Resulta T (10)=168.38170 C.
Efectuamos un segundo disparo con α=20. Repitiendo el proceso anterior obtenemos
ahora T (10)=285.90190 C. Se demuestra que para ecuaciones lineales como (6), el resultado T (L) depende linealmente del parámetro de disparo α. Para comprobar este punto
numéricamente, seguir los siguientes pasos:
• Interpolar linealmente para obtener el valor α (12.6907),que conduce a T (10)=200 0 C.
• Repetir la integración y comprobar el resultado.
• Obtener la temperatura de la barra en el punto x=6m interpolando linealmente en
la tabla de valores obtenida para T .
2
Aplicación: Reacción quı́mica con difusión
Para ilustrar la utilización del método en un caso más general, consideramos como ejemplo
una reacción quı́mica irreversible de primer orden que tiene lugar en un sistema con una
sola fase y que contiene únicamente el reactivo A y el producto B de la reacción. La
difusión del reactivo A y la reacción quı́mica se producen simultaneamente, lo que conduce
a la ecuación diferencial de segundo orden:
k
d 2 CA
=
CA
2
dz
DAB
(7)
donde CA es la concentración del reactivo A, z la distancia a lo largo de la mezcla de
reacción, k es la constante de velocidad de la reacción, y DAB el coeficiente de difusión
2
mutuo. Las condiciones de contorno para este problema son:
CA (z = 0) = CA0 ;
dCA
= 0 para z = L
dz
(8)
Como en el ejemplo anterior, la ecuación diferencial de segundo orden puede expresarse
como un sistema de dos ecuaciones de primer orden:
dCA
= u
dz
k
du
=
CA
dz
DAB
(9)
donde se desconoce la condición inicial para la segunda ecuación. Como en el ejemplo
anterior, buscaremos la condición inicial que conduzca a una solución que verifique la
condición de contorno u(L) = 0, y para ello vamos a emplear un método iterativo que
proporciona el valor óptimo del parámetro α = u(z = 0). Sea u(z; α) la solución del
sistema (9) obtenida resolviendo dicho sistema con la condición inicial u(0) = α. El valor
de α buscado es por tanto la solución de la ecuación
u(L; α) = 0
(10)
y para resolver esta ecuación algebraica empleamos el método de Newton-Raphson:
αk+1 = αk −
u(L; αk )
u0 (L; αk )
(11)
donde
du(L; α)
(12)
dα
Para obtener el valor de u0 , introducimos las derivadas de la funciones CA y u con respecto
al parámetro α:
dv
∂CA
, w(z) =
= u0 (z; α)
(13)
v(z) =
∂α
dz
Derivando con respecto a α las ecuaciones (9) y teniendo en cuenta que z es independiente
de α, resulta el sistema de cuatro ecuaciones diferenciales:
u0 (L; α) =
dCA
dz
du
dz
dv
dz
dw
dz
= u
k
CA
DAB
=
= w
k
v
DAB
=
(14)
3
que resolveremos con las condiciones iniciales:
CA (z = 0) = CA0 ; u(z = 0) = α ; v(z = 0) = 0 ; w(z = 0) = 1 ,
(15)
donde la condición inicial de la tercera ecuación es consecuencia de imponer que la
condición fı́sica del problema (CA (z = 0) = CA0 ) se verifique para cualquier α. La
cuarta condicion inicial proviene de que u(z = 0) = α, de manera que
∂u =1
(16)
w(z = 0) =
∂α z=0
resulta de derivar u con respecto a α en z = 0. Una vez resuelto este sistema de ecuaciones
diferenciales, disponemos de las cantidades necesarias para aplicar el método de NewtonRaphson: u(L; α) y u0 (L; α) = w(L)
El esquema del procedimiento es por tanto:
1. Resolución del sistema (14) con un valor inicial α1 del parámetro de disparo. La
solución nos proporciona la fución w(z)
2. Obtención de un segundo valor α2 empleando el método de Newton-Raphson (ecuación
(11)) con u0 (L; α1 ) = w(L).
3. Proceso iterativo en el que se emplean los valores sucesivos obtenidos en (11) como
condiciones iniciales del sistema (14).
Este procedimiento se ha programado en el programa difusi.m que llama a la subrutina
disparo.m que evalúa a su vez los segundos miembros de las ecuaciones (14). Para ejecutarlo basta escribir:
difusi
El programa nos muestra gráficamente la solución CA (z) en cada paso del procedimiento
iterativo. Es necesario pulsar ”return” para que el procedimiento continúe.
Como ejercicio, modificar el programa disparo.m para que resuelva el problema
análogo:
d 2 CA
k
=
C2
(17)
2
dz
DAB A
con las mismas constantes que el ejemplo anterior. Modificar asimismo el programa
difusi para que llame al nuevo programa. Salvar los programas modificados con nombres
distintos.
4
3
La ecuación de Schrödinger para el oscilador armónico.
3.1
Planteamiento del problema.
La resolución de la ecuación de Schrödinger para un sistema monodimensional es un
problema similar a los que hemos discutido anteriormente. En concreto, la vibración de
dos núcleos o fragmentos moleculares A-B puede describirse mediante la ecuación (en
unidades atómicas):
1 d2 ψ 1 2
−
+ kx ψ(x) = E ψ(x)
(18)
2µ dx2
2
donde:
• x es la separación internuclear, x = 0 es la posición de equilibrio, x < 0 indica que
los núcleos están más próximos que en la posición de equilibrio x > 0 indica que
están más alejados.
• µ es la masa reducida:
1
1
1
=
+
µ
MA MB
(19)
• k es la constante de fuerza del oscilador armónico que empleamos como modelo para
este problema. Este oscilador tiene un frecuencia
s
1
k
ω
ν=
=
(20)
2π µ
2π
• E es la energı́a de vibración.
• ψ(x), solución de (18), es la función de onda.
Las soluciones del problema deben verificar la condición de normalización
Z +∞
ψ 2 dx = 1
(21)
−∞
Para satisfacer esta condición es preciso que ψ se anule cuando |x| se hace muy grande,
lo que significa que los núcleos permanecen enlazados.
3.2
Resolución numérica.
Para simplificar el problema, es conveniente introducir las magnitudes escaladas
2E
E
=
πν
ω
r
µν
x
X =
2π
=
5
(22)
que permiten escribir la ecuacion (18) en una forma independiente de los parámetros del
problema:
d2 ψ
+ ( − X 2 ) ψ(X) = 0
(23)
dX 2
La resolución del problema sigue los mismos pasos que en los ejemplos anteriores.
Comenzamos escribiendo la ecuación diferencial de segundo orden como un sistema de
ecuaciones diferenciales de primer orden:
dψ
= z
dX
dz
= (X 2 − ) ψ
dX
que resolvemos con la condición inicial
ψ(X0 ) = 0, z(X0 ) = α
(24)
(25)
donde el punto X0 (X0 < 0) está suficientemente alejado de la posición de equilibrio. El
valor de α se selecciona para lograr que se cumpla la condición de normalización (21).
En la práctica disponemos del programa oscila.m, que resuelve las ecuaciones (24)
para una determinada energı́a
3.3
Cuantización de la energı́a.
Comprobamos que sólo algunos valores de la energı́a conducen a soluciones que se anulan
a largas distancias y por tanto conducen a soluciones aceptables. Para ello ejecutamos
el programa oscila con energı́as 0.2hν y 0.5hν. El programa nos preguntará si queremos
iterar sobre el valor de α. De momento contestar a la pregunta con un cero, para indicar
que no queremos iterar:
oscila
energia en unidades de omega 0.2
iteracion 0/1 0
El programa resuelve las ecuaciones (24) con el valor de α inicial y nos presenta gráficamente
la solución ψ(X) Vemos que para E=0.2ω la función diverge para valores grandes de X,
mientras que para E=0.5ω la función se anula cuando x crece y por tanto puede normalizarse. Repetir el proceso para energı́as 1.0ω, 1.5ω y 2.5ω y decidir qué valores son
posibles.
3.4
Normalización.
La aplicación del método de disparo a este problema se basa en iterar sobre α hasta que
se cumple la condición (21). El programa oscila calcula la integral empleando la regla
6
trapezoidal y utiliza el método de Newton-Raphson para buscar el valor de α. Repetimos
el cálculo anterior para obtener la solución normalizada:
oscila
energia en unidades de omega 0.5
iteracion 0/1 1
Repetir el proceso para E=1.5ω y 2.5ω y contestar la pregunta 4 de la hoja de resultados.
7
4
Resultados.
NOMBRE Y APELLIDOS:
1. Escribir el valor obtenido para la temperatura en x = 6 en el ejemplo de la barra
de la sección 1.
2. Obtener el valor de CA (L) en el problema modificado con dependencia cuadrática
en CA .
3. Decidir si son posibles los valores 1.0ω, 1.5ω y 2.5ω para el oscilador armónico.
¿Cuál de las dos ecuaciones siguiente: En = nω o En = (n + 1/2)ω (n=0,1,2,...)
verifican las energı́as de este sistema?
8
4. ¿En cuántos puntos se anulan las funciones obtenidas en el apartado 3.4? Relacionar
el número de ceros de la función con el valor del número cuántico n de la cuestión
anterior.
5. Describir brevemente el método de diferencias finitas para resolver problemas de
contorno.
9