Universidad De Guadalajara

Centro Universitario De Ciencias Exactas e Ingenierı́as

División de Ciencias Básicas

Método de diferencias finitas con funciones de base radial

aplicado a un modelo espacial para la distribución de microalgas
en un sistema marino

Tesis que para obtener el tı́tulo de

Licenciado en matemáticas
presenta
Francisco Seneca Bryant Morales

Director: Dr. Juan Antonio Licea Salazar

Guadalajara, Jalisco. Enero 2023.

 I

Los seres humanos, institucionalmente, hemos ido conformando el mundo

para que sea matemático [...] El libro de la naturaleza no viene escrito en
caracteres matemáticos: son los cientı́ficos los que lo escriben usando
caracteres matemáticos.
Carlos Madrid, conferencia: ¿Qué es la filosofı́a de las matemáticas?
Agradecimientos

Resulta complicado para una persona que aprecia tanto la presencia del otro en su vida y que cree que
los demás dan razón a su existencia escribir unos agradecimientos que estén al nivel de la importancia
que tienen esas personas que me han marcado (y más si lo va a escribir por primera vez). Escribir, por
ejemplo, ((gracias, amigo(a) por estar allı́, te quiero mucho)) resulta insuficiente para describir el hecho de
que esa persona es, en buena medida, responsable de que hayas pasado algunos de los mejores momentos
de tu vida y que probablemente, si no hubiera estado allı́, tu paso por aquella institución educativa serı́a un
mero trámite sin sentido que incluso podrı́a quedar inconcluso. Sin embargo, es necesario escribirlo por dos
razones: por tradición (toda tesis lleva una) y porque a esas personas se les tiene que agradecer de todas las
formas posibles (la escrita es una y que mejor que un((trabajo especia))). Escribiré los agradecimientos a las
personas de los ámbitos que considero más importantes en mi vida, a saber: familia, profesores y amigos (a
partir de secundaria); el resto de gente también es valiosa, no obstante son relaciones esporádicas o no han
tenido la misma profundidad, además de que este agradecimiento se volverı́a gigante. Habrá descripciones
que serán más largas para unas personas que para otras, esto no significa que se estime más a una que otras;
simplemente no se me ocurren palabras para escribir todo lo que quiero expresar, o no me quiero repetir, o
una de ellas me dio permiso de conocerla/convivir más.

Familia

Realmente solo existen 4 personas que estén relacionadas genéticamente conmigo y pueda afirmar que
marcaron mi vida: Mi madre (Marı́a de los Ángeles), mi padre (Francis Edward), mi abuela materna (Marı́a
de la Luz) y mi hermana (Sabina Amor). Los primero dos, claramente, me trajeron a la vida y me cuidaron,
pero también formaron a una persona con gusto por enseñar lo que sabe, las ciencias, artes y filosofı́a, que
trata de comprender al resto y ser cercano a los discriminados y subalternos; a la tercera por cuidarme cuando
yo no era consciente de mi existencia, sin esperar nada a cambio, y por dotar de sentido a mi vida cuando
yo tuve que cuidarla; y la cuarta es, simplemente, la persona más importante, no sé qué harı́a si me faltara.

II
 III

Profesores
Un profesor se supone que da herramientas y te guı́a a través de cierta rama del conocimiento. Los hay
malos, regulares, buenos y que son un parteaguas de tu vida. De los últimos solo he tenido cuatro, dos en
secundaria y dos en universidad. Los de secundaria son el profesor Tovar (artes plásticas) y el profesor Piña
(matemáticas I), el primero me dio otra visión de las artes y la docencia, y el segundo es el culpable principal
de que ahora esté escribiendo los agradecimientos de una tesis para la licenciatura en matemáticas: me hizo
enojar tanto que tuve que estudiar matemáticas de forma obsesiva, descubriendo ası́ que tenı́a una facilidad
y gusto por ellas; los de universidad son el profesor Casillas (Cálculo I y II) y el profesor Licea (Teorı́a de
ecuaciones diferenciales II y director de esta tesis), de ellos agradezco que dieron clases excelentes, fueron
comprensivos pero disciplinados, y fueron muy cercanos y respetuosos en el trato: el primero con amistad
y discusiones interesantes, el segundo dirigiendo esta tesis y apoyando lo más posible en cualquier asunto
relacionado con la escritura de textos académicos.

Amigos
En secundaria, preparatoria y universidad he tenido las relaciones más fuertes. Algunas son presentes,
otras ya son esporádicas y otras solo están como recuerdos bastante lindos.

Quizás pueda decir que el tercer año de secundaria es el mejor año de mi vida en términos escolares:
iba contento a la escuela porque sabı́a que tanto la mayorı́a de las clases como las pláticas con mis
amigos serı́an espectaculares. Aquı́ agradezco a Edgar León e Iván Vergara por tener las pláticas más
jocosas y forjar gran parte de mi humor, y a Grecia Rodrı́guez, Lilia Álvarez y Kate Adamary por
ser las chicas con las que mejor me lleve en su tiempo y las primeras personas con las que forjé una
amistad profunda.

En la preparatoria desde primero tuve grandes amigos, entre ellos Gabriel, Nilton Palomera, Alfredo
Flores y Martı́n Arguelles, quien es el mejor amigo que podrı́a tener. Particularmente quiero agradecer
a este último, ya que la convivencia ha sido tan cercana, respetuosa, incondicional, inteligente y sana
que su existencia me produce esperanza y admiración (más conociendo su historia). Más tarde se
unieron al grupo lindas personas como Diana, Karla Cibrı́an (y su hermana Kenia, aunque no estuviera
en la preparatoria), Aylin Rodriguez y Rosalba Hernández. Aunque hayan llegado o desarrollado la
amistad después, también son importantı́simas para mı́.

En la universidad, agradezco a dos grupos de personas: a los que cayeron en primer semestre pero la
amistad fue tan fuerte que aún los recuerdo, y mis amigos de toda la universidad. Del primer grupo,
están Gamaliel y Keliev; del segundo, Irma Yolanda, Renata Godı́nez, Miriam Rolón, Gustavo Álvarez
y Stephanie Esquivel. Esta última tiene un papel en mi vida análogo al de Martı́n en la preparatoria.
IV

Ahora, la persona más cercana y que más amo en estos momentos: Vika Ogas. Eres una mujer increı́ble,
inteligente y trabajadora, a la cual estimo y admiro mucho. Eres la persona que más me ha enseñado y me
ha hecho reı́r en el último año y medio. Nos hemos apoyado entre los dos y hemos visto el desarrollo para
bien (y a veces para mal) de la relación y de cada uno. Momentos desde hermosos hasta extraños contigo
se vuelven siempre un tema de conversación, y obviamente tengo muchı́simas ganas de seguir pasando ese
tipo de cosas contigo. Solo sé que quiero pasar más tiempo contigo.
Por último, agradecer a esas personas que explı́citamente no vienen en este agradecimiento pero me
ayudaron en cualquier cosa o me hicieron reir en algún moemento.
Índice general

1. Introducción 1

2. Marco teórico 3
2.1. Sistemas de reacción-difusión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1.1. Análisis de inestabilidad de Turing . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2. Interpolación con funciones de Base Radial . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1. Conceptos básicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.2. Invertibilidad de matriz de la interpolación. . . . . . . . . . . . . . . . . . . . . . . 13

3. Interacción plancton-nutriente 18
3.1. Ecosistemas marinos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.1.1. Cadena trófica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2. Modelo espacial para la interacción plancton-nutriente . . . . . . . . . . . . . . . . . . . . 21
3.2.1. Análisis de inestabilidad de Turing. . . . . . . . . . . . . . . . . . . . . . . . . . . 25

4. Solución numérica 26
4.1. Introducción al método DF-FBR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.2. DF-FBR por stencils . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.2.1. Solución de la parte espacial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.2.2. Solución de la parte temporal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

5. Resultados 31
5.1. Simulaciones numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.1.1. Primera simulación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.1.2. Segunda simulación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

6. Conclusiones 38

V
VI ÍNDICE GENERAL

A. Método de diferencias finitas 39

B. Código en Python 43

Referencias. 50
Capı́tulo 1

Introducción

En cualquier sistema marino con cierta profundidad, en las partes más altas, donde la luz solar alcanza a
penetrar, encontramos una dinámica peculiar entre el plancton y los nutrientes. Resulta importante el estudio
de lo anterior ya que el plancton es la base de la cadena trófica (marina) y tiene utilidades humanas valio-
sas, como la producción de bio-diesel; mientras que los nutrientes es uno de los principales componentes
abióticos de cualquier ecosistema marino y su concentración en este está relacionado directamente con la
concentración de CO2 en la atmósfera [27, 31]. En la literatura cientı́fica, uno de los modelos matemáticos
que ha tenido mayor eficacia para predecir el comportamiento de la interacción plancton-nutriente son los
sistemas de reacción-difusión (ver [2, 6, 10, 31]).
Los sistemas de reacción-difusión son modelos matemáticos que describen cómo una o más sustancias
distribuidas en un espacio determinado cambian bajo la influencia de dos procesos: la difusión (movimiento
local) y la reacción (crecimiento, interacciones, cambios de estado) [24]. El estudio formal de estos sistemas
empezó en la década de los treinta con trabajos de Kolmogorov, Petrovskii y Piscounov [18], y Fisher
[12]. Sin embargo, el potencial de estos sistemas lo encontró A. Turing en 1952, al notar que la reacción
quı́mica y la difusión molecular, procesos homogeneizadores por separados, pueden dar lugar a patrones
estacionarios cuando se les combina de forma especial [30]. Este último fenómeno se puede encontrar en
diversas áreas de la biologı́a y quı́mica, en particular, donde se involucre la dispersión e interacción de
individuos, células o especies quı́micas dentro de una región: crecimiento de poblaciones, la regeneración
de tejidos y la formación de tumores [13].
Encontrar una solución explı́cita para un sistema de reacción-difusión puede ser extremadamente com-
plicado o imposible (esto es lo que ocurre con mayor frecuencia), por lo que la solución es recomendable
(o debe) darse a través de métodos numéricos. Existen métodos ((clásicos)) para esto, como el de diferencias
finitas o elemento finito, sin embargo se optará por usar un método basado en el primero mencionado, pero
usando unas funciones especiales llamadas funciones de base radial (FBR), el cual ha sido desarrollado en

1
2 CAPÍTULO 1. INTRODUCCIÓN

los últimos años.

El uso de FBR para interpolar datos tiene sus inicios en el año 1971, donde el geofı́sico Rolland Hardy
utilizó una de estas funciones (multicuádrica) para interpolar datos geográficos obteniendo muy buenos
resultados [14]. Para interpolar datos que procedan de una función univariable, elegir entre usar FBR o poli-
nomios no es crucial en términos de eficacia o dificultad de plantear; sin embargo, para espacios euclidianos
de dimensión mayor o igual a dos, interpolar con FBR resultaba una buena técnica y más sencilla que la
polinómica.
A finales del siglo XX las FBR empezaron utilizarse para la resolución de ecuaciones diferenciales (ED).
El primer ejemplo de esto es el método de colocación simétrica propuesto por Kansa en 1990 [16, 17], el
cual se basaba en proponer como solución global de una ED un interpolador por FBR. Desde entonces ha
sido utilizado para resolver problemas de reacción-difusión, resultando una gran herramienta debido a su
precisión y a que no depende de mallas. Sin emabargo, no parece haber aplicaciones de este método al
modelo que nos compete.
El presente trabajo está organizado como sigue: en el capı́tulo 2 se presenta los sistemas de reacción-
difusión y el análisis de inestabilidad de Turing, ası́ como la una introducción a la interpolación con FBR que
constituyen el marco teórico necesario para el desarrollo de este proyecto de tesis; en el capı́tulo 3 se presenta
el modelo matemático que se propone para describir la interacción plancton-nutriente en un ecosistema
marino; en el capı́tulo 4 se construye un método numérico DF-FBR y este es empleado para realizar algunas
simulaciones numéricas del modelo plancton-nutriente con variación espacial; en los capı́tulos 5 y 6 se
presentan los resultados y conclusiones, respectivamente.
Capı́tulo 2

Marco teórico

En este capı́tulo se realiza una introducción a los conceptos base para el desarrollo de este proyecto que
son los sitemas de reacción-difusión y el análisis de inestabilidad de Turing de dichos sistemas.

2.1. Sistemas de reacción-difusión

Los sistemas de reacción-difusión son modelos matemáticos que han sido utilizados para representar
cierto tipo de fenómenos fı́sicos y quı́micos en donde dos o más ((sustancias))1 distribuidas en un espacio
dado cambian e interactúan bajo la influencia de dos procesos; a saber, la difusión (movimiento local) y la
reacción (crecimiento, interacciones, cambios de estado) [24].
Por un lado, los procesos de difusión se presentan en aquellos fenómenos en donde las partı́culas de
varias sustancias se mueven como grupo de acuerdo con la trayectoria irregular que tenga cada una de ellas.
Por ejemplo, la distribución de calor a través de un objeto homogéneo; el modelo matemático más sencillo
(y famoso) que representa el fenómeno anterior es la ecuación de calor

∂u
= D∆u, (2.1)
∂t

donde D > 0 y es llamado coeficiente de difusión, y ∆ es el operador laplaciano [28].

En general, para modelar un sistema de difusión de concentraciones de sustancias en n dimensiones,
necesitamos ecuaciones de la forma

∂u
= ∇ · (D(u, x)∇u), (2.2)
∂t
1
Se utilizará este término para referirse desde sustancias en sı́ hasta especı́es, objetos, etc. En cuanto se llegue a un caso particular
se especificará qué se está modelando exactamente.

3
4 CAPÍTULO 2. MARCO TEÓRICO

donde u = u(x, t), con x ∈ Rn y t ∈ R+ , es la concentración de la sustancia a través del tiempo. Notemos
que en (2.2) el término de difusión es una función que depende de la concentración y del espacio que ocupa;
sin embargo, solo nos concentraremos en términos de difusión constantes como en (2.1).
Por otro lado, estas sustancias pueden tener una reacción a lo largo del tiempo, ya sea de manera es-
pontánea o por la interacción con otras sustancias. Este tipo de fenómenos se presenta principalmente en la
dinámica de poblaciones de algunas especies y en reacciones quı́micas. Ejemplos de estos sistemas que son
utilizados para modelar los anteriores fenómenos son la ecuación de conservación (2.3), la cual se basa solo
en los nacimientos y muertes de cada cierto tiempo, y la ecuación de crecimiento logı́stico (2.4), que está
basado en la taza de nacimientos per cápita [21]:
dP
= bP − dP, (2.3)
dt  
dP P
= rP 1 − , (2.4)
dt K
donde b, d, r y K son constantes positivas y P (t) es el tamaño de la población en el tiempo t.
Notemos que podemos expresar el lado derecho de las ecuaciones anteriores como una función f (P (t)) =
f (P ). A este lo llamaremos término de reacción y marca la taza de crecimiento de la población estudiada.
Ası́, un sistema de reacción básico se puede expresar de la siguiente forma
dP
= f (P ). (2.5)
dt
Para modelar la interacción de un conjunto de n sustancias en un espacio acotado Ω ⊂ Rn , se utiliza el
sistema de reacción-difusión
∂u1
= Du1 ∆u1 + f1 (u1 , u2 , . . . , un ),
∂t
∂u2
= Du2 ∆u2 + f2 (u1 , u2 , . . . , un ),
∂t (2.6)
..
.
∂un
= Dun ∆un + fn (u1 , u2 , . . . , un ),
∂t
donde las ui : Ω × R+ −→ R son funciones que representan el ((movimiento)) de la i-ésima sustancia en la
región Ω a través del tiempo, las fi : Rn −→ R son los términos de reacción de cada sustancia en relación
al resto, y las Di son las respectivas constantes de difusión (positivas).
Para completar el sistema, se imponen las condiciones

ui (x, 0) = gi (x), x ∈ Ω, (2.7)

∂ui
= Gi (x, t), x ∈ ∂Ω, (2.8)
∂~n
2.1. SISTEMAS DE REACCIÓN-DIFUSIÓN 5

donde (2.7) es la condición inicial con la que se especifica la cantidad de la sustancia ui en tiempo t = 0, y
(2.8) es la condición de frontera de Neumann, que fijan la cantidad de flujo de ui en dirección normal a la
frontera de la región Ω.

2.1.1. Análisis de inestabilidad de Turing

Análisis de inestabilidad de Turing es el proceso por el cual se encuentran las condiciones para que
un punto crı́tico asintóticamente estable de un sistema de reacción se convierta en inestable al agregarle
difusión.
Se realizará el análisis de inestabilidad de Turing a un sistema de reacción-difusión con dos sustancias.
Sea
∂u
= Du ∆u + f (u, v),
∂t (2.9)
∂v
= Dv ∆v + g(u, v),
∂t
con las condiciones inciales

u(x, y, 0) = u0 (x, y) y v(x, y, 0) = v0 (x, y), (2.10)

y condiciones de frontera
∂u ∂v
=0 y = 0. (2.11)
∂~n ∂~n
Observación 1. Las funciones u y v son sustancias que tienden a ocupar toda una region Ω ⊂ R2 por la
difusión, f y g son términos de reacción (interacción) y Du , Dv son coeficientes de difusión positivos. Es
importante que las condiciones de frontera (2.11) sean homogéneas ya que impiden el flujo de las sustancias
ui hacia afuera de la región Ω.

Dicho análisis consiste de dos partes. Dada una solución (ū, v̄) de equilibrio del sistema (2.9) sin difu-
sión.

1. Establecer las condiciones para que (ū, v̄) sea localmente asintóticamente estable, las cuales son

Tr(J(ū, v̄)) < 0 y det(J(ū, v̄)) > 0, (2.12)

f
u fv

donde J = gu gv , es la matriz jacobiana.

2. Establecer las condiciones para que (ū, v̄) se vuelva inestable al incorporar la difusión. Para esto,
primero linealizamos el sistema:
" # " # " #" #
ut u Du 0 ∆u
= J(ū, v̄) + . (2.13)
vt v 0 Dv ∆v
6 CAPÍTULO 2. MARCO TEÓRICO

La solución del sistema linealizado sigue dependiendo de las variables independientes x, y y t por
contener el término de difusión. Ası́ que, supondremos que la soluciones se pueden escribir de la
siguiente forma

u(x, y, t) = ak (t)eik·x ,
(2.14)
v(x, y, t) = bk (t)eik·x ,

que es equivalente a la expresión matricial

" #
u
= ck (t)eik·x , (2.15)
v

con ck (t) = ck , k y x vectores. Sustituyendo lo anterior en el sistema linealizado, llegarı́amos a

∂
ck eik·x = J(ū, v̄)ck eik·x + D∆(ck eik·x ),
∂t
∂
eik·x eik·x eik·x
 2 (2.16)
ck =  J(ū, v̄)ck −  |k| Dck ,
  

∂t
d
ck = Bck ,
dt
donde D = D0u D0v y B = J(ū, v̄) − r2 D, con r = |k|.
 

La inestabilidad se consigue si Tr B > 0 o det B < 0, para cualquier k ∈ R2 . Sin embargo, la

primera desigualdad no se cumple: del análisis de estabilidad sabemos que Tr J < 0, ası́

Si k = (0, 0), entonces Tr B = Tr J < 0,

Si k 6= (0, 0), entonces Tr B = Tr J − r2 Tr D < 0 (r2 , Tr D > 0).

En consecuencia, la inestabilidad se dará cuando det B < 0, es decir,

r4 det D − (fu Dv + gv Du )r2 + det J < 0. (2.17)

Analizaremos cuándo (2.17) se cumple. Primero, notemos que la expresión de la izquierda puede verse
como un polinomio de segundo grado en la variable s = r2 ,

p(s) = (det D)s2 − (fu Dv + gv Du )s + det J. (2.18)

Dicho polinomio representa gráficamente una parábola cóncava hacia arriba, debido a que det D = Du Dv >
0. Ası́ que, si existe algún intervalo en donde p(s) < 0 deberá ser una vecindad del punto s que cumpla
p0 (s) = 2sDu Dv − (fu Dv + gv Du ) = 0, el cual es
 
1 fu gv
s= + . (2.19)
2 Du Dv
2.1. SISTEMAS DE REACCIÓN-DIFUSIÓN 7

Sustituyendo (2.19) en (2.18) y simplificando expresiones, llegamos a que la inestabilidad se podrı́a presen-
tar cuando

gv 2 4 det J
 
fu
+ > ,
Du Dv Du Dv (2.20)
(Dv fu + Du gv )2 − 4Du Dv det J > 0.

Observación 2. Dado que s = r2 , es necesario que se cumpla que el número

s  
∗ 1 fu gv
r = + , (2.21)
2 Du Dv

llamado número de onda crı́tico, exista. Es decir, que

Dv fu + Du gv > 0. (2.22)

En conclusión, de las desigualdades (2.12) y (2.20) obtenemos las condiciones que se deben de cumplir
para que exista la posibilidad de formar patrones. Estas son

fu + gv < 0,
fu gv + fv gv > 0,
(2.23)
Dv fu + Du gv > 0,
(Dv fu + Du gv )2 − 4Du Dv (fu gv + fv gv ) > 0.

Al conjunto de parámetros que cumplan con (2.23) es llamado espacio de Turing del sistema [32].

Observación 3. Es importante aclarar que, aun cumpliéndose (2.23), no necesariamente se deben formar
patrones en todo el espacio Ω, ya que dichas condiciones son locales. Con esto, también se debe exigir que
las condiciones iniciales (2.10) sean de la forma

u0 (x, y) = ū + ε1 ru (x, y),

(2.24)
v0 (x, y) = v̄ + ε2 rv (x, y),

donde ε1 y ε2 son números suficientemente pequeños.

Ejemplo 1 (Modelo de Schnakenberg). En 1979, Schnackenberg modela un fenómeno de reacción trimo-

lecular con las funciones

F (A, B) = k1 − k2 A + k3 A2 B,
(2.25)
G(A, B) = k4 − k3 A2 B,
8 CAPÍTULO 2. MARCO TEÓRICO

donde A es un quı́mico activador y B uno inhibidor. En el modelo, el término A2 B es la producción de A en

presencia de B en la primera ecuación, y el consumo de B en presencia de A en la segunda. El término −A
representa la degradación del mismo en la primera ecuación, produciendo un aumento de B en la segunda
[26, 19].
Al agregar difusión al sistema de reacción con las funciones anteriores y adimensionalizando, obtenemos
el modelo de Schnackenberg:

∂u
= ∆u + γ(a − u + u2 v),
∂t (2.26)
∂v
= d∆v + γ(b − u2 v),
∂t
donde u es la densidad del quı́mico activador y v es la densidad del quı́mico inhibidor, las constantes a, b > 0
Dv
representan la producción de u y v, respectivamente [13], d = D u
es la proporción de los coeficientes de
difusión y γ es un parámetro de escala del sistema que puede tener diferentes interpretaciones2 [22].
b
El único punto crı́tico de (2.26) es (ū, v̄) = (a + b, (a+b)2 ), ası́

" #
b−a
a+b (a + b)2
J(ū, v̄) = 2b
. (2.27)
− a+b −(a + b)2

Las condiciones para la estabilidad asintótica son

b − a − (a + b)3 < 0,
(2.28)
−(a + b)(b − a) + 2b(a + b)2 = (a + b)2 > 0.

Las condiciones para la inestabilidad son

(b − a)Dv − (a + b)3 Du > 0,

p (2.29)
(b − a)Dv − (a + b)3 Du  > 2(a + b)2 Du Dv .
 

Teniendo en cuenta que a, b, Du y Dv son reales positivos, las relaciones (2.28)-(2.29) pueden simplifi-
carse y llegar al siguiente espacio de Turing:

b − a − (a + b)3 < 0,
p (2.30)
(b − a)Dv − (a + b)3 Du > 2(a + b)2 Du Dv .

Observación 4. Es necesario que Du < Dv para que exista posibilidad de formación de patrones, es decir,
el inhibidor v se debe propagar más rápido que el activador u.
2
Como representar la fuerza relativa de los términos de reacción o una disminución de la relación del coeficiente de difusión d
(al aumentar γ).
2.1. SISTEMAS DE REACCIÓN-DIFUSIÓN 9

El modelo de Schnakenberg es de los sistemas de reacción-difusión más famosos y utilizados dentro de

la literatura cientı́fica ya que modela desde reacciones quı́micas hasta la formación y crecimiento de huesos
[13].
Este modelo tiene una particularidad interesante. La forma y tamaño del dominio afecta significativa-
mente la configuración de los patrones, independientemente de que se le impongan las mismas condiciones
iniciales y de frontera. Por ejemplo, diferentes autores han realizado simulaciones numéricas en donde el
dominio original crece con una tasa constante (ver 2.1), y también pruebas con figuras totalmente indepen-
dientes (como un cuadrado unitario en la Figura 2.2 o un rectángulo ((delgado)) en la Figura 2.3), y en ambas
se observa el fenómeno al principio mencionado.

Figura 2.1: Imagen extraı́da de [13]. Simulación numérica con elemento finito (parte espacial) y Euler ha-
cia atrás (temporal) utilizando los siguientes datos: el dominio inicial fue un cuadrado unitario, la tasa de
crecimiento uniforme es de 10−3 , a = 0.1, b = 0.9, γ = 10 y d = 10.
10 CAPÍTULO 2. MARCO TEÓRICO

(a) Primeros dos momentos (b) Últimos dos momentos

Figura 2.2: Simulación mediante el método de diferencias finitas (parte espacial) y Euler hacia atrás (tem-
poral) con los siguientes datos: γ = 100, a = 0.1305, b = 0.7695, D1 = 0.05 y D2 = 1.

2.2. Interpolación con funciones de Base Radial

Las funciones de base radial son un tipo de funciones que están intimamente relacionadas con la mag-
nitud de los vectores de su dominio. Estos habı́an sido utilizados prácticamente desde el inicio del cálculo,
obviamente sin conocer el potencial de estos (solo como una función más). No fue hasta 1971 que Ro-
lland Hardy en Multiquadric equations of topography and other irregular surfaces [14] utiliza una de estas
funciones (multicuádrica) para interpolar datos topográficos, obteniendo resultados bastante precisos.

Al momento no existe una sistematización exhaustiva sobre la teorı́a de funciones de base radial y sus
aplicaciones, sin embargo sı́ existen trabajos que lo han tratado de forma relativamente profunda, como
[4, 8], en los cuales nos vamos a basar.
2.2. INTERPOLACIÓN CON FUNCIONES DE BASE RADIAL 11

Figura 2.3: Se utilizan los mismos métodos y datos que la Figura 2.2, salvo que el dominio es un rectángulo
de 3 × 0.2.

2.2.1. Conceptos básicos

El problema de interpolar un conjunto de puntos es un clásico en el análisis numérico: dado un conjunto
de datos {(x1 , y1 ), . . . , (xn , yn )}, donde xi ∈ Rd y yi ∈ R con i = 1, . . . , n, encontrar una función
f : Rd −→ R tal que f (xi ) = yi . Este problema queda conceptualizado en la siguiente definición.

Definición 1. Supongamos que Ω ⊂ Rd contiene al menos N puntos y sea V ⊂ C(Ω) un espacio vectorial
N -dimensional. Entonces se dice que V es un espacio de Haar de dimensión N sobre Ω si para cualesquiera
puntos distintos x1 , . . . , xN ∈ Ω y cualesquiera u1 , . . . uN ∈ R existe una única función s ∈ V tal que
s(xj ) = uj , j = 1, . . . , N .

La interpolación polinomial es comúnmente empleada para aproximar funciones (véase [5, 25]) debido
a que es una herramienta sencilla. Esta técnica aproxima de forma ((eficaz)) funciones continuas (teorema
de aproximación de Weierstrass) y siempre existe un (único) polinomio de grado N − 1 que cumpla con la
tarea, es decir, V = {a0 + a1 x + · · · + an xN −1 : ai ∈ R, i = 0 . . . N } es un espacio de Haar de dimensión
N sobre cualquier subconjunto de N puntos Ω ⊂ R.
12 CAPÍTULO 2. MARCO TEÓRICO

Sin embargo, al tratar de interpolar puntos en RN , con N ≥ 2, la interpolación polinomial se vuelve

una tarea complicada y además, es imposible fijar un espacio de Haar N dimensional. Esto es debido al
siguiente resultado sobre la existencia de espacios de Haar.

Teorema 1 (Mairhuber-Curtis). Supongamos que d ≥ 2 y que Ω ⊂ Rd contiene un punto interior. Entonces

no existe ningún espacio de Haar sobre Ω de dimensión N ≥ 2.

En una dimensión, escoger entre interpolación polinomial o mediante FBR no resulta primordial, en
dimensiones mayores definitivamente es más sencillo y práctico utilizar la segunda herramienta. Ası́ que a
continuación empezaremos un estudio introductorio sobre las FBR.

Definición 2. Una función g : Rm −→ R es un radial si el valor de la función solo depende de la magnitud

del argumento, es decir, existe φ : [0, ∞) −→ R y || · || una norma tal que g(x) = φ(||x||) = φ(r). A la
función φ se le llama función de base radial (FBR).

(a) Radial g(x, y) = exp(−||x − y||2 ). (b) FBR φ(r) = exp(−r2 ).

Figura 2.4: Función radial y su respectiva FBR.

Ası́, si φ es una FBR y {(x1 , y1 ), . . . , (xn , yn )} un conjunto de datos, una función ξ es un interpolador
por funciones de base radial si es de la forma
n
X
ξ(x) = ci φ(||x − xi ||), x ∈ Rd , (2.31)
i=1

donde ci ∈ R y ξ(xj ) = yj , i. e.,

n
X
ci φ(||xj − xi ||) = yj , j = 1, . . . n. (2.32)
i=1
2.2. INTERPOLACIÓN CON FUNCIONES DE BASE RADIAL 13

Notemos que, para nuestro problema, las incógnitas son los escalares ci . Estos se pueden conocer desa-
rrollando (2.32) como el sistema de ecuaciones
    
φ(||x1 − x1 ||) · · · φ(||xn − x1 ||) c1 y1
 .
.. . .. .
..
 .   . 
  ..  =  ..  ,

     (2.33)
φ(||x1 − xn ||) · · · φ(||xn − xn ||) cn yn

que podemos simplificar por

ΦC = Y. (2.34)

Debido a las propiedades de la norma, Φ es una matriz cuadrada de n × n y simétrica. Además, para
que (2.33) tenga solución única, necesitamos que Φ sea no singular, lo que también nos darı́a la solución al
problema de interpolación.

2.2.2. Invertibilidad de matriz de la interpolación.

En la sección anterior vimos que el problema de interpolación se resume en resolver un sistema de ecua-
ciones lineales, y que este tendrı́a solución si la matriz es no singular. En general, las matrices Φ que surgen
de la interpolación con funciones de base radial no son invertibles, pero bajo ciertas condiciones adicionales,
se puede lograr que el problema de interpolación tenga solución única. No existe una sistematización sobre
la invertibilidad de todas las matrices de interpolación asociadas a una FBR, sin embargo esto se puede
abordar desde los conceptos de (((condicionalmente) definida positiva)).

Definición 3. Sea A una matriz real y simétrica de N × N . Se dice que A es semidefinida positiva si para
b ∈ RN se cumple que bT Ab ≤ 0, es decir,
N X
X N
bi bj Aij ≥ 0. (2.35)
i=1 j=1

Si la desigualdad es estricta para todo b 6= 0, se dice que A es definida positiva.

Definición 4. Sea φ : Rn → C una función continua. Se dice que φ es semidefinida positiva si para todo
N ∈ N y todos los conjuntos de puntos distintos por pares {x1 , . . . , xN } ⊂ Rn y b ∈ CN , se cumple que
N X
X N
bi bj φ(xi − xj ) ≥ 0. (2.36)
i=1 j=1

Si la desigualdad es estricta para todo b 6= 0, se dice que φ es definida positiva.

Enunciamos algunas proposiciones y ejemplos, de los cuales se puede revisar la demostración y/o desa-
rrollo en [7].
14 CAPÍTULO 2. MARCO TEÓRICO

Proposición 2. Toda combinación lineal de constantes no negativas con funciones semidefinidas positivas,
es una función semidefinida positiva. Si en la combinación lineal existe una función definida positiva y su
respectiva constante es positiva, entonces la combinación lineal también es definida positiva.

Proposición 3. El producto de funciones definidas positivas también es una función definida positiva.
T
Ejemplo 2. La función φ(x) = eix y , con y ∈ Rn fijo, es definida positiva ya que
 2
N X N XN 
T
X
bi bj φ(xi − xj ) =  bi ei(xi ) y  ≥ 0. (2.37)
 
 
i=1 j=1 i=1

Proposición 4. Si A es una matriz definida positiva, entonces A es invertible.

Corolario 5. Si φ es una función de base radial definida positiva (o combinación lineal de ellas), su matriz
de interpolación Φ es invertible.

Algunas FBR generan una matriz de interpolación Φ definida positiva. Ejemplos de estas son la gaus-
2 1
siana φ1 (r) = e−(cr) o la multicuádrica inversa φ2 (r) = (r2 + c2 )− 2 , c > 0. Existen FBR usadas en
multitud de aplicaciones las cuales no generan matrices definidas positivas. Para esto es necesario aumentar
el interpolador de FBD con una parte polinomial, lo cual nos lleva al concepto de matriz condicionalment
definida positiva.
Dado que no toda matriz A real y simétrica de N × N cumplirá la desigualdad (2.35) para todo vector b,
podemos relajar esta condición pidiendo que cumpla dicha desigualdad pero solo con los vectores b ∈ RN
que cumplan con que
N
X
bj = 0. (2.38)
i=1
Esta restricción es conocida como condición de ortogonalidad, y toda matriz A que cumpla con todo lo
anterior es llamada condicionalmente semidefinida positiva de orden 1. El aplicar la condición de ortogo-
nalidad no es arbitrario: si A es una matrı́z condicionalmente definida positiva y P = (1, . . . , 1)T ∈ RN ,
entonces el sistema " #" # " #
A P b Y
T
= , (2.39)
P 0 d 0
tiene solución única.
Notemos que si una matriz A condicionalmente definida positiva, cumple con la condición de ortogo-
nalidad para un número infinito de vectores: si b 6= 0 es un vector que cumple con esta condición, entonces
cb, con c ∈ R, también la cumple. Esto puede verse como que el conjunto formado por cada entrada del
vector b cumple con ser una combinación lineal no trivial que anula cualquier base del espacio vectorial
real de polinomios de grado 0 (constantes). Esta última idea es la que genera la definición completa de una
matriz/función condicionalmente definida positiva de orden m.
2.2. INTERPOLACIÓN CON FUNCIONES DE BASE RADIAL 15

Definición 5. Sea A real y simétrica de N × N . Se dice que A es condicionalmente semidefinida positiva

de orden m en Rn si se cumple que
N X
X N
bi bj Aij ≥ 0, (2.40)
i=1 j=1

para todos los conjuntos de puntos distintos por pares {x1 , . . . , xn } ⊂ Rn , siempre que
n
X
bk p(xk ) = 0, (2.41)
k=1

para cualquier polinomio p de valores complejos y grado menor o igual a m − 1. Si la desigualdad es estricta
para todo b 6= 0, se dice que A es condicionalmente definida positiva.

Definición 6. Sea φ : Rn → C una función continua. Se dice que φ es condicionalmente semidefinida

positiva de orden m en Rn si se cumple que
N X
X N
bi bj φ(xi − xj ) ≥ 0, (2.42)
i=1 j=1

para todos los conjuntos de puntos distintos por pares {x1 , . . . , xn } ⊂ Rn , siempre que
n
X
bk p(xk ) = 0, (2.43)
k=1

para cualquier polinomio p de valores complejos y grado menor o igual a m − 1. Si la desigualdad es estricta
para todo b 6= 0, se dice que φ es condicionalmente definida positiva.

Notemos que las condiciones de ortogonalidad se pueden verificar para un número finito de ecuaciones,
en particular, para el conjunto de polinomios que genera a Pm−1 (Rn ), es decir, se cumpla que
n
X
bk pl (xk ) = 0, 1 ≤ l ≤ Q, (2.44)
k=1

donde Q = dim(Pm−1 (Rn )) y gen({p1 , . . . , pQ }) = Pm−1 (Rn ).

La invertibilidad estará asociada a los anteriores conceptos de la siguiente manera.

Definición 7. Un conjunto X = {x1 , . . . , xN } ⊂ Rn se llama m-unisolvente, si el único polinomio de

grado a lo más m que interpola los datos cero sobre X es el polinomio cero.

Ejemplo 3. El conjunto {1, 2, 3} ⊂ R es 1-unisolvente y 2-unisolvente, sin embargo no es m-unisolvente

para m ≥ 3 ya que el polinomio f (x) = (x − 1)(x − 2)(x − 3) tiene estas raı́ces (ver Figura 2.5).
16 CAPÍTULO 2. MARCO TEÓRICO

Figura 2.5: El conjunto {1, 2, 3} no es m-unisolvente para m ≥ 3.

Proposición 6. Sea A una matriz simétrica real de N × N condicionalmente definida positiva de orden m
y {x1 , . . . , xN } forman un conjunto (m − 1)-unisolvente. Entonces el sistema
" #" # " #
A P b Y
T
= , (2.45)
P 0 d 0
tiene solución única.
Ası́, si tenemos una FBR φ condicionalmente definida positiva de orden m, podemos generar una matriz
de interpolación invertible siempre y cuando al interpolador se le agregue un polinomio de grado m y las
condicones de ortogonalidad. Lo cual queda materializado en los siguientes párrafos.
Sea Pdm−1 [x] el conjunto de todos los polinomios p : Rd −→ R de grado a lo más m − 1. Si {pj }m
j=1 es
b

una base para Pdm−1 [x], al interpolador ξ lo podemos completar de la siguiente manera.
n
X m
b
X
ξ(x) = ci φ(||x − xi ||) + λj pj (x), x ∈ Rd , ci , λj ∈ R; (2.46)
i=1 j=1

con la restricción
n
X
λi pj (xi ) = 0, j = 1, . . . , m.
b (2.47)
i=1
En este caso se tendrı́a que resolver el sistema de ecuaciones en donde ci y λj son las incógnitas. El
sistema matricial equivalente a (2.46) y (2.47) es
" #" # " #
Φ P C Y
= . (2.48)
PT 0 λ 0
2.2. INTERPOLACIÓN CON FUNCIONES DE BASE RADIAL 17

A continuación se anexa el cuadro 2.1 (extraı́do y traducido de [8]) en donde se muestran las FBR más
utilizadas. Todas funciones del cuadro son (condionalmente) positivas definidas. Cabe destacar que, para
efectos prácticos, la totalidad de estas dependen de un parámetro real positivo llamado parámetro de forma.
El estudio de cuál parámetro es el óptimo en cada FBR es todo un campo de investigación y en ocasiones se
requiere particularizar al problema tratado. Es puesto por el usuario y afecta la eficiencia del interpolador (y
en consecuencia de cualquier método que esté basado en FBR), sin embargo no se ahondará mucho en este.

Función de base radial φ(r) Parámetro Orden

Gaussiana e (cr)2 c>0 0
Spline poliharmónico r2k−1 k∈N m=k
r2k log(r) k∈N m=k+1
√
Multicuádrica r2 + c2 c>0 1
Multicuádrica inversa √ 1 c>0 0
r2 +c2
1
Cuádrica inversa r2 +c2
c>0 0

Cuadro 2.1: Funciones de base radial comunes.

Capı́tulo 3

Interacción plancton-nutriente

En este capı́tulo se presenta la interacción del fitoplancton, zooplancton y distintos tipos de nutrientes
dentro de de un ecosistema marino. Además, se mostrarán algunos modelos básicos (clásicos) que son
utilizados para representar distintos fenómenos biológicos asociados a la interacción entre especies y su
entorno.

3.1. Ecosistemas marinos1

En lagos y lagunas, ası́ como mares y océanos, podemos encontrar una zona en la cual la luz penetra de
tal forma que facilita que ciertos organismos realicen la fotosı́ntesis. A dicha zona se le llama zona eufótica y
su profundidad varı́a dependiendo de la intensidad de la luz, la turbidez del agua, sustancias externas, etc. La
zona en donde no llega a penetrar la luz (menos del 1 %), recibe el nombre de zona afótica y, generalmente,
es la zona más profunda de un sistema acuático (salvo que sea poco profundo como un manantial o charco).
Debido a las condiciones de ambas zonas, en la fótica se encuentra la mayor parte de la biota2 acuática,
aproximadamente un 90 % [31]. La organismos que conviven en la zona eufótica se dividen principalmente
en dos tipos:

El plancton, que son organismos flotantes cuyos movimientos dependen más o menos de las corrientes
(ver Figura 3.1a).

El necton, que son todos los organismos flotantes que pueden navegar a voluntad. Por ejemplo, los
peces, anfibios, grandes insectos nadadores, etc. (ver Figura 3.1b).

En este trabajo nos interesan los primeros. El plancton también puede clasificarse en dos grandes grupos:
1
La información mostrada en esta sección fue obtenida principalmente de Ecologı́a de Eugene Odum [23] y Essentials of
oceanography de Trujillo y Thurman [29].
2
Conjunto de organismos vivos.

18
3.1. ECOSISTEMAS MARINOS 19

(a) Plancton (b) Necton

Figura 3.1: Imágenes obtenidas de https://www.ecologiaverde.com.

Si el organismo planctónico puede llevar a cabo el proceso de fotosı́ntesis y ası́ autoalimentarse,

lo llamaremos fitoplancton. Estos llegan a medir aproximadamente 100 micras o menos, lo cual es
aproximadamente el grosor de un cabello (ver Figura 3.2a).

Si el organismo planctónico no puede producirse su propia comida y tiene que aprovecharse de los ali-
mentos producidos por otros, lleva el nombre de zooplancton. Estos llegan a medir aproximadamente
2 milimetros (ver Figura 3.2b).

Dado tamaño de el zooplancton, este se alimenta principalmente de fitoplancton, por lo que se puede
proponer un modelo de dinámica de poblaciones. Sin embargo, además de este factor, existen ciertas sus-
tancias, como el nitrógeno, fósforo, silicio y el azufre, llamadas nutrientes, que actúan como componentes
abióticos en los ecosistemas marinos [31]. Por lo que incluir nutrientes a dichos modelos lo mejorarı́a de
manera notable.

3.1.1. Cadena trófica

Daremos un pequeño resumen de cómo es que el fitoplacton interactúa con los diversos tipos de nutrien-
tes que se pueden encontrar en cualquier ecosistema marino. Esta información es recabada de [27].
Como vimos anterirormente, el fitoplancton es, prácticamente, una planta que realiza el proceso de
fotosı́ntesis. Ası́ como las plantas terrestres, estas necesitan ((fertilizantes)) para fabricar sus proteinas, ácidos
nucléicos y otras partes de la célula que necesitan para sobrevivir y reproducirse; los nutrientes juegan este
papel. Los principales nutrientes que la mayorı́a de especies de fitoplacton necesitan son el nitrógeno y el
fósforo (y algunos otros tipos en mucha menor cantidad como el sı́licio y algunos metales, a saber el hierro,
el zinc o el cobalto).
20 CAPÍTULO 3. INTERACCIÓN PLANCTON-NUTRIENTE

(a) Fitoplancton (b) Zooplancton

Figura 3.2: Representación gráfica de algunas especies de planctons. Gráficos obtenidos de Essentials of
oceanography de Trujillo y Thurman [29].

Los nutrientes proceden de forma ((natural)) de la meteorización de las rocas y el suelo, y de la conversión
del nitrógeno atmosférico en formas biológicamente útiles. De forma especı́fica, el fósforo también procede
del uso de detergentes y del tratamiento de aguas residuales; el nitrógeno de la intensa actividad agrı́cola,
del uso excesivo de fertilizantes y del arado de la tierra.
El fitoplancton crece muy rápido, pero también muere apresuradamente: aproximadamente tiene una
esperanza de vidad de un dı́a. Este es consumido (generalmente ya muerto) por bacterias y zooplancton,
que convierten la materia orgánica en CO2 y en los nutrientes básicos. El CO2 es expulsado hacia el aire
o es reutilizado junto con los otros nutrientes liberados para nuevos procesos de fotosı́ntesis. De ocurrir lo
anterior, no habrı́a cambios significativos en la concentración de CO2 en la atmósfera. A este proceso se le
conoce como remineralización (ver Figura 3.3). De lo contrario, si el fitoplancton se hunde lo suficiente, el
CO2 y los nutrientes son almacenados en el fondo del sistema acuático lo cual provoca una disminución de
CO2 en la atmósfera3 .

3
Esto puede llegar a ser una bomba biológica, ya que los nutrientes se vuelven sedimentos y el CO2 es liberado a la atmósfera
hasta que la circulación oceánica eleve el agua del fondo hasta la superficie, lo cual puede llegar a tardar hasta un milenio.
3.2. MODELO ESPACIAL PARA LA INTERACCIÓN PLANCTON-NUTRIENTE 21

Figura 3.3: Proceso de remineralización.

3.2. Modelo espacial para la interacción plancton-nutriente

Utilicemos P y N para representar la cantidad de fitoplancton y nutriente, respectivamente. En el capı́tu-

lo anterior vimos que un modelo de reacción-difusión (2 × 2) para la representar la interacción de estos dos
22 CAPÍTULO 3. INTERACCIÓN PLANCTON-NUTRIENTE

componentes debe tener la forma

∂N
= D1 ∆N + f (P, N ),
∂t (3.1)
∂P
= D2 ∆P + g(P, N ),
∂t
con sus respectivas condiciones iniciales y de frontera.
Notemos que el término de interés es el de reacción, a saber, las funciones f y g.
El modelo más sencillo y que sirve como base de cualquier sistema que represente dinámica de po-
blaciones es la ecuación de conservación. Si P representa la población de cualquier especie, entonces el
modelo antes mencionado tiene la forma
dP
= bP − dP, (3.2)
dt
donde b, d son constantes positivas que representan la cantidad de nacimientos y muertes, respectivamente,
cada cierto tiempo.
Es claro que el modelo anterior no es realista (en un tiempo relativamente extenso) ni representa con fi-
delidad cualquier fenómeno poblacional: la solución es P (t) = P0 e(b−d)t , con P0 = P (0) > 0 la población
inicial, lo cual significarı́a que el comportamiento serı́a exponencial, tanto en el crecimiento cuando b > d,
como en el decrecimiento cuando b < d (ver Figura 3.4). Es imposible que una población de cualquier
tipo crezca o decrezca de esa forma dado que siempre existen factores como depredadores, espacio, comida
limitada, etc. Por lo que es primordial agregar complejidad al término reactivo.
Un modelo que es central en cualquier descripción de parasitismo o depredación y expresa de forma
efectiva la influencia del comportamiento de los enemigos naturales como individuos sobre la dinámica de
poblaciones, son los de respuesta funcional [11].
En 1959, Holling en su artı́culo Some characteristics of simple types of predation and parasitism [15]
describe tres tipos de respuesta funcional, a saber [3, 11]:

Tipo I: el depredador destruye un porcentaje fijo de la densidad de la presa en función del incremento
de densidad poblacional de la misma, esto es, existe un aumento lineal de la tasa de ataque del de-
predador respecto a la densidad de la presa, hasta llegar a un punto a patir del cual la máxima tasa de
ataque permanece constante. La matanza y, por ende, la depredación cesa cuando el depredador llega
a un nivel de saciedad (ver Figura 3.5a).

Tipo II: el número de individuos de la presa atacada y destruida por el depredador se incrementa
pero con una tasa decreciente, en función del incremento de la densidad de la presa. Aparece otro
parámetro, llamado ((tiempo de manipulación)), que representa el tiempo dedicado a perseguir, domi-
nar, consumir y/o digerir las presas, y a prepararse para la siguiente búsqueda. Como el tiempo de
3.2. MODELO ESPACIAL PARA LA INTERACCIÓN PLANCTON-NUTRIENTE 23

Figura 3.4: Comportamiento de la solución de la ecuacion de conservación con P0 = 1. En negro cuando

b = d + 0.25 y en azul cuando b = d − 0.25.

búsqueda se reduce (debido al tiempo de manipulación), este tipo de respuesta resulta en un aumento
desacelerado a medida que aumentan las presas consumidad, hasta llegar a una ası́ntota en la cual se
expresa la máxima tasa de ataque. Ası́, los tiempos de manipulación largos conducen a bajas tasas de
ataque y viceversa (ver Figura 3.5b).

Tipo III: el depredador aprende a concentrarse sobre la presa cuando la población de la presa se
incrementa. El número de presas consumidas por unidad de tiempo se acelera con el aumento de la
densidad de presas, hasta que el tiempo de manipulación comienza a limitar su consumo. Por lo que,
a bajas densidades de presa, una mortalidad denso-dependiente directa. Una vez superado el punto de
inflexión, esta respuesta tiende a parecerse a la tipo II (ver Figura 3.5c).

La respuesta funcional que más se apega al fenómenos a estudiar, y que es utilizado más en la literatura
cientı́fica, es el tipo II. Esta viene representado por funciones de la forma
v
f (u, v) = −α1 u,
k+v
v (3.3)
g(u, v) = −α2 u,
k+v
donde u y v son las las funciones que denotan la cantidad de presa y parásito (o depredador), respectivamen-
te, y los parámetros α1 y α2 (α1 ≥ α2 ) marcan una tasa de absorción/consumo y crecimiento/conversión
para ambas poblaciones, y k es una constante.
24 CAPÍTULO 3. INTERACCIÓN PLANCTON-NUTRIENTE

(a) RF Holling Tipo I (b) RF Holling Tipo II (c) RF Holling Tipo III

Figura 3.5: Tres respuestas funcionales de Holling.

Para el fitoplancton P y los nutrientes N , el sistema de reacción-difusión tendrı́a la forma

∂N
= D1 ∆N + f ∗ (N ) − H1 (N, P ),
∂t (3.4)
∂P ∗
= D2 ∆P + g (P ) + H1 (N, P ),
∂t
donde las funciones f ∗ y g ∗ denotan un comportamiento análogo al de la ecuación de conservación vistas al
principio de la sección, y H1 , H2 son funciones de Holling tipo II.
Para proponer las funciones f ∗ y g ∗ del sistema (3.4), note que por un lado, dado que los nutrientes no se
reproducen pero sı́ se sedimentan, la función f ∗ debe ser de la forma f ∗ (N ) = I − bN , donde I es el ı́ndice
de entrada de nutrientes que fluyen en el agua y b es la razón de supresión natural del nutriente (en forma
de sedimentación); en cambio, el fitoplancton sı́ tiene un ı́ndice de muertes, pero la reproducción quedará
marcada principalmente por su relación con los nutrientes, por lo que g ∗ es de la forma g ∗ (P ) = −µP , con
µ la taza de muertes del fitoplancton.
Por lo tanto, el sistema de reacción-difusión propuesto en este trabajo es
∂N N
= D1 ∆N + I − bN − α P,
∂t k+N
(3.5)
∂P N
= D2 ∆P + β P − µP,
∂t k+N
con las condiciones inciales
N (x, y, 0) = (N , P ) + ε1 N0 (x, y)
(3.6)
P (x, y, 0) = (N , P ) + ε2 P0 (x, y),

donde (N , P ) es el punto de equilibrio y ε1 , ε2 son constantes positivas ((pequeñas)), y de frontera

∂N
=0
∂~n (3.7)
∂P
= 0.
∂~n
3.2. MODELO ESPACIAL PARA LA INTERACCIÓN PLANCTON-NUTRIENTE 25

Observación 5. En el modelo, los parámetros α y β (α ≥ β) son la tasa de absorción de nutrientes por

la población de fitoplancton y la tasa de conversión de nutrientes para el crecimiento de la población de
fitoplancton, respectivamente, y K es una constante de saturación media.

3.2.1. Análisis de inestabilidad de Turing.

Se consiguen dos puntos asintóticamente estables:
   
I µk β(µ(bk + I) − βI)
(N̄1 , P̄1 ) = ,0 y (N̄2 , P̄2 ) = , . (3.8)
b β−µ αµ(µ − β)

Debido a la complicación que podrı́a provocar utilizar las anteriores expresiones, el análisis de inestabi-
lidad se hará de forma general suponiendo que las soluciones de equilibrio son (N̄ , P̄ ).
Calculando el jacobiano,
" k
#
−b − α (N̄ +k) 2 P̄ −α N̄N̄+k
J(N̄ , P̄ ) = k
, (3.9)
β (N̄ +k)2 P̄ β N̄N̄+k − µ

de donde obtenemos que

β N̄ (N̄ + k) − αk P̄
Tr J(N̄ , P̄ ) = − (b + µ),
(N̄ + k)2
(3.10)
αµk P̄ − bβ N̄ (N̄ + k)
det J(N̄ , P̄ ) = + bµ.
(N̄ + k)2

Ası́, el espacio de Turing del sistema (3.5) son los parámetros que cumplen con las siguientes desigual-
dades.

β N̄ (N̄ + k) − αk P̄ < (N̄ + k)2 (b + µ),

bβ N̄ (N̄ + k) − αµk P̄ < (N̄ + k)2 bµ,
(3.11)
D1 β N̄ (N̄ + k) − D2 αk P̄ > (N̄ + k)2 (D2 b + D1 µ),
[D1 β N̄ (N̄ + k) − D2 αk P̄ − (N̄ + k)2 (D2 b + D1 µ)]2 > 4D1 D2 (N̄ + k)2 R,

donde R = βµ(N̄ + k)2 + αµk P̄ − bβ N̄ (N̄ + k).

Capı́tulo 4

Solución numérica

En este capı́tulo se va a dar un bosquejo del método DF-FBR para resolver sistemas de EDP con valores
en la frontera de la forma

Y (u) = f (x), x ∈ Ω,
(4.1)
B(u) = g(x), x ∈ ∂Ω,

donde Ω ⊂ R2 es el dominio, f, g : R2 −→ R y Y es un operador diferencial y B de frontera.

4.1. Introducción al método DF-FBR

El primero en usar FBR para dar solución a un sistema de ecuaciones diferenciales parciales fue Edward
Kansa para tratar de resolver problemas relacionados con dinámica de fluı́dos [16, 17]. Simplemente propone
un interpolador como solución a la ED y resuelve un sistema de ecuaciones.
El método propuesto por Kansa, llamado de ((colocación asimétrica)) (CA), es bastante simple. Sea φ
una FBR condicionalmente definida positiva y considerermos el problema con valores en la frontera (4.12).
Discreticemos la región Ω en N puntos, denotemos esto por Ω̂ = {x1 , . . . , xN }. El método consiste en
proponer la solución del sistema como un interpolante del tipo

n
X m
b
X
ξ(x) = ci φ(||x − xi ||) + λj pj (x), x ∈ Rd , ci , λj ∈ R, (4.2)
i=1 j=1

donde m
b es el orden de φ.

26
4.2. DF-FBR POR STENCILS 27

Ası́, para que ξ(x) sea solución del sistema debe de resolverse el siguiente sistema

n
X m
b
X
ci Y (φ(||x − xi ||)) + λj Y (pj (x)) = f (xi ), i ∈ Ω tal que xi ∈ int(Ω),
i=1 j=1
n
X m
b
X
ci B(φ(||x − xi ||)) + λj B(pj (x)) = g(xi ), i ∈ Ω tal que xi ∈ ∂Ω, (4.3)
i=1 j=1
N
X
λj pj (x) = 0, 1 ≤ j ≤ m.
b
j=1

Este es el método más simple para la solución de ED usando FBR; sin embargo tiene algunas desven-
tajas importantes como que la matriz resultante del sistema anterior es llena y mal condicionada, y que la
invertibilidad de esta sigue siendo un problema abierto [1]. Existen varios método más modernos y potentes,
como los métodos DF-FBR.

4.2. DF-FBR por stencils

Existen varios esquemas del método DF-FBR, pero en este trabajo utilizaremos uno basado en el trabajo
[33] junto con la implementación de nodos fantasma.
Para una región Ω se puede hacer una discretización Ω̂ = {x1 , . . . , x}, con M nodos interiores (N −
M son exteriores). Además, a esta discretización podemos agregar nodos fantasmas1 . Para cada vecindad
(discreta) de x ∈ Ω̂ definamos un stencil Jn (x) como un conjunto conformado por los ı́ndices de los
elementos de dicha vecindad, donde n es la cantidad de ı́ndices2 (ver Figura 4.1). Al igual que Ω̂, cualquier
subconjunto de este puede separarse en puntos interiores y exteriores; ası́, un stencil puede separarse en
ı́ndices de puntos interiores y de frontera, los cuales vamos a denotar por In (x) y Kn (x), respectivamente3 .
El conjunto In (x) tendrá m < n elementos y Kn (x) tendrá n − m.
Recordemos que la idea principal del método de diferencias finitas es la discretización de la derivada de
forma local, esto se obtiene a través de encotrar unos ((pesos)). Para trasladar esta idea con FBR, supondremos
que para cada x ∈ Ω̂ y su respectivo stencil Jn (x), podemos interploar (localmente) con la función
X
ξ(x) = aj φ(||x − xj ||) + bT (x − x) + c, (4.4)
j∈Jn (x)

1
Esto para garantizar la discretización de algunos operadores de frontera y marcar la dirección del vector normal.
2
Tı́picamente n < 10.
3
También el conjunto Kn (x) puede separarse conforme al tipo de condiciones (Dirichlet o Neumann) que se les aplique, sin
embargo en este trabajo no hará falta.
28 CAPÍTULO 4. SOLUCIÓN NUMÉRICA

Figura 4.1: Región cuadrangular en donde se muestra cómo estarı́a formado un stencil de 5 elementos con
centro en xk , tomando en cuenta la norma euclidiana.

1 P
donde φ es una es una FBR y x = n j∈Jn (x) xj . Notemos que la expresión anterior puede reescribirse
como
ξ(x) = Ψ (x)e, (4.5)
con e = (a; b; c)4 y Ψ (x) = φ(||x − x1 ||), . . . , φ(||x − xn ||), (x − x)T , 1 .

Para calcular e es necesario construir una matriz de interpolación y saber el valor del interpolador en
ciertos puntos especı́ficos, para esto utilizaremos las condiciones de interpolación en cada stencil. Para cada
x ∈ Ω̂ y su respectivo Jn (x), se tiene lo siguiente:

Para los puntos interiores tendrı́amos la condición ξ(xi ) = Ψ (x)e, es decir,

X
aj φ(||xi − xj ||) + bT (xi − x) + c = ξi , i ∈ In (x), (4.6)
j∈Jn (x)

4
El sı́mbolo ((;)) significa concatenación de vectores columna.
4.2. DF-FBR POR STENCILS 29

donde ξi = ξ(xi ).

Es apropiado exigir exactitud del interpolador para cualquier función lineal, por lo que también se
deben imponer las siguientes condiciones:
X X
aj = 0 y aj (xj − x) = 0. (4.7)
j∈Jn (x) j∈Jn (x)

Además, esto asegura una matriz de interpolación simétrica5 .


∂ξ 
Para los puntos en la frontera con condiciones de Neumann tendrı́amos que se debe cumplir ∂~
n x =
k
∂Ψ 

n x e, es decir,
∂~ k

X ∂
aj [φ(||x − x||)]xk + bT ~n = ξ k , k ∈ Kn (x), (4.8)
∂~n
j∈Jn (x)

∂ξ 
donde ξ k = ∂~
n x .
k

Observación 6. En este punto ya se ha supuesto que el operador de frontera es el de derivada en dirección

∂
normal B = ∂~ n . No se pierde ninguna generalización y estamos particularizando el esquema a nuestro
problema. También, ya trabajaremos con Y = ∆.
0 = [ξ ; 0; ξ ] ∈ Rn+3 , el sistema de interpolación local está
Definiendo ξloc = [ξi ], ξ loc = [ξ k ] y ξloc loc loc
dado por
0
Ge = ξloc , (4.9)

donde G es una matriz de tamaño (n + 3) × (n + 3), con las primeras m + 3 filas representadas por las
ecuaciones (4.6) y (4.7), y las últimas n − m por la ecuación (4.8).

4.2.1. Solución de la parte espacial

Ahora, si tenemos la siguiente EDP
∆(ξ) = f, (4.10)

donde ∆ es el operador laplaciano y f una función conocida, para los M nodos interiores tendrı́amos que la
ecuación (4.10) tendrı́a la forma

∆(Ψ (x))|xi e = f (xi ), i = 1, . . . , M. (4.11)

5
También se tienen beneficios adicionales dependiendo de la FBR. Por ejemplo, como la matriz de interpolación de la FBR
multicuádrica siempre es invertible [9], se seguirá preservando la invertibilidad.
30 CAPÍTULO 4. SOLUCIÓN NUMÉRICA

De la ecuación (4.9) se llega a que e = G−1 ξloc

0 , por lo que (4.11) es equivalente a

∆(Ψ (x))|xi G−1 ξloc

0
= f (xi ), i = 1, . . . , M. (4.12)

La ecuación (4.12) representa la forma final del método DF-FBR para los M nodos internos.

4.2.2. Solución de la parte temporal

Para la parte temporal, utilizaremos el siguiente esquema: tomaremos a ≤ t ≤ b y tamaños de paso h;
si un y vn son las soluciones evaluadas en tn = t0 + nh, entonces

un+1 = un + hD1 ∆un+1 + hf (un , vn ),

(4.13)
vn+1 = vn + hD2 ∆vn+1 + hg(un , vn ).
Del cual resulta el sistema
(I − hD1 ∆)un+1 = un + hf (un , vn ),
(4.14)
(I − hD2 ∆)vn+1 = vn + hg(un , vn ),
donde I es el operador identidad.
Sin embargo, los operadores de la parte izquierda de (4.14) se discretizarán de manera especial utilizando
el método DF-FBR: para cada paso de tiempo, tendrı́amos el sistema

B1 un+1 = un + hf (un , vn ),
(4.15)
B2 vn+1 = vn + hg(un , vn ),
donde
Bi = E − Di hAI + ABC , i = 1, 2. (4.16)
Ası́ mismo, Bi = (I − hDi ∆), E = IInt + IBou y AI = AInt + ABou . Las configuraciones de las anteriores
matrices son las siguientes.

IInt es una matriz diagonal con 1 en la intersección de la diagonal con la fila respectiva al ı́ndice de
los nodos interiores y cero en el resto de la matriz.

IBou es una matriz diagonal con 1 en la intersección de la diagonal con la fila respectiva al ı́ndice de
los nodos fantasma y cero en el resto de la matriz.

AInt es la matriz extendida correspondiente a la discretización del laplaciano en los nodos interiores.

ABou es la matriz extendida correspondiente a la discretización del laplaciano en los nodos fantasma.

ABC es la matriz extendida correspondiente a la discretización de la derivada normal en los nodos

frontera.
Capı́tulo 5

Resultados

En este apartado se presentan algunas simulaciones numéricas del modelo para interacción plancton-
nutriente introducido en el capı́tulo 3. Para esto, se obtendrán soluciones numéricas mediante el método
DF-FBR introducido en el capı́tulo 4 en torno a soluciones de equilibrio positivas, es decir, para las que
N̄1 > 0, P̄1 > 0. Se reportará tanto la solución numérica como un gráfico que muestre el comportamiento
de la cantidad de los factores del sistema.

Se usarán dos conjuntos de parámetros basados en los trabajos de Chatterjee y Samares, y Upadhyay
et al. [6, 31].

1. I = 0.9, b = 0.3, α = 0.9, β = 0.8, µ = 0.5, k = 0.6, D1 = 0.5, D2 = 0.009.

2. I = 0.2, b = 0.8, α = 1.5, β = 0.5, µ = 0.08, k = 0.6, D1 = 0.5, D2 = 0.009.

Se destaca que estos parámetros cumplen con las desigualdades del espacio de Turing para nuestro
modelo. Los puntos crı́ticos son (N̄1 , P̄1 ) ≈ (0.138461, 5.087179) y (N̄1 , P̄1 ) ≈ (0.114285, 0.452381).

Se implementarán FBR multicuádricas con parámetro de forma 5 y stencils de tamaño 9.

Para el primer conjunto de parámetros se tomará el intervalo de tiempo [0, 400] con una partición de
en 10000 puntos, ası́ h1 ≈ 0.0.040004. Para el segundo [0, 1000] con 10000, por lo que h2 ≈ 0.10001.

Haremos la prueba del método tomando como dominio Ω = [100, 100]2 , al cual vamos discretizar con
2120 puntos utilizando la función ((min energy nodes)) del paquete rbf en python (ver Figura 5.1).

Observación 7. Se utilizará factorización LU para las matrices B1 y B2 mostradas en el esquema de la

solución numérica vista en la subsección 4.2.2.

31
32 CAPÍTULO 5. RESULTADOS

Figura 5.1: Discretización del dominio.

Observación 8. Se utilizará regla trapezoidal en dos dimensiones para gráficar la densidad del fitoplancton
y nutriente a través del tiempo.

Las condiciones iniciales están dadas por el punto de estabilidad más una ((pequeña)) acumulación central
modelada por la función

f (x, y) = 0.05 exp(−0.0025((x − 50)2 + (y − 50)2 )). (5.1)

5.1. Simulaciones numéricas

A continuación se anexan las dos simulaciones númericas.

5.1.1. Primera simulación

Con el primer conjunto de parámetros tenemos lo siguiente.

En la Figura 5.2 se representa la solución numérica del sistema a través de gráficas de contorno.
A la izquierda encontramos la acumulación de nutrientes, y a la derecha la de fitoplacton. De toda la
partición del tiempo, se toman 4 momentos que muestren un comportamiento interesante. En el primer
5.1. SIMULACIONES NUMÉRICAS 33

par de imágenes (t = 0) notamos los cúmulos de nutrientes y fitoplancton en el centro del dominio,
descrita por la condición inicial antes mencionada. En el segundo par de imágenes (t = 160), la
cantidad de nutriente y fitoplancton empieza a converger a valores cercanos a los de estabilidad, y se
nota una acumulación de fitoplancton y falta de nutrientes en la zona central del dominio. A partir
del tercer par de imágenes (t = 172) ya existe una ((convivencia heterogénea)) entre fitoplancton y
nutriente, aunque se sigue notando el comportamiento análogo al anterior en el centro. En las últimas
imágenes (t = 400) se tiene completamente la convivencia heterogénea. Cabe destacar que del tercer
al cuarto momento no existe un cambio significativo entre cantidad de los nutrientes y el fitoplancton.

En la Figura 5.3 se muestra la densidad del nutriente a través del tiempo; en la Figura 5.4 del fitoplac-
ton. Vemos que ya en t = 100 (incluso antes) la densidad de ambos componentes no tiene cambios
significativos. Esto indica que, aunque la convivencia heterogénea se da aproximadamente en t = 172,
desde t = 100 no existe cambio significativo de la cantidad total del nutriente y el fitoplancton.

5.1.2. Segunda simulación

Para el segundo conjunto de parámetros tenemos lo siguiente.

En la Figura 5.5 se representa la solución numérica del sistema a través de gráficas de contorno.
Notamos un comportamiento similar al de la primer simulación (Figura 5.2), salvo que la convivencia
heterogénea se presenta más tarde, a partir de t = 650 (aproximadamente).

En la Figura 5.6 se muestra la densidad del nutriente a través del tiempo; en la Figura 5.7 del fitoplac-
ton. También existe un comportamiento similar al caso anterior (Figura 5.6 y Figura 5.7), salvo que la
ausencia de cambios significativos se empieza a notar en t = 200.
34 CAPÍTULO 5. RESULTADOS

Figura 5.2: Comportamiento del nutriente (izquierda) y el fitoplancton (derecha) en t = 0, 160, 172, 400.
5.1. SIMULACIONES NUMÉRICAS 35

Figura 5.3: Densidad del nutriente a través del tiempo para el primer conjunto de datos parámetros.

Figura 5.4: Densidad del fitoplancton a través del tiempo para el primer conjunto de parámetros.
36 CAPÍTULO 5. RESULTADOS

Figura 5.5: Comportamiento del nutriente (izquierda) y el fitoplancton (derecha) en t = 0, 600, 650, 1000.
5.1. SIMULACIONES NUMÉRICAS 37

Figura 5.6: Densidad del nutriente a través del tiempo para el segundo conjunto de parámetros.

Figura 5.7: Densidad del fitoplancton a través del tiempo para el segundo conjunto de parámetros.
Capı́tulo 6

Conclusiones

En este trabajo se ha logrado lo siguiente.

Se ha dado una introducción a temas varios tales como sistemas de reacción-difusión y la inestabi-
lidad de Turing, e interpoladores por funciones de base radial y su uso combinado con métodos de
diferencias finitas.

Se ha propuesto un modelo espacial para describir la interacción entre el placton y nutriente, y me-
diante un análisis de inestabilidad se determinó su espacio de Turing.

Se utilizó el método DF-FBR para dar solución numérica al sistema. Se obtuvieron un par de simula-
ciones con valores de parámetros reportados en la literatura existente sobre este tipo de problema.

Se anexaron gráficas de la solución y las densidades de fitoplancton y nutriente.

Las evidencias dadas en este trabajo coinciden con las esperadas teóricamente y, además son análogas
a las presentadas en artı́culos en donde se utilizaban otros métodos numéricos, especı́ficamente aquellos de
donde se obtuvieron las constantes (ver [6, 31]), incluso sin tomar en cuenta un factor tan importante como la
influencia del zooplancton. Ası́, podemos aportar evidencia sobre la utilidad del modelo plancton-nutriente
propuesto y la efectividad del método numérico y sus ventajas (libre de malla).
En cuanto a trabajo futuro podemos proponer las siguientes lı́neas de investigación.

Agregar más aspectos o factores al sistema, tales como la influencia del zooplacton y/o necton, la
temperatura, la contaminación, etc.

Uso del método DF-FBR en este u otros sistemas variando y estudiando el comportamiento de la
solución al cambiar la FBR, el parámetro de forma y el tamaño de los stencil.

Probar este y demás modelos con datos obtenidos de mediciones experimentales.

38
Apéndice A

Método de diferencias finitas

En este apéndice se dará un pequeño repaso a los métodos de diferencias finitas (DF). No existe un
modelo general para un sistema de ecuaciones tanto ordinarias o parciales ya que el método se construye
a partir de aproximaciones discretas de cada derivada y estas son distintas por cada orden. Por lo que aquı́
solo ejemplificaremos su uso con la ecuación de Poisson debido a que utiliza el mismo operador diferencial
que el modelo de esta tesis, el laplaciano. Si desea, puede revisar literatura básica sobre el tema, tales como
en [5] y [25].
En esencia, el método de diferencias finitas se basa en generar unos ((pesos)) que conformarán una matriz
que aproxime los operadores de derivada y frontera en puntos de una malla que después se invertirá. Es decir,
si ui es la aproximación de la solución de una ecuación diferencial en un punto xi , obtendrı́amos al final un
sistema del tipo P U = fi , donde P es la matriz de pesos, U es un vector de aproximaciones y fi son los
valores que toma la parte (no) homogénea de la ecuación diferencial en los puntos xi .
Vamos a esquematizar el método DF para un PVF con la ecuación de Poisson en una variable sujeta
condiciones de Dirichlet. Ası́ resultará más sencillo extender el método a más variables.

−y 00 (x) = f (x), y(0) = 0, y(0) = 0, (A.1)

El método consiste en seleccionar un N ∈ N para subdividir el intervalo [0, 1] en N + 1 intervalos equidis-

tantes, generando una sucesión de puntos xi = ih, i = 0, 1, ..., N, N + 1, con h = 1/(N + 1). Ası́, x0 = 0,
xN +1 = 1 y los puntos medios de esta subdivisión tendrán la forma xi = ih, i = 1, 2, ..., N .
La primera (si hubiera) y la segunda derivada estarı́an reemplazadas por las fórmulas del punto medio de
tres puntos, para la primera, y la del punto medio, para la segunda, i.e., si y(x) es la solución de una EDO,

39
40 APÉNDICE A. MÉTODO DE DIFERENCIAS FINITAS

entonces
1 h2
y 0 (xi ) = [y(xi + h) − y(xi − h)] − y (3) (ξ1 )
2h 6
1 h2 (3)
= [y(xi+1 ) − y(xi−1 )] − y (ξ1 )
2h 6 (A.2)
1 h2
y 00 (xi ) = 2 [y(xi − h) − 2y(xi ) + y(xi + h)] − y (4) (ξ2 )
h 12
1 h2 (4)
= 2 [y(xi−1 ) − 2y(xi ) + y(xi+1 )] − y (ξ2 ),
h 12
para algunos ξ1 , ξ2 ∈ [xi−1 , xi+1 ].

Notemos que si y (3) y y (4) son continuas en [xi−1 , xi+1 ], las aproximaciones de (A.2) son O(h2 ). A
estas aproximaciones las llamaremos fórmulas de diferencias centradas.

Sustituyendo (A.2) en la EDO (A.1), llegamos a la ecuación

1 h2 (4)
−
[y(x i−1 ) − 2y(x i ) + y(x i+1 )] = f (x i ) − y (ξ2 ). (A.3)
h2 12
Incluyendo las condiciones de frontera, haciendo wi una aproximación de y(xi ) en (A.3) y omitiendo el
término de error, llegamos al sistema de ecuaciones

w0 = 0 wN +1 = 0

y
1
− [wi−1 − 2wi + wi+1 ] = f (xi ), (A.4)
h2
para cada i = 1, . . . , N .

Si reescribimos (A.4) de la forma

−wi−1 + 2wi − wi+1 = h2 f (xi ),

llegamos al sistema de ecuaciones lineales N × N

Aw = b,

donde A es la siguiente matriz tridiagonal

 
−12 0
 −1 . . . . . .
 

A= , (A.5)
 
.. ..

 . . −1 

0 −1 2
 41

además,    
w1 h2 f (x1 )
h2 f (x2 )
   
 w2  
.. ..
   
w= b=
   


 .  . 

w h2 f (x
N −1 )
 
 N −1  
wN h2 f (xN )
Para extender el método a dos dimensiones, tomemeos la ecuación de Poisson con condiciones de fron-
tera de Dirichlet definida sobre el rectángulo Ω = {(x, y) ∈ R2 |a ≤ x ≤ b, c ≤ y ≤ d; a, b, c, d ∈ R},

∆u = f (x, y), ∀(x, y) ∈ int(Ω),

(A.6)
u = g(x, y), ∀(x, y) ∈ ∂Ω,

empezaremos tomando n, m ∈ Z+ para definir los tamaños de paso h = (b − a)/n y k = (d − c)/m y

luego hacer una división en partes iguales de los intervalos [a, b] (eje x) y [c, d] (eje y), respectivamente. Ası́,
tendrı́amos la sucesión de números xi = a + ih, con i = 0, 1, ..., n y yj = c + jk, con j = 0, 1, ..., m.
Para cada punto (xi , yj ) en el interior de la malla (i = 1, ..., n − 1 y j = 1, ..., m − 1), podemos usar
la serie de Taylor en la variable x y y alrededor de xi y yi , respectivamente, para generar las fórmulas de
diferencias centradas.
∂2u u(xi+1 , yj ) − 2u(xi , yj ) + u(xi−1 , yj ) h2 ∂ 4 u
(x i , yj ) = − (ξi , yj )
∂x2 h2 12 ∂x4
(A.7)
∂2u u(xi , yj+1 ) − 2u(xi , yj ) + u(xi , yj−1 ) k 2 ∂ 4 u
(x i , yj ) = − (xi , ηj ),
∂y 2 k2 12 ∂y 4
donde ξi ∈ (xi−1 , xi+1 ) y ηj ∈ (yj−1 , yj+1 ).

Las condiciones de frontera estarı́an dadas por las relaciones

u(x0 , yj ) = g(x0 , yj ) y u(xn , yj ) = g(xn , yj );

(A.8)
u(xi , y0 ) = g(xi , y0 ) y u(xi , ym ) = g(xi , ym ),

para cada j = 0, . . . , m e i = 1, . . . , n − 1.

Definiendo wi,j como una aproximación de u(xi , yj ), omitiendo los términos de error de (A.7) y con-
siderando las condiciones (A.8), llegamos al esquema de método de diferencias finitas para la ecuación de
Poisson (A.6) es
"  #
h 2
 2
h
2 + 1 wi,j − (wi+1,j − wi−1,j ) − (wi,j+1 + wi,j−1 ) = −h2 f (xi , yi ), (A.9)
k k
42 APÉNDICE A. MÉTODO DE DIFERENCIAS FINITAS

para cada i = 1, 2, . . . , n − 1 y j = 1, 2, . . . , m − 1, y

w0,j = g(x0 , yj ) y wn,j = g(xn , yj );

(A.10)
wi,0 = g(xi , y0 ) y wi,m = g(xi , ym ),

para cada j = 0, . . . , m e i = 1, . . . , n − 1.
Apéndice B

Código en Python

A continuación anexamos el código (python 3.8) utilizado en este trabajo. No se agrega la importación
de los paquetes más conocidos.

from rbf.sputils import expand_rows

from rbf.pde.fd import weight_matrix
from rbf.pde.nodes import min_energy_nodes

#Definiendo el dominio y el tipo de condiciones de frontera.

Lx = 100
Ly = 100
vert = np.array([[0.0, 0.0], [Lx, 0.0], [Lx, Ly], [0.0, Ly]])
smp = np.array([[0, 1], [1, 2], [2, 3], [3, 0]])
boundary_groups = {’free’: [0, 1, 2, 3]}

spacing = 2.0 #Separación aproximada entre nodos.

n = 9 #Tamaño de stencils

eps=5 #Parámetro de forma.

phi = ’mq’ #FBR

order = 1 #Orden del polinomio adicional.

43
44 APÉNDICE B. CÓDIGO EN PYTHON

#Función generadora de puntos.

def rho(x):
return 0.1 + np.maximum(abs(x[:, 0] - Lx/2), abs(x[:, 1] - Ly/2))

#Generando los nodos

nodes, groups, normals = min_energy_nodes(
spacing,
(vert, smp),
rho=rho,
boundary_groups=boundary_groups,
boundary_groups_with_ghosts=[’free’],
include_vertices=True)

#Generando las matrices de peso con las FBR

A_interior = weight_matrix(
x=nodes[groups[’interior’]],
p=nodes,
n=n,
diffs=[[2, 0], [0, 2]],
phi=phi,
order=order,
eps=eps)

Eye_interior = weight_matrix(
x=nodes[groups[’interior’]],
p=nodes,
n=n,
diffs=[[0, 0]],
phi=phi,
order=order,
eps=eps)

A_boundary = weight_matrix(
x=nodes[groups[’boundary:free’]],
 45

p=nodes,
n=n,
diffs=[[2, 0], [0, 2]],
phi=phi,
order=order,
eps=eps)

Eye_boundary = weight_matrix(
x=nodes[groups[’boundary:free’]],
p=nodes,
n=n,
diffs=[[0, 0]],
phi=phi,
order=order,
eps=eps)

A_free = weight_matrix(
x=nodes[groups[’boundary:free’]],
p=nodes,
n=n,
diffs=[[1, 0], [0, 1]],
coeffs=[normals[groups[’boundary:free’], 0],
normals[groups[’boundary:free’], 1]],
phi=phi,
order=order,
eps=eps)

N = nodes.shape[0] #Cantidad de nodos.

A_int_expanded = expand_rows(A_interior, groups[’interior’], N)

A_int_expanded += expand_rows(A_boundary, groups[’ghosts:free’], N)

Eye_int_expanded = expand_rows(Eye_interior, groups[’interior’], N)

Eye_int_expanded += expand_rows(Eye_boundary, groups[’ghosts:free’], N)
46 APÉNDICE B. CÓDIGO EN PYTHON

A_BC_expanded = expand_rows(A_free, groups[’boundary:free’], N)

#Coeficientes para el modelo

D1 = 0.5
D2 = 0.009
I = 0.9
b = 0.3
alpha = 0.9
beta = 0.8
mu = 0.15
k = 0.6

#Funciones de reacción.
def f(N,P):
return I-b*N-alpha*(N/(k+N))*P

def g(N,P):
return beta*(N/(k+N))*P-mu*P

#Puntos de estabilidad.
N_bar = (mu*k)/(beta - mu)
P_bar = beta*(mu*(b*k+I)-beta*I)/(alpha*mu*(mu - beta))

#Condiciones iniciales.
def u_0(x):
return N_bar + 0.05*np.exp(-0.0025*((x[:,0]-Lx/2)**2+(x[:,1]-Ly/2)**2))

def v_0(x):
return P_bar + 0.05*np.exp(-0.0025*((x[:,0]-Lx/2)**2+(x[:,1]-Ly/2)**2))

# Método DF-FBR
T = 500 # Tiempo máximo.
Nt = 2000 # Cantidad de puntos en el tiempo.
t = np.linspace(0,T,Nt)
dt = t[1]-t[0] #Tamaño de paso.
 47

u = np.zeros((N,Nt))
v = np.zeros((N,Nt))

#Matrices de aproximación al operador (I-D*dt*Laplaciano)

B_u = Eye_int_expanded - D1*dt*A_int_expanded + A_BC_expanded
B_v = Eye_int_expanded - D2*dt*A_int_expanded + A_BC_expanded

#Factorización LU.
lu=splu(csc_matrix(B_u))
lv=splu(csc_matrix(B_v))

#Solución en la condición inicial.

u[groups[’interior’],0] = u_0(nodes[groups[’interior’]])
u[groups[’ghosts:free’],0] = u_0(nodes[groups[’boundary:free’]])
v[groups[’interior’],0] = v_0(nodes[groups[’interior’]])
v[groups[’ghosts:free’],0] = v_0(nodes[groups[’boundary:free’]])

#Iteración en el tiempo.
for n in range(Nt-1):
rhsu = u[:, n] + dt * f(u[:, n], v[:, n])
rhsu[groups[’boundary:free’]] = 0.0
rhsv = v[:, n] + dt * g(u[:, n], v[:, n])
rhsv[groups[’boundary:free’]] = 0.0
u[:, n + 1] = lu.solve(rhsu)
v[:, n + 1] = lv.solve(rhsv)
Referencias

[1] Pablo Alexei, ¡No Más Mallas! Solución Numérica de EDP’s con funciones de base radial, Congreso,
Octubre 2015, 6to Aquelarre Matemático en la Facultad de Ciencias de la UNAM.

[2] Ovide Arino, Khalid Boushaba y Ahmed Boussouar, A mathematical model of the dynamics of the
phytoplankton-nutrient system, Nonlinear analysis: real world applications 1 (2000), no 1, 69–87.

[3] M. H. Badii, J. Landeros, H. Rodrı́guez, E. Cerna, J. Valenzuela y Y. Ochoa, Algunos Aspectos de

Depredación, Daena: International Journal of Good Conscience 8 (2013), no 1, 148–158.

[4] Martin D Buhmann, Radial basis functions: theory and implementations, vol. 12, Cambridge university
press, 2003.

[5] Richard L. Burden, Douglas J. Faires y Annette M. Burden, Análisis numérico, décima ed., CENGAGE
Learning, México, 2017.

[6] Anal Chatterjee y Samares Pal, Plankton nutrient interaction model with effect of toxin in presence
of modified traditional holling type II functional response, Systems Science & Control Engineering 4
(2016), no 1, 20–30.

[7] Stefano De Marchi y E Perracchiono, Lectures on radial basis functions, Preprint (2018).

[8] Wilna Du Toit, Radial Basis Function Interpolation, Tesis de Licenciatura, University of Stellenbosch,
Sudáfrica, Marzo 2008.

[9] Gregory E Fasshauer, Meshfree approximation methods with MATLAB, vol. 6, World Scientific, 2007.

[10] Zhilan Feng, Rongsong Liu y Donald L DeAngelis, Plant–herbivore interactions mediated by plant
toxicity, Theoretical Population Biology 73 (2008), no 3, 449–459.

[11] Valeria Fernández-Arhex y Juan C Corley, La respuesta funcional: una revisión y guı́a experimental,
Ecologı́a austral 14 (2004), no 1, 83–93.

48
REFERENCIAS 49

[12] R. Aylmer Fisher, The wave of advance of advantageous genes, Annals of eugenics 7 (1937), no 4,
355–369.

[13] Libardo A. González, Juan C. Vanegas y Diego A. Garzón, Formación de patrones en sistemas de
reacción-difusión en dominios crecientes, Revista Internacional de Métodos Numéricos 25 (2009),
no 2, 145–161.

[14] Rolland L Hardy, Multiquadric equations of topography and other irregular surfaces, Journal of
geophysical research 76 (1971), no 8, 1905–1915.

[15] Crawford S Holling, Some characteristics of simple types of predation and parasitism, The Canadian
Entomologist 91 (1959), no 7, 385–398.

[16] Edward J Kansa, Multiquadrics—A scattered data approximation scheme with applications to compu-
tational fluid-dynamics—I surface approximations and partial derivative estimates, Computers & Mat-
hematics with applications 19 (1990), no 8-9, 127–145.

[17] , Multiquadrics—A scattered data approximation scheme with applications to computational

fluid-dynamics—II solutions to parabolic, hyperbolic and elliptic partial differential equations, Com-
puters & mathematics with applications 19 (1990), no 8-9, 147–161.

[18] A. Kolmogorov, I. Petrovskii y N. Piscounov, A Study of the Diffusion Equation with Increase in the
Amount of Substance, and its Application to a Biological Problem, Selected Works of A. N. Kolmo-
gorov. Volume I: Mathematics and Mechanics (V. M. Tikhomirov, ed.), Kluwer Academic Publishers,
Netherlands, 1991, págs. 242–270, ISBN 90-277-2796-1.

[19] Aldo Durán Ledesma, Patrones de turing en sistemas biológicos, Tesis de Doctorado, Universidad
Autónoma Metropolitana - Iztapalapa, 2012.

[20] Anotida Madzvamuse, Time-stepping schemes for moving grid finite elements applied to reaction–
diffusion systems on fixed and growing domains, Journal of computational physics 214 (2006), no 1,
239–263.

[21] J. D. Murray, Mathematical biology: I. An introduction, tercera ed., Interdisciplinary Applied Mathe-
matics, vol. 17, Springer, Estados Unidos de America, 2002, ISBN 0-387-95223-3.

[22] , Mathematical biology II: Spatial models and biomedical applications, tercera ed., In-
terdisciplinary Applied Mathematics, vol. 18, Springer, Estados Unidos de America, 2003,
ISBN 0-387-95228-4.

[23] E. Odum, Ecologı́a, tercera ed., Nueva editorial interamericana, Distrito Federal, 1971.
50 REFERENCIAS

[24] Rosa E. Peña y Jesus A. Espı́nola, Ecuaciones de Reacción-Difusión, Proyecto por el Centro de Inves-
tigación en Matemáticas, México, Diciembre 2012.

[25] Timothy Sauer, Análisis numérico, segunda ed., Pearson Educación, México, 2013.

[26] J. Schnakenberg, Simple chemical reaction systems with limit cycle behaviour, Journal of theoretical
biology 81 (1979), no 3, 389–400.

[27] Lucinda Spokes, El fitoplancton y los nutrientes de los océanos, Proyecto ESPERE-ENC, Oc-
tubre 2013, http://klimat.czn.uj.edu.pl/enid/2__Nutrientes_en_el_oc_ano/
-_El_fitoplancton_y_los_nutrientes_3vx.html, Revisión cientı́fica: Prof. Tim Jic-
kells y Dr. Keith Weston.

[28] Aleksandr N. Tijonov y Aleksandr A. Samarsky, Ecuaciones de la fı́sica matemática, segunda ed.,
Editorial MIR, Moscú, 1980.

[29] Alan P. Trujillo y Harold V. Thurman, Essentials of oceanography, doceava ed., Pearson, Boston, 2018,
ISBN 978-0-134-07354-5.

[30] A. Turing, The Chemical Basis of Morphogenesis, Philosophical Transactions of the Royal Society of
London. Series B, Biological Sciences 237 (1952), no 641, 37–72, ISSN 00804622, http://www.
jstor.org/stable/92463.

[31] Ranjit Kumar Upadhyay, Sarita Kumari, Pramod Kumar y Vikas Rai, Spatial distribu-
tion of microalgae in marine systems: A reaction–diffusion model, Ecological Comple-
xity 39 (2019), 100771, ISSN 1476-945X, https://www.sciencedirect.com/science/
article/pii/S1476945X1930025X.

[32] Juan Vanegas, Nancy Landinez y Diego Garzón, Análisis de la inestabilidad de Turing en modelos
biológicos, Dyna 76 (2009), no 158, 123–134.

[33] Riccardo Zamolo y Enrico Nobile, Solution of incompressible fluid flow problems with heat transfer
by means of an efficient RBF-FD meshless approach, Numerical Heat Transfer, Part B: Fundamentals
75 (2019), no 1, 19–42.