MAESTRÍA EN MATEMÁTICA APLICADA

Proyecto final
Andrés Felipe Ortiz Florez

Modelo Matemático
Profesor: Diego Armando Rueda-Gómez

FACULTAD DE CIENCIAS
UNIVERSIDAD INDUSTRIAL DE SANTANDER
9 de diciembre de 2023
Índice
1. Introducción 1

2. Esquema numérico 2
2.1. Planteamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.2. Conservación de propiedades físicas . . . . . . . . . . . . . . . . . . . . . . . 2
2.3. Existencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3. Implementación 10

4. Bibliografía 11
1. Introducción
En este informe, se explora un modelo matemático de competencia interespecífica con qui-
miorrepulsión, para el cual se propone un esquema numérico mixto con diferencias finitas
en tiempo y elementos finitos en espacio. Para este esquema, se prueba la conservación de
propiedades físicas, así como la existencia y unicidad de la solución en cada paso de tiem-
po. Además, se llevaron a cabo simulaciones para analizar la dinámica del modelo y validar
numéricamente los órdenes de convergencia esperados.

Descripción del Modelo Matemático: El sistema de EDP propuesto describe la evolución

temporal de las densidades de dos especies competidoras, u y v, así como la densidad de una
señal química, w:

ut = div(du∇u + χu∇w) + µ1 u(1 − u − a1 v),

vt = dv∆v + µ2 v(1 − v − a2 u),

τ wt = dw∆w − λw + αv,

para todo x ∈ Ω y t > 0, con Ω un dominio acotado con frontera suave; sujeto a las condiciones
de frontera:

∂u ∂v ∂w
= = = 0, x ∈ ∂Ω
∂ν ∂ν ∂ν

∂
donde ∂ν
representa la derivada normal con respecto a la frontera ∂Ω, con condiciones ini-
ciales:

u(x, 0) = u0 (x),

v(x, 0) = v0 (x),

w(x, 0) = w0 (x).

Además, los parámetros du , χ, µ1 , a1 , dv , µ2 , a2 , τ, dw , λ, α son constantes positivas.

1
2. Esquema numérico

2.1. Planteamiento

Para el esquema numérico, consideramos X = U × V × W , el espacio de solución del esque-

ma original. Tomamos sus respectivos espacios finito-dimensionales Xm = Um × Vm × Wm ,
los cuales se obtienen tomando una base de Schwartz (xk )∞
k=1 de X y definiendo Xn =

spam{x1 , . . . , xm }.

Inicialización: (u0 , v0 , w0 ) = (u(0), v(0), w(0)).

n-ésimo paso de tiempo: Dado (un−1 , vn−1 , wn−1 ) ∈ Xm hallar (un , vn , wn ) ∈ Xm tal que

 (δ u , ū) + du (∇un , ∇ū) + χ (un ∇wn , ∇ū) − µ1 (un , ū) + µ1 ([un ]2 , ū) + a1 µ1 (un vn , ū) = 0, ∀ū ∈ Un ,
 t n


(δt vn , v̄) + dv (∇vn , ∇v̄) − µ2 (vn , v̄) + µ2 ([vn ]2 , v̄) + a2 µ2 (un vn , v̄) = 0, ∀v̄ ∈ Vn ,


 τ (δ w , w̄) + d (∇w , ∇w̄) + (w , w̄) − α(v , w̄) = 0, ∀w̄ ∈ W .

t n w n n n n
(1)

2.2. Conservación de propiedades físicas

Para el esquema propuesto, la conservación de las propiedades físicas consiste en mostrar que
si los datos iniciales u0 , v0 , y w0 son positivos, entonces para cada paso de tiempo n, se tiene
que un , vn , y wn también son positivos. Para lograr esto, consideramos el siguiente esquema
modificado:

 (δ u , ū) + du (∇un , ∇ū) + χ ((un )+ ∇wn , ∇ū) − µ1 (un , ū) + µ1 ((un )2+ , ū) + a1 µ1 (un vn , ū) = 0, ∀ū ∈ Un ,
 t n


(δt vn , v̄) + dv (∇vn , ∇v̄) − µ2 (vn , v̄) + µ2 ((vn )2+ , v̄) + a2 µ2 ((un )+ vn , v̄) = 0, ∀v̄ ∈ Vn ,


 τ (δ w , w̄) + d (∇w , ∇w̄) + (w , w̄) − α(v , w̄) = 0, ∀w̄ ∈ W ,

t n w n n n n
(2)

luego, demostramos que toda solución de (2) es positiva, ya que toda solución positiva de
(2) es también solución de (1); además vamos asumir que 1 > ∆t máx{µ1 , µ2 }. Sea n > 0

2
y supongamos que un−1 , vn−1 , y wn−1 son positivos, veamos que un , vn , y wn también son
positivos.
Positividad de vn : En (2) tomamos v̄ = rh ((vn )− ), donde rh es el interpolador nodal de
elementos finitos P1 y aplicamos mass-lumping a todos los términos de la ecuación de v,
excepto a dv (∇vn , ∇rh ((vn )− )) para obtener:

0 = (δt vn , rh ((vn )− ))h + dv (∇vn , ∇rh ((vn )− )) − µ2 (vn , rh ((vn )− ))h

(3)
+ µ2 ((vn )2+ , rh ((vn )− ))h + a2 µ2 ((un )+ vn , rh ((vn )− ))h .

Primero estimamos (δt vn , rh ((vn )− ))h :

Z
h
(δt vn , rh ((vn )− )) = rh ((δt vn ) (rh ((vn )− )))
Ω
 n
v − v n−1
Z 
= rh · rh (vn )−
Ω ∆t
Z
1 

= rh (vn rh ((vn )− )) − rh vn−1 rh (vn )− ,
∆t Ω

como rh conserva el signo, vn−1 ≥ 0 y (vn )− ≤ 0 entonces −rh vn−1 rh (vn )− ≥ 0, por lo
tanto
1
(δt vn , rh ((vn )− ))h ≥ (vn , rh (vn )− )h . (4)
∆t

Ahora estimamos dv ∇vn , ∇rh (vn )− ): De la inentidad
Z Z Z
∇rh [vn ]+ ∇rh [vn ]− + ∥∇rh (vn )∥ ∥22






∇vn ∇rh [vn ]− = ∇rh [vn ]− ∇rh [vn ]− +
Ω Ω Ω

R 2
y de Ω
∇rh [vn ]− · ∇rh [vn ]− = ∇rh [vn ]− 2
≥ 0 obtenemos

∇vn , ∇rh [vn ]− ≥ ∇rh [vn ]+ , ∇rh [vn ]− ,

3
note que
Z n+1
X Z
∇rh [vn ]+ ∇rh [vn ]− = [vn ]+ (xi ) [vn ]− (xj ) ∇ϕi ∇ϕi ,
Ω i,j=0 Ω

si i = j entonces [vn ]+ (xi ) [vn ]− (xj ) = 0 y si i ̸= j entonces [vn ]+ (xi ) [vn ]− (xj ) ≤ 0 y
R
Ω
∇ϕi ∇ϕi ≤ 0, por lo tanto


dv ∇vn , ∇rh (vn )− ) ≥ 0. (5)

h
Ahora estimamos µ2 [vn ]2+ , rh (vn )− :

" n+1
n+1 X
Z X #

h
[vn ]2+ , rh (vn )− = [vn ]2+ (xj )(vn )− (xi ) ϕi (xj ) ϕj (x),
Ω j=0 i=0

note que si i ̸= j entonces ϕi (xj ) = 0 y si i = j entonces [vn ]2+ (xj )(vn )− (xi ) = 0, por lo tanto

h
µ2 [vn ]2+ , rh (vn )− = 0. (6)

Estimemos a2 µ2 ((un )+ vn , rh ((vn )− ))h :

Z n+1
!
X
((un )+ vn , rh ((vn )− ))h = rh u+ vn (vn )− (xi ) ϕi (x)
Ω i=0
n+1
Z X
= [(un )+ (xj )vn (xi )(vn )− (xi ) ϕi (xj )] ϕj (x),
Ω i,j=0

si i ̸= j entonces ϕi (xj ) = 0 y si i = j entonces (un )+ (xj )vn (xi )(vn )− (xi ) ϕi (xj ) ≥ 0, por lo
tanto

((un )+ vn , rh ((vn )− ))h ≥ 0. (7)

4
De (3)-(7) y 1 > ∆t máx{µ1 , µ2 } se tiene que (vn , rh (vn )− )h ≤ 0, y como
Z

h
vn , rh [vn ]− = rh [vn ]2− ,
Ω

entonces rh [vn ]2− = 0, lo que implica que [vn ]− = 0. En consecuencia

vn ≥ 0. (8)

Positividad de un : En (2) tomamos ū = rh ((un )− ), donde rh es el interpolador nodal

de elementos finitos P1 y aplicamos mass-lumping a todos los términos de la ecuación de u,
excepto a du (∇un , ∇rh ((un )− )) para obtener:

0 = (δt un , rh ((un )− ))h + du (∇un , ∇rh ((un )− )) + χ ((un )+ ∇wn , ∇rh ((un )− ))h
(9)
− µ1 (un , rh ((un )− ))h + µ1 ((un )2+ , rh ((un )− ))h + a1 µ1 (un vn , rh ((un )− ))h .

De manera análoga obtenemos

1
(δt un , rh ((un )− ))h ≥ (un , rh (un )− )h ,
∆t
du (∇un , ∇rh ((un )− )) ≥ 0,
(10)
µ1 ((un )2+ , rh h
((un )− )) = 0,

a1 µ1 (un vn , rh ((un )− ))h ≥ 0.

Note que
" n+1 #
Z n+1 X
((un )+ ∇wn , ∇rh ((un )− ))h = ∇wn (xj ) (un )+ (xj ) ∇(un )− (xi ) ϕi (xj ) ϕj (x),
Ω i,j=0

5
si i = j entonces (un )+ (xj ) ∇(un )− (xi ) = 0 y si i ̸= j entonces ϕi (xj ) = 0, por lo tanto

((un )+ ∇wn , ∇rh ((un )− ))h = 0. (11)

De (9)-(11) y 1 > ∆t máx{µ1 , µ2 } se tiene que (un , rh (un )− )h ≤ 0, nuevamente, esto implica
que
un ≥ 0. (12)

Positividad de wn : En (2) tomamos w̄ = rh ((wn )− ) y aplicamos mass lumping para obtener

τ (δt wn , rh ((wn )− ))h + dw (∇wn , ∇rh ((wn )− )) + (wn , rh ((wn )− ))h − α(vn , rh ((wn )− )) = 0.
(13)

En este caso también se tiene que

1
(δt wn , rh ((wn )− ))h ≥ (wn , rh (wn )− )h ,
∆t (14)
dw (∇wn , ∇rh ((wn )− )) ≥ 0.

Además, como α > 0, vn ≥ 0 y rh ((wn )− ) ≤ 0, entonces

−α(vn , rh ((wn )− )) ≥ 0. (15)

De (13)-(15) se tiene que (wn , rh (wn )− )h ≤ 0, nuevamente esto implica que

wn ≥ 0. (16)

6
2.3. Existencia

Para probar que en cada paso de tiempo existe solución de (1), usamos el teorema de punto
fijo de Leray-Schauder. Para esto definimos el operador R : Xm → Xm por R(ūn , v̄n , w̄n ) =
(un , vn , wn ) de manera que (un , vn , wn ) resuelva los siguientes problemas lineales desacoplados

 (δ u , ū) + du (∇un , ∇ū) + χ (ūn ∇w̄n , ∇ū) − µ1 (ūn , ū) + µ1 ([ūn ]2 , ū) + a1 µ1 (ūn v̄n , ū) = 0, ∀ū ∈ Un ,
 t n


(δt vn , v̄) + dv (∇vn , ∇v̄) − µ2 (v̄n , v̄) + µ2 ([v̄n ]2 , v̄) + a2 µ2 (ūn v̄n , v̄) = 0, ∀v̄ ∈ Vn ,


 τ (δ w , w̄) + d (∇w , ∇w̄) + (w̄ , w̄) − α(v̄ , w̄) = 0, ∀w̄ ∈ W .

t n w n n n n
(17)
En este caso como estamos en espacios de dimensión finita solo debemos probar que R esta
bien definido, los puntos fijos de λR están limitados para 0 < λ ≤ 1 y que R es continuo.

1. R está bien definido. Para mostrar que R está bien definido, debemos probar que
dado (ūn , v̄n , w̄n ) existe (un , vn , wn ) solución de (17). Para mostrar la existencia de un ,
definimos la forma bilineal au : Um × Um → R por

1
au (un , ū) = (un , ū) + du (∇un , ū),
∆t

y el funcional Fu : Um → R por

1
Fu (ū) = (un−1 , ū) − χ (ūn ∇w̄n , ∇ū) + µ1 (ūn , ū) − µ1 ([ūn ]2 , ū) − a1 µ1 (ūn v̄n , ū).
∆t

Se puede verificar que au es bilineal, continua y coercitiva, y que Fu es lineal y continuo.

Por lo tanto, del lema de Lax-Milgram se sigue que existe una única un que satisface
(17). Para probar la existencia de vn y wn también usamos el lema de Lax-migram
definiendo las formas bilineales av : Vm × Vm → R y aw : Wm × Wm → R por

1
av (vn , v̄) = (vn , v̄) + dv (∇vn , v̄),
∆t
τ
aw (wn , w̄) = (wn , w̄) + dw (∇wn , w̄),
∆t

7
 respectivamente, y los funcionales Fv : Vm → R y Fw : Wm → R por

1
(vn−1 , v̄) − a2 µ2 (ūn v̄n , v̄) − µ2 (v̄n , v̄) + µ2 ([v̄n ]2 , v̄) ,


Fu (ū) =
∆t
τ
Fw (w̄) = (wn−1 , w̄) − (w̄n , w̄) + α(v̄n , w̄),
∆t

respectivamente.

2. Ahora, probemos que todos los posibles puntos fijos de λR, con λ ∈ [0, 1], están acota-
dos. Sea (un , vn , wn ) un punto fijo de λR, entonces

λ 1
(un−1 , ū) = (un , ū) + du (∇un , ∇ū) + χλ (un ∇wn , ∇ū) − µ1 λ(un , ū) + µ1 λ([un ]2 , ū)
∆t ∆t
+ λa1 µ1 (un vn , ū),
λ 1
(vn−1 , v̄) = (vn , v̄) + dv (∇vn , ∇v̄) − µ2 λ(vn , v̄) + µ2 λ([vn ]2 , v̄) + a2 µ2 λ(un vn , v̄),
∆t ∆t
λτ τ
(wn−1 , w̄) = (wn , w̄) + dw (∇wn , ∇w̄) + λ(wn , w̄) − αλ(vn , w̄),
∆t ∆t
(18)

para todo (ū, v̄, w̄) ∈ Xm . Para este sistema se conservan las propiedades de positividad,
por lo tanto podemos asumir que un , vn y wn son no negativos. Tomando (ū, v̄, w̄) =
(un , vn , wn ) en la ecuación (18) obtenemos
 
λ 1
(un−1 , un ) = − µ1 λ ∥un ∥22 + du ∥∇un ∥22 + χλ (un ∇wn , ∇un ) + µ1 λ([un ]2 , un )
∆t ∆t
+ λa1 µ1 (un vn , un ),
 
λ 1
(vn−1 , vn ) = − µ2 λ ∥vn ∥22 + du ∥∇vn ∥22 + µ2 λ([vn ]2 , vn ) + a2 µ2 λ(un vn , vn ),
∆t ∆t
λτ  τ 
(wn−1 , w̄) = + λ ∥wn ∥22 + dw ∥∇wn ∥22 − αλ(vn , wn ).
∆t ∆t

Como un , vn y wn son no negativos, entonces los términos µ1 λ([un ]2 , un ), λa1 µ1 (un vn , un ),

8
 µ2 λ([vn ]2 , vn ), a2 µ2 λ(un vn , vn ) y −αλ(vn , wn ) son no negativos, por lo tanto
 
λ 1
(un−1 , un ) ≥ − µ1 λ ∥un ∥22 + du ∥∇un ∥22 ,
∆t ∆t
 
λ 1
(vn−1 , vn ) ≥ − µ2 λ ∥vn ∥22 + dv ∥∇vn ∥22 ,
∆t ∆t
λτ  τ 
(wn−1 , w̄) ≥ + λ ∥wn ∥22 + dw ∥∇wn ∥22 ,
∆t ∆t

de las desigualdades de Young y Holder se sigue que

Cϵ λ
∥un−1 ∥22 ≥ c∥un ∥2H 1 ,
∆t
Cϵ λ
∥vn−1 ∥22 ≥ c∥vn ∥2H 1 ,
∆t
Cϵλτ
∥wn−1 ∥22 ≥ c∥wn ∥2H 1 ,
∆t
 1−ϵ
donde c = mı́n ∆t
− µ1 λ, 1−ϵ
∆t
τ
− µ2 λ, ∆t + λ , du , dv , dw y ϵ es lo suficientemente pe-
1−ϵ
queño para que ∆t
−µ1 λ, 1−ϵ
∆t
−µ2 λ y τ
∆t
+λ sean no negativos. Como (un−1 , vn−1 , wn−1 )
es dado entonces existe una constante C > 0 tal que

C ≥ ∥un ∥2H 1 + ∥vn ∥2H 1 + ∥wn ∥2H 1 .

Por lo tanto los puntos fijos de λR están limitados.


3. R es continuo. Sea (ūl , v̄ l , w̄l ) l∈N
una sucesión tal que

(ūl , v̄ l , w̄l ) → (ū, v̄, w̄). (19)

  l l l
Por lo tanto, (ūl , v̄ l , w̄l )
l∈N
está acotada, y en consecuencia (u , v , w ) = R(ūl , v̄ l , w̄l ) l∈N
 r r r
está acotada. Luego, existe una subsucesión R(ūl , v̄ l , w̄l ) r∈N tal que

r r r
R(ūl , v̄ l , w̄l ) → (u′ , v ′ , w′ ). (20)

Así, a partir de (19) y (20), un paso estándar al límite cuando r → +∞ en (17) permite

9
 r r r
deducir que R(ūl , v̄ l , w̄l ) = (u′ , v ′ , w′ ). Por lo tanto, cualquier subsucesión convergente

de R(ūl , v̄ l , w̄l ) l∈N converge a R(ū, v̄, w̄). A partir de la unicidad de R(ū, v̄, w̄), se
concluye que R(ūl , v̄ l , w̄l ) → R(ū, v̄, w̄). Así, R es continua.

3. Implementación
Dado que el esquema (1) es no lineal aproximamos su solución usando el método de Picard.
Dado (un−1 , vn−1 , wn−1 ) ∈ Xm (asumiendo (u0n , vn0 , wn0 ) = (uk−1 k−1 k−1
n , vn , wn ) en el primer paso

de iteración) hallar (ukn , vnk , wnk ) tal que







 δt ukn , ū + du ∇ukn , ∇ū + χ ukn ∇wnk−1 , ∇ū − µ1 (ukn , ū) + µ1 ([unk−1 ]2 , ū) + a1 µ1 (ukn vnk−1 , ū) = 0,



δt vnk , v̄ + dv ∇vnk , ∇v̄ − µ2 (vnk , v̄) + µ2 ([vnk−1 ]2 , v̄) + a2 µ2 (uk−1 k
n vn , v̄) = 0,


 τ δ wk , w̄
+ d ∇wk , ∇w̄
+ (wk , w̄) − α(v k−1 , w̄) = 0,

t n w n n n
(21)
todo (ū, v̄, w̄) ∈ Xm , hasta que se cumpla
para   el siguiente criterio de parada:
∥ukn −uk−1
n ∥2 ∥vnk −vnk−1 ∥2 ∥wnk −wnk−1 ∥2
máx , vk−1 , ≤ tol.
∥uk−1
n ∥2 ∥ n ∥2 ∥wnk−1 ∥2

10
4. Bibliografía
Allaire, G. (2007). Numerical analysis and optimization: an introduction to mathematical
modelling and numerical simulation. OUP Oxford.
Guillén-González, F., Rodríguez-Bellido, M. A., Rueda-Gómez, D. A. (2020). Study of a
chemo-repulsion model with quadratic production. Part II: analysis of an unconditionally
energy-stable fully discrete scheme. Computers Mathematics with Applications, 80(5), 636-
652.

11