Aulas do Curso de Modelagem matemática IV da FGV-EMAp

Path: Modelagem-Matematica-IV/Planilhas Sage / Lista de Exercícios 1-Soluções.ipynb

Views: ²⁴¹²
License: GPL3
Image: default

Kernel: SageMath 8.1

1ª Lista de Exercícios

Cinética Química

Neste exercícios vamos construir modelos, resolvê-los analitica e numéricamente, e discutir a estabilidade dos seus equilíbrios.

I. Considere a equação de produção e decaimento introduzida abaixo. Em biologia celular, é comum a produção por tempo limitado, de uma molécular como resposta a um estimulo nervoso ou hormonal.

\frac{dx}{dt}= I - \gamma x

Modifique o modelo de Produção e decaimento de forma que $I$ seja positivo apenas entre os tempos $t=4$ e $t=10$ . Dica: Use a função de Heaviside mostrada abaixo.
Como você escreveria o modelo resultante em notação matemática?
Resolva o modelo e discuta o seu equilíbrio.
Verifique se o equilíbrio do modelo depende da sua condição inicial, x(0)

In [1]:

%display typeset

In [2]:

heaviside?

In [3]:

plot(heaviside(x-4),(x,0,10))

In [4]:

plot(heaviside(x-4)-heaviside(x-10),(x,0,15))

In [5]:

def fun(t,x):
    I = heaviside(t-4)-heaviside(t-10)
    return [I - 0.3*x[0]]
T = ode_solver()
T.algorithm = "rk8pd"
T.function = fun
T.ode_solve(y_0=[2],t_span=[0,15], num_points=100)
#T.ode_solve?

In [6]:

T.plot_solution(i=0,interpolate=False, gridlines=true)

II. Uma substância com meia-vida de um dia é produzida a uma taxa constante de $10μMh^{−1}$ . Suponha que a sua concentração seja denotada por $C(t)$ .

Use uma equação diferencial para descrever a dinâmica deste processo. Usando o Sage encontre a solução deste modelo analíticamente.
Se um inibidor for aplicado em $t=0$ , de forma que a substância não seja mais produzida. Encontre a solução $C(t)$ e a use para mostrar the $C(t)→0$ quando $t→∞$
Suponha agora que começando de $C(0)=C_0$ , aplica-se uma droga que inibe o decaimento da substância completamente, sem afetar a sua taxa de produção. Mostre que a substância se acumulará a uma taxa linear $C(t)=C_0+kt$ e encontre o valor de $k$ .

In [7]:

var('C t k C0 I')
C=function('C')(t)
dCdt = diff(C,t)== I-k*C
sol = desolve(dCdt,C,ivar=t, ics=[0,C0])
show(sol)

In [8]:

simplify(4*k*exp(-k*24)/k)

In [9]:

f(t) = sol(I=0,k=1/50,C0=4)
plot (f,0,24)

III. Imagine dois tanques de reação $T$ e $U$ , que se comunicam com fluxos $k_1=k_2$ A substância $A$ é introduzida no tanque $T$ a uma taxa constante $I$ e a substância $B$ é introduzida no tanque $U$ a um taxa constante $J$ . A substância $B$ decai a uma taxa $\gamma_2$ apenas no tanque $T$ enquanto que a substância $A$ decai a uma taxa $\gamma_1$ , apenas no tanque $U$ .

Escreva um sistema de equações diferenciais descrevendo a evolução temporal das duas substâncias em cada um dos Tanques
Resolva o sistema analíticamente, usando o Sage. e discuta os equilíbrios e sua estabilidade.
Apresente a evolução da solução em um campo vetorial
Resolva numericamente o sistema no Sage
Assuma que $I=sen(t)$ , encontre a nova solução do sistema e seu(s) equilibrio(s)
Assuma que a produção $J$ de $B$ é inibida pela concentração de $A$ no tanque $U$ , e vice-versa: $A$ produção $I$ de $A$ é inibida pela Concentração de $B$ no tanque $T$ . Encontre a solução do sistema e analise os seus equilíbrios.

In [10]:

forget()
var('C_at C_au C_bt C_bu I J k_1 k_2 gamma_1 gamma_2 t')
C_at = function('C_at')(t)
C_au = function('C_au')(t)
C_bt = function('C_bt')(t)
C_bu = function('C_bu')(t)
dcatdt =  I-k_1*C_at + k_2*C_au
dcaudt =  k_1*C_at - gamma_1*C_au - k_2*C_au
dcbtdt =  k_2*C_bu - (k_1+gamma_2)*C_bt
dcbudt =  J-k_2*C_bu +k_1*C_bt

#assume(k_2^2+2*k_1*k_2+2*gamma_1*k_2+k_1^2-2*gamma_1*k_1+gamma_1^2>0)
#assume(k_1>0)
#assume(k_2>0)
#assume(gamma_1>0)
#assume(gamma_2>0)
#k_1=k2=0.1
#gamma_1=gamma_2=0.01
#I=J=1

sol = desolve_system([dcatdt,dcaudt,dcbtdt,dcbudt],[C_at,C_au,C_bt,C_bu],ivar=t, ics=[0,1,1,1,1])
show(sol)

In [11]:

forget()
assumptions()

In [12]:

var('C_at C_au C_bt C_bu I J k_1 k_2 gamma_1 gamma_2')
ss = solve([I-k_1*C_at + k_2*C_au,k_1*C_at - gamma_1*C_au - k_2*C_au, k_2*C_bu - (k_1+gamma_2)*C_bt, J-k_2*C_bu +k_1*C_bt],C_at,C_au,C_bt,C_bu)
ss

In [13]:

jac = jacobian([I-k_1*C_at, k_1*C_at - gamma_1*C_au, k_2*C_bu - gamma_2*C_bt, J-k_2*C_bu],[C_at,C_au,C_bt,C_bu])
jac

In [14]:

jac.eigenvalues()

In [15]:

J=1
I=1
solAT = sol[0].rhs()
solAU = sol[1].rhs()
solBT = sol[2].rhs()
solBU = sol[3].rhs()

plot([solAT(I=1,k_1=0.5,gamma_1=1,gamma_2=1,k_2=1),solAU(I=1,gamma_1=1)],(t,0,10))#, color='red')

In [16]:

plot([solBT(J=1,k_2=1,gamma_2=2), solBU(J=1,k_2=1,gamma_2=3,k_1=0.5)],(t,0,10))

In [17]:

solBT,solBU(I=1,k_1=0.5)

In [18]:

def fun3(t,y, params):
    C_at,C_au,C_bt,C_bu = y
    k_1,k_2,gamma_1,gamma_2,I,J = params
    I = sin(t)
    return[I-k_1*C_at + k_2*C_au,k_1*C_at - gamma_1*C_au - k_2*C_au, k_2*C_bu - (k_1+gamma_2)*C_bt, J-k_2*C_bu +k_1*C_bt]
    
show(fun3(0,[10,10,10,10],(1,1,1,1,.5,.5)))

T = ode_solver()
T.algorithm = "rk8pd"
T.function = fun3
T.ode_solve(y_0=[10,10,10,10],t_span=[0,20],params=(1,1,1,1.5,2,2), num_points=100)

In [19]:

def plot_sol(sol):
    cat = list_plot([(i[0],i[1][0]) for i in sol],color='red', plotjoined=True, legend_label=r'$C_{AT}$', alpha=.8)
    cau = list_plot([(i[0],i[1][1]) for i in sol],color='blue', plotjoined=True, legend_label=r'$C_{AU}$', alpha=.8, axes_labels=["t","x"],)# marker='v')#, gridlines=True)
    cbt = list_plot([(i[0],i[1][2]) for i in sol],color='green', plotjoined=True, legend_label=r'$C_{BT}$', alpha=.8)
    cbu = list_plot([(i[0],i[1][3]) for i in sol],color='black', plotjoined=True, legend_label=r'$C_{BU}$', alpha=.8)#, gridlines=true)
    
    show(cat+cau+cbt+cbu)
plot_sol(T.solution)

Análise Dimensional e Adimensionalização

IV. Considere o seguinte modelo de crescimento populacional: $dN/dt=rN$ , onde $N(0)=N_0>0$ e $r>0$ . Neste modelo $r$ é a taxa de crescimento e $N_0$ é a população inicial. Lembre-se que a solução deste modelo é $N(t)=N_0 e^{rt}$ .

Re-escale este modelo em unidades da população inicial, ou seja, defina $y(t)=N(t)/N_0$ . Qual a equação resultante e quais as condições iniciais correspondentes?
Quais as unidades de $r$ ?
Qual o "tempo de duplicação" desta população? isto é o tempo $τ$ para o qual $y(τ)=yN_0$ .
Mostre que é possivel definir um tempo adimensional $s$ tal que o modelo se transforme em $dy/ds=y$ , $y(0)=1$ .

Solução:

Como $N_0$ é uma constante, temos que $y(t) = \frac{N(t)}{N_0}$ representa apenas uma mudança de escala do modelo. Temos, portanto, que $N(t) = N_0y(t)$ e a equação original se transforma em:

\frac{d(N_0y(t))}{dt} = r (N_0y(t))

N_0 \frac{dy(t)}{dt} = rN_0y(t) .

Se $N_0 \neq 0$ , temos

\frac{dy}{dt} = ry(t) .

Como $y(t) = \frac{N(t)}{N_0}$ , segue que

y(0) = \frac{N(0)}{N_0} = \frac{N_0}{N_0} = 1

e, consequentemente, $y(t) = e^{rt}$ .

Para que o modelo seja coerente, é necessário que a expressão $rN(t)$ tenha unidade $\frac{M}{T}$ . Como $[N(t)] = M$ , segue que $[r] = T^{-1}$ .

Para encontrarmos o "tempo de duplicação" da população, fazemos:

y(\tau) = 2y_0

y(\tau) = 2

e^{r\tau} = 2

r\tau = \log{2}

\tau = \frac{ \log 2}{r} .

V - Considere a seguinte equação para o crescimento de uma única espécie de organismo:

$\frac{dP}{dt}=\frac{v P}{K+P}-\mu P$

Interprete o que estas equações estão dizendo
defina $x=P / Q$ e $s=t/\tau$ onde $Q$ , $\tau$ são escalas a serem escolhidas. Converta a equação para uma forma adimensional em termos destas novas escalas.
Qual seria uma escolha razoável para $Q$ ? e para $\tau$ ?

Solução:

In [20]:

var('P t v mu K Q tau s alpha M')
f(P) = v*(P/(K+P)) - mu*P
f(P)

In [21]:

fd = expand(f(P)).substitute(P=M,K=M)
fd

Como a expressão acima deve ter dimensão $M/T$ , segue que $\mu$ tem unidade $T^-1$ e $v$ tem unidade $M/T$ .

Defina $x=P / Q$ e $s=t/\tau$ onde $Q$ , $\tau$ são escalas a serem escolhidas. Converta a equação para uma forma adimensional em termos destas novas escalas.

In [22]:

eq = P/t == v*P/(K+P) - mu*P
eq

In [23]:

eq = eq.subs(P=x*Q, t=tau*s)
eq

Agora multiplicamos ambos os lados da equação por $τ/Q$ :

In [24]:

eq = eq.multiply_both_sides(tau/Q)
show(eq)
show(html("Expandindo o lado direito...<br>"))
eq.expand()

Expandindo o lado direito...

Agora precisamos escolher as escalas de $Q$ e $t$ . Eles precisam ser definidos como operações entre parâmetros existentes. Note que na equação original $v$ tem unidade $M/t$ para manter o equilíbrio dimensional. como $Q$ tem unidades $M$ , ou massa, as alternativas para sua escala são: $Q=K$ ou $Q=v/μ$ . Já $\tau$ tem unidades de tempo então pode assumir as seguintes escalas: $K/v$ ou $1/\mu$

In [25]:

eq.expand().subs(Q=K, tau=1/mu).simplify()

Agora podemos notar que podemos combinar conjunto de parâmetros delimitado por colchetes abaixo, em um único parâmetro α:

\left[\frac{v}{K \mu}\right] \frac{x}{x+1}

\alpha = \frac{v}{K \mu}

In [26]:

eq.expand().subs(Q=K, tau=1/mu).simplify().subs(v=K*mu*alpha)

Finalmente, obtemos uma versão adimensional simplificada do modelo: $\frac{dx}{ds}=\frac{\alpha x}{K(x+1)} -x$

VI. A dinâmica da lagarta do pinheiro pode ser descrita pelo modelo proposto por Ludwig, Jones e Holling. Este inseto se reproduz e é predado por pássaros.

$\frac{dB}{dt} = r_B B \left(1-\frac{B}{K_B}\right) - \beta \frac{B^2}{\alpha^2 + B^2}$

Explique o significado dos termos desta equação
Re-escreva esta equação em forma adimensional. Há duas escolhas de escalas para a densidade de da lagarta e duas para o tempo.

Solução:

In [27]:

var('t B r_B K_B beta alpha')

In [28]:

B = function('B')(t)
dBdt = diff(B,t) == r_B*B(1-(B/K_B)) - beta((B**2)/(alpha**2 + beta**2))
show(dBdt)

/usr/lib/python2.7/dist-packages/IPython/core/interactiveshell.py:2882: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...)
See http://trac.sagemath.org/5930 for details.
  exec(code_obj, self.user_global_ns, self.user_ns)

O termo $r_B B \left(1-\frac{B}{K_B}\right)$ representa um crescimento logístico da população de lagartas ( $B$ ) com taxa intrínseca de crescimento igual a $r_B$ e capacidade de suporte $K_B$ , como vimos em sala. Já o termo $\beta \frac{B^2}{\alpha^2 + B^2}$ reflete o efeito da predação realizada por pássaros à população de lagartas. Esse termo nos diz que existe um limite superior para a taxa de mortalidade provocada pela predação. Isso se deve ao fato de que o número de predadores é limitado por seu comportamento territorial e, além disso, eles possuem baixa capacidade de busca, o que significa dizer que uma grande população de lagartas não necessariamente implica em uma maior ocorrência de predação por parte dos pássaros. Evidentemente, o limite superior da predação é igual a $\beta$ .

Da equação acima, sabemos que $[\beta] = \frac{M}{T}$ , $[r_B] = T^{-1}$ e $[\alpha] = M$ .

Vamos introduzir novos parâmetros $\mu = \frac{B}{\alpha}$ e $t' = \frac{\beta t}{\alpha}$ . Perceba que isso implica que $\mu$ e $t'$ são adimensionais. A equação original se transforma em:

\frac{d(\alpha \mu)}{d(\frac{\alpha t'}{\beta})} = r_B(1 - \frac{\alpha \mu}{K_B}) - \frac{\alpha \beta \mu}{\alpha^2(1+\mu^2)}

e, consequentemente,

\frac{d\mu}{d t'} = \frac{\alpha r_B}{\beta}(1-\frac{\alpha \mu}{K_B}) - \frac{\mu}{1+\mu^2} .

Agora, fazemos

R = \frac{\alpha r_B}{\beta}, \ Q = \frac{K_B}{\alpha}

e reescrevemos a equação como:

\frac{d\mu}{dt'} = R(1-\frac{\mu}{Q}) - \frac{\mu}{1+\mu^2} .

In [ ]: