1034.sagews

Suppose my initial traffic density is given by $\rho(x,0)=e^{-(x-2)^2}$.  I'll call this initial function $f(x)$, so $\rho(x,0)=f(x)$.

var('x,t')
f(x)=e^(-(x-2)^2)
plot(f(x), (x,0,20), axes_labels=['$x$',r'$\rho(x,0)$'])

Looking at the plot, we see a high-density area centered around $x=2$ on the roadway at time 0.

Furthermore, suppose that the partial differential equation modeling my traffic density is $3\rho_x+\rho_t=0$.  Then from what we did in lecture, $c=3$.

c=3

The method of characteristics helps us understand that for this type of partial differential equation, $\rho(x,t)$ will have the same values along the line $x-ct=D$, where $D$ is some constant.  Let's draw a bunch of those lines on the $x,t$ plane.  First, we'll solve for $t$ in terms of $x$, to get $t=\frac{x}{c}+E$ for a constant $E$ (we know that $E=-\frac{D}{c}$, but it is still just a constant).

Along each of the lines in the plot below, the function $\rho(x,t)$ is constant.  In other words, these are level curves or contour lines of $\rho(x,t)$.

characteristics=sum(plot(x/c+E, (x,0,20)) for E in [-5..5])
characteristics

To solve for $\rho(x,t)$, we need to find the traffic density, given a position $x$ and a time $t$.  Suppose we were given the position $x=10$ and time $t=3$.

characteristics+point([12,3],color='red', pointsize=50)

Since we know that $\rho(x,t)$ is the same along any level curve above, we know that value $\rho(12,3)$ is the same as the value of $\rho(3,0)$.  In other words, $\rho$ has the same value for both points shown below, since they are on the same level curve.

characteristics+point([[12,3], [3,0]],color='red', pointsize=50)

We know what $\rho(3,0)$ is, since $\rho(x,0)=f(x)$.  So $\rho(12,3)=\rho(3,0)=f(3)=e^{-(3-2)^2}=e^{-1}\approx .367879$.

f(3).n()

0.367879441171442

In fact, since the slope of each contour line is $1/c$, if I'm at the point $(x,t)$, then the point $(x-ct,0)$ is also on the same contour line (since then the slope would be $\frac{t-0}{x-(x-ct)}=\frac{1}{c}$).  Therefore, $\rho(x,t)=\rho(x-ct,0)=f(x-ct)=e^{-{(3t - x + 2)}^{2}}$.

rho(x,t)=f(x-c*t)
rho(x,t)

e^(-(3*t - x + 2)^2)

Now we can find the densities at specific places and times easily, like a distance of 6 at time $t=1.9$.

rho(6,1.9)

0.0555762126114832

Now that we know what $\rho(x,t)$ is, let's plot it in 3d.

plot3d(rho(x,t), (x,0,20), (t,0,12),plot_points=100)

We can also check our contour lines for $\rho(x,t)$.

contour_plot(rho(x,t), (x,0,20), (t,0,12),fill=False,cmap='winter')

We can also plot the density $\rho(x,t)$ for various values of $t$ to see how the traffic density "moves" with time.

html.table([[plot(rho(x,t), (x,0,20),axes_labels=['$x$',r'$\rho(x,%d)$'%t]) for t in r] for r in [[0,1],[2,3],[4,5]]])

Or we can animate it.

animate([plot(rho(x,t), (x,0,20))+text("time=%d"%t,(10,.5)) for t in [0..5]]).show() # Wait a few seconds for this one

Or we could make it an interact.

@interact
def plot_rho(t=(0,5)):
    plot(rho(x,t), (x,0,20),axes_labels=['$x$',r'$\rho(x,%s)$'%t]).show()