In [1]:

```
from sympy import * #Import symbolic python
import numpy as np #Import the numerical library
import matplotlib.pyplot as plt #Import the plotting library
%config InlineBackend.figure_format = "retina" #This makes graphs crisper
plt.rcParams['figure.figsize'] = 4, 3 #Change the dimension of plots
x, t = symbols('x t') #Define some variables
```

The differential equation is $$\frac{dP}{dt}=kP$$

This is an easily separable equation, but let's have sympy solve it.

In [2]:

```
P = Function('P') #Create the function P(t)
k = Symbol('k') #Define a constant
eq1 = Eq(diff(P(t),t),k*P(t));eq1 #Create and display an equation
```

Out[2]:

$\displaystyle \frac{d}{d t} P{\left(t \right)} = k P{\left(t \right)}$

In [3]:

```
sol1 = dsolve(eq1);sol1 #Solve and display the DE
```

Out[3]:

$\displaystyle P{\left(t \right)} = C_{1} e^{k t}$

The logistic growth equation is given by: $$\frac{dP}{dt}=b P\left(1-\frac{P}{K}\right)$$

Notice: if $P=0$ the growth rate of the population is 0, and if $P=K$ the growth rate is also 0. $K$ is called the carrying capicity, in which the population cannot grow larger than this value. Solve for $P$.

In [4]:

```
b, K = symbols('b K')
P = Function('P')
```

In [5]:

```
eq2 = Eq(diff(P(t),t),b*P(t)*(1-P(t)/K))
```

In [6]:

```
sol2 = dsolve(eq2); sol2
```

Out[6]:

$\displaystyle P{\left(t \right)} = \frac{K e^{C_{1} K + b t}}{e^{C_{1} K + b t} - 1}$

dsolve certainly gave us a solution for $P(t)$, but it can be simplified (by hand) into a nicer expression

$$P(t) = \frac{K}{1+ae^{-bt}}$$In [ ]:

```
```

In [7]:

```
t = symbols('t')
Q = Function('Q')
eq3 = Eq(diff(Q(t),t),200-Q(t)/25);eq3
```

Out[7]:

$\displaystyle \frac{d}{d t} Q{\left(t \right)} = 200 - \frac{Q{\left(t \right)}}{25}$

In [8]:

```
sol3 = dsolve(eq3, ics = {Q(0):1000});sol3 #Solve the equation with initial condition Q(0)=1000
```

Out[8]:

$\displaystyle Q{\left(t \right)} = 5000 - 4000 e^{- \frac{t}{25}}$

In [9]:

```
sol3.subs(t,10).evalf() #Evaluate Q(10) using the subs command
```

Out[9]:

$\displaystyle Q{\left(10 \right)} = 2318.71981585744$

The concentration will be $C(10)=\frac{Q(10)}{500}$:

In [10]:

```
2318.72/500
```

Out[10]:

4.63744

What we know is true, is that the *rate of change in quantity* is equal to *the rate of change of quantity coming in* minus *the rate of change in quantity going out*, or

and since the rate of the quantity of ammonia coming in is $conc_{in}\times rate_{in}$ and simlar to the amount going out, we have

$$\frac{dQ}{dt}=C_{in}\cdot r_{in}-C_{out}\cdot r_{out}$$This leads us to the equation $$\frac{dQ}{dt}=5g/L\cdot 10L/m - \frac{Q(t)}{1200+5t}\cdot 5L/m$$

or $$Q^{\prime}=50-\frac{5Q}{1200+5t}$$

We can simplify this to $$Q^{\prime}=50-\frac{Q}{240+t}$$

Finally, we can solve this using the initial condition of $Q(0)=24000g$ and the tank will become full after 160 minutes.

In [11]:

```
Q = Function('Q')
sol4 = dsolve(Eq(diff(Q(t),t),50-Q(t)/(240+t)), ics = {Q(0):24000});sol4
```

Out[11]:

$\displaystyle Q{\left(t \right)} = \frac{25 t^{2} + 12000 t + 5760000}{t + 240}$

In [12]:

```
sol3.subs(t, 160)
```

Out[12]:

$\displaystyle Q{\left(160 \right)} = 5000 - \frac{4000}{e^{\frac{32}{5}}}$

There will be 20800 grams of ammonia in the 2000 L giving a concentration of 10.4g/L

In [13]:

```
plot(sol4.rhs, (t,0,160),title='Quantity',xlabel='time',ylabel='grams')
```

Out[13]:

<sympy.plotting.plot.Plot at 0x7fc4083c46d0>

In [14]:

```
plot(sol4.rhs/(1200+5*t),(t,0,160),title='Concentration',xlabel='time',ylabel='g/L')
```

Out[14]:

<sympy.plotting.plot.Plot at 0x7fc4083eb280>

If $X(t)$ is the quantity in tank A and $Y(t)$ is the quantity in tank B, this gives the equations:

$$ \begin{cases} X^\prime=& 0-20\cdot \frac{X}{500} \\ Y^\prime=& 20\cdot \frac{X}{500}-20\cdot\frac{Y}{800} \\ \end{cases} $$In [15]:

```
X = Function('X') #First solve for X
eqX = Eq(diff(X(t),t),-20*X(t)/500)
solX = dsolve(eqX, ics = {X(0): 4000});solX
```

Out[15]:

$\displaystyle X{\left(t \right)} = 4000 e^{- \frac{t}{25}}$

In [16]:

```
Y = Function('Y') #Substitute X in and solve for Y
eqY = Eq(diff(Y(t),t),20*solX.rhs/500-20*Y(t)/800)
solY = dsolve(eqY, ics ={Y(0): 1600});solY
```

Out[16]:

$\displaystyle Y{\left(t \right)} = - \frac{32000 e^{- \frac{t}{25}}}{3} + \frac{36800 e^{- \frac{t}{40}}}{3}$

In [17]:

```
plot(solX.rhs,solY.rhs,(t,0,100),ylabel="Quantity (oz)")
```

Out[17]:

<sympy.plotting.plot.Plot at 0x7fc4083c1e20>

In [18]:

```
plot((solX.rhs)/500,(solY.rhs)/800,(t,0,100),ylabel="Concentration (oz/gal)")
```

Out[18]:

<sympy.plotting.plot.Plot at 0x7fc43977e250>

First load sympy and matplotlib, and create the variables $t$, $P$, and the function $Q(t)$.

The DE to find the yearly payment is
$$\frac{dQ}{dt}=0.065Q-P$$

where $P$ is the yearly payment. The initial value is $Q(0)=\$500000$ and we also have $Q(30)=0$.

First, solve for $Q(t)$ in terms of $t$ and $P$.

In [19]:

```
P = symbols('P')
sol6 = dsolve(diff(Q(t),t)-0.065*Q(t)+P, ics = {Q(0): 500000});sol6
```

Out[19]:

$\displaystyle Q{\left(t \right)} = 15.3846153846154 P + \left(500000.0 - 15.3846153846154 P\right) e^{0.065 t}$

We can solve for the yearly payment by substituting $t=30$ and $Q=0$, or $Q(30)=0$:

In [20]:

```
Pmt=solve((sol6.rhs).subs(t,30),P)[0];Pmt
```

Out[20]:

$\displaystyle 37890.8913947774$

Suppose we can increase the yearly payment by 5% each year.

First, solve for the yearly payment with increases of 5% each year. (Pv is the varying payment.) That is solve:

$$\frac{dP}{dt}=0.05\cdot P$$where $P_0 = \text{Pmt}$.

In [21]:

```
Pv = Function('Pv')
Pmt2 = dsolve(diff(Pv(t),t)-0.05*Pv(t), ics = {Pv(0):Pmt}); Pmt2
```

Out[21]:

$\displaystyle \operatorname{Pv}{\left(t \right)} = 37890.8913947774 e^{0.05 t}$

Solve the original DE using this varying payment in place of the fixed yearly payment.

In [22]:

```
sol62 = dsolve(diff(Q(t),t)-0.065*Q(t)+Pmt2.rhs, ics = {Q(0): 500000});sol62
```

Out[22]:

$\displaystyle Q{\left(t \right)} = 2526059.42631849 e^{0.05 t} - 2026059.42631849 e^{0.065 t}$

Finally, solve this equation for $t$ to find how long it will take to pay off the mortgage.

In [23]:

```
time = solve(sol62.rhs,t)[0];time
```

Out[23]:

$\displaystyle 14.7045208178974$

In [24]:

```
plot((sol6.rhs).subs(P,Pmt),sol62.rhs,(t,0,30),ylim=(0,500000))
```

Out[24]:

<sympy.plotting.plot.Plot at 0x7fc42863b550>

In [25]:

```
eq7 = Eq(diff(Q(t),t),0.03*0.5 - Q(t)/2000*0.5);eq7
```

Out[25]:

$\displaystyle \frac{d}{d t} Q{\left(t \right)} = 0.015 - 0.00025 Q{\left(t \right)}$

In [26]:

```
sol7 = dsolve(eq7, ics = {Q(0): 0});sol7
```

Out[26]:

$\displaystyle Q{\left(t \right)} = 60.0 - 60.0 e^{- 0.00025 t}$

In [27]:

```
solve(sol7.rhs - 0.00015*2000,t)
```

Out[27]:

[20.0501672941771]

It will take a little over 20 minutes to reach 0.00015 or about 150 parts per million.

In [ ]:

```
```