Phys 263 Notebook for 17 May 12

The goal for the day is to work through as much of the homework as we can.
Your job is to evaluate each cell, think through why you expected the answer
Mathematica gave (or why it didn't do what you expected) and typically to play
with the cell, e.g. changing parameters.
You may find it useful to create your own notebook, copying cells out this one.
And of course you can come back later and finish reading through as you do the homework.

Morin 4.35

asks about the stadard issue 3 spring / 2 mass problem where the masses are driven
by an external force which just happens to go at 2ω with twice as much force on one of them

Here is the drawing function which created the above:

state[x_,y_]:= (L=12;Show[{ParametricPlot[{0+(t*x),Sin[2 Pi 6 t]},{t,0,1},Axes->True,PlotStyle->Red,PlotRange->{{0,L},{-1,1}}],Graphics[Rectangle[{x-1/4,-1/4},{x+1/4,1/4}]],
ParametricPlot[{x+t*(y-x),Sin[2 Pi 6 t]},{t,0,1},PlotStyle->Blue],
Graphics[Rectangle[{y-1/4,-1/4},{y+1/4,1/4}]],
ParametricPlot[{y+t*(L-y),Sin[2 Pi 6 t]},{t,0,1},PlotStyle->Orange]}])

Here are the equations of motion, including damping and forcing:

eq1 = x1''[t]+ 2 γ x1'[t] +2  ω0^2 x1[t] - ω0^2 x2[t] == F1 Cos[ω t]
eq2 = x2''[t]+ 2 γ x2'[t] +2  ω0^2 x2[t] - ω0^2 x1[t] == F2 Cos[ω t]

Pick some parameters:

p={ω0->1,γ->0,F1->1,F2->2,ω->2}

And ask for the solution (you might choose to Simplify):

generalsolution =DSolve[{eq1,eq2}/.p,{x1[t],x2[t]},t]

motion = {x1[t],x2[t]}/.generalsolution[[1]]//Simplify

Maybe it will look simpler when we pick boundary conditions, say starting at zero at rest:

s=DSolve[{eq1,eq2,x1'[0]==0,x2'[0]==0,x1[0]==0,x2[0]==0}/.p,{x1[t],x2[t]},t]//Simplify

Which part of this solution solution is the "particular" solution, i.e. the piece which comes from the driving
force and has to be there in every solution?  And which piece is the "homogeneous" solution, i.e. the modes
that would be there even without a driving force?
Your task: find boundary conditions which excite only the "particular" solution.
For the homework you want to work "by hand", taking as an ansatz a motion which wiggles
with frequency 2ω and seeing what works.
Can you think of a way to solve this system for a general force?  Try solving for the motion
when the force is the *same* on both particles, and when the forces are *opposite*.
What happens when the driving frequency is more general?
To explore, you may find it useful to just solve things numerically, as below:

s = NDSolve[{eq1,eq2,x1'[0]==0,x2'[0]==0,x1[0]==0,x2[0]==0}/.p,{x1[t],x2[t]},{t,0,100}]

Either way, you'll want to see the movie:

y1[t_] = (x1[t]/.s[[1]])
y2[t_] = (x2[t]/.s[[1]])
Animate[state[4+y1[t],8+y2[t]],{t,0,100}]

BTM 9.6.4

on page 276 asks about the quadratic form (i.e. quadratic function of the variables)

f[x_,y_] = x + y - x y - y^2

To find the critical point (there is only one):

eqnx = D[f[x,y],x]==0
eqny = D[f[x,y],y] ==0

{x0,y0} = {x,y}/.Solve[{eqnx,eqny}][[1]]

Note that means we can find the matrix of second derivatives at this point (actually, the matrix 
of second derivatives is the *same* for any point):

(hessian = D[f[x,y],{{x,y},2}])//MatrixForm

And we can write the form in this way:

g[x_,y_] =(1/2) {x-x0,y-y0}.hessian.{x-x0,y-y0}

Simplify[f[x,y]-g[x,y]]

The character of the critical point is determined by the Hessian:

Eigensystem[hessian]

Are the eigenvalues both positive (in which case we have a minimum), both negative
(max) or opposite in sign (saddlepoint)?

So now you won't be surprise to see the following:

Plot3D[f[x,y],{x,-2,0},{y,0,2}]

BTM 9.6.5 (unassigned, but still good to know about)

BTM 9.6.5 asks us to solve the eigenvalue problem for a general 2x2 hermitian matrix, and verify
that the determinant is in fact the product of eigenvalues, and the trace is in fact the sum of eigenvalues.
Defining a general matrix

(matrix = {{α,β+I γ},{β-I γ,δ}})//MatrixForm

we can ask for the eigenvalues, i.e. the roots of some quadratic

{λ1,λ2} =Eigenvalues[matrix]

Check out the sum and product

λ1 + λ2//Simplify

Tr[matrix]

λ1*λ2 //Simplify

Det[matrix]

Thus if a matrix has negative determinant we know the product of eigenvalues is negative,
i.e. one is positive and one is negative, and of the matrix is a Hessian, we have a saddle point.
To distinguish between maxima and minima (both of which have positive det) we have to look 
at the trace: positive trace means both eigenvalues are positive, i.e. a local minimum, and negative trace means max.

BTM 9.6.6

ask you to use the insight of BTM 9.6.5 as you find and analyze the critical points of the functions

f[x_,y_] = (x-1)3 - 3 x - 4 y^2

f[x_,y_] = x^4 + y^4 - 2 x^2(1-y^2)+2y^2

To get you going, here is the function in Figure 9.4 (page 276):

f[x_,y_] = x y (4 + x) - (2x+y)^2

p1=ContourPlot[f[x,y],{x,-4,4},{y,-3,6}]

We can look at this in 3D as well.

Plot3D[f[x,y],{x,-4,4},{y,-3,6}]

To find the critical points (now there are several)

eqx = D[f[x,y],x]==0
eqy = D[f[x,y],y]==0

s=Solve[{eqx,eqy}]

Make a list of the points and then plot them as red dots:

criticalpoints = {x,y}/.%

p2=ListPlot[criticalpoints,PlotStyle->{PointSize[.03],Red}]

Show[p1,p2]

To add infcrmation on the eigendirections, we construct the Hessian at each point:

(h=D[f[x,y],{{x,y},2}])//MatrixForm

The direction vectors at one of the critical points are:

pt = {x,y}/.s[[2]]

eig=Eigenvectors[h/.s[[2]]]

To draw the vectors we need more functions:

Needs["VectorFieldPlots`"]

p3=ListVectorFieldPlot[{{pt,eig[[1]]},{pt,eig[[2]]}}]

Show[p1,p2,p3]

Your job: add p4 and p5, i.e. vectors for the other two points,
and then produce analagous plots for the other target functions.

BTM 9.7.1

asks you to check that these functions are properly normalized, i.e. their absolute squares
integrate to 1 over the domain (here the interval from 0 to L).

ϕ[n_,x_] = Sqrt[2/L]Sin[n π x/L]

The problem only asks about n=1 and n=2, but we can be general (and so can you, by
writing the sin(θ) as (1/2I)(Exp[I θ]-Exp[-I θ]) ):

Integrate[ϕ[n,x]^2,{x,0,L}]

Simplify[%,Assumptions->Element[n,Integers]]

Likewise we check orthogonality for modes 1 and 2

Integrate[ϕ[1,x]*ϕ[2,x],{x,0,L}]

Or in general:

Integrate[ϕ[n,x]*ϕ[m,x],{x,0,L}]

Simplify[%,Assumptions->{Element[m,Integers],Element[n,Integers]}]

(Looks like Mathematica also assumes m not equal to n.)

BTM 9.7.2

but more often than not we choose the basis of periodic functions:

ϕ[n_,x_] = (1/Sqrt[L]) Exp[2 π I n x/L]

Check that we normalized them correctly:

Integrate[ϕ[n,x]*ϕ[-n,x],{x,0,L}]

Your task: verify they are orthogonal for general integers m and n, just as we did in the section above.

Fourier Series

Choose an interval, e.g. -1 to 1:

x1=-1;x2=1;L=x2-x1;

Define normalized periodic basis functions:

ϕ[n_,x_] = Exp[2 π I n x/L]/Sqrt[L]

Check that they are orthogonal etc.:

overlap = Integrate[ϕ[1,x]*ϕ[2,x],{x,x1,x2}]

overlap = Integrate[ϕ[n,x]*ϕ[-m,x],{x,x1,x2}]

overlap/.{m->3,n->1}

Limit[overlap,m->n]

Pick a function (evaluate one or the other of these)

f[x_] = If[x<0,1+x,1-x]

f[x_] = x^2

Plot[f[x],{x,x1,x2}]

Compute the Fourier series:

a[n_] = Integrate[f[x]*ϕ[-n,x],{x,x1,x2}]

Note that the n=0 mode has to be treated as a limit:

Limit[a[n],n->0]

a[0]=Limit[a[n],n->0]

Simplify[a[n],Assumptions->{Element[n,Integers]}]

Define the n-th approximation by summing the modes from -n to n:

g[n_,x_] := Sum[a[j]ϕ[j,x],{j,-n,n}]

Plot[{f[x],g[5,x]},{x,x1,x2},PlotStyle->{Red,Blue}]

The value of at x=1 is1 :

f[1]

On the other hand the Fourier series evaluated at that point is the sum of terms:

table=Table[a[n]ϕ[n,1]+a[-n]ϕ[-n,1],{n,0,5}]

Note we are double counting the n = 0, but we can fix it with:

table[[1]]=table[[1]]/2

The sum of all these numbers is f(1), i.e. 1, or equivalently, the sum of the following is π^2/4:

table*Pi^2/4

Since π^2/4 - π^2/6 = π^2/12, we've just computed the following sum:

Sum[1/n^2,{n,1,Infinity}]

BTM 9.7.3 etc.

Asks for the same game to be applied to:

Clear[h,L]

f[x_] = If[x<L/2,2x h/ L,2 h (L-x)/L]

Before plotting, take a look at f(x) in your mind's eye: what does it look like for x<L/2?
what is the value f(L/2)?  what does it look like for x>L/2?   Now peek:

p={h->1,L->1}

Plot[f[x]/.p,{x,0,L}/.p]

Plot[f[x]/.p,{x,0,L/.p}]

And go compute the coefficients:

ϕ[n_,x_] = (1/Sqrt[L]) Exp[2 π I n x/L]

fm[m_] = Integrate[f[x]*ϕ[-m,x],{x,0,L},Assumptions->L>0]

Note we can simplify this further when m's are known to be integers:

fm[m_] = Simplify[fm[m],Assumptions->Element[m,Integers]]

Note that this expression still isn't right for m=0, so we fill in the gap by hand:

fm[0] = Integrate[f[x]*ϕ[0,x],{x,0,L},Assumptions->L>0]

And here is how to express the partial sums:

g[n_,x_] := Sum[fm[m]*ϕ[m,x],{m,-n,n}]

And here is the approximation with just the three contributions from m=-1,0,1.
Plot and then repeat with more modes.

g[1,x]

Plot[{f[x],g[1,x]}/.p,{x,0,L/.p}]

Finally, noting that f[L/2] = h, we learn how to sum another infinite series, just along
the lines of Sum[1/n^2] above.

BTM 9.7.5

Drive an LRC circuit with this v[t]:

p = {R->60, L->1, c->2/10^6,V0->200,T->1/100}

v[t_] = (V0/T Mod[t,T]/.p)

Plot[v[t],{t,0,3T/.p}]

First ask for the numerical solution:

s=NDSolve[{L Q''[t]+R Q'[t] + Q[t]/c ==v[t],Q[0]==0,Q'[0]==0}/.p,Q[t],{t,0,.4}][[1]]

q[t_] = Q[t]/.s

And take a look.  q[t] itself looks kind of normal, but the current q'[t] has an
interesting hiccup at the time the applied voltage jumps:

Plot[q[t],{t,.3,.32}]

Plot[q'[t],{t,.3,.32}]

To solve by Fourier series, resolve the applied voltage in a sum of exp()'s.
Then recall that for one pure frequency ω, i.e. v=vn Exp[I ωn t]  the solution for q(t) is:

qn[t_] = vn/ (-ωn^2 L + I ωn R + 1/c) Exp[I ωn t]

Your task: resolve v[t] in to a sum Sum[vn Exp[I ωn t]], i.e. find the Fourier coeffiecients vn.
Then q[t] can be found as Sum[qn[t]].  Taking one time derivative we then get the current,
which according to the back of BTM should be:

curr[m_,t_] := Re[Sum[100 I Exp[200 n π I t]/(60 n π+I(200 n^2 π^2 - 2500)),{n,1,m}]
+ Sum[100 I Exp[200 (-n) π I t]/(60 (-n) π+I(200 n^2 π^2 - 2500)),{n,1,m}]]

Note that after a dozen or so modes even that little kink starts to become apparent:

Plot[{curr[12,t],q'[t]},{t,.2,.22},PlotStyle->{Blue,Red}]

Solution: Expanding as usual in

ϕ[n_,t_] = 1/Sqrt[T] Exp[2 π I n t/T]

We paused to check the orthonormality of the basis:

overlap=Integrate[ϕ[-n,t]*ϕ[m,t],{t,0,T}]

Simplify[overlap,Assumptions->Element[{m,n},Integers]]

Limit[overlap,m->n]

Then we expand:

f[m_] = Integrate[ϕ[-m,t]* V0 t/T,{t,0,T},Assumptions->{Element[m,Integers],T>0}]

Note this gives a bogus result for m=0, when we want

f[0] = Integrate[ϕ[0,t]* V0 t/T,{t,0,T}]

Now we can solve for the current:  the vn are the f's we've just computed, while ωn is 2 π n /T:

q[n_,t_] =f[n] / (-(2π n /T)^2 L + I(2π  n/T) R + 1/c) Exp[I (2π n t /T)]

Note this general expression is messed up for n=0 again.  But that n=0 mode drops out when we
take a derivative to find the current:

qndot =D[q[n,t],t]

qndot/.p