Lab 16: Vectors and Vector Products

MAT 229, Spring 2021

Vector review

A vector is an object that has magnitude and direction.

Vectors in the plane

Can be represented in component form (Cartesian coordinates) as (Overscript[v, ->])=⟨a,b⟩.

Can be represented in magnitude, angle form (polar coordinates) as |(Overscript[v, ->])|=(Sqrt[(a^2)+(b^2)]),  and tan(θ)=(b/a) where θ is the angle between the positive x-axis and the vector.

Vectors in space

Can be represented in component form (Cartesian coordinates) as (Overscript[v, ->])=⟨a,b,c⟩.

Mathematica

Mathematica can work with vectors. Execute the following commands as you encounter them.

Vector components

To write a vector in component form use curly braces to delimit the components. For example to define the vector (Overscript[v, ->])=⟨1,2.45,-3/17⟩ in Mathematica, enter

v={1,2.45,-3/17}

Scalar multiplication and vector addition

These two arithmetic operations in Mathematica are just as you would expect.

2v

v+{1,2,3}

Dot product review

There are two ways to view the dot product of two vectors:

If (Overscript[u, ->])=⟨a,b,c⟩ and (Overscript[v, ->])=⟨α,β,γ⟩, then (Overscript[u, ->])·(Overscript[v, ->])=a×α+b×β+c×γ.

If θ is the angle between (Overscript[u, ->]) and (Overscript[v, ->]), then (Overscript[u, ->])·(Overscript[v, ->])=|(Overscript[u, ->])||(Overscript[v, ->])|cos(θ).

Dot product

To compute the dot product of two vectors in Mathematica, just place a period between them. For example, to compute (Overscript[v, ->])·(Overscript[v, ->]) where (Overscript[v, ->]) is defined above, enter

v.v

Questions

What is the angle between ⟨1,2⟩ and ⟨-5,3⟩?

Are the two vectors ⟨1,0,-1⟩ and ⟨2,-3,2⟩ perpendicular or not?

Vector length

To compute the length of a vector in Mathematica use the fact that (Overscript[v, ->])·(Overscript[v, ->])=((|(Overscript[v, ->])|)^2). For example, to compute the length of the vector defined above enter

Sqrt[v.v]

Unit vectors

Definition

A unit vector is a vector with magnitude 1. Its magnitude is fixed, but its direction can be any direction.

Questions

Let (Overscript[u, ->])=⟨1,2⟩

Is (Overscript[u, ->]) a unit vector?

Find a scalar λ such that λ (Overscript[u, ->]) is a unit vector.

Find a general method for taking a vector and finding a unit vector pointing in the same direction.

Does this work for all vectors?

Questions

What are the unit vectors parallel to the tangent line to the curve y=2sin(x) at the point (π/6,1)?

What are the unit vectors perpendicular to the tangent line to the curve y=2sin(x) at x=π/6?

Application: Kepler’s First Law:

A planet revolves around the sun in an elliptical orbit with the sun at one focus.

Andy Long, Spring, 2019

This is from Chapter 13 of Stewart (end of section 13.4): to understand why planets travel in ellipses (Kepler’s first law)

I chose the following at (somewhat) random: you want a little “eccentricity” to your ellipse (not too much like a circle). h is related to the angular momentum. c is the length of the constant vector we obtain in the derivation, and determines the eccentricity -- that is, how much the ellipse varies from a circle (c=0).

h=1;
c=10;
eccentricity=.6; (* Equals c/(GM) in reality -- so dependent as well on universal gravitation and the mass of the sun. *)
d=h^2/c;
r[θ_]:=eccentricity d/(1+eccentricity Cos[θ])
x[θ_]:=r[θ]Cos[θ]
y[θ_]:=r[θ]Sin[θ]
p1=ParametricPlot[{x[θ],y[θ]},{θ,0,2π}];
p2 = ListPlot[{{0,0}},PlotMarkers->{""}];
Show[p1,p2]

Most of the sausage-making is done in the initialization cells found hidden here-below....

If you plot r(θ) as θ changes uniformly, you’ll see that the “planet” is moving at the wrong speeds: it should be faster when it’s near the sun, and slower when it’s far away. That’s because the angular change in the plot is constant, which is not the case for planetary motion:

The velocity vector runs along the curve traced out by r:

th0=2Pi/3;
dydx[t_]:=y'[t]/x'[t]
tangentLine[t_]:=y[th0]+dydx[th0](x[t]-x[th0])
p3=ParametricPlot[{x[t],tangentLine[t]},{t,th0-Pi/12,th0+Pi/12},PlotLabel->StringJoin["θ =",ToString[th0,InputForm]],
PlotRange->{{-.2,.1},{-.1,.1}}];
p4=ParametricPlot[t{x[th0],y[th0]},{t,0,1}];
p5=ParametricPlot[
t(
r'[th0]/r[th0]{x[th0],y[th0]}+{-y[th0],x[th0]}
)/r[th0]^2
+
{x[th0],y[th0]},
{t,0,1},PlotStyle->Orange];
Show[p3,p1,p4,p5]

To get the speed right, we would need to solve a “differential equation” to find θ as a function of time t -- that is, θ(t). (Make sure that you take DEs one day!):

From (UnderBar[h])=(r^2)((OverHat[u]) x (OverHat[u])’):

Since the length of (OverHat[u]) can’t change, the only change possible is in direction -- perpendicular to the direction of (OverHat[u]) (i.e., increasing or decreasing θ, with no change in the length “r” -- which is one in this case). Hence, computing the norm of (UnderBar[h]), we get

h=(r^2)| (OverHat[u])’ | .

The change  (OverHat[u])’ will be ((r dθ)/dt), but the length r is 1 (this is a unit vector). Therefore
 | (OverHat[u])’ | = ((dθ)/dt)
Finally,
 ((dθ)/dt)=(h/(r^2))=(h/((r[θ])^2))
 (which I couldn’t solve analytically). It’s separable, so  ∫((r[θ])^2)dθ =  h∫dt = h t + C.
 
 By the way, this proves Kepler’s equal areas  formula, too: if we turn these into definite integrals, then it says for a given time span T,
 
 (∫_0^T)((r[θ(t)])^2)dθ(t) =  h(∫_0^T)dt = h T 
 
 which says the area traced out is ((h T)/2) -- equal areas are traced out in equal times. 
 
 I tried using Mathematica to compute this integral, but it just spun:
 
Integrate[(eccentricity d/(1+eccentricity Cos[th]))^2, {th,0,θ}]

So we’ll do it numerically (approximately): the end time was chosen to give us values of θ from 0 to 2π (a complete period of the planet’s orbit) :

end=.0442;(* this works for the parameter values I chose! *)
s=NDSolve[{θ'[t]==h/r[θ[t]]^2,θ[0]==0},θ,{t,0,end}];

What we obtain is a really good numerical solution for θ(t). 

Now we want to plot the orbit as a function of time, including the sun’s location (indicated by the dodecahedron), and the planet’s position (indicated by the icosahedron). These markers are chosen in honor of Kepler, i.e. 
https://upload.wikimedia.org/wikipedia/commons/1/19/Kepler-solar-system-1.png

Manipulate[
Module[{point=θ[t]/.s},
point=point[[1]];
Show[
p1,
p2,
ListPlot[{{r[point]Cos[point],r[point]Sin[point]}},PlotMarkers->{}]
]
],{t,0,0.044}
]

Now I’ll add in a position vector, velocity, and acceleration, and we’ll have it all!

Manipulate[
Module[{rvector,vvector,secondVvector, avector, 
angle,th0=θ[t]/.s,dthdt,
norma, normv,normr,normsecond,
p4, p5, p6,p7},
(* θ is a function of time, which is created in the initialization cells. So time will increase steadily, constantly,but θ will not: the rate of change of θ with respect to time is not constant. *)
th0=th0[[1]];
rvector={x[th0],y[th0]}; (* The position vector *)
normr=Norm[rvector]; (* The length of the position vector -- distance from the sun *)
dthdt=h/(normr^2);
(* Chain rule and derivatives of vectors: v = dr/dt = dr/dθ dθ/dt *)
vvector={x'[th0],y'[th0]} dthdt;
normv=Norm[vvector];(* The length of the velocity vector -- speed of the planet *)
avector=
{x''[th0],y''[th0]}dthdt^2
+
vvector (-2  Sin[th0]/d)/normr;(* The acceleration vector *)
norma=Norm[avector];(* The length of the acceleration vector, proportional to the force of the sun on the planet *)
angle=ArcCos[Dot[vvector,rvector]/normv/normr]; (* just used as a label *)
p5=ParametricPlot[t vvector+rvector,{t,0,.005},PlotStyle->{Red,Thick},
PlotRange->{{-.25,.1},{-.1,.1}},PlotLabel->angle];
p6=ParametricPlot[t avector+rvector,{t,0,.00003},PlotStyle->{Black,Thick}];
p4=ParametricPlot[t rvector,{t,0,1}];
(* Looks to me like I could replace {x[th0],y[th0]} in the following with rvector: *)
secondVvector=r'[th0]rvector/r[th0]^3+{-y[th0],x[th0]}/r[th0]^2;
normsecond=Norm[secondVvector];
p7=ParametricPlot[t secondVvector+{x[th0],y[th0]},{t,0,.005},PlotStyle->{Orange,Thick}];
Show[p5,p1,p2,p4,p6,
ListPlot[{{r[th0]Cos[th0],r[th0]Sin[th0]}},PlotMarkers->{}]
]
],
{t,0,end}
]

Questions:

At which points in the motion is the position vector perpendicular to the velocity vector?

Where is the size of the acceleration an extreme?

In whose honor do I choose the dodecahedron and icosahedron for the point masses?