La Transformada de Laplace y sus aplicaciones a la resolución de ecuaciones diferenciales
Abstract
En este trabajo nos centraremos en el estudio de la Transformada de Laplace. Presentaremos conceptos básicos, como su definición, teoremas de existencia y unicidad y las principales propiedades, viendo después como se aplican estas para la resolución de ecuaciones diferenciales ordinarias y ecuaciones diferenciales parciales. Por último, veremos como se utiliza la Transformada de Laplace para resolver un problema realista como es la dispersión de un contaminante en la atmósfera.
Full text
Traballo Fin de Grao La Transformada de Laplace y sus aplicaciones a la resolución de ecuaciones diferenciales Alba Rey Santiago 2021/2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA
GRAO DE MATEMÁTICAS Traballo Fin de Grao La Transformada de Laplace y sus aplicaciones a la resolución de ecuaciones diferenciales Alba Rey Santiago 2021/2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA
Trabajo propuesto Área de Conocimiento: Análisis matemático. Título: La Transformada de Laplace y sus aplicaciones a la resolución de ecuaciones diferenciales. Breve descripción del contenido Se trata de estudiar las propiedades más importantes de la transformada de Laplace así como sus aplicaciones a la resolución de ecuaciones diferenciales. En particular, se hará especial hincapié en el análisis matemático y resolución de ecuaciones en derivadas parciales relacionadas con problemas medioambientales. Recomendaciones iii
Índice Resumen vii Introducción ix 1. Denición y propiedades básicas. 1 1.1. Denición........................................ 1 1.2. Existenciayunicidad.................................. 3 1.3. Propiedadesbásicas .................................. 5 1.4. Transformadainversa.................................. 9 2. Aplicación a la resolución de EDOs 13 3. Aplicación a la resolución de EDPs 17 4. Aplicación a un problema de dispersión de un contaminante. 25 4.1. Modelo Gaussian Plume ................................ 26 4.2. Solución analítica del modelo Gaussian Plume .................... 29 A. Códigos de Matlab desarrollados 33 Bibliografía 43 v
Resumen En este trabajo nos centraremos en el estudio de la Transformada de Laplace. Presentaremos conceptos básicos, como su denición, teoremas de existencia y unicidad y las principales propiedades, viendo después como se aplican estas para la resolución de ecuaciones diferenciales ordinarias y ecuaciones diferenciales parciales. Por último, veremos como se utiliza la Transformada de Laplace para resolver un problema realista como es la dispersión de un contaminante en la atmósfera. Abstract In this memory we will focus on the study of the Laplace Transform. We will present basic concepts, as its denition, existence and uniqueness theorems and the main properties, then we will see how these are applied to solve ordinary dierential equations and partial dierential equations. Finally, we will see how the Laplace Transform is used to solve a realistic problem such as the dispersion of a pollutant in the atmosphere. vii
4 1. Denición y propiedades básicas. Por la hipótesis de continuidad a trozos de f(t) , tenemos que la primera integral del segundo miembro existe. Vamos a ver que existe la segunda: Zb t0 e−ptf(t)dt≤Zb t0|e−ptf(t)|dt ≤MZb t0 e−ptectdt =MZb t0 e(c−p)tdt =−M"e(c−p)t p−c#b t0 =Me(c−p)b−e(c−p)t0 c−p. Y aplicando ahora el límite tenemos l´ım b→∞ Zb t0|e−ptf(t)|dt ≤l´ım b→∞ Me(c−p)b−e(c−p)t0 c−p=Me−(p−c)t0 p−c cuando p>c . Como la parte de la derecha de la desigualdad converge, la integral Z∞ t0 e−ptf(t)dt converge por comparación. Por lo tanto, tenemos así que la transformada de Laplace existe cuando p>c , como queríamos demostrar. Teorema 1.4 (Unicidad o teorema de Lerch) . Si F(p) = Z∞ 0 e−ptf(t)dt, p > γ, (1.1) se cumple para una función continua f(t) , donde γ es una constante, no existe otra función continua que satisfaga dicha ecuación. Antes de proceder con la demostración del resultado, es conveniente enunciar el siguiente lema, cuya demostración se puede encontrar en [2, página 23] Lema 1.5. Sea ψ(x) una función continua en [0,1] y Z1 0 xn−1ψ(x)dx = 0, n = 1,2, ... . Entonces ψ(x)=0, en 0≤x≤1. Veamos entonces la demostración del Teorema 1.4 Demostración. Suponemos que g(t) es otra función continua que satisface (1.1) y denimos h(t) = f(t)−g(t) , que por ser la diferencia de dos funciones continuas también es continua. Entonces Z∞ 0 e−sth(t)dt = 0, s ≥γ. (1.2)
1.3. Propiedades básicas 5 Sea s=γ+n , donde n es cualquier número entero positivo. Luego Z∞ 0 e−(γ+n)th(t)dt =Z∞ 0 e−nt(e−γth(t))dt =e−nt Zt 0 e−γuh(u)du∞ 0 +nZ∞ 0 e−nt Zt 0 e−γuh(u)dudt =nZ∞ 0 e−nt Zt 0 e−γuh(u)dudt, y de (1.2) se sigue que Z∞ 0 e−ntϕ(t)dt = 0, (1.3) donde: ϕ(t) = Zt 0 e−γth(t)dt. (1.4) Ahora, consideremos en (1.3) el cambio de variable x=e−t ( t=ln(1/x) , dx =−e−tdt ⇒dt = −x−1dx ): Z∞ 0 e−ntϕ(t)dt =Z0 1 xnϕ(ln(1/x)) −x−1dx =Z1 0 xn−1ϕ(ln(1/x)) dx, esto es: Z1 0 xn−1ψ(x)dx = 0, n ∈N, donde ψ(x) = ϕ[ln(1/x)] . Además, es claro que ϕ(x) es continua en el intervalo cerrado [0,1] , tomando: ψ(0) = l´ım t→∞ ϕ(t)y ψ(1) = ϕ(0) = 0. Se sigue entonces del lema 1.5 que ψ(x)≡0 en 0≤x≤1 y, por lo tanto, ϕ(t) = Zt 0 e−γth(t)dt = 0,∀t≥0. Como e−γth(t) es continua cuando t≥0 , se sigue de la igualdad anterior que e−γth(t) = 0 , ∀t≥0 , es decir h(t)=0 , cuando t≥0 . En consecuencia, g(t) = f(t) cuando t≥0 como queríamos probar. Observación 1.6 . En nuestro caso será más conveniente hablar de inyectividad, en vez de unicidad, ya que no puede haber más de una función continua f(t) que tenga la misma transformada, es decir, la Transformada de Laplace es inyectiva. 1.3. Propiedades básicas Antes de enumerar las propiedades vamos a dar la siguiente denición, que será imprescindible para el segundo teorema de traslación:
6 1. Denición y propiedades básicas. Denición 1.7 (Función escalón unitario) . La función U(t−a) se dene como sigue: U(t−a) = 0,0≤t < a, 1, t ≥a. Proposición 1.8. Las principales propiedades de la transformada de Laplace son las siguientes: i . Linealidad: L{αf(t) + βg(t)}(p) = αF(p) + βG(p). ii . Primer teorema de traslación: L{eatf(t)}(p) = F(p−a). iii . Segundo teorema de traslación: Si F(s) = L{f(t)} y a > 0 , entonces L{f(t−a)U(t−a)}(p) = e−apF(p). iv . Transformada de la derivada: L{f′(t)}(p) = pF(p)−f(0) v . Transformada de la derivada de orden n : L{f(n)(t)}(p) = pnL{f(t)}−pn−1f(0) −pn−2f′(0) −pn−3f′′(0) −... −f(n−1)(0) vi . Transformada de una integral: LZt 0 f(t)dt(p) = F(p) p. vii . Derivada de la transformada: L{(−1)ntnf(t)}(p) = F(n)(p). viii . Integral de la transformada: Lf(t) t(p) = Z∞ p F(p)dp. ix . Convolución: L{(f∗g)(t)}(p) = LZt 0 f(t−τ)g(τ)dτ(p) = F(p)G(p). Se dene la convolución como la integral del producto de ambas funciones después de desplazar una de ellas una distancia t .
1.3. Propiedades básicas 7 Demostración. i . Sean f y g dos funciones para las que existe la Transformada de Laplace y α, β ∈R . Entonces: L{αf(t) + βg(t)}(p) = Z∞ 0 e−st(αf(t) + βg(t))dt =αZ∞ 0 e−stf(t)dt +βZ∞ 0 e−stg(t)dt =αL{f(t)}(p) + βL{g(t)}(p) =αF(p) + βG(p). ii . La demostración es inmediata a partir de la denición: L{eatf(t)}(p) = Z∞ 0 e−pteatf(t)dt =Z∞ 0 e(a−p)tf(t)dt =F(p−a). iii . Expresamos Z∞ 0 e−ptf(t−a)U(t−a)dt como la suma de dos integrales y tenemos en cuenta la denición 1.6: L{f(t−a)U(t−a)}=Za 0 e−ptf(t−a)U(t−a)dt +Z∞ a e−ptf(t−a)U(t−a)dt =Z∞ a e−ptf(t−a)dt. Ahora hacemos el siguiente cambio de variable, v=t-a, dv=dt y entonces L{f(t−a)U(t−a)}=Z∞ 0 e−p(v+a)f(v)dv =e−ap Z∞ 0 e−pvf(v)dv =e−apL{f(t)}. iv . Si aplicamos integración por partes a la denición: L{f′(t)}(p) = Z∞ 0 e−ptf(t)dt =e−ptf(t)∞ 0−Z∞ 0 (d dte−pt)f(t)dt =−f(0) + pZ∞ 0 e−ptf(t)dt =pF(p)−f(0). v . Para esta demostración partimos de la propiedad anterior: L{f′(t)}(p) = pF(p)−f(0)
8 1. Denición y propiedades básicas. y aplicando inducción llegamos al resultado que queríamos: L{f(n)(t)}(p) = pnL{f(t)}−pn−1f(0) −pn−2f′(0) −pn−3f′′(0) −... −f(n−1)(0) vi . Para esta demostración vamos a aplicar la propiedad III, la transformada de la derivada: L{f(t)}(p) = Ld dt Zt 0 f(t)dt(p) =pLZt 0 f(t)dt(p)−Z0 0 f(t)dt =pLZt 0 f(t)dt(p). Despejando obtenemos: L{f(t)}(p) p=F(p) p=LZt 0 f(t)dt(p). vii . Aplicando la denición e intercambiando el orden de integración obtenemos lo siguiente: Z∞ p F(p)dp =Z∞ pZ∞ 0 e−ptf(t)dtdp =Z∞ 0Z∞ p e−ptf(t)dpdt =Z∞ 01 −te−pt∞ p f(t)dt =Z∞ 0 e−pt f(t) tdt =Lf(t) t(p). viii . Procediendo del mismo modo que en el apartado anterior tenemos: LZx 0 f(x−t)g(t)dt=Z∞ 0 e−px Zx 0 f(x−t)g(t)dtdx =Z∞ 0Zx 0 e−pxf(x−t)g(t)dtdx =Z∞ 0Z∞ t e−pxf(x−t)dxg(t)dt. Ahora, realizando el cambio de variable x=s+t , se tendrá que: Z∞ 0Z∞ t e−pxf(x−t)dxg(t)dt =Z∞ 0Z∞ 0 e−p(s+t)f(s)dsg(t)dt =Z∞ 0Zx 0 e−p(s+t)f(s)g(t)dsdt =Z∞ 0 e−psf(s)dsZ∞ 0 e−ptg(t)dt =F(p)G(p).
1.4. Transformada inversa 9 1.4. Transformada inversa Denición 1.9 (Transformada inversa de Laplace) . Llamamos transformada inversa de Laplace de F(p) a la siguiente función: f(t) = L−1{F(p)}. El siguiente resultado fue obtenido del artículo [1]. Teorema 1.10 (Existencia de la Transformada Inversa de Laplace) . Sea f: [0,∞)→R una función continua, tal que e−ptf(t)∈L1[0,∞) para algún p∈R y existan las derivadas laterales de f y sean nitas en [0,∞) , esto quiere decir que Z∞ 0|e−ptf(t)|dt < ∞ y que los siguientes límites, f′ +(x0) = l´ım h→0+ f(x0+h)−f(x0) h y f′ −(x0) = l´ım h→0− f(x0+h)−f(x0) h , existen y son nitos en [0,∞) , entonces podemos armar que existe la Transformada inversa de Laplace de una función f(t) . Ejemplo 1.11. Veamos algunos ejemplos: L−1np p2+4o= cos(2x) L−1n2 p2+4o= sin(2x) L−1n1 p2o=x Observación: En el caso de la transformada inversa de Laplace de una función F(p) puede no haber unicidad, es posible que L{f(t)}=L{g(t)} y, sin embargo, f=g . Para nuestro próposito no nos va a preocupar esto, ya que si f y g son continuas a trozos y de orden exponencial en [0,∞) , y si L{f(t)}=L{g(t)} , las funciones f y g son esencialmente iguales, esto quiere decir, que pueden diferir sólo en los puntos de discontinuidad. Proposición 1.12. Las principales propiedades de la transformada inversa de Laplace son las siguientes: i . Linealidad: L−1{αF(p) + βG(p)}=αL−1{F(p)}+βL−1{G(p)}=αf(t) + βg(t). ii . Forma inversa del primer teorema de traslación: L−1{F(p−a)}=L−1{F(p)|p→p−a}=eatf(t). donde p→p−a signica que se sustituye p por p−a en F(p) .
10 1. Denición y propiedades básicas. iii . Forma inversa del teorema de convolución. Tiene sentido aplicarla cuando f(t) y g(t) son dos funciones continuas a trozos en [0,∞) y de orden exponencial: L−1{F(p)G(p)}= (f∗g)(t) = Zt 0 f(t−τ)g(τ)dτ Demostración. i . Sean F y G dos funciones para las que existe la transformada inversa de Laplace y α, β ∈R . Entonces: L−1{αF(p) + βG(p)}=L−1αZ∞ 0 e−stf(t)dt +βZ∞ 0 e−stg(t)dt =αL−1{F(p)}+βL−1{G(p)} =αf(t) + βg(t). ii . Simplemente aplicando la denición tenemos: L−1{F(p−a)}=L−1Z∞ 0 e(a−p)tf(t)dt =L−1Z∞ 0 e−pteatf(t)dt =L−1{L{eatf(t)}} =eatf(t). iii . Sean F(p) = L{f(t)}(p) = Z∞ 0 e−ptf(t)dt y G(p) = L{g(τ)}(p) = Z∞ 0 e−pτ g(τ)dτ, entonces: L−1{F(p)G(p)}=L−1Z∞ 0 e−ptf(t)dtZ∞ 0 e−pτ g(τ)dτ =L−1Z∞ 0Zx 0 e−p(t+τ)f(t)g(τ)dtdτ =L−1Z∞ 0 g(τ)dτ Z∞ 0 e−p(t+τ)f(t)dt. Mantenemos ja la τ y realizamos el cambio de variable x=τ+t : L−1Z∞ 0 g(τ)dτ Z∞ 0 e−p(t+τ)f(t)dt=L−1Z∞ 0 g(τ)dτ Z∞ τ e−pxf(x−τ)dx
1.4. Transformada inversa 11 Puesto que f y g son continuas a trozos en [0,∞) y de orden exponencial podemos intercambiar el orden de integración: L−1Z∞ 0 g(τ)dτ Z∞ τ e−pxf(x−τ)dx=L−1Z∞ 0 e−pxdx Zx 0 g(τ)f(x−τ)dτ =L−1Z∞ 0 e−px Zx 0 g(τ)f(x−τ)dτ =L−1LZx 0 g(τ)f(x−τ)dτ =L−1{L{(f∗g)(t)}} = (f∗g)(t).
12 1. Denición y propiedades básicas.
Capítulo 2 Aplicación a la resolución de EDOs Un área importante para la aplicación de la transformada de Laplace es la resolución de ecuaciones diferenciales lineales con coecientes constantes. Este tipo de ecuación diferencial se puede reducir a una ecuación algebraica en la variable transformada. Para entenderlo, veamos el problema del valor inicial: an dny dtn+an−1 dn−1y dtn−1+... +a1 dy dt +a0y=g(t), y(0) = y0, y′(0) = y1, ..., y(n−1)(0) = yn−1, (2.1) donde ai, i = 0,1, .., n y y0, y1, .., yn−1 son constantes. Aplicando la transformada de Laplace a dicha ecuación diferencial y teniendo en cuenta la linealidad de la misma tenemos: anLdny dtn+an−1Ldn−1y dtn−1+... +a1Ldy dt +a0L{y}=L{g(t)}. (2.2) Ahora, aplicando la propiedad IV de la transformada de Laplace vista en el capítulo anterior, la ecuación (2.2) equivale a: an[pnY(p)−p(n−1)y(0) −.. −y(n−1)(0)] + an−1[pn−1Y(p)−p(n−2)y(0) −.. −y(n−2)(0)] + .. +a0Y(p) = G(p). (2.3) Despejando nos quedaría lo siguiente: [anpn−an−1p(n−1) +.. +a0]Y(p) = an[pn−1y(0) + .. +yn−1] +an−1[pn−2y0+.. +yn−2] + .. +G(p), donde Y(p) = L{y(t)} y G(p) = L{g(t)} . Despejamos Y(p) en (2.3) y ya podemos obtener y(t) mediante la transformada inversa: y(t) = L−1{Y(p)}. 13
20 3. Aplicación a la resolución de EDPs Ejemplo 3.3. Vamos a ilustrar lo anterior con el siguiente problema: x∂y ∂x +∂y ∂t +ay =bx2, x > 0, t > 0, a, b constantes, y(0, t)=0, y(x, 0+)=1. (3.6) Sea L{y(x, t)}=Y(x, p) y tomando la transformada de Laplace en ambos lados de la igualdad tenemos xYx(x, p) + pY (x, p)−y(x, 0+) + aY (x, p) = bx2 p, y utilizando una de las condiciones iniciales nos queda xdY dx + (p+a)Y=bx2 p o de otro modo, dY dx +(p+a) xY=bx p(p > 0). Resolvemos la EDO de primer orden utilizando factores de integración y tenemos Y(x, p) = bx2 p(p+a+ 2) +cx−(p+a)(x > 0, p > −a). Tomando ahora la transformada de la condición de frontera inicial nos queda Y(0, p) = L{y(0, t)}= 0, llegando así a que c= 0 . Por lo tanto, Y(x, p) = bx2 p(p+a+ 2), y por último, invirtiendo nos queda, y(x, t) = bx2 a+ 2(1 −e−(a+2)t). En la gura 3.1 mostramos una representación gráca de la solución tomando a= 1 y b= 1 en el rectángulo {(x, t)∈R2:x∈[0,10], t ∈[0,30]} .
21 (a) Líneas de contorno asociadas a la solución (b) Grafo asociado a la solución Figura 3.1: Representación gráca de la solución obtenida y(x, t) = bx2 a+2(1 −e−(a+2)t) en el rectángulo (x, t)∈[0,10] ×[0,30] . Ahora, como comentamos al principio del capítulo, nos vamos a centrar en la resolución de la ecuación del calor y la ecuación de ondas. En primer lugar vamos a mencionar el siguiente teorema, que nos será útil a continuación: Teorema 3.4. Se tendrá que: (i) L−1{e−a√p}=a 2√πt3e−a2/4t(a > 0). (ii) L−1e−a√p √p=1 √πte−a2/4t(a > 0). Ejemplo 3.5. En este caso, vamos a resolver la ecuación del calor (3.2) para c=1, con unas condiciones iniciales dadas, es decir, ∂2u ∂x2=∂u ∂t , x > 0, t > 0, para (i) u(x, 0+)=0, (ii) u(x, t) = f(t), t > 0, (iii) l´ım x→∞ u(x, t)=0.
22 3. Aplicación a la resolución de EDPs En primer lugar aplicamos la transformada en los dos lados de la igualdad L∂2u ∂x2=L∂u ∂t ⇒d2 dx2U(x, p) = pU(x, p)−u(x, 0+). Aplicando la condición inicial (i) y despejando nos queda una EDO homogénea d2 dx2U(x, p)−pU(x, p)=0 cuya solución viene dada por U(x, p) = c1e√px +c2e−√px y en vista de la condición (iii) nos queda U(x, p) = c2e−√px. Ahora, utilizando (ii) tenemos U(0, p) = L{u(0, t)}=L{f(t)}=F(p), por lo que c2=F(p) y U(x, p) = F(p)e−√px. Haciendo uso del Teorema (3.3) que vimos anteriormente y de la convolución, nos queda u(x, t) = Zt 0 x 2√πτ3e−x2/4τf(t−τ)dτ, y por último, haciendo la substitución σ2=x2/4τ , llegamos a nuestra solución u(x, t) = 2 √πZ∞ x/2√t e−σ2ft−x2 4σ2dσ. En la gura 3.2 mostramos una representación gráca de la solución tomando f(t) = sen(t) en el rectángulo {(x, t)∈R2:x∈[0,10], t ∈[0,5]} . (a) Líneas de contorno asociadas a la solución (b) Grafo asociado a la solución Figura 3.2: Representación gráca de la solución obtenida u(x, t) = 2 √πR∞ x/2√te−σ2sin t−x2 4σ2dσ en el rectángulo (x, t)∈[0,10] ×[0,5] .
23 Ejemplo 3.6. Vamos a resolver la ecuación de ondas (3.3) con unas condiciones iniciales concretas. ∂2y ∂t2=a2∂2y ∂x2, x > 0, t > 0, para (i) y(x, 0+) = 0, x > 0, (ii) yt(x, 0+) = 0, x > 0, (iii) y(0, t) = f(t) (f(0) = 0), (iv) l´ımx→∞ y(x, t)=0. Primero aplicamos la transformada en ambos lados de la igualdad L∂2y ∂t2=La2∂2y ∂x2⇒p2Y(x, p)−py(x, 0+)−∂ ∂ty(x, 0+) = a2d2Y dx2, que aplicando la condición inicial (i) y despejando nos queda como d2Y dx2−p2 a2Y= 0. Resolviendo la EDO homogénea tenemos Y(x, p) = c1e(p/a)x+c2e−(p/a)x. De la condición (iv) podemos obtener que c1= 0 , por lo tanto Y(x, p) = c2e−(p/a)x, y por la condición (iii), Y(0, p) = L{y(0, t)}=L{f(t)}=F(p) , por lo que c2=F(p) , y entonces Y(x, p) = F(p)e−(p/a)x. Y por último, invertimos haciendo uso del segundo teorema de traslación, y nos queda lo siguiente y(x, t) = ux a(t)ft−x a o lo que es lo mismo, y(x, t) = ft−x a, t ≥x a, 0, t < x a. Como podemos ver, la cuerda permanece en reposo hasta que el tiempo es t=x a , y a partir de aquí realiza el mismo movimiento que el nal en x= 0 , con un retardo de tiempo de x a .
24 3. Aplicación a la resolución de EDPs En la gura 3.3 mostramos una representación gráca de la solución tomando f(t) = sin(t) en el rectángulo {(x, t)∈R2:x∈[0,10], t ∈[0,10]} . (a) Líneas de contorno asociadas a la solución (b) Grafo asociado a la solución Figura 3.3: Representación gráca de la solución obtenida y(x, t) en el rectángulo (x, t)∈[0,10]× [0,10] .
Capítulo 4 Aplicación a un problema de dispersión de un contaminante. En este capítulo abordaremos el estudio de un problema realista asociado a la dispersión de contaminantes en la atmósfera, y para ello nos centraremos en el artículo [7]. En primer lugar vamos a contextualizar un poco dicho problema. Nos vamos a centrar en el estudio de la descripción matemática del transporte de contaminantes en la atmósfera, lo cual recibe el nombre de modelado de dispersión atmosférica. El término dispersión en este contexto se utiliza para describir la combinación de difusión, por el movimiento turbulento de remolinos; y advección, por el efecto del viento; que se da en el aire que está cerca de la supercie de la Tierra. El motivo de que nos interese especialmente este problema es que la concentración de un contaminante liberado en el aire puede describirse mediante la ecuación de advección-difusión, que es una ecuación diferencial parcial de segundo orden de tipo parabólico, como la ecuación del calor 3.2 que veíamos en el capítulo anterior. La solución exacta más simple para el problema de la dispersión atmosférica recibe el nombre de Gaussian plume , que corresponde a una fuente puntual continua que emite contaminantes en un viento unidireccional que sopla en un dominio innito. Los modelos de Gaussian plume se han aplicado en el estudio de emisiones de grandes operaciones industriales, así como en la liberación de cenizas de erupciones volcánicas, dispersión de semillas, polen e insectos o en la propagación de olores de las instalaciones ganaderas. A parte de estos casos, también hay un gran interés en aplicaciones relacionadas con la liberación de contaminantes nucleares y biológicos. 25
26 4. Aplicación a un problema de dispersión de un contaminante. Figura 4.1: Representación gráca del modelo Gaussian plume (Imagen obtenida del artículo [7]) Como podemos ver en la gura 4.2, una pluma contaminante es emitida desde una fuente puntual continua, con la dirección del viento alineada con el eje x. Los perles de concentración se dan en dos ubicaciones a favor del viento, en vertical en rojo, y en horizontal en azul, y la forma gaussiana de las secciones transversales de la pluma se muestran con respecto a la línea central, y es por esto que el modelo recibe el nombre de Gaussian plume . 4.1. Modelo Gaussian Plume Ahora que ya estamos en contexto, comenzaremos derivando la solución de Gaussian plume para la ecuación de advección-difusión, estudiando sus propiedades matemáticas y sacando conclusiones sobre su utilidad y sus limitaciones. El modelo se ilustra utilizando una versión simplicada de un escenario de emisiones industriales real en el que los contaminantes transportados por el aire se liberan de una operación de fundición. En nuestro caso nos centraremos en el transporte de un solo contaminante, cuya concentración de masa (o densidad) en el lugar x = (x, y, z)∈R3 y tiempo t≥0 pueden describirse mediante una función uniforme C(x, t) . La ley de conservación de la masa para C puede expresarse en forma diferencial como ∂C ∂t +∇· J=S, (4.1) donde S(x, t) es un término fuente o sumidero y la función vectorial J representa el ujo de masa de contaminante debido a los efectos combinados de difusión y advección. La difusión del ujo surge del movimiento de remolinos turbulentos en la atmósfera, y se puede suponer que la
4.1. Modelo Gaussian Plume 27 difusión atmosférica sigue la ley de Fick, que establece que el ujo de difusión es proporcional al gradiente de concentración o JD=−K∇C . El signo negativo nos dice que el contaminante uye desde las regiones de alta concentración a regiones de baja concentración, y el coeciente de difusión K(x) = diag(Kx, Ky, Kz) es una matriz diagonal cuyas entradas son las difusividades de remolinos turbulentos que en general son funciones de posición. La segunda contribución al ujo se debe a la simple advección lineal del viento, que puede expresarse como JA=Cu , donde u es la velocidad del viento. Juntando las dos, obtenemos el ujo total J= JD+ JA=Cu −K∇C , y después de sustituir en la ecuación de la conservación de la masa (4.1) obtenemos la ecuación tridimensional de advección-difusión ∂C ∂t +∇·(Cu) = ∇·(K∇C) + S. (4.2) El siguiente paso es hacer una serie de hipótesis que daremos por ciertas con el n de obtener una solución analítica: H1. El contaminante se emite a un ritmo constante Q desde un solo punto fuente x = (0,0, H) situado a una altura H sobre la supercie del suelo. El término fuente puede escribirse de la siguiente manera S(x) = Qδ(x)δ(y)δ(z−H), donde δ(·) es la función delta de Dirac. H2. La velocidad del viento es constante y está alineada con el eje positivo x, de modo que u = (u, 0,0) para alguna constante u≥0 . Más adelante, tomaremos una hipótesis menos restrictiva que nos permita un campo de viento uniforme que varía en el tiempo u =u(t) , con |u| ≥ 0 . H3. La solución es el estado estacionario, lo cual es razonable ya que la velocidad del viento y los demás parámetros son independientes del tiempo y la escala de tiempo que nos interesa no es suciente amplia. H4. Las difusividades de los remolinos son funciones de la distancia x a favor del viento, y la difusión es isotrópica, por lo tanto, Kx(x) = Ky(x) = Kz(x) =: K(x) . H5. La velocidad del viento es lo sucientemente alta como para que la difusión en la dirección x sea mucho más pequeña que la advección, así el término Kx∂2 xC puede despreciarse. H6. Las variaciones en la topografía son insignicantes, por lo que tomaremos la supercie del suelo como el plano z= 0 . H7. El contaminante no penetra en el suelo.
28 4. Aplicación a un problema de dispersión de un contaminante. Una vez visto esto, haciendo uso de las hipótesis H1-H6, la ecuación 4.2 se reduce a u∂C ∂x =K∂2C ∂y2+K∂2C ∂z2+Qδ(x)δ(y)δ(z−H), (4.2a) y solo nos preocupa la solución para los valores de x, z ∈[0,∞) e y∈(−∞,∞) . Para obtener un problema bien planteado, necesitamos un conjunto apropiado de condiciones de contorno, que serán las siguientes: C(0, y, z)=0, C(∞, y, z)=0, C(x, ±∞, z) = 0, C(x, y, ∞)=0. (4.2b) La primera condición es una consecuencia del viento unidireccional y la hipótesis de que no hay fuentes contaminantes para x < 0 . Las condiciones restantes en el innito son consecuentes con el requisito de que la masa total de contaminante debe permanecer nita. De acuerdo con la hipótesis H7, el ujo vertical en el suelo debe desaparecer, lo que nos da la última condición de contorno K∂C ∂z (x, y, 0) = 0. (4.2c) Con las condiciones 4.2a, 4.2b y 4.2c, tenemos un problema bien planteado para un estado estable de la concentración de contaminante C(x, y, z) . Una formulación equivalente de este teorema podemos encontrarla eliminando el término fuente de la PDE e introduciendo en su lugar una función delta en las condiciones de contorno: u∂C ∂x =K∂2C ∂y2+K∂2C ∂z2, (4.3a) C(0, y, z) = Q uδ(y)δ(z−H), (4.3b) C(∞, y, z)=0, C(x, ±∞, z) = 0, C(x, y, ∞)=0, (4.3c) K∂C ∂z (x, y, 0) = 0. (4.3d) La demostración de la equivalencia entre los problemas (4.2a)-(4.2c) y (4.3a)-(4.3d), para x > 0 , la podemos encontrar en [6, p.59]. En nuestro caso utilizaremos esta segunda forma de las ecuaciones para derivar la solución analítica.
4.2. Solución analítica del modelo Gaussian Plume 29 4.2. Solución analítica del modelo Gaussian Plume En la mayoría de libros y artículos se presenta la solución de Gaussian plume omitiendo detalles importantes de la derivación y de las hipótesis consideradas. Por ello, en nuestro caso lo vamos a resolver con suciente detalle para que el problema pueda resolverse generalizando a otras situaciones más complicadas. Los coecientes de difusión de remolinos en la capa límite de la atmósfera son funciones fuertes de la distancia a favor del viento, y varían con las condiciones climáticas y el tiempo, por lo que en la práctica son difíciles de determinar, y es por esto que suele reemplazarse x con la nueva variable independiente r=1 uZx 0 K(ξ)dξ (4.4) cuyas unidades de medida son los [m2] . Este cambio de variables elimina los coecientes K que teníamos antes, lo que nos deja el siguiente problema de coecientes constantes para c(r, y, z) := C(x, y, z) : ∂c ∂r =∂2c ∂y2+∂2c ∂z2 (4.5) Las condiciones de contorno para c son idénticas a las de C de 4.3b, 4.3c y 4.3d, con la diferencia de que x se reemplaza por r . Ahora aplicamos el método de separación de variables a 4.5, suponiendo que la dependencia de la solución de y y z se puede separar de la siguiente manera c(r, y, z) = Q ua(r, y)·b(r, z). (4.6) Así obtenemos dos problemas de dimensión reducida que tienen forma de ecuaciones de difusión de dos dimensiones: ∂a ∂r =∂2a ∂y2para 0≤r < ∞y−∞ < y < ∞, (4.7a) a(0, y) = δ(y), a(∞, y)=0, a(r, ±∞)=0, (4.7b) y
36 A. Códigos de Matlab desarrollados % Declaramos la función asociada a las condiciones de frontera syms t; f=@(t) cos(t).*sin(t); % Parte 2.- Codigo asociado a la representacion grafica % ----------------------------------------------------- % Declaramos la funcion que luego tenemos que integrar syms x sigma; fun=@(x,t,sigma) (2/sqrt(pi))*exp(-sigma.^2).*f(t-x^2./(4*sigma.^2)); % Generamos el mallado para la representacion grafica tt=linspace(Tmin,Tmax,ndivT*(Tmax-Tmin)); xx=linspace(Xmin,Xmax,ndivX*(Xmax-Xmin)); [TT,XX] = meshgrid(tt,xx); nval=size(TT); usol=zeros(nval); % Evaluamos la funcion en los nodos del mallado for i=1:nval(1) % Bucle en tiempo for j=1:nval(2) % Bucle en espacio % Comprobamos si el tiempo es cero if(TT(i,j)==0) usol(i,j)=0.0; else % Calculamos el valor de la integral en el nodo usol(i,j)=integral(@(sigma) ... fun(XX(i,j),TT(i,j),sigma),XX(i,j)/(2*sqrt(TT(i,j))),1.e4); end end end % Realizamos la representacion grafica figure contourf(XX,TT,usol); title('Solución de la ecuación del calor') xlabel('Variable espacial') ylabel('Variable temporal')
37 colorbar pause figure [M,c] = contour3(XX,TT,usol,100); c.LineWidth = 3; title('Solución de la ecuación del calor') xlabel('Variable espacial') ylabel('Variable temporal') colorbar Código para la gráca 3.3. % ********************************************************** % Script de Matlab para graficar la solución de la ecuación % de ondas del ejemplo 5.9 del libro: % % Joel L. Schiff. The Laplace transform. Undergraduate % Texts in Mathematics. Springer- Verlag, New York, 1999. % Theory and applications % % Alba Rey Santiago, 17/06/2021 % ********************************************************** clear all format long e % Parte 1.- Parametros modificables por el usuario % ------------------------------------------------ % Definimos algunos parametros que emplearemos Tmin=0; % Tiempo inicial Tmax=10; % Tiempo final de la simulacion Xmin=0; % x inicial Xmax=10; % x final ndivT=100; % Numero de divisiones de la unidad para el mallado temporal ndivX=100; % Numero de divisiones de la unidad para el mallado espacial a=1.0; % Termino asociado a la ecuacion de ondas
38 A. Códigos de Matlab desarrollados % Declaramos la función asociada a las condiciones de frontera %syms t; f=@(t) sin(t); % Parte 2.- Codigo asociado a la representacion grafica % ----------------------------------------------------- % Generamos el mallado para la representacion grafica tt=linspace(Tmin,Tmax,ndivT*(Tmax-Tmin)); xx=linspace(Xmin,Xmax,ndivX*(Xmax-Xmin)); [TT,XX] = meshgrid(tt,xx); nval=size(TT); usol=zeros(nval); % Evaluamos la funcion en los nodos del mallado for i=1:nval(1) % Bucle en tiempo for j=1:nval(2) % Bucle en espacio % Comprobamos si TT(i,j)>=XX(i,j)/a if(TT(i,j)>=XX(i,j)/a) usol(i,j)=f(TT(i,j)-XX(i,j)/a); else usol(i,j)=0.0; end end end % Realizamos la representacion grafica figure contourf(XX,TT,usol); title('Solución de la ecuación del ondas') xlabel('Variable espacial') ylabel('Variable temporal') colorbar pause figure [M,c] = contour3(XX,TT,usol,100);
39 c.LineWidth = 3; title('Solución de la ecuación del ondas') xlabel('Variable espacial') ylabel('Variable temporal') colorbar Código para la gráca 4.2. % ********************************************************** % Script de Matlab para graficar la solucion en el plano % z=0 del modelo Gaussian Plume deducido en % % J. Stockie. The mathematics of atmospheric dispersion % modeling. SIAM Review, 53(2):349372, 2011. % % Alba Rey Santiago, 29/06/2022 % ********************************************************** clear all format long e % Parte 1.- Parametros modificables por el usuario % ------------------------------------------------ % Definimos algunos parametros que emplearemos Xmin=0.0; % x inicial Xmax=10; % x final Ymin=-10.0; % y inicial Ymax=10; % y final ndivX=100; % Numero de divisiones de la unidad para el mallado espacial ndivY=100; % Numero de divisiones de la unidad para el mallado espacial % Parametros necesarios Q=10.0; K=1.0; H=1.0; u=2.0;
40 A. Códigos de Matlab desarrollados % Declaramos la funcion asociada a la solucion f=@(x,y) Q*exp(-u*(y.^2+H^2)./(4*K*x))./(2*pi*K*x); % Parte 2.- Codigo asociado a la representacion grafica % ----------------------------------------------------- % Generamos el mallado para la representacion grafica xx=linspace(Xmin,Xmax,ndivX*(Xmax-Xmin)); yy=linspace(Ymin,Ymax,ndivY*(Ymax-Ymin)); [XX,YY] = meshgrid(xx,yy); nval=size(XX); usol=zeros(nval); % Evaluamos la funcion en los nodos del mallado for i=1:nval(1) % Bucle en tiempo for j=1:nval(2) % Bucle en espacio % Comprobamos si el tiempo es cero if(XX(i,j)==0) usol(i,j)=0.0; else % Calculamos el valor de la integral en el nodo usol(i,j)=f(XX(i,j),YY(i,j)); end end end % Realizamos la representacion grafica figure contourf(XX,YY,usol); title('Solución del modelo Gaussiano') xlabel('Variable espacial x') ylabel('Variable espacial y') colorbar pause figure
41 [M,c] = contour3(XX,YY,usol,100); c.LineWidth = 3; title('Solución del modelo Gaussiano') xlabel('Variable espacial x') ylabel('Variable espacial y') colorbar
42 A. Códigos de Matlab desarrollados
Bibliografía [1] Alberto Cabada. An alternative explicit expression of the kernel of the one dimensional heat equation with dirichlet conditions. Applied Mathematics Letters , 89(4):97102, 2019. [2] A.M. Cohen. Numerical Methods for Laplace Transform Inversion . Numerical Methods and Algorithms. Springer US, 2007. [3] Brian Davies. Integral Transforms and Their Applications . Texts in Applied Mathematics 41. Springer-Verlag New York, 3 edition, 2002. [4] Joel L. Schi. The Laplace transform . Undergraduate Texts in Mathematics. Springer-Verlag, New York, 1999. Theory and applications. [5] G. F. Simmons. Dierential equations with applications and historical notes . Textbooks in Mathematics. CRC Press, Boca Raton, FL, 2017. Third edition. [6] Ivar Stakgold. Boundary value problems of mathematical physics. Vol. II . The Macmillan Company, New York; Collier-Macmillan Ltd., London, 1968. [7] J. Stockie. The mathematics of atmospheric dispersion modeling. SIAM Review , 53(2):349 372, 2011. [8] D.G. Zill and F.S. Fragoso. Ecuaciones diferenciales con aplicaciones de modelado . Thomson, 2007. 43