Feb 10 homework

For these first two, solve, graph, discuss limits, ... whatever can be done.
In each case, describe what method(s) worked.  
All of these are from Boyce/DiPrima 9th, chap 2.

prelude

directionField[dydx_,xmin_,xmax_,ymin_,ymax_]:=
VectorPlot[
{1/Sqrt[1+((dydx[x,y])^2)],dydx[x,y]/Sqrt[1+((dydx[x,y])^2)]},
{x,xmin,xmax},{y,ymin,ymax},
VectorStyle->"Segment",VectorScale->{0.02,Automatic,Automatic}
]

Solve   t (dy/dt)+ 2 y = (t^2)-t+1   with y=1/2 at t=1, for t > 0.

First, here's the picture.

dydt[t_,y_]=-2y/t+t-1+1/t;
p1=directionField[dydt,0.1,2,0,2]

dydt[1.0,0.5]

Using integrating factor method :
	t dy + 2 y dt = ((t^2)-t+1) dt
A(t) t dy + 2 A(t) y dt = A(t) ((t^2)-t+1)
then
	d(F(t) y) = F(t)dy + (dF/dt)y dt
so	F(t)=t A(t)
((dF(t))/dt)=2A(t) = 2 ((F(t))/t)
which gives
	(dF/F)=2 (dt/t)
or
	ln(F)=2ln(t)+C
F(t) = C (t^2)
A(t) = C t
and the orginal equation becomes (once the C's cancel)
	(t^2)dy + 2 y t dt = d((t^2)y)=((t^3)-(t^2)+t) dt
or
	(t^2)y = ((t^4)/4)-((t^3)/3)+((t^2)/2)+C    (for some other C)
or
	y(t) = ((t^2)/4)-(t/3)+(1/2)+(C/(t^2))
Since we want to have an initial condition of y=0.5 at t=1, we need C to satisfy
	(1/2)=(1/4)-(1/3)+(1/2)+C
or
	C = (1/2)-(1/4)+(1/3)-(1/2)=(1/12)

p2=Plot[((t^2)/4)-(t/3)+(1/2)+(((1/12))/(t^2)),{t,0.03,2},PlotStyle->Red]

p3=ListPlot[{{1,1/2}},AxesLabel->{t,y},PlotStyle->PointSize[0.02]]

Show[p3,p1,p2]

Solve  x (dy/dx)=(Sqrt[1-(y^2)])

This one is separable :
	(dy/(Sqrt[1-(y^2)]))=(dx/x)
The right side is just ln(x); the left side falls to a trig substitution :
	y = sin(u)    ;     dy = cos(u) du
	(dy/(Sqrt[1-(y^2)]))=(( cos(u) du)/(Sqrt[1-(sin^2)(u)]))=((cos(u) du)/(cos(u)))=  du
Integrating the original equation gives
	u = ln(x) + C
	y(x) = sin(ln(x)+C)

Here's what it looks like for several values of C and x > 0.  
(It oscillates very fast near the origin, and isn't defined for x=0.)

Plot[
Table[Sin[Log[x]+C],{C,.5,2.5,.5}],
{x,0.01,6},AxesLabel->{x,y}]

A ball with mass 0.15 kg is thrown upward with initial velocity 20 m/sec from the roof of a building 20 m high.

The equation of motion is
	m(dv/dt)= - m g - AirFriction
where the friction term is always in the opposite direction from the velocity.

For each of the following 3 cases, find (i) the maximum height, (ii) the time when the ball hits the ground (assuming it doesn't hit the building), and (iii) plot height vs time.  (All three on the same graph, please.)

a) Air friction = 0

(dv/dt)=-g
	(dx/dt)=v = -g t + (Subscript[v, 0])
	x(t)=-(1/2)g (t^2)+(Subscript[v, 0])t + (Subscript[x, 0])
or
	((x(t))/(Subscript[x, 0]))=- (1/2)(((g (Subscript[x, 0]))/((Subscript[v, 0])^2)))((((t/((Subscript[x, 0])/(Subscript[v, 0])))))^2) + (t/(((Subscript[x, 0])/(Subscript[v, 0]))))+1
or
	w(τ) = - γ τ + 1
	s(τ) = - (1/2)γ(τ^2) + τ + 1
The maximum height happens when the velocity is 0, or τ = 1/γ = 2.04

At that time the height is
	(Subscript[s, max])=- (1/2)γ(((1/γ))^2) + (1/γ) + 1=1+(1/(2γ))= 2.02
	(Subscript[x, max])= (20 m) * 2.02 = 40.4 m
	
It hits the ground when s=0, which (using the quadratic formulat) is at
	(Subscript[τ, ground])=((1+(Sqrt[1+2γ]))/γ) = 4.9125
	(Subscript[t, ground])=(1 sec) (Subscript[τ, ground]) = 4.9125 sec
	
I've defined characterstic dimensions and dimensionless variables :
	(Subscript[t, 0])=(Subscript[x, 0])/(Subscript[v, 0])=(20 m)/(20 m/sec)=characteristic time
(Subscript[x, 0])=20m = initial distance
(Subscript[v, 0])=20 m/sec = initial velocity
τ = t/(Subscript[t, 0])=t/(1 sec) = dimensionless time
s = x/(Subscript[x, 0])=x/(20 m) = dimensionless distance
w = v/(Subscript[v, 0])= (ds/(dτ))= dimensionless velocity
γ = (((g (Subscript[x, 0]))/((Subscript[v, 0])^2))) = (g (((Subscript[t, 0])^2)/(Subscript[x, 0])))  = (9.8 m/(sec^2))((((1 sec))^2)/(20 m))=0.49 =dimensionless gravity
Note that the the original equation is
	(dw/(dτ))=-γ - (AirFriction)((((Subscript[t, 0])^2)/(m (Subscript[x, 0]))))

γ=0.49
pg1=Plot[-(1/2)γ (τ^2)+τ+1,{τ,0,5}]

b) Air friction = abs (v)/(m/sec) (1/30) Newtons

In dimensionless units, the friction term is
	Friction = | w | * (1/30) Newtons / ( (Subscript[m, 0])(Subscript[x, 0])/((Subscript[t, 0])^2)) = |w| α
	(Subscript[m, 0])= 0.15 kg = characteristic mass
	α = ((Newton/30))(((0.15kg)(20 m))/(((1 sec))^2))= 0.100
and the equation of motion is
	 (dw/(dτ))=-γ -α |w|
Separating variables and integrating gives us two equations, 
one for positive velocity and one for negative
w>0 : 
∫ (dw/(-γ - α w))=-((1/α))ln(-γ-α w)=τ + C
-γ-α w = K (e^(-α τ))
w(τ) = -(1/α)(γ + K (e^(-α τ)))
            = -(1/α)(γ + (-γ-α)(e^(-α τ)))        so that w=1 at τ=0
            = - (γ/α)+((γ/α)+1)(e^(-α τ))
	Does this have the right behavior for α->0 ?
		w ~ - (γ/α)+((γ/α)+1)(1-ατ + O((τ^2)))
    ~ 1 - γ τ -α τ + O((τ^2))                    which looks OK.

Solving for w=0 gives
	0 =  - (γ/α)+((γ/α)+1)(e^(-α τ))
	(Subscript[τ, w=0]) = -((1/α))ln(((γ/α)/((γ/α)+1))) = ((1/α))ln(((γ+α)/γ))

And the plot looks like this.

α = 0.1;
tauZero = (1/α) Log[(γ+α)/γ];
Plot[-(γ/α)+((γ/α)+1)Exp[-α τ],{τ,0,tauZero},AxesLabel->{τ,w}]

tauZero

To find the height, we need one more integral.
	(ds/(dτ))=w= - (γ/α)+((γ/α)+1)(e^(-α τ))
	s = ∫dτ {- (γ/α)+((γ/α)+1)(e^(-α τ))}=-(γ/α)τ +((1/(-α)))((γ/α)+1)(e^(-α τ)) +C
where the constant should be chosen so that s=1 at τ=0.
	1 = ((1/(-α)))((γ/α)+1) + C
	C = 1 + ((1/α))((γ/α)+1)
so
	s(τ) = -(γ/α)τ +((1/(-α)))((γ/α)+1)(e^(-α τ)) +1 + ((1/α))((γ/α)+1)        for 0 ≤ τ ≤ (Subscript[τ, w=0])

and the peak height is at

sLinearFric1[τ_]:= -(γ/α)τ +((1/(-α)))((γ/α)+1)Exp[-α τ] +1 + ((1/α))((γ/α)+1)

sLinearFric1[tauZero]

which in meters is

sLinearFric1[tauZero] (20 meters)

which gives us a plot like this :

pp1=Plot[sLinearFric1[τ],{τ,0,tauZero},AxesLabel->{τ,s}]

The equation when w < 0 is the same but with α->-α , and the height and time adjusted to match the first version at the peak.
w < 0
	(dw/(dτ))=-γ +α w
	∫ (dw/(-γ + α w))=((1/α))ln(-γ+α w)=τ + C
w(τ) = ((1/α))(γ + K (e^(α τ)))
            = ((1/α))(γ -γ (e^(α (τ - tauZero)))      so that w=0 at τ=tauZero
s(τ) = ∫ dτ { ((1/α))(γ -γ (e^(α (τ - tauZero))))}
          = (γ/α)τ-(γ/(α^2)) (e^(α (τ - tauZero)))+ C

with C chosen so that at tauZero this agrees with the earlier result.
	C =  sLinearFric1[tauZero] - ((γ/α)tauZero-(γ/(α^2)))

sLinearFric2[tau_]:=(γ/α)tau - (γ/(α^2))Exp[α(tau-tauZero)]+{ sLinearFric1[tauZero] - ((γ/α)tauZero-(γ/(α^2)))}

pp2=Plot[sLinearFric2[tau],{tau,tauZero,5}]

... which looks OK.

Finding when this hits the ground without resorting to numerical methods is getting a bit iffy though; we have to solve at equation with τ in both a linear and exponential term.

But we can always ask Mathematica for a numeric solution :

Solve[sLinearFric2[t]==0,t]

Plotting both of these positive and negative velocity cases gives

sLinearFric=
UnitStep[tauZero-tau] sLinearFric1[tau]+
UnitStep[tau-tauZero] sLinearFric2[tau];

pg2=Plot[sLinearFric,{tau,0,5}]

Here are the positive and negative velocity solutions superimposed.  The friction term is small, so they're fairly similar.

Plot[{sLinearFric1[t],sLinearFric2[t]},{t,0,5}]

And here are the first two solutions.  The one with friction doesn't go as high, and hits the ground sooner.

Show[pg1,pg2]

c) Air friction = (v^2)/(((m/sec))^2)(1/1325) Newtons

where a "Newton" is the unit of force in the MKS system,

I didn't get to this one; the approach is pretty much the same as (b), but with a different integral.

Epidemic

Assume a given population can be divided into sick and healthy, with y the proportion of the sick.  Since people getting infected are proportional to both the number of sick people (y) and the number of healthy people (1-y), the rate of infection can be modeled as
	(dy/dt)=α y (1-y)

where α is a constant.  Let (Subscript[y, 0])be the initial proportion of sickies.

a) Find the equilibrium points for the system, and determine whether each is stable, unstable, or semistable.

alpha = 4;
dydt[t_,y_]= alpha y (1 - y);
p1=directionField[dydt,-1,2,-1,2]

We find the stable points by setting (dy/dt)=0, which gives y=0 and y=1.
The direction plot shows these two "interesting" values for y; in either case, there is no further change.

To see whether each is stable or unstable, we look at the second derivative and its sign.
	(((d^2)y)/(dt^2))= α - 2 α y

At y=0, this is positive, which means its unstable.  The direction curves suggest that: as you move away towards higher t, the lines move away from y=0.

At y=1, this is negative, which means it is stable.  Again, the direction curves suggest that.

b) Solve the initial value problem and verify that the conclusions you reached in a) are correct.  Show that y(t)->1 as t-> ∞ in this model (i.e. everyone gets sick).

(dy/dt)=α y (1-y)
This is separable.
	∫ (dy/(y(1-y)))=∫α dt= α t + C
Now all we need to do is solve the integral on the left, which can be done with partial fractions.
	(1/(y(1-y)))= (a/y)+(b/(1-y))=((a(1-y) + b y)/(y(1-y)))=((a + y (b -a))/(y(1-y))) 
which works if a=b=1.  So
	∫ (dy/(y(1-y))) = ∫(dy/y)+∫(dy/(1-y))=ln(y) -ln(1-y)
and the whole solution is
	ln(y) - ln(1-y) = ln((y)/(1-y)) = α t + C
	((y(t))/(1-y(t))) = K (e^(α t)) 
where the initial condition implies that
	K = (Subscript[y, 0])/(1-(Subscript[y, 0]))
Multiplying through and solving for y(t) gives
	y(t) = (((Subscript[y, 0])(e^(α t)))/(1 - (Subscript[y, 0])+ (Subscript[y, 0])(e^(α t))))
For any initial (Subscript[y, 0]), in the limit as t->∞, this approaches 1.

alpha=4;
y0 = 0.02;
Plot[ 
((y0 Exp[alpha t])/(1 - y0+ y0 Exp[alpha t])),
{t,0,2.5},
PlotRange->{{0,2.5},{0,1}}]