# Direction fields¶

In [1]:
t, w, x, y, n = var('t w x y n')
palette = [(215/255, 0/255, 132/255), (255/255, 1/255, 73/255), (255/255, 121/255, 1/255), (255/255, 210/255, 0/255)]
cool_palette = [(0/255, 150/255, 173/255), (0/255, 200/255, 146/255)]


## Filling a water tank¶

### Scenario¶

You're pumping water into the bottom of a cylindrical tank. You want to know $y(t)$, the water level after $t$ seconds of pumping. The tank starts out empty. The pump pushes water into the tank with a constant pressure. As the water level rises, the pressure working against the pump grows proportionally.

### Example model¶

Let's suppose the flow rate is proportional to the pressure difference between the tank and the pump. (This assumption is based on the "Hagenâ€“Poiseuille law," which describes flow through a long, narrow pipe.) If the pump pressure is $100$ (in kilopascals), the tank pressure is $10y$ (in kilopascals per meter), and the constant of proportionality is $\tfrac{1}{500}$ (in meters per minute per kilopascal), the flow rate is related to the amount of water in the tank by the equation

$$y' = \tfrac{1}{500}(100 - 10y).$$

In [122]:
water_direction = plot_slope_field((1/500) * (100 - 10*y), (t, 0, 150), (y, 0, 12))
water_solutions = [plot(exp(-t/50) * 4*k + 10 * (1 - exp(-t/50)), (t, 0, 150), thickness = 3, rgbcolor = palette[k]) for k in range(4)]
water_labels = ["$t$ : seconds", "$y(t)$ : meters"]
show(water_direction, axes_labels = water_labels)
show(water_direction + water_solutions[0], axes_labels = water_labels)


### Other initial water levels¶

What if the tank starts out with some water in it?

In [123]:
show(water_direction + sum(water_solutions), axes_labels = water_labels)


## Wearing out truck tires¶

### Scenario¶

You order a batch of new tires for a very large fleet of trucks. You want to know $n(t)$, the number of tires that are still intact after $t$ distance units of driving.

### Example model¶

$$n' = -5t^4 n$$

In [2]:
tires_direction = plot_slope_field(-5*t^4 * n, (t, 0, 2), (n, 0, 4000))
tires_solutions = [plot(1000*(4-k) * exp(-t^5), (t, 0, 2), thickness = 3, rgbcolor = palette[k]) for k in range(4)]
tires_labels = ["$t$ : distance units", "$n(t)$ : tires"]
show(tires_direction, axes_labels = tires_labels)
show(tires_direction + tires_solutions[0], axes_labels = tires_labels)


### Other initial numbers of tires¶

What if you start out with fewer tires?

In [3]:
show(tires_direction + sum(tires_solutions), axes_labels = tires_labels)


## Swinging a pendulum¶

### Scenario¶

You attach a heavy weight to the ceiling using a light string. You pull the weight 60° back from vertical, keeping the string taut, and let it go. You want to know $\theta(t)$, the angle of the string after $t$ seconds of swinging.

### Example model¶

$$\theta'' = -\sin \theta$$

In [114]:
energies = [1 - cos(k*pi/12) for k in reversed(range(1,5))]
pendulum_solutions = [plot(asin(sqrt(en/2) * jacobi_sn(t, en/2)), (t, elliptic_kc(en/2), 9*elliptic_kc(en/2)), thickness = 3, rgbcolor = color) for (en, color) in zip(energies, palette)]
pendulum_labels = ["$t$ : seconds", "$\\theta(t)$ : radians"]
show(pendulum_solutions[0], axes_labels = pendulum_labels)