Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download
563 views
Tweet
Share
Share

Modelo matemático para problemas de contaminación

Ecuación Streeter-Phelps

Tuberías de Aguas Residuales 
### Ecuación de Streeter - Phelps
L0, k1, k2, D0 = var('L0, k1, k2, D0')

D(t) = L0*k1/(k2 - k1)*(e^(-k1*t) - e^(-k2*t)) - D0*e^(-k2*t)
print("La Ecuación Streeter-Phelps estudiada es:")
print('D(t) = ')
show(D(t))

## Donde:
#  D(t) = deficit de oxígeno en el tiempo t (mg de O2/L)
#  L0 = Demanda bioquímicade oxígeno inicial
#  k1 = tasa de desoxigenación
#  k2 = tasa de ocigenación
#  D0 = Déficit de oxígeno inicial
#  L0, k1, k2 y D0 = constantes
#  t = tiempo
La Ecuación Streeter-Phelps estudiada es: D(t) =
−L0k1(e(−k1t)−e(−k2t))k1−k2−D0e(−k2t)\displaystyle -\frac{L_{0} k_{1} {\left(e^{\left(-k_{1} t\right)} - e^{\left(-k_{2} t\right)}\right)}}{k_{1} - k_{2}} - D_{0} e^{\left(-k_{2} t\right)}

Derivada de la función D(t)

# Derivamos la Ecuacion D(t):
dDdt(t) = diff(D(t), t);
print("La derivada de la función D(t), es: ")
print("dDdt(t) = ")
show(dDdt(t))
La derivada de la función D(t), es: dDdt(t) =
D0k2e(−k2t)+(k1e(−k1t)−k2e(−k2t))L0k1k1−k2\displaystyle D_{0} k_{2} e^{\left(-k_{2} t\right)} + \frac{{\left(k_{1} e^{\left(-k_{1} t\right)} - k_{2} e^{\left(-k_{2} t\right)}\right)} L_{0} k_{1}}{k_{1} - k_{2}}

Reesolvemos para los valores del inciso a)

valores del inciso a):

L0=−0.5L0 = -0.5
k1=−0.3k1 = -0.3
k2=1k2 = 1
D0=−5D0 = -5
## Puntos críticos para el inciso a):
L0 = -0.5
k1 = -0.3
k2 = 1
D0 = -5
D(t) = L0*k1/(k2 - k1)*(e^(-k1*t) - e^(-k2*t)) - D0*e^(-k2*t);
dDdt(t) = diff(D(t), t);

# Igualamos la derivado  a cero (dDdt(t) = 0) y resolvemos para "t":
puntos_criticos = solve(dDdt(t) == 0, t);
print("Los tiempo críticos son: ")
for pts in puntos_criticos:
    show(pts)
print("\nEl único tiempo real, es: ")
print("\nt = ")
t_min = puntos_criticos[-1].rhs()
show(n(t_min))
Los tiempo críticos son:
t=10 log⁡(19⋅127011391213e(213i π))\displaystyle t = 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}} e^{\left(\frac{2}{13} i \, \pi\right)}\right)
t=10 log⁡(19⋅127011391213e(413i π))\displaystyle t = 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}} e^{\left(\frac{4}{13} i \, \pi\right)}\right)
t=10 log⁡(19⋅127011391213e(613i π))\displaystyle t = 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}} e^{\left(\frac{6}{13} i \, \pi\right)}\right)
t=10 log⁡(19⋅127011391213e(813i π))\displaystyle t = 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}} e^{\left(\frac{8}{13} i \, \pi\right)}\right)
t=10 log⁡(19⋅127011391213e(1013i π))\displaystyle t = 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}} e^{\left(\frac{10}{13} i \, \pi\right)}\right)
t=10 log⁡(19⋅127011391213e(1213i π))\displaystyle t = 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}} e^{\left(\frac{12}{13} i \, \pi\right)}\right)
t=−12013i π+10 log⁡(19⋅127011391213)\displaystyle t = -\frac{120}{13} i \, \pi + 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}}\right)
t=−10013i π+10 log⁡(19⋅127011391213)\displaystyle t = -\frac{100}{13} i \, \pi + 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}}\right)
t=−8013i π+10 log⁡(19⋅127011391213)\displaystyle t = -\frac{80}{13} i \, \pi + 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}}\right)
t=−6013i π+10 log⁡(19⋅127011391213)\displaystyle t = -\frac{60}{13} i \, \pi + 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}}\right)
t=−4013i π+10 log⁡(19⋅127011391213)\displaystyle t = -\frac{40}{13} i \, \pi + 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}}\right)
t=−2013i π+10 log⁡(19⋅127011391213)\displaystyle t = -\frac{20}{13} i \, \pi + 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}}\right)
t=10 log⁡(19⋅127011391213)\displaystyle t = 10 \, \log\left(\frac{1}{9} \cdot 1270^{\frac{1}{13}} 9^{\frac{12}{13}}\right)
El único tiempo real, es: t =
3.80734430932032\displaystyle 3.80734430932032

Gráfica de D(t) para los valores del inciso a): 

# Gráfica de la Ecuación de Streeter - Phelps:
L0 = -0.5
k1 = -0.3
k2 = 1
D0 = -5
D(t) = L0*k1/(k2 - k1)*(e^(-k1*t) - e^(-k2*t)) - D0*e^(-k2*t);

## Funcion D(t):
D_plot = plot(D(t),(-10,30),ymax=100, ymin=-1)

# Valores mínimo:
t_min = puntos_criticos[-1].rhs()
D_min = D(t_min)
print("Min on interval [-10, 30]: f({0}) = {1}".format(n(t_min), n(D_min)))

# Gráfica del punto mínimo:
min_point = point((t_min, D(t_min)), color='red', size=50)
show(D_plot + min_point, figsize=(4, 4))
Min on interval [-10, 30]: f(3.80734430932032) = 0.470049771963661

Reesolvemos para los valores del inciso b)

valores del inciso b):

L0=5L0 = 5
k1=2.6k1 = 2.6
k2=1.2k2 = 1.2
D0=5D0 = 5
## Puntos críticos para el inciso b):
L0 = 5
k1 = 2.6
k2 = 1.2
D0 = 5
D(t) = L0*k1/(k2 - k1)*(e^(-k1*t) - e^(-k2*t)) - D0*e^(-k2*t);
dDdt(t) = diff(D(t), t);

# Igualamos la derivado  a cero (dDdt(t) = 0) y resolvemos para "t":
puntos_criticos = solve(dDdt(t) == 0, t);
print("Los tiempo críticos son: ")
for pts in puntos_criticos:
    show(pts)
print("\nEl único tiempo real, es: ")
print("\nt = ")
t_min = puntos_criticos[-1].rhs()
show(n(t_min))
Los tiempo críticos son:
t=5 log⁡(136⋅169173667e(27i π))\displaystyle t = 5 \, \log\left(\frac{1}{36} \cdot 169^{\frac{1}{7}} 36^{\frac{6}{7}} e^{\left(\frac{2}{7} i \, \pi\right)}\right)
t=5 log⁡(136⋅169173667e(47i π))\displaystyle t = 5 \, \log\left(\frac{1}{36} \cdot 169^{\frac{1}{7}} 36^{\frac{6}{7}} e^{\left(\frac{4}{7} i \, \pi\right)}\right)
t=5 log⁡(136⋅169173667e(67i π))\displaystyle t = 5 \, \log\left(\frac{1}{36} \cdot 169^{\frac{1}{7}} 36^{\frac{6}{7}} e^{\left(\frac{6}{7} i \, \pi\right)}\right)
t=−307i π+5 log⁡(136⋅169173667)\displaystyle t = -\frac{30}{7} i \, \pi + 5 \, \log\left(\frac{1}{36} \cdot 169^{\frac{1}{7}} 36^{\frac{6}{7}}\right)
t=−207i π+5 log⁡(136⋅169173667)\displaystyle t = -\frac{20}{7} i \, \pi + 5 \, \log\left(\frac{1}{36} \cdot 169^{\frac{1}{7}} 36^{\frac{6}{7}}\right)
t=−107i π+5 log⁡(136⋅169173667)\displaystyle t = -\frac{10}{7} i \, \pi + 5 \, \log\left(\frac{1}{36} \cdot 169^{\frac{1}{7}} 36^{\frac{6}{7}}\right)
t=5 log⁡(136⋅169173667)\displaystyle t = 5 \, \log\left(\frac{1}{36} \cdot 169^{\frac{1}{7}} 36^{\frac{6}{7}}\right)
El único tiempo real, es: t =
1.10455698319069\displaystyle 1.10455698319069

Gráfica de D(t) para los valores del inciso b): 

# Gráfica de la Ecuación de Streeter - Phelps:
L0 = 5
k1 = 2.6
k2 = 1.2
D0 = 5
D(t) = L0*k1/(k2 - k1)*(e^(-k1*t) - e^(-k2*t)) - D0*e^(-k2*t);

## Funcion D(t):
D_plot = plot(D(t),(-5,20),ymax=5, ymin=-20)

# Valores mínimo:
puntos_criticos = solve(dDdt(t) == 0, t);
t_min = puntos_criticos[-1].rhs()
D_min = D(t_min)
print("Min on interval [-5, 20]: f({0}) = {1}".format(n(t_min), n(D_min)))

# Gráfica del punto mínimo:
min_point = point((t_min, D(t_min)), color='red', size=50)
show(D_plot + min_point, figsize=(4, 4))
Min on interval [-5, 20]: f(1.10455698319069) = 0.613104211144102

Reesolvemos para los valores del inciso c)

valores del inciso c):

L0=5L0 = 5
k1=2.4k1 = 2.4
k2=0.3k2 = 0.3
D0=−3D0 = -3
## Puntos críticos para el inciso c):
L0 = 5
k1 = 2.4
k2 = 0.3
D0 = -3
D(t) = L0*k1/(k2 - k1)*(e^(-k1*t) - e^(-k2*t)) - D0*e^(-k2*t);
dDdt(t) = diff(D(t), t);

# Igualamos la derivado  a cero (dDdt(t) = 0) y resolvemos para "t":
puntos_criticos = solve(dDdt(t) == 0, t);
print("Los tiempos críticos son: ")
for pts in puntos_criticos:
    show(pts)

###
f = e^(1/5*t) == n(1/61*61^(11/12)*5^(1/12)*sqrt(2))*e^(1/40*t)
sol_t = solve([f],t)
show(sol_t)
###
f1 = e^(1/5*t)/e^(1/40*t) == (8353829)/(7276132)
###
f2 = log(e^(7/40*t)) == log((8353829)/(7276132))
###
f3 = t == 40/7*log((8353829)/(7276132))
show(simplify(f3))
###
f4 = t == n(40/7*log((8353829)/(7276132)))
print("\nEl único tiempo real, es: ")
show(simplify(f4))
Los tiempos críticos son:
e(15 t)=1122⋅6111122(51123+i⋅5112)e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = \frac{1}{122} \cdot 61^{\frac{11}{12}} \sqrt{2} {\left(5^{\frac{1}{12}} \sqrt{3} + i \cdot 5^{\frac{1}{12}}\right)} e^{\left(\frac{1}{40} \, t\right)}
e(15 t)=1122⋅6111122(i⋅51123+5112)e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = \frac{1}{122} \cdot 61^{\frac{11}{12}} \sqrt{2} {\left(i \cdot 5^{\frac{1}{12}} \sqrt{3} + 5^{\frac{1}{12}}\right)} e^{\left(\frac{1}{40} \, t\right)}
e(15 t)=161i⋅61111251122e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = \frac{1}{61} i \cdot 61^{\frac{11}{12}} 5^{\frac{1}{12}} \sqrt{2} e^{\left(\frac{1}{40} \, t\right)}
e(15 t)=1122⋅6111122(i⋅51123−5112)e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = \frac{1}{122} \cdot 61^{\frac{11}{12}} \sqrt{2} {\left(i \cdot 5^{\frac{1}{12}} \sqrt{3} - 5^{\frac{1}{12}}\right)} e^{\left(\frac{1}{40} \, t\right)}
e(15 t)=−1122⋅6111122(51123−i⋅5112)e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = -\frac{1}{122} \cdot 61^{\frac{11}{12}} \sqrt{2} {\left(5^{\frac{1}{12}} \sqrt{3} - i \cdot 5^{\frac{1}{12}}\right)} e^{\left(\frac{1}{40} \, t\right)}
e(15 t)=−161⋅61111251122e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = -\frac{1}{61} \cdot 61^{\frac{11}{12}} 5^{\frac{1}{12}} \sqrt{2} e^{\left(\frac{1}{40} \, t\right)}
e(15 t)=−1122⋅6111122(51123+i⋅5112)e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = -\frac{1}{122} \cdot 61^{\frac{11}{12}} \sqrt{2} {\left(5^{\frac{1}{12}} \sqrt{3} + i \cdot 5^{\frac{1}{12}}\right)} e^{\left(\frac{1}{40} \, t\right)}
e(15 t)=−1122⋅6111122(i⋅51123+5112)e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = -\frac{1}{122} \cdot 61^{\frac{11}{12}} \sqrt{2} {\left(i \cdot 5^{\frac{1}{12}} \sqrt{3} + 5^{\frac{1}{12}}\right)} e^{\left(\frac{1}{40} \, t\right)}
e(15 t)=−161i⋅61111251122e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = -\frac{1}{61} i \cdot 61^{\frac{11}{12}} 5^{\frac{1}{12}} \sqrt{2} e^{\left(\frac{1}{40} \, t\right)}
e(15 t)=−1122⋅6111122(i⋅51123−5112)e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = -\frac{1}{122} \cdot 61^{\frac{11}{12}} \sqrt{2} {\left(i \cdot 5^{\frac{1}{12}} \sqrt{3} - 5^{\frac{1}{12}}\right)} e^{\left(\frac{1}{40} \, t\right)}
e(15 t)=1122⋅6111122(51123−i⋅5112)e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = \frac{1}{122} \cdot 61^{\frac{11}{12}} \sqrt{2} {\left(5^{\frac{1}{12}} \sqrt{3} - i \cdot 5^{\frac{1}{12}}\right)} e^{\left(\frac{1}{40} \, t\right)}
e(15 t)=161⋅61111251122e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = \frac{1}{61} \cdot 61^{\frac{11}{12}} 5^{\frac{1}{12}} \sqrt{2} e^{\left(\frac{1}{40} \, t\right)}
[e(15 t)=83538297276132 e(140 t)\displaystyle e^{\left(\frac{1}{5} \, t\right)} = \frac{8353829}{7276132} \, e^{\left(\frac{1}{40} \, t\right)}]
t=407 log⁡(83538297276132)\displaystyle t = \frac{40}{7} \, \log\left(\frac{8353829}{7276132}\right)
El único tiempo real, es:
t=0.789260538867\displaystyle t = 0.789260538867

Gráfica de D(t) para los valores del inciso c): 

# Gráfica de la Ecuación de Streeter - Phelps:
L0 = 5
k1 = 2.4
k2 = 0.3
D0 = -3
D(t) = L0*k1/(k2 - k1)*(e^(-k1*t) - e^(-k2*t)) - D0*e^(-k2*t);

## Funcion D(t):
D_plot = plot(D(t),(-5,20),ymax=10, ymin=-20)

# Valores mínimo:
t_min = n(40/7*log((8353829)/(7276132)))
D_min = D(t_min)
print("Min on interval [-5, 20]: f({0}) = {1}".format(n(t_min), n(D_min)))

# Gráfica del punto mínimo:
min_point = point((t_min, D(t_min)), color='red', size=50)
show(D_plot + min_point, figsize=(4, 4))
Min on interval [-5, 20]: f(0.789260538866877) = 6.01739331162348

Reesolvemos para los valores del inciso d)

valores del inciso c):

L0=−1.4L0 = -1.4
k1=−0.4k1 = -0.4
k2=0.3k2 = 0.3
D0=−3D0 = -3
## Puntos críticos para el inciso d):
L0 = -1.4
k1 = -0.4
k2 = 0.3
D0 = -3
D(t) = L0*k1/(k2 - k1)*(e^(-k1*t) - e^(-k2*t)) - D0*e^(-k2*t);
dDdt(t) = diff(D(t), t);

# Igualamos la derivado  a cero (dDdt(t) = 0) y resolvemos para "t":
puntos_criticos = solve(dDdt(t) == 0, t);
print("Los tiempos críticos son: ")
for pts in puntos_criticos:
    show(pts)
print("\nEl único tiempo real, es: ")
print("\nt = ")
t_min = puntos_criticos[-1].rhs()
show(n(t_min))
Los tiempos críticos son:
t=10 log⁡(116⋅33171667e(27i π))\displaystyle t = 10 \, \log\left(\frac{1}{16} \cdot 33^{\frac{1}{7}} 16^{\frac{6}{7}} e^{\left(\frac{2}{7} i \, \pi\right)}\right)
t=10 log⁡(116⋅33171667e(47i π))\displaystyle t = 10 \, \log\left(\frac{1}{16} \cdot 33^{\frac{1}{7}} 16^{\frac{6}{7}} e^{\left(\frac{4}{7} i \, \pi\right)}\right)
t=10 log⁡(116⋅33171667e(67i π))\displaystyle t = 10 \, \log\left(\frac{1}{16} \cdot 33^{\frac{1}{7}} 16^{\frac{6}{7}} e^{\left(\frac{6}{7} i \, \pi\right)}\right)
t=−607i π+10 log⁡(116⋅33171667)\displaystyle t = -\frac{60}{7} i \, \pi + 10 \, \log\left(\frac{1}{16} \cdot 33^{\frac{1}{7}} 16^{\frac{6}{7}}\right)
t=−407i π+10 log⁡(116⋅33171667)\displaystyle t = -\frac{40}{7} i \, \pi + 10 \, \log\left(\frac{1}{16} \cdot 33^{\frac{1}{7}} 16^{\frac{6}{7}}\right)
t=−207i π+10 log⁡(116⋅33171667)\displaystyle t = -\frac{20}{7} i \, \pi + 10 \, \log\left(\frac{1}{16} \cdot 33^{\frac{1}{7}} 16^{\frac{6}{7}}\right)
t=10 log⁡(116⋅33171667)\displaystyle t = 10 \, \log\left(\frac{1}{16} \cdot 33^{\frac{1}{7}} 16^{\frac{6}{7}}\right)
El único tiempo real, es: t =
1.03416977032385\displaystyle 1.03416977032385

Gráfica de D(t) para los valores del inciso d): 

# Gráfica de la Ecuación de Streeter - Phelps:
L0 = -1.4
k1 = -0.4
k2 = 0.3
D0 = -3
D(t) = L0*k1/(k2 - k1)*(e^(-k1*t) - e^(-k2*t)) - D0*e^(-k2*t);

## Funcion D(t):
D_plot = plot(D(t),(-20,20),ymax=100, ymin=-20)

# Valores mínimo:
t_min = n(40/7*log((8353829)/(7276132)))
D_min = D(t_min)
print("Min on interval [-20, 20]: f({0}) = {1}".format(n(t_min), n(D_min)))

# Gráfica del punto mínimo:
min_point = point((t_min, D(t_min)), color='red', size=50)
show(D_plot + min_point, figsize=(4, 4))
Min on interval [-20, 20]: f(0.789260538866877) = 2.83314562606401