Aplicaciones de la ecuación de reacción-difusión

Aplicaciones de la Ecuaci´on de Reacci´on-Difusi´on Universidad de Sevilla Facultad de F´ısica Departamento de An´alisis Matem´atico Trabajo Fin de Grado Grado en F´ısica Curso 2022-2023 Dirigido por : Mar´ıa Anguiano Moreno Fdo : Mario P´erez Diego Sevilla, Mayo de 2023 A mis padres Manuel y Nieves, mi familia y a todos mis amigos que me han ayudado durante este camino a llegar hoy aqu´ı . Resumen En este documento se presenta un estudio de la ecuaci´on de Reacci´on-Difusi´on, tanto matem´aticamente c´omo sus aplicaciones en f´ısica, qu´ımica, biolog´ıa, etc.. En la primera parte veremos la parte te´orica de dicha ecuaci´on y trabajando en los espacios de funciones y gracias al Teorema de Hille-Yosida podremos garantizar la existencia y unicidad de sus soluciones. Despu´es, en la segunda parte veremos aplicaciones y ejemplos de esta ecuaci´on en distintos contextos. Aparte la resolveremos para algunos casos de inter´es. Abstract In this work, we present a study of the Reaction-Diffusion equation, mathematically and its applications in physics, chemistry, biology, etc... In the first part, we will study the theoretical aspect of this equation, working in function spaces and using the Hille- Yosida theorem we will be able to ensure the existence and uniqueness of its solutions. Then, in the second part, we will explore applications and examples in different contexts. Additionally, we will solve it for some interesting cases. . ´ Indice 1. Introducci´on 1 2. La ecuaci´on de Reacci´on-Difusi´on 3 2.1. Los espacios Lp: Definici´on y propiedades elementales ........... 3 2.2. Espacios de Hilbert .............................. 12 2.3. Teorema de Hille-Yosida ........................... 13 2.3.1. El caso Autoadjunto ......................... 21 2.4. Espacios de Sobolev .............................. 22 2.5. Ecuaci´on de Reacci´on-Difusi´on: Existencia y unicidad ........... 24 2.5.1. Ecuaci´on del calor .......................... 24 2.5.2. Ecuaci´on general del calor ...................... 28 3. Aplicaciones de la ecuaci´on de reacci´on-difusi´on 30 3.1. Ejemplos de reacciones t´ıpicas ........................ 31 3.1.1. Din´amica de poblaciones ....................... 31 3.1.2. Presa-depredador ........................... 33 3.1.3. Competici´on .............................. 35 3.1.4. Simbiosis ............................... 38 3.1.5. Reacciones qu´ımicas ......................... 40 3.1.6. Otras aplicaciones ........................... 41 3.2. Ecuaci´on de difusi´on ............................. 42 3.2.1. Ley de Flick .............................. 42 3.2.2. Soluci´on fundamental de la ecuaci´on de difusi´on homog´enea . . . 43 3.2.3. Soluci´on general de la ecuaci´on de difusi´on no homog´enea ..... 46 3.2.4. Camino aleatorio y Movimiento Browniano ............. 48 5 4. Resumen y conclusiones 51 Bibliograf´ıa 52 . Cap´ıtulo 1 Introducci´on En la naturaleza hay muchos sistemas diferentes de reacci´on-difusi´on, algunos manejables, otros no tanto, algunos simples y otros que presentan un comportamiento muy complejo. El sistema de reacci´on-difusi´on m´as simple es la propia ecuaci´on de difusi´on: ut= ∆u(1.1) d´onde u(x, t) puede representar muchas cosas, c´omo la densidad de una poblaci´on, la temperatura,etc... en cierto instante de tiempo ten un punto x∈Rn. ∆ denota el operador laplaciano con respecto a las variables espaciales x. Si tenemos (1.1) en todo x∈Rn, con la condici´on inicial apropiada u(x,0) = u0(x) (1.2) Entonces podemos resolver la ecuaci´on de reacci´on difusi´on usando m´etodos de transformadas de Fourier. Sin embargo, si tenemos (1.1) en un cierto subconjunto Ω ⊂Rn, tendremos que imponer ciertas condiciones de contorno de usobre ∂Ω, el contorno de Ω. El estudio de las soluciones de (1.1) es bien conocido y estudiado, y existen varios m´etodos para ellos, c´omo la trasformadas de Fourier, el desarrollo en autofunciones, separaci´on de variables,etc... que permiten expresar la soluci´on en t´erminos de cosas conocidas. 1 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 8 Usando (2.8): Z|αf(x) + βg(x)|pdx < Z2p−1(|αf(x)|p+|βg(x)|p)dx Por lo que obtenemos Z|αf(x) + βg(x)|pdx ≤2p−1|α|p·Z|f(x)|pdx +|β|p·Z|g(x)|pdx Se puede ver claramente que αf(x)+βg(x)∈Lpya que las constantes |α|y|β|son siempre mayores o iguales que cero y que f(x), g(x)∈Lp, por lo que cumplen individualmente (2.1). Esto demuestra que efectivamente Lpes un espacio vectorial: Z|αf(x) + βg(x)|pdx < ∞ → αf(x) + βg(x)∈Lp⇔Lpes un espacio vectorial □ Teorema 6. .- ∥ ∥Lpes una norma ∀1≤p≤ ∞. Demostraci´on.Para que sea norma tienen que cumplirse tres condiciones ∀f, g ∈Lpy k∈R: 1. ∥f(x)∥Lp≥0 si f(x)= 0 y∥f(x)∥Lp= 0 ⇒f(x) = 0 2. ∥k·f(x)∥=|k|·∥f(x)∥ 3. ∥f+g∥Lp≤ ∥f∥Lp+∥g∥Lp(Desigualdad de Minkowski) Las propiedades 1 y 2 son triviales de demostrar usando la definici´on de norma (2.2). Entonces lo que queda es probar la desigualdad de Minkowski: Los casos p= 1 y p=∞ son evidentes. Vamos a centrarnos en el caso 1 <p<∞: ∥f+g∥p Lp=Z|f+g|p−1|f+g| ≤ Z|f+g|p−1|f|+Z|f+g|p−1|g| Podemos ver que |f+g|p−1∈Lq. |f+g|p−1∈Lq⇒Z|f+g|q(p−1) =Z|f+g|p<∞ Esto es claro usando la definici´on de los exponentes conjugados (q(1 −p) = p): f+g∈Lp⇒ |f+g|p−1∈Lq Tenemos, por tanto, que se cumple la desigualdad de H¨older en los t´erminos de la derecha: Z|f+g|p−1|f| ≤ ∥f∥Lp·∥ |f+g|p−1∥Lq=∥f∥Lp·∥f+g∥p−1 Lp 8 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 9 Z|f+g|p−1|g| ≤ ∥g∥Lp·∥ |f+g|p−1∥Lq=∥g∥Lp·∥f+g∥p−1 Lp Si sumamos ambas desigualdades: ∥f+g∥p Lp≤ ∥f∥Lp·∥f+g∥p−1 Lp+∥g∥Lp·∥f+g∥p−1 Lp Dividiendo por ∥f+g∥p−1 Lpla inecuaci´on obtenemos la desigualdad de Minkowski ∥f+g∥Lp≤ ∥f∥Lp+∥g∥Lp(2.9) y por tanto tenemos probada que ∥ ∥Lpes efectivamente una norma. □ Teorema 7. Lpes un espacio de Banach (4) para todo 1≤p≤ ∞. Demostraci´on. 1) Supongamos primero que p=∞. Sea fnuna sucesi´on de Cauchy en L∞. Dado un entero k≥1 existe un Nktal que ∥fm−fn∥L∞≤1 k∀m, n ≥Nk Y existe Ekde medida cero tal que usando la definici´on de ∥ ∥L∞en (2.3) |fm(x)−fn(x)| ≤ 1 k∀x∈Ω\Ek∀m, n ≥Nk(2.10) Escribiendo E=[ k Ek(Ees de medida cero), se observa que, ∀x∈Ω\E, la sucesi´on fn(x) es de Cauchy en R. Al pasar al l´ımite en (2.1) cuando m→ ∞ se obtiene |f(x)−fn(x)| ≤ 1 k∀x∈Ω\Ek∀n≥Nk As´ı f∈L∞y∥f−fn∥L∞≤1 k∀n≥Nk. Por tanto ∥f−fn∥L∞→0. 4Un espacio de Banach es un espacio vectorial V, normado y completo sobre RoCcon una norma ||·|| tal que toda sucesi´on de Cauchy tiene un l´ımite en V. 9 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 10 2) Supongamos ahora que 1≤p < ∞. Sea fnuna sucesi´on de Cauchy en Lp. Para la demostraci´on es suficiente con probar que existe una subsucesi´on es convergente en Lp, porque una sucesi´on de Cauchy que tiene una subsucesi´on convergente es siempre convergente Se extrae una subsucesi´on (fnk) tal que ∥fnk+1 −fnk∥Lp≤1 2k∀k≥1 [Para tener esto, elegimos n1tal que ∥fm−fn∥Lp≤1 2∀m, n ≥n1, luego elegimos n2≥n1 tal que ∥fm−fn∥Lp≤1 22∀n≥n2, etc...]. Demostraremos que fnkconverge en Lp. Para una mayor simplificaci´on de la notaci´on escribiremos fken lugar de fnk, de forma que se tiene ∥fk+1 −fk∥Lp≤1 2k∀k≥1 (2.11) Definimos gn(x) = n X k=1 |fk+1 −fk|(2.12) Haciendo un peque˜no c´alculo usando la Desigualdad de Minkowski ∥gn(x)∥Lp=     n X k=1 |fk+1 −fk|    Lp≤ n X k=1 ∥fk+1 −fk∥Lp= n X k=1 1 2k Usando un resultado conocido de series geom´etricas (5), obtenemos que ∥gn(x)∥Lp≤ n X k=1 1 2k=1/2 (1 −1/2n) 1−1/2=1−1 2n≤1 Por lo que ∥gn(x)∥Lp≤1 Esto implica que la funci´on g(x) definida por l´ım n→∞gn(x) = g(x) pertenece a Lp. Para ver esto, se observa que |gn|p≤ |gn+1|p, luego gnes decreciente en Lp. Adem´as. tenemos el l´ımite puntual l´ım n→∞|gn(x)|p= n X k=1 |fk+1 −fk| p =|g(x)|p por lo cual, las condiciones del teorema de la convergencia mon´otona de Lebesgue son satisfechas (es decir, gnes creciente en Lpyg(x) = l´ım n→∞gn(x) ). 5La suma en´esima de una serie geom´etrica de razon res n X k=1 rk=r(1 −rn) 1−rsii |r|<1 10 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 11 Adem´as (6) ZΩ|g|pdx = l´ım n→∞ZΩ|gn|pdx < 1⇒g∈Lp(Ω) (2.13) Por otro lado, para n, m ≥2 tenemos |fm(x)−fn(x)| ≤ |fm(x)−fm−1(x)|+|fm−1(x)−fm−2(x)|+···|fn+1(x)−fn(x)| ≤ g(x)−gn−1(x) con g(x) = +∞ X k=1 |fk+1(x)−fk(x)|. C´omo fn(x) es de Cauchy, converge a un l´ımite finito llamado f(x). Entonces , pasando al l´ımite en la anterior desigualdad: |f(x)−fn(x)| ≤ g(x)∀n≥2 d´onde g(x) es un sumatorio con ∞t´erminos positivos. Por tanto f(x)≤g(x) + fn(x)⇒f∈Lp donde se ha demostrado que g(x)∈Lpanteriormente y fn(x)∈Lppor hip´otesis. Finalmente aplicamos el Teorema de la convergencia dominada (T.C.D) de Lebesgue fn→fpuntualmente |fn−f|p≤g(x)p            T.C.D =⇒ ∥fn−f∥Lp→0 donde g(x)p∈Lp(Ω) ya que ZΩ gpdx < +∞pues g∈Lp. Teorema 8. Sean (fn)una sucesi´on en Lpyf∈Lp, tales que ∥fn−f∥Lp→0. Entonces existe una subsucesi´on (fnk)tal que (a) fnk→f(X)c.t.p en Ω (b) |fnk| ≤ h(x)∀ky c.t.p en Ω,con h∈Lp 6Aqu´ı se ha usado que ∥gn∥Lp≤1⇒ ∥gn∥p Lp=RΩ|gn|pdx ≤1 11 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 12 Demostraci´on.El resultado es evidente para p=∞. Supongamos entonces el caso 1≤p≤ ∞. Como la sucesi´on (fn) es de Cauchy, se puede repetir la demostraci´on del Teorema 7y extraer una subsucesi´on (fnk) verificando (2.11). Desarrollando c´omo en la demostraci´on del Teorema 7se puede ver que fnkconverge c.t.p a un l´ımite, que designaremos por ¯ f(x)(7). Adem´as se tiene, gracias a (2.13) que |¯ f−fnk| ≤ g(x)∀kc.t.p en Ω,con g∈Lp(2.14) De d´onde resulta que ¯ f∈Lpy que fnk→¯ fen Lppor el teorema de Lebesgue. Por consiguiente se pueden demostrar las proposiciones (a) y(b) haciendo en cada caso lo siguiente (a) f(x) = ¯ f(x) (b) h(x) = |¯ f(x)|+g(x) 2.2. Espacios de Hilbert Los espacios de Hilbert son un concepto que permiten la generalizaci´on de los m´etodos usados en ´ Algebra lineal y C´alculo para espacios eucl´ıdeos (dimension finita) a espacios que pueden llegar a ser de dimensi´on infinita. Normalmente los espacios de Hilbert son conocidos como espacios de funciones y se definen de la siguiente forma. Definici´on. Sea H un espacio vectorial. Un producto escalar (u, v) es una forma bilineal de H×Hen R, sim´etrica, definida positiva, esto es que (u, v)=(v, u)∀u, v ∈H(Simetr´ıa) (u, v)≥0∀u, v ∈H(Positividad) (u, u)>0∀u= 0 (Positividad) (u, u)= 0 ∀u= 0 (Definida) Definici´on. Un espacio de Hilbert es un espacio vectorial H dotado de un producto escalar (u, v) y que es completo para la norma (u, u)1/2. 7A priori es necesario distinguir entre fy¯ f: sabemos que fn→fen Lpy que fnk→¯ fc.t.p en Ω 12 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 13 A continuaci´on, H designar´a siempre un espacio de Hilbert. Un ejemplo b´asico es L2(Ω), que dotado del producto escalar (u, v) = ZΩ u(x)v(x)dx (2.15) es efectivamente un espacio de Hilbert (8). Recordemos que todo producto escalar cumple la desigualdad de Cauchy-Schwarz(9): |(u, v)| ≤ (u, u)1/2(v, v)1/2∀u, v ∈H(2.16) Notemos que |u|= (u, u)1/2es una norma (10) tambi´en. [En efecto,|u+v|2=|u|2+|v|2+ 2(u, v)≤ |u|2+|v|2+ 2|u||v|]. La demostraci´on de la desigualdad de Cauchy-Schwarz es sencilla habi´endo usando la ya demostrada desigualdad de H¨older, ya que es un caso particular de esta para p=q= 2 Z|u(x)·v(x)| ≤ ∥u(x)∥L2∥v(x)∥L2→(u, v)≤Z|u|21/2Z|v|21/2 Por lo que se demuestra sustituyendo por (2.15) (u, v)=(u, u)1/2(v, v)1/2 2.3. Teorema de Hille-Yosida Definici´on. Sea A : D(A)⊂H→Hun operador lineal no acotado. Se dice que Aes mon´otono(11) si: (Av, v)≥0∀v∈D(A) yAes m´aximal mon´otono si adem´as R(I+A) = H, es decir, que: ∀f∈H, ∃u∈D(A) tal que u+Au=f Vamos a ver ahora una serie de definiciones y propiedades que necesitaremos para la demostraci´on del teorema de Hille-Yosida. 8el espacio de Sobolev H1en el cap´ıtulo 2.4 es un espacio de Hilbert ”modelado”sobre L2. 9Para establecer la desigualdad de Cauchy-Schwarz no se ha usado la hip´otesis (u, u)>0 si u= 0. 10La norma asociada a un producto escalar, a menudo, se notar´a c´omo | | en vez de ∥ ∥. 11Tambi´en se puede expresar c´omo que A es acretivo o disipativo. 13 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 14 Definiciones.-Sea Aun operador maximal mon´otono. Se dice que ∀λ > 0 Jλ= (I+λA)−1yAλ=1 λ(I−Jλ) d´onde Jλes la resolvente de AyAλes la regularizaci´on Yosida de A Proposici´on 1. Sea A un operador mon´otono. Se cumplen las siguientes propiedades: a) Aλv=A(Jλv)∀v∈Hy∀λ > 0 b) Aλv=Jλ(Av)∀v∈D(A)y∀λ > 0 c) |Aλ| ≤ |Av| ∀v∈D(A)y∀λ > 0 d) l´ım λ→0Jλv=v∀v∈H e) l´ım λ→0Aλv=Av∀v∈D(A) f) (Aλ, v)≥0∀v∈Hy∀λ > 0 g) |Aλv| ≤ 1 λ|v| ∀v∈Hy∀λ > 0 h) |Jλv| ≤ |v| ∀v∈Hy∀λ > 0 Demostraci´on.La demostraci´on de estas propiedades se puede encontrar en el cap´ıtulo VII del Ha¨ım Br´ezis [1]. Teorema 9. (Hille-Yosida).-Sea A un operador maximal mon´otono en un espacio de Hilbert H. Entonces ∀u0∈D(A)existe una ´unica funci´on u∈C1([0,+∞[ ; H)∩C([0,+∞[ ; D(A)) tal que            du dt +Au = 0 en [0,+∞[ u(0) = uo(Condici´on inicial) (2.17) Y adem´as se verifica: |u(t)| ≤ |u0|y du(t) dt =|Au(t)| ≤ |Au0| ∀t≥0 14 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 15 El inter´es principal del estudio del Teorema de Hille-Yosida se debe a que la resoluci´on del problema de evoluci´on (2.17) se reduce a comprobar que Aes maximal mon´otono, esto es, estudiar la ecuaci´on estacionaria u+λAu =f. Demostraci´on.La vamos a descomponer en 6 partes. 1°Parte: Unicidad.- Sean uy ˜udos soluciones de (2.17). Tenemos que: d dt(u−˜u) = −A(u−˜u)→d dt(u−˜u), u −˜u=−(A(u−˜u), u −˜u)≤0 Usando un resultado que nos dice que podemos escribir d dt |φ|2= 2 dφ dt , φ(12) tenemos 1 2 d dt|u(t)−˜u(t)|2=d dt(u(t)−˜u(t)), u(t)−˜u(t) por lo que d dt|u(t)−˜u(t)|2≤0 Tenemos entonces que la funci´on |u(t)−˜u(t)|es decreciente en el intervalo [0,+∞[. Como tanto uc´omo ˜uson soluciones, cumplen la misma condici´on inicial y se tiene que |u(0) −˜u(0)|= 0. Tenemos entonces una funci´on |u(t)−˜u(t)|con valor inicial 0, que es decreciente y positiva en el intervalo [0,+∞[, por lo que se deduce que |u(t)−˜u(t)|= 0 ∀t≥0 de aqu´ı se comprueba la unicidad de las soluciones. Para demostrar la existencia de u, se sustituye Apor su regularizaci´on Yosida Aλ, se establecen diversas estimaciones independientes de λy se pasa al l´ımite cuando λ→0. Sea uλla soluci´on del problema            duλ dt +Aλuλ= 0 en [0,+∞[ uλ(0) = uo∈D(A) (2.18) Se ve que uλexiste gracias al teorema de Cauchy-Lipschitz-Picard (13) aplicado a un operador F=−Aλ. 12Notar que esto solo ocurre si φ∈C1([0,+∞[ ; H) y que, entonces, |φ|2∈C1([0,+∞[ ; R) 13Para la formulaci´on y demostraci´on completa ver el Teorema VII.3 del H.Br´ezis [1], resumidamente es un teroema que garantiza la existencia y unicidad para cierto tipo de EDO. 15 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 16 2°Parte:.-Se cumple la estimaci´on  duλ(t) dt =|Aλuλ(t)| ≤ |Au0| ∀t≥0∀λ≥0(2.19) Esta desigualdad es consecuencia inmediata del Lema 2.Sea w∈C1([0,+∞[; H)una funci´on que verifica dw dt +Aλw= 0 en [0,+∞[.(2.20) Entonces las funciones |w(t)|y dw dt (t)=|Aλw(t)|son decrecientes sobre [0,+∞[ Demostraci´on.Se tiene que dw dt , w+ (Aλw, w) = 0. Pero usando la Proposici´on 1fse tiene que (Aλw, w)≥0 y por tanto 1 2 d dt|w|2≤0 Por otra parte, como Aλes un operador lineal acotado, se deduce de (2.20) que wes C∞y que d dt dw dt +Aλdw dt = 0 (2.21) Se aplica entonces lo anterior a dw dt . 3°Parte.- Demostraremos que para todo t≥0, uλ(t) converge, cuando λ→0, a un l´ımite notado por u(t); adem´as esta convergencia es uniforme en t sobre cada intervalo acotado [0, T]. En efecto sean λ, µ > 0. Se tiene que duλ dt −duµ dt +Aλuλ−Aµuµ= 0 y por consiguiente 1 2 d dt|uλ−uµ|2+ (Aλuλ−Aµuµ, uλ−uµ) = 0 (2.22) Si desarrolamos el segundo t´ermino de (2.22) y usamos las definiciones de JλyAλy la Proposici´on 1b          (Aλuλ−Aµuµ, uλ−uµ) = (Aλuλ−Aµuµ, uλ−Jλuλ+Jλuλ−Jµuµ+−uµ) = (Aλuλ−Aµuµ, Aλuλ−µAµuµ)+(A(Jλuλ−Jµuµ), Jλuλ−Jµuµ) ≥(Aλuλ−Aµuµ, λAλuλ−µAµuµ) (2.23) 16 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 17 Si usamos la desigualdad de Cauchy-Schwarz (2.16) en la expresi´on anterior obtenemos (Aλuλ−Aµuµ, λAλuλ−µAµuµ)≤ |Aλuλ−Aµuµ|·|λAλuλ−µAµuµ| Que por la desigualdad triangular |Aλuλ−Aµuµ|·|λAλuλ−µAµuµ| ≤ (|Aλuλ|+|−Aµuµ|) (|λAλuλ|+|−µAµuµ|) = =λ|Aλuλ|2+ (λ+µ)|Aλuλ||Aµuµ|+µ|Aµuµ|2 Usando la estimaci´on de (2.19) para λyµ λ|Aλuλ|2+ (λ+µ)|Aλuλ||Aµuµ|+µ|Aµuµ|2≤2(λ+µ)|Auo|2 se obtiene que (Aλuλ−Aµuµ, λAλuλ−µAµuµ)≤2(λ+µ)|Auo|2(2.24) Se tiene entonces el siguiente sistemas de ecuaciones formando por (2.22),(2.23)y(2.24)                      1 2 d dt|uλ−uµ|2+ (Aλuλ−Aµuµ, uλ−uµ)=0 (Aλuλ−Aµuµ, uλ−uµ)≥(Aλuλ−Aµuµ, λAλuλ−µAµuµ) (Aλuλ−Aµuµ, λAλuλ−µAµuµ)≤2(λ+µ)|Auo|2 (2.25) Combin´andolas se puede obtener la siguiente expresi´on 1 2 d dt|uλ−uµ|2≤2(λ+µ)|Auo|2(2.26) y al integrar |uλ(t)−uµ(t)|2≤4(λ+µ)t|Auo|2 i.e |uλ(t)−uµ(t)| ≤ 2p(λ+µ)t|Auo|(2.27) Resulta que, para cada t≥0, uλes de Cauchy, y por lo tanto, es convergente cuando λ→0 a un l´ımite que notamos por u(t). Si en (2.27) pasamos al l´ımite cuando µ→0 se obtiene |uλ(t)−u(t)| ≤ 2√λt|Auo| Consiguientemente, la convergencia es uniforme en tsobre cada intervalo acotado [0, T] yu∈C([0,+∞[, H) 17 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 24 2.5. Ecuaci´on de Reacci´on-Difusi´on: Existencia y unicidad Vamos a empezar probando la existencia y unicidad de la soluci´on de la ecuaci´on de trasmisi´on del calor sin fuentes ni p´erdidas, que es un caso particular de la ecuaci´on de reacci´on difusi´on en la que el coeficiente D, de difusi´on, es k, la constante de difusividad t´ermica, que consideraremos igual a 1 por simplicidad matem´atica, y f= 0. Despu´es veremos una serie de teoremas para el caso de a˜nadir una fuente f, que es el caso de la ecuaci´on general del calor, o ecuaci´on inhomog´enea. 2.5.1. Ecuaci´on del calor Notaci´on. Sea Ω ⊂RNun abierto de frontera Γ. Se nota Q= Ω ×]0,+∞[.Σ = Γ ×]0,+∞[ Σ es la frontera lateral del cilindro Q. Consideremos el problema siguiente. Hallar una funci´on u(x, t) : Ω ×[0,+∞[→Rque ∂u ∂t = ∆uen Q(2.37) u= 0 sobre Σ (2.38) u(x,0) = u0(x) en Ω (2.39) 24 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 25 d´onde ∆ = N X i=1 ∂2 ∂xi2designa el Laplaciano respecto de las variables espaciales,t es la variable temporal yu0(x)es una funci´on dada. La ecuaci´on (2.37) se le denomina ecuaci´on del calor, pues modela la distribuci´on de temperatura uen el dominio Ω en un instante t. La ecuaci´on del calor y sus variantes intervienen en numerosos fen´omenos de difusi´on. La ecuaci´on del calor es el ejemplo m´as sencillo de ecuaci´on parab´olica (18). La ecuaci´on (2.38) es la condici´on de contorno de Dirichlet, que se puede sustituir por la condici´on de Neumann ∂u ∂n = 0 sobre Σ (2.40) o por la condici´on de contorno tipo Robin a∂u ∂n +bu = 0 sobre Σ (2.41) D´onde nes el vector unitario de la normal exterior a Σ La condici´on (2.38) expresa que el borde de Γ de Ω se mantiene la temperatura 0; la condici´on (2.40) expresa que el flujo de calor a trav´es de Γ es 0. La ecuaci´on (2.39) es la condici´on inicial o dato de Cauchy. Resolveremos el problema (2.37) (2.38) (2.39) considerando u(x, t)c´omo una funci´on definida en [0,+∞[ con valores en un espacio H, donde H es un espacio de funciones que s´olo dependen de x; por ejemplo H=L2(Ω), o H=H1 0(Ω), etc... As´ı pues, la notaci´on u(t) designar´a un elemento de H, es decir, la funci´on x→u(x, t) con tfijo. Este planteamiento permite obtener una soluci´on del problema (2.37) (2.38) (2.39) muy f´acilmente, combinando el teorema de Hille-Yosida con los Espacios de Sobolev. Supondremos que Ω es de clase C∞con Γ acotada. Teorema 11. Supongamos que u0∈ L2(Ω). Entonces existe una ´unica funci´on que verifica (2.37) (2.38) (2.39) y u∈C0,∞; L2(Ω)∩C(]0,∞; H2(Ω) ∩H1 0(Ω),(2.42) u∈C1(]0,∞; L2(Ω).(2.43) 18La clasificaci´on tradicional de las EDP es en tres tipos: ((el´ıpticas)),((parab´olicas))y((hiperb´olicas)) 25 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 26 Adem´as u∈C∞(¯ Ω×[ε, ∞[) ∀ε > 0.(2.44) Por ´ultimo u∈L2(0,∞;H1 0(Ω)) y se verifica 1 2|u(T)|2 L2(Ω) +ZT 0|∇u(t)|2 L2(Ω)dt =1 2|u0|2 L2(Ω) ∀T>0 (19) (2.45) Demostraci´on.Se aplica el teorema de Hille-Yosida en el espacio H=L2(Ω ). Para ello se introduce el operador no acotado A:D(A)⊂H→Hdefinido por    D(A) = H2(Ω) ∩H1 0(Ω) Au=−∆u Es importante observar que la condici´on de contorno (2.38) se incluye a la definici´on del dominio de A. Se puede comprobar que Aes maximal mon´otono y autoadjunto. Entonces podemos usar el teorema 10 y ver que existe una soluci´on ´unica para (2.37) (2.38) (2.39) (2.42) y (2.43). i)A es mon´otono. En efecto, si u∈D(A) se verifica (Au, u)L2=Z(−∆u)u=Z|∇u|2⩾0. ii)A es maximal mon´otono. Basta con demostrar que R(I + A) = H = L2. Ahora bien, para toda f∈L2existe u∈H2∩H1 0soluci´on ´unica de la ecuaci´on u−∆u+f(20). iii)A es autoadjunto. Gracias a la Proposici´on 2es suficiente comprobar que A es sim´etrico. Si u, v ∈D(A) se tiene (Au, v)L2=Z(−∆u)v=Z∇u∇v (u, Av)L2=Zu(−∆v) = Z∇u∇v 19Precisemos las notaciones |u(T)|2 L2(Ω) =ZΩ|u(x, T)|2dx, |∇u(t)|2 L2(Ω) = N X t=1 ZΩ ∂u ∂xi (x, t) 2 dx. 20La demostraci´on de esto se puede ver en el Teorema IX.25 de [1]. 26 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 27 y por tanto (Au, v)=(u, Av), es decir Aes sim´etrico y por tanto autoadjunto. Hay un resultado (21) que nos dice que D(A′)⊂H2(Ω) con inyecci´on continua (22) ; m´as exactamente, se tiene D (A′) = u∈H2l(Ω); u= ∆u=. . . ∆l−1u= 0 sobre Γ Por el teorema (9) tenemos que la soluci´on ude (2.37) (2.38) (2.39) pertenece a Ck(]0,+∞; D Al ∀k, ∀l y entonces u∈Ck(]0,+∞; H2l(Ω)∀k, ∀l. Resulta (23) que u∈Ck( ]0,+∞;Ck(¯ Ω)∀k Demostremos (2.45) ; formalmente, podemos multiplicar (2.37) por uy se integra sobre Ω×]0,T. Pero, hay que tener cuidado, pues u(t) es diferenciable en ]0,+∞[ pero no en [0,+∞[. Consideremos una funci´on φ(t) = 1 2|u(t)|2 L2(Ω) :φes de clase C1en ]0,+∞[ .Tenemos que usando (2.44) obtenemos φ′(t) = u(t),du dt(t)L2 = (u, ∆u)L2=−ZΩ|∇u|2. Por consiguiente, si 0 < ε < T<∞, tenemos φ(T) −φ(ε) = ZT ε φ′(t)dt=−ZT ε|∇u(t)|2 L2(Ω) dt. Cuando ε→0, φ(ε)→1 2|u0|2 L2(Ω) y se obtiene por tanto (2.45). Haciendo hip´otesis suplementarias sobre u0, la funci´on use hace m´as regular en torno a t= 0 (recordemos que por el teorema 11 se sabe que u∈Cω(¯ Ω) ×[ε, ∞[∀ε > 0) 21Ver Teorema IX.25 en H.Brezis [1]. 22Si X⊂Ycon inyecci´on continua, entonces X⊂Yy adem´as ∥ ∥Y≤ ∥ ∥X. 23Gracias al Corolario IX.15 del H.Brezis, [1]) 27 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 28 2.5.2. Ecuaci´on general del calor Hemos probado la existencia y unicidad de la ecuaci´on del calor, ahora vamos a ver de nuevo la existencia y unicidad de sus soluciones a˜nadiendo un t´ermino f. En el contexto de la ecuaci´on general del calor, esto se puede ver como a˜nadir una fuente a nuestro sistema. Definimos QT= Ω ×]0, T[ ΣT= Γ ×]0, T[ D´onde Γ es la frontera lateral del cilindro Q, parecido al visto en la figura 3.1. Consideremos el problema          ∂u ∂t = ∆u+f(x, t) en QT u= 0 sobre ΣT(24) u(x,0) = u0(x) en Ω (2.46) Supongamos tambi´en, que Ω es acotado de clase C∞. Teorema 12. (Regularidad L2).- Dadas f∈L2(QT)yu0∈H1 0(Ω), existe una ´unica soluci´on ude (2.46) tal que u∈C[0,T]; H1 0(Ω)∩L20,T; H2(Ω) ∩H1 0(Ω) y ∂u ∂t ∈L20, T ;L2(Ω) Demostraci´on.Para la demostraci´on ver Lions-Magenes [3] o H.Br´ezis-G.Tronel [2] Con m´as generalidad, en los espacios Lpse tiene el Teorema 13. (Regularidad Lp)- Dadas f∈Lp(Ω×]0,T[ ) con 1 < p < ∞y u0= 0(25)existe una ´unica soluci´on de (2.46) tal que u, ∂u ∂t ,∂u ∂xi ,∂2u ∂xi∂xj∈Lp(Ω×]0,T[) ∀i, j 24Se podr´ıa haber dado una condici´on de contorno no homog´enea, u(x, t) = g(x, t) sobre Γ+]0, T[, pero para simplificar, nos restringimos al caso en que g= 0. 25Para una mayor simplificaci´on 28 Cap´ıtulo 2 F´ısica US, Curso 2022-2023 29 Teorema 14. (Regularidad H¨olderiana).-Sea 0< α < 1. Supongamos que f∈Cα,α/2(¯ Ω×[0, T ]) (29)y que u0∈C2+α(¯ Ω) verifican las condiciones de compatibilidad naturales u0= 0 en Γy−∆u0=f(x, 0) en Γ Entonces (2.46) tiene soluci´on ´unica tal que u, ∂u ∂t ,∂u ∂xi ,∂2u ∂xi∂xj∈Cα,α/2(¯ Ω×[0, T ]) ∀i, j Demostraci´on.Las demostraciones de los teoremas 13 y14 son complicadas (excepto el caso p= 2 del teorema 13 ). Para la demostraci´on de Lpver Grisvard [6] (Secci´on 9) yStroock-Varadhan [7]. Para el teorema 14 ver Ladyzhenskaya-Solonnikov-Uraltseva [8]. Nota 5. En este trabajo se ha visto la existencia y unicidad para un t´ermino de reacci´on f(x, t), pero tal y c´omo hemos visto en la introducci´on cabr´ıa plantearse un problema en que la fuente dependiera de la propia funci´on u, sus derivadas o incluso el Laplaciano. Estas demostraciones de existencia y unicidad no son nada sencillas ni cortas y se exige que la funci´on fcumpla ciertas propiedades para garantizar su soluci´on ´unica en los diferentes casos (32). 31Es decir, |f(x1, t1)−f(x2, t2)| ≤ C|x1−x2|2+|t1−t2|2α/2∀x1, x2, t1, t2 32Se puede ver el caso en que F(u, t) = F(u) + h(t) en [16] o incluso casos m´as complicados. 29 Cap´ıtulo 3 Aplicaciones de la ecuaci´on de reacci´on-difusi´on Una ecuaci´on de reacci´on-difusi´on, c´omo su propio nombre indica, la forman un t´ermino llamado de reacci´on y otro de difusi´on de la siguiente forma ∂u ∂t =D∆u+f(u) d´onde u=u(x, t) puede verse en algunos casos c´omo la densidad/concentraci´on de una substancia, poblaci´on... en una posici´on x∈Ω∈Rnpara un cierto tiempo t(con Ω un conjunto abierto). ∆ es el denominado operador Laplaciano. Entonces, el primer t´ermino de la derecha es el que describe la llamada ((Difusi´on)), con Duna constante llamada coeficiente de difusi´on. El segundo t´ermino de la derecha, f, es el t´ermino de ((Reacci´on)), es una funci´on suave f:R→Ry describe procesos que realmente cambian a u, es decir, algo que le ocurre a ´el (nacimiento, muerte, reacci´on qu´ımica...) y no simplemente una difusi´on en el espacio. Es tambi´en posible que el t´ermino no dependa ´unicamente de u, sino que pueda depender de ∇uo incluso tener una dependencia expl´ıcita de las coordenadas espaciales xo del tiempo t. Eso es ut=D∆u+f(u,∇u, x, t) En vez de una ecuaci´on escalar, se pueden introducir los sistemas de ecuaciones de reacci´on-difusi´on, que son del tipo    ∂u ∂t −D∆u =f(u) + Condiciones de contorno y dato inicial 30 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 31 D´onde u(x, t) es un vector de mcomponentes, Des una matriz diagonal m×myf: Rm→Rmes una aplicaci´on no lineal. Estas ecuaciones modelan fen´omenos que aparecen en campos muy variados: qu´ımica, neurofisiolog´ıa, biolog´ıa, epidemiolog´ıa, gen´etica de poblaciones, etc... 3.1. Ejemplos de reacciones t´ıpicas En esta parte, vamos a considerar ecuaciones de reacci´on t´ıpicas. Si consideramos que estamos ignorando el t´ermino de difusi´on, entonces las ecuaciones en derivadas parciales se transformar´an en ecuaciones difrenciales ordinarias (ODEs) du dt = ˙u=f(u) 3.1.1. Din´amica de poblaciones La ecuaci´on de reacci´on-difusi´ontiene su aplicaci´on en din´amica de poblaciones, ya que es usada para describir el crecimiento de estas en el espacio. Vamos a ver unos conceptos generales de din´amica de poblaciones. Generalmente, los estados estacionarios posibles, i.e ˙u= 0, y su estabilidad son de mucho inter´es, ya que corresponden a poblaciones que no cambian de tama˜no con el tiempo. Crecimiento exponencial f(u) = au con a∈Runa constante (tasa de crecimiento) y u(0) = u0la condici´on inicial. Se puede resolver f´acilemente la EDO obteni´endose el sistema de crecimiento exponencial o Malthusiano(1) u(t) = u0eat que se basa en la idea de Malthus de que todas las formas de vida son propensas al crecimiento exponencial, i.e la poblaci´on es proporcional a la velocidad con la que esta misma crece , cuando los recursos son abundantes, pero que esta misma est´a limitado por los recursos disponibles. 1En honor a Thomas Robert Malthus, pionero en el estudio de din´amica de poblaciones 31 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 32 Crecimiento log´ıstico f(u) = au 1−u K Es una mejora del crecimiento exponencial ya que a˜nade una capacidad K∈Rc´omo limitador del crecimiento, por lo que este no crecer´a infinitamente sino que se mantendr´a constante cuando u(t) = K. Podemos resolver la EDO con c.i u(0) = u0obteni´endose u(t) = K K u0−1e−at + 1 −→ l´ım t→∞u(t) = K Efecto Allee f(u) = au n K0−11−n K Este modelo parte del crecimiento log´ıstico pero a˜nade que cuando la poblaci´on es muy baja, esta se morir´a y desaparecer´a. Este fen´omeno aparece debido que, al ser tan baja la poblaci´on, su tasa de supervivencia y/o su tasa reproductiva desciende debido a la falta de individuos y la defensa del grupo ante los depredadores y por tanto a˜nade a la ecuaci´on un t´ermino adicional u K0−1. Este efecto es m´as estudiado en animales. Los sistemas anteriormente mencionados han sido descritos para el crecimiento de una ´unica poblaci´on. Por tanto, si consideramos varias poblaciones que interaccionan entre s´ı, tendremos la din´amica de varias poblaciones, un caso de mucho inter´es. El prototipo de interacci´on de dos poblaciones se puede formular generalmente c´omo un sistema de ODEs. ˙x=f(x, y) ˙y=g(x, y) Pudi´endose extender f´acilmente a sistemas de m´as de 2 poblaciones. Hay tres casos t´ıpicos de interacci´on entre dos poblaciones. 32 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 33 3.1.2. Presa-depredador Vamos a ver el origen de los sistemas presa-depredador, el sistema Lotka-Volterra. El sistema de Lotka-Volterra es de un gran inter´es hist´orico, pero, sin embargo, no es un buen modelo matem´atico ya que no es estable estructuralmente. Este modelo se propuso gracias a Umberto D’Ancona alrededor de los 1920s, un bi´ologo marino, que hizo un an´alisis estad´ıstico del pescado que era vendido en el mercado de Venecia. La pesca fue suspendida por la primera guerra mundial entre 1914 y 1918 y D’Ancona observ´o que esto coincidi´o con el aumento de la frecuencia relativa de algunas especies frente al descenso de otras. El sistema se puede describir c´omo un sistema de dos ODEs, siendo N(t) el n´umero o densidad de presas y P(t) el de depredadores (del ingl´es ((Predators))) ˙ N=rN −cNP (3.1a) ˙ P=bNP −mP (3.1b) Vamos a explicar brevemente el significado f´ısico de las constantes y t´erminos: rN : Las presas crecer´an exponencialmente en la ausencia de depredadores. −cNP : La presas disminuir´an debido a la presencia de depredadores. bNP : Los depredadores aumentar´an debido a la presencia de presas. −mP : Los depredadores disminuir´an exponencialmente en la ausencia de presas. Hacemos el siguiente cambio de variable x=b mN , y =c rP(3.2) Podemos reescribir las ecuaciones (3.1), quedando ˙x=r(1 −y)x(3.3a) ˙y=m(x−1)y(3.3b) El efecto de esta transformaci´on es reescalar el equilibrio positivo para que ocurra siempre en (x, y) = (1,1). Podemos definir, usando las ecuaciones (3.3), el siguiente cociente dy dx =m(x−1)y r(1 −y)x(3.4) 33 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 40 3.1.5. Reacciones qu´ımicas En muchos ´ambitos, las reacciones qu´ımicas son importantes, en este apartado vamos a ver unas ideas b´asica de c´omo transformar dichas reacciones en ODEs. Ley de acci´on de masas Consideramos una reacci´on irreversible simple A+Bk −→ C D´onde kes la llamada velocidad de reacci´on. Suponemos lo siguiente. El cambio del productor en el tiempo se debe al n´umero de colisiones entre las mol´eculas A y B, multiplicada por la probabilidad de que esa reacci´on ocurra (r2) en caso de la colisi´on(4). Definimos a= [A], b= [B] y c= [C] c´omo las concentraciones de las mol´eculas A,B, y C. El t´ermino r1ab∆taproximar´a el n´umero de colisiones que ocurren en un tiempo ∆t. Por lo que el cambio de en el tiempo viene dado por ∆c=abk∆t con k=r1r2, si hacemos ∆t→0 obtenemos dc dt = ˙c(t) = kab la conocida ((Ley de Acci´on de Masas))(5). Ahora vamos a ver el caso de las Reacciones reversibles, tal que A+Bk+ ⇌ k− C Con los t´erminos de reacci´on k+yk−. Asumimos que el proceso reversible (separaci´on de C en A y B) es proporcional a c. El sistema de ODEs queda: dc dt =k+ab −k−c da dt =k−c−k+ab db dt =k−c−k+ab. 4Por ejemplo, que haya suficiente energ´ıa cin´etica disponible para el inicio de la reacci´on 5Es en realidad un modelo matem´atico, no una ley. Hay que considerar tambi´en que las concentraciones son estrictamente positivas 40 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 41 Sistema general de reacciones Vamos a generalizar el caso anterior para un Sistema general de reacciones, es decir, para un n´umero rde reacciones entre sreactivos Ui(con concentraciones ui), i= 1,...s, que interaccionan a la vez: s X i=1 lijUi kj −→ s X i=1 rijUi, j = 1, . . . r. D´onde lij, rij son los coeficientes estequiom´etricos, que describen las ganancia y p´erdida del n´umero de mol´eculas Uien una reacci´on j;kj(t) es el correspondiente velocidad de reacci´on, que puede llegar a depender del tiempo debido a factores c´omo la temperatura. Definimos por tanto la funci´on ratio gj(t, u) = kj(t) s Y n=1 (un)lnj que se corresponde con la velocidad de reacci´on jy es usada para la formulaci´on de las ODEs para las uj. Se obtiene c´omo resultado final el sistema ˙ui= r X j=1 (rij −lij)gj(t, u(t)), i = 1, . . . , s, 3.1.6. Otras aplicaciones Hay otras muchas aplicaciones m´as de la ecuaci´on de reacci´on, por ejemplo: Si f(u) = u(1 −u), obtenemos la ecuaci´on de Fisher o de Kolmogorov-Petrovsky- Piscounov, utilizada para describir la expansi´on de las poblaciones biol´ogicas, [18], [19]. Si f(u) = u(1 −u2), obtenemos la ecuaci´on Newell-Whitehead-Segel para describir la Inestabilidad de Benard, [21]. Si f(u) = u(1 −u)(u−) y 0 << 1, obtenemos la ecuaci´on Zeldovich que aparece en la teor´ıa de la combusti´on, [20]. Estas son algunas de las muchas aplicaciones que hay de esta ecuaci´on. 41 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 42 3.2. Ecuaci´on de difusi´on Imagin´emonos un caso cotidiano, un terr´on de az´ucar, que al echarlo sobre un vaso de agua o caf´e, las mol´eculas de sacarosa se difunden por todo el l´ıquido. Esto y otros ejemplos nos muestran la distribuci´on espacial de mol´eculas no debe de ser homog´enea, debe existir una diferencia o gradiente de concentraci´on entre dos puntos del medio. Esta es una idea de el origen de la ecuaci´on de difusi´on. 3.2.1. Ley de Flick La forma m´as sencilla de entender la difusi´on es usando la denominada Ley de Flick o de conservaci´on de masas. Empezamos con u(x, t)∈R3, una densidad de poblaci´on o concentraci´on de una substancia dentro de un contenedor. Hay un flujo, denotado por J(x, t)∈R3, es decir, un vector que apunta en la direcci´on general de movimiento y |J(x, t)|es proporcional al n´umero de part´ıculas que van en esa direcci´on por unidad de tiempo. Elejimos un volumen cualquiera Ω con contorno Γ. Si no hay reacciones, el ´unico factor que influencia un cambio en la densidad en Ω es el flujo a trav´es de Γ. Figura 3.6: Volumen test Esto quiere decir que d dt ZΩ u(x, t)dV =−ZΓ J(x, t)dS (3.12) d´onde dV denota la integral de volumen en R3(6), dS la superficie de integraci´on (R3−1). 6Es sencillo extenderlo a Rn 42 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 43 El teorema de la divergencia de Gauss dice ZΓ J(x, t)dS =ZΩ (∇·J(x, t)) dV (3.13) Por lo que podemos reescribir (3.12) ZΩd dtu+∇·J(x, t)dV = 0 C´omo la integral se cumple para todo los vol´umenes test Ω, podemos quedarnos con d dtu+∇·J(x, t) = 0 (3.14) la conocida c´omo Ley de Flick (relaci´on entre la primera derivada temporal de la densidad y su flujo cuando asumimos la conservaci´on de masa). La Segunda Ley de Flick dice que el flujo tiene la direcci´on negativa del gradiente de la distribuci´on de part´ıculas, i.e J=−D∇u(3.15) Metiendo (3.15) en (3.14) se obtiene la ecuaci´on de difusi´on. ut=D∇u tambi´en conocida c´omo la ecuaci´on del calor. 3.2.2. Soluci´on fundamental de la ecuaci´on de difusi´on homog´enea En esta secci´on vamos a considerar la ecuaci´on de difusi´on en 1D a trav´es del eje x, i.e −∞ <x<∞yt≥0. S´olo necesitaremos una condici´on inicial, no har´an falta condiciones de contorno. Es decir nuestro problema es el siguiente ut=Duxx para −∞ <x<∞, t > 0 (3.16) u(x, 0) = ϕ(x)(3.17) En un primer paso, vamos a resolver el problema para una funci´on especial ϕ(x) y despu´es, derivaremos la soluci´on general. Para ello, vamos a usar, unas llamadas propiedades de invarianza de la ecuaci´on (3.24): 43 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 44 Proposici´on 3. Sea una soluci´on de (3.24), se cumple que a) La traslaci´on u(x−d, t)para cada soluci´on u(x, t)es tambi´en una soluci´on. b) Cualquier derivada (ux, uxx, ut) de una soluci´on es tambi´en soluci´on. c) Cualquier combinaci´on lineal de soluciones de (3.24) es soluci´on tambi´en. d) Cada integral de una soluci´on es tambi´en soluci´on. Si tenemos una funci´on S(x, t)que es soluci´on de (3.24), por lo que existe tambi´en S(x−d, t)por a) y v(x, t) = Z+∞ −∞ S(x−y, t)g(y)dy para cada funci´on g(y)si la integral converge de manera adecuada. e) Cada deformaci´on u(√ax, at)con a > 0de una soluci´on es tambi´en soluci´on. Se puede probar usando la desigualdad de la cadena. Buscamos una soluci´on particular, Q(x, t) de (3.24) que cumpla esta condici´on inicial particular (7) Q(x, 0) =    1 parax > 0 0 para x < 0 (3.18) Determinamos Q en tres pasos. Paso 1: Gracias a la propiedad e) sabemos que x→√ax,t→at, dejando Qt=DQxx y la condici´on inicial (3.18) invariante. Esto es solo posible si la dependencia de Qcon xytes del tipo x/√t. La deformaci´on lleva x √ta√ax √at =x √t. Por lo que buscamos una Q(x, t) de la siguiente forma Q(x, t) = g(p) para p=x √4Dt (3.19) D´onde ges una funci´on que tenemos que determinar y el factor √4Des un factor adicional que nos servir´a m´as tarde. Paso 2: Podemos formular una ODE para gusando (3.24) y (3.19): Qt=∂g ∂p ∂p ∂t =−x 2t√4Dtg′(p) Qx=∂g ∂p ∂p ∂x =1 √4Dtg′(p) Qxx =∂Qx ∂p ∂p ∂x =1 4Dtg′′(p) 7Esta condici´on inicial no var´ıa ante e) 44 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 45 Por lo que 0 = Qt−DQxx =−x 2t√4Dtg′(p)−1 4tg′′(p) = 1 t−p 2g′(p)−1 4g′′(p) La ODE queda: g′′(p)+2pg′(p) = 0 (3.20) Resolvi´endola para g′obtenemos g′(p) = C1e−R2p dp =C1e−p2y por tanto Q(x, t) = g(p) = C1Ze−p2dp +C2 Paso 3: Para los l´ımites de la integral elegimos Q(x, t) = g(p) = C1Zx/√4Dt 0 e−p2dp +C2 los cu´ales son v´alidos para t > 0. Teniendo en cuenta la condici´on inicial para Q(x, t), expresada c´omo un l´ımite t↘0 (8), obtenemos(9: Para x > 0 l´ım t↘0Q(x, t) = C1Z∞ 0 e−p2dp +C2=C1 √π 2+C2= 1 Para x > 0 l´ım t↘0Q(x, t) = C1Z−∞ 0 e−p2dp +C2=−C1 √π 2+C2= 0 por lo que se deduce que C1=1 √π,C2=1 2y Q(x, t) = 1 2+1 √πZx/√4Dt 0 e−p2dp para t > 0 satisface todas las condiciones que necesitamos. Por la propiedad b) tenemos que, S=∂Q ∂x es una soluci´on de (3.24) tambi´en. Para una funci´on difrenciable arbitraria ϕcon ϕ(x) = 0 para un largo |y|, definimos u(x, t) = Z+∞ −∞ S(x−y, t)ϕ(y)dy para t > 0 (3.21) De acuerdo con d),ues tambi´en soluci´on de (3.24). Usando integraci´on por partes u(x, t) = Z+∞ −∞ ∂Q ∂x (x−y, t)ϕ(y)dy =−Z+∞ −∞ ∂ ∂x[Q(x−y, t)]ϕ(y)dy = −Z+∞ −∞ Q(x−y, t)ϕ′(y)dy − :0 [Q(x−y, t)ϕ(y)]y=+∞ y=−∞ 8Esta notaci´on significa que ttiende a 0 de forma mon´otona decreciente 9Se ha hecho uso de la integral de Gauss: Z+∞ −∞ ea(x+b)2=rπ a 45 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 46 Por tanto u(x, 0) = Z+∞ −∞ Q(x−y, 0)ϕ′(y)dy =Zx −∞ ϕ′(y)dy =ϕ(x) Por lo que (3.21) se corresponde con la f´ormula de la funci´on deseada, d´onde S=∂Q ∂x =1 √4πDte−x2 4Dt para t > 0 i.e u(x, t) = 1 √4πDt Z+∞ −∞ e−(x−y)2 4Dt ϕ(y)dy (3.22) Ses la llamada ((Soluci´on Fundamental))o((funci´on de Green)). En la mayor´ıa de casos, es imposible resolver la integral en (3.22) anal´ıticamente, pero para algunos valores especiales de ϕ(x) se puede escribir de forma bonita usando la llamada ((funci´on error))(muy usada en estad´ıstica): Erf(x) = 2 √πZx 0 e−z2dz (3.23) Las propiedades principales de la funci´on error son Erf(0) = 0 l´ım x→∞Erf(x) = 1 3.2.3. Soluci´on general de la ecuaci´on de difusi´on no homog´enea En este apartado, vamos a ver la soluci´on de la ecuaci´on de difusi´on no homog´enea. Es decir queremos resolver el siguiente problema ut=Duxx +f(x, t)para −∞ <x<∞, t > 0 (3.24) u(x, 0) = ϕ(x)(3.25) d´onde ϕ(x) describe la distribuci´on espacial inicial y f(x, t) es una fuente o sumidero adicional. La soluci´on de este problema es: u(x, t) = Z+∞ −∞ S(x−y, t)ϕ(y)dy +Zt 0Z∞ −∞ S(x−y, t −s)f(y, s)dyds (3.26) Una idea para deducir esta f´ormula es partir del siguiente problema de valores iniciales para una ODE: du dt +Au(t) = f(t)A=Cte (3.27) 46 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 47 Este problema admite tiene una soluci´on del tipo u(t) = e−At +Zt 0 e(s−t)Af(s)ds (3.28) Ahora imaginemos que tenemos u(t)∈Rn,ϕ∈Rn,A∈Rn×n, es decir (3.27) corresponde a un sistema de nODEs acopladas. Para f(t) = 0 la soluci´on se escribe c´omo u(t) = S(t)ϕ, d´onde S(t) = e−At. Sin embargo, para f(t)= 0, multiplicamos (3.27) por un factor integrante S(−t) = eAt: d dt[S(−t)u(t)] = S(−t)du(t) dt +S(−t)Au(t) = S(−t)f(t) que si integramos desde 0 a t S(−t)u(t)−S(0)u(0) | {z } ϕ =Zt 0 S(−s)f(s)ds Si multiplicamos por S(t) u(t) = S(t)ϕ |{z} Soluci´on Ecuaci´on Homog´enea +Zt 0 S(t−s)f(s)ds Esta f´ormula es muy similar a (3.26), hay un t´ermino correspondiente a la soluci´on del problema homog´eneo y un segundo t´ermino que da cuenta de la inhomogeneidad. Para ver la soluci´on de la ecuaci´on inhomog´enea usamos S, ya definida en (3.22). Se le conoce tambi´en por el nombre de ((operador fuente)), ya que transforma cada funci´on ϕen una nueva funci´on. Lo ´unico que queda por demostrar es que efectivamente (3.26) es soluci´on de (3.27) para las condiciones iniciales. C´omo el t´ermino que incluye a ϕya ha sido examinado, hacemos ϕ= 0 por simplicidad pare comprobar el segundo t´ermino. Derivando (3.26) obtenemos: ∂u ∂t =∂ ∂t Zt 0Z∞ −∞ S(x−y, t −s)f(y, s)dyds =Zt 0Z∞ −∞ ∂S ∂t (x−y, t −s)f(y, s)dyds + l´ım s→tZ∞ −∞ S(x−y, t −s)f(y, s)dy, Teniendo en cuenta la singularidad de S(x−y, t−s) para t−s= 0. C´omo S(x−y, t−s) satisface la ecuaci´on de difusi´on, tenemos ∂u ∂t =Zt 0Z∞ −∞ D∂2S ∂x2(x−y, t −s)f(y, s)dyds + l´ım ε→0Z∞ −∞ S(x−y, ε)f(y, t)dy =D∂2 ∂x2Zt 0Z∞ −∞ S(x−y, t −s)f(y, s)dyds +f(x, t) = D∂2u ∂x2+f(x, t), 47 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 48 que correponde con la ecuaci´on inhomog´enea. Vamos por ´ultimo a comprobar la condici´on inicial, hacemos t→0 en (3.26) y obtenemos l´ım t→0u(x, t) = ϕ(x) + Z0 0··· =ϕ(x) Concluimos por tanto, que (3.26) es efectivamente la soluci´on deseada. 3.2.4. Camino aleatorio y Movimiento Browniano La ecuaci´on de difusi´on se puede derivar tambi´en del llamado camino aleatorio o movimiento browniano. Aqu´ı vamos a considerar un caso sencillo de 1D. Para el movimiento browniano: por unidad de tiempo τ, las part´ıculas se mueven a izquierda o derecha con una longuitud de paso ∆x, empezando en un xcualquiera. La direcci´on en la que la part´ıcula se mueva es determinada de forma aleatoria, no hay conexi´on entre los pasos (asumimos que la probabilidad de saltar a la izquierda y la derecha son iguales λl=λr=1 2). Llamamos ξ∈(−∆x, ∆x) al movimiento en el intervalo de tiempo [(k−1)∆t], k∆t, 1 ≤n y sin p´erdida de generalidad elegimos el punto de inicio en x= 0. Despu´es de un tiempo t=n∆t, el reocorrido completo de la part´ıcula ser´a xn= n X k=1 ξk. Si las part´ıcula hace rpasos por unidad de tiempo, entonces cada paso ∆xnecesitar´a ∆t=1 τde unidad de tiempo. La posici´on de la part´ıcula, xn, en un tiempo t=n∆tpuede interpretarse c´omo una variable aleatoria. Denominamos u(x, t) = P(xn=x) a la probabilidad de la part´ıcula de estar en una posici´on xen un tiempo t=n·∆t. Esta puede ser dada expl´ıcitamente por una distribuci´on binomial. Las probabilidades de las dos direcciones son iguales, por lo que podemos escribir x=m·∆x, d´onde m∈Z. Llamamos ral n´umero de pasos a la derecha y nla la izquierda dados por la part´ıcula. Obviamente n=nl+nr. Asumiendo que la part´ıcula 48 Cap´ıtulo 3 F´ısica US, Curso 2022-2023 49 se encuentra en m·∆xtras npasos, otra condici´on extra ser´a nr−nl=m. Junt´andolas se obtiene 2nr=n+m⇔nr=n+m 2 Entonces podemos aplicar ya la probabilidad correspondiente a la distribuci´on binomial P(xn=m·∆x) =  n nr 1 2nr ·1 2n−nr =1 2n · n n+m 2 =1 2n· n n+m 2 . El valor esperado < xn>y la varianza σ(xn) pueden ser determinada < xn>=*n X k=1 ξk+= n X k=1 ⟨ξk⟩= n X k=1 ((−∆x)·P(ξk=−∆x) | {z } =1 2 +(∆x)·P(ξk= ∆x) | {z } =1 2 )=0, V(xn) = V n X k=1 ξk!= n X k=1 V(ξk) = n X k=1 ((−∆x)2·P(ξk=−∆x) | {z } −1 2 +(∆x)2·P(ξk= ∆x) | {z } −1 2 ) = (∆x)2·n= (∆x)2·n t·t= (∆x)2·r·t. Obtenemos entonces que pV(xn)∼√t. Esto significa que se espera que la part´ıcula est´e en el punto inicial, debido a la simetr´ıa que existe entre la probabilidad de ir en las dos direcciones, sin embargo, a medida que pasa el tiempo t, la probabilidad de que la part´ıcula se encuentre en el punto inicial decrece. Ahora, para abordar el caso browniano, vamos a pasar del modelo discreto al continuo, haciendo tender la longuitud del paso ∆xa cero y el n´umero de pasos por unidad de tiempo r→ ∞, de manera que l´ım ∆x→0, r→∞(∆x)2r= 2D, con D= 0. El n´umero de part´ıculas en el intervalo [x, x + ∆x] cuando ha pasado un tiempo test´a descrita por u(x, t). La correspondiente ecuaci´on discreta es: u(x, t + ∆t) = u(x, t) + λru(x−∆x, t)−λru(x, t) + λlu(x+ ∆x, t)−λlu(x, t) (3.29) Podemos usar un desarrollo de Taylor de u(x, t) u(x, t + ∆t) = u(x, t) + ∂u ∂t ∆t+1 2 ∂2u ∂t2∆t2+. . . (3.30) u(x±∆x, t) = u(x, t)±∂u ∂x∆x+1 2 ∂2u ∂x2∆x2±. . . (3.31) Usando (3.30) y (3.31) en (3.29) y que λl=λr=1 2: 49