Physics 229& 100 Homework #6
Name: Nasser Abbasi

1. Propagation of waves through a thin film in 1 dimension. Waves propagating to the right are represented by ^( k x) terms and waves propagating to the left are represented by ^(-  k x)  terms. k is the magnitude of the wave vector which is (2π)/λ= ( ω)/c , where ω is the angular frequency , and  c is the speed of the wave in the external media. Inside the film, the wave speed is slower by a factor of n ( in optics, this is the index of refraction) , so the k vector kf=(n ω)/c. ω is determined by the source of the waves and is the same inside and outside the film.  

[Graphics:HTMLFiles/index_6.gif]

a) On the left hand side, the incident wave has unit amplitude, and the reflected wave has amplitude R , so the total field is ^( k x)+ R ^(- k x). On the right hand side, there is only the transmitted wave with amplitude T. Inside the film, waves can go both directions so the field is some combination of left and right going waves, as indicated in the figure. The film has thickness d.
The reflection and transmission coefficients are determined by the boundary conditions at the interfaces of the film, which are that the field amplitude is continuous and the derivative of the field is continuous; the same type of logic holds for acoustic waves, light waves, and quantum waves. Calculate the  amplitude transmission coefficient T in terms of k, kf and d.

Remove["Global`*"]

fieldLeft = ^( k x) + R ^(- k x) ;

fieldCenter = A ^( k_f x) + B ^(- k_f x) ;

fieldRight = T ^( k x) ;

Now set the 4 equations using the boundary conditions given

eq1 = (fieldLeft == fieldCenter)/.x→0

1 + R == A + B

eq2 = (fieldCenter == fieldRight)/.x→d

B ^(- d k_f) + A ^( d k_f) == ^( d k) T

eq3 = (D[fieldLeft, x] == D[fieldCenter, x])/.x→0//Simplify

k R + A k_f == k + B k_f

eq4 = (D[fieldCenter, x] == D[fieldRight, x])/.x→d//Simplify

^(- d k_f) (-B + A ^(2  d k_f)) k_f == ^( d k) k T

So we have now 4 equations, here they are

{eq1, eq2, eq3, eq4}//TableForm

1+R==A+B
B ^(- d k_f) + A ^( d k_f) == ^( d k) T
k R + A k_f == k + B k_f
^(- d k_f) (-B + A ^(2  d k_f)) k_f == ^( d k) k T

Now I need to solve the above to find T

eq5 = Flatten[eq3/.Solve[{eq1, eq2}, R]]

{k (-1 + A - A ^(2  d k_f) + ^( d k +  d k_f) T) + A k_f == k + B k_f}

sol = Flatten[Solve[{eq4, eq5[[1]]}, {A, B}]]

eq6 = eq2/.sol ;

sol2 = Flatten[Solve[eq6, T]] ;

myT = T/.sol2 //FullSimplify ;

T = myT

-(4 ^(- d (k - k_f)) k k_f)/(^(2  d k_f) (k - k_f)^2 - (k + k_f)^2)

b) T is the amplitude transmission coefficient; the intensity of the transmitted wave or the probability of transmission is measured by the energy transmission coefficient TE= T T*.  Find a compact expression for this quantity in terms of k, kf and d.

T_e = T ComplexExpand[Conjugate[T]]//FullSimplify

-(16 k^2 k_f^2)/(2 Cos[2 d k_f] (k^2 - k_f^2)^2 - 2 (k^4 + 6 k^2 k_f^2 + k_f^4))

c) The k vector inside the film kf is related to the k vector in the external media by kf=n k, where n is the optical index, or the ratio of the wave speeds. Express TE in terms of k , d and n.

T_e = T_e/.k_f→n k

-(16 k^4 n^2)/(-2 (k^4 + 6 k^4 n^2 + k^4 n^4) + 2 (k^2 - k^2 n^2)^2 Cos[2 d k n])

c) TE and n are dimensionless ratios. The only  terms in the above expression that have dimensions are k and d, but they appear in the dimensionless combination k d = (2π d)/λ  . Let d→dr λ , where dr is a "reduced"  film thickness measured in units of the wavelength,  and express TE in terms of n and dr

T_e = T_e/.k→ (2π)/λ ;

T_e = T_e/.d→d_r λ //Simplify

(8 n^2)/(1 + 6 n^2 + n^4 - (-1 + n^2)^2 Cos[4 n π d_r])

d) Let n→1.5 ( typical for glass) and make a plot as a function of dr. Where do the maxima and minima occur? Comment on the physical significance of the plot.

T_e = T_e/.n->3/2

18/(313/16 - 25/16 Cos[6 π d_r])

Plot[T_e, {d_r, 0, 8}, PlotRange→All, AxesLabel→ {d , "Intensity of transmitted wave"}] ; r

[Graphics:HTMLFiles/index_45.gif]

2.  Using the discussion of polar coordinates as a guide, compute the expression for the laplacian in spherical coordinates. This will require expressing , for example, ∂/∂x = ∂r/∂x∂/∂r+∂θ/∂x∂/∂θ+∂φ/∂x ∂/∂φ . Use makeop to make a ddx operator .  Compare your answer with the expression obtained from the VectorAnalysis package using Laplacian.

Remove["Global`*"] ;

Needs["Calculus`VectorAnalysis`"]

Needs["MMHTools`DForm`"]

(j = JacobianMatrix[Spherical[r, θ, φ]])//MatrixForm

(jT = Transpose[j])//MatrixForm

(jInvT = Inverse[jT])//MatrixForm

Now each row in the inverse of the transpose of the jacobian matrix represent partial derivative of x or y or z w.r.t. to the spherical coordinates. Use makeop to make partial operators as follows

ddxop = makeop[jInvT[[1]], {r, θ, φ}] ;

ddyop = makeop[jInvT[[2]], {r, θ, φ}] ;

ddzop = makeop[jInvT[[3]], {r, θ, φ}] ;

Now that we have have constructued the first partial derivatives operators, we apply the operator again to obtain the second partial derivatives, and then add them to make the Laplacian

myAnswer = ddxop[ddxop[g[r, θ, φ]]] + ddyop[ddyop[g[r, θ, φ]]] + ddzop[ddzop[g[r, θ, φ]]]//Simplify ;

DForm[myAnswer]

(2 r (∂g)/(∂r) + Cot[θ] (∂g)/(∂θ) + r^2 ∂^2g/(∂r^2) + ∂^2g/(∂θ^2) + Csc[θ]^2 ∂^2g/(∂φ^2))/r^2

Now find the Laplacian in spherical using Mathematica and then compare the answers

(mmaAnswer = Laplacian[g[r, θ, φ], Spherical[r, θ, φ]])//DForm

myAnswer - mmaAnswer//Simplify

0

The above verifies the answer is the same.

3. Construct the matrix which maps the components of a Cartesian vector into the components of a vector in Spherical coordinates. Use JacobianMatrix and ScaleFactors from the VectorAnalysis package. Use the mapping to convert the Cartesian vector  Overscript[z,^] = {0,0,1}  into spherical coordinates.

Remove["Global`*"] ;

Needs["Calculus`VectorAnalysis`"]

Needs["MMHTools`DForm`"]

(j = JacobianMatrix[Spherical[r, θ, φ]])//MatrixForm

sc = ScaleFactors[Spherical[r, θ, φ]]

{1, r, r Sin[θ]}

(sc = DiagonalMatrix[sc])//MatrixForm

( {{1, 0, 0}, {0, r, 0}, {0, 0, r Sin[θ]}} )

(A = j . Inverse[sc])//MatrixForm

( {{Cos[φ] Sin[θ], Cos[θ] Cos[φ], -Sin[φ]}, {Sin[θ] Sin[φ], Cos[θ] Sin[φ], Cos[φ]}, {Cos[θ], -Sin[θ], 0}} )

{rCoord, θCoord, φCoord} = Transpose[A/.θ→ 0] . {0, 0, 1}

{1, 0, 0}

Now Verify above answer using Mathematica

{rCoord, θCoord, φCoord} = CoordinatesFromCartesian[{0, 0, 1}, Spherical]

{1, 0, 0}

4 a)  In Spherical coordinates, calculate ∇^2( r ^n ), and ∇^2 ( r^( 2) Y_2^1(θ,φ))  (see SphericalHarmonicY in Help).

Remove["Global`*"] ;

Needs["Calculus`VectorAnalysis`"]

Needs["MMHTools`DForm`"]

Laplacian[r^n, Spherical[r, θ, φ]]//DForm

n (1 + n) r^(-2 + n)

Laplacian[r^2 SphericalHarmonicY[1, 2, θ, φ], Spherical[r, θ, φ]]//DForm

0

b) In Spherical coordinates, calculate ∇·(r ^n Overscript[r,^]) , where Overscript[r,^] is the unit vector in the r direction.

rhat = {1, 0, 0}

{1, 0, 0}

v = r^n rhat

{r^n, 0, 0}

Div[v, Spherical[r, θ, φ]]

(2 + n) r^(-1 + n)

c) In Cylindrical coordinates, calculate ∇×Overscript[φ,^] , where Overscript[φ,^] is the unit vector in the φ direction.

tauHat = {0, 1, 0}

{0, 1, 0}

Curl[tauHat, Cylindrical[r, θ, z]]

{0, 0, 1/r}

5.  Prove that the shortest path between 2 points in the plane is a straight line.  The length of a curve defined by y[x] is ∫ (1 + y '[x]^2)^(1/2) dx . Find the (2nd order) differential equation that determines the minimum value of the integral. Solve it and show that it represents a line. Also set each of the conserved quantities to a constant, and show that the solution of the resulting (1st order) differential equation also gives a straight line.

Remove["Global`*"]

F[x, y[x], y '[x]] := (1 + (y '[x])^2)^(1/2) ;

Write Euler equation which is d/dx(dF/dy ')-dF/dy

ode = ∂_x (∂_y '[x] F[x, y[x], y '[x]]) - ∂_y[x] F[x, y[x], y '[x]] //Simplify

y^′′[x]/(1 + y^′[x]^2)^(3/2)

DSolve[ode == 0, y[x], x]

{{y[x] →C[1] + x C[2]}}

We see from the above that the solution is a stright line.

First find what is dF/dy

D[ F[x, y[x], y '[x]], y[x] ]

0

Since dF/dy=0, then this means d/dx(dF/dy ')=0 or in otherwords, dF/dy 'is constant, say k. hence conserved.

We solve this first order ODE

ode = D  [ F[x, y[x], y '[x]], y '[x] ] == k

y^′[x]/(1 + y^′[x]^2)^(1/2) == k

DSolve[ode, y[x], x]

{{y[x] → -(k x)/(1 - k^2)^(1/2) + C[1]}, {y[x] → (k x)/(1 - k^2)^(1/2) + C[1]}}

Hence we see from the above that the solution y(x) also gives straight line

ode = 1 + (y '[x])^2 == k

1 + y^′[x]^2 == k

DSolve[ode, y[x], x]

{{y[x] → -(-1 + k)^(1/2) x + C[1]}, {y[x] → (-1 + k)^(1/2) x + C[1]}}

We see the the solution also gives a stright line.
Hence we have 2 conserved quantities, and each generate first order ODE and each give solution as constant line.

6. A mechanics problem that is found in most text books is the double pendulum, i.e. 2 massless rods of lengths L1 and L2  with masses m1 and m2 at the lower ends arranged so that L1 is attached from the ceiling and L2 is attached at the end of L1. Their positions are described in terms of angular deviation from the vertical of each rod, θ1[t] and θ2[t].

a) Draw a picture of this situation, including the relevant variables, and insert the drawing into this notebook. PowerPoint is probably the easiest to use; PowerPoint drawings can be simply copied and pasted into Mathematica. Any other drawing program can be used, but you will have to convert vector graphics to bitmaps. Alt PrintScreen will make a bitmap image of the monitor which can be edited and cropped in Paint.

[Graphics:HTMLFiles/index_143.gif]

b) Assume the double pendulum moves in a plane. Explicitly write the 2D position vectors p1 and p2 of the two masses. Find the velocity vectors v1 and v2. Find expressions for the kinetic energy KE and the gravitational potential energy PE.

Remove["Global`*"]

First write the coordinates of each bob

x1 = L1 Sin[θ1[t]] ;

y1 = -L1 Cos[θ1[t]] ;

x2 = L1 Sin[θ1[t]] + L2 Sin[θ2[t]] ;

y2 = -L1 Cos[θ1[t]] - L2 Cos[θ2[t]] ;

Now write the vectors for each mass position

p1 = {x1, y1}

{L1 Sin[θ1[t]], -L1 Cos[θ1[t]]}

p2 = {x2, y2}

{L1 Sin[θ1[t]] + L2 Sin[θ2[t]], -L1 Cos[θ1[t]] - L2 Cos[θ2[t]]}

v1 = D[p1, t]

{L1 Cos[θ1[t]] θ1^′[t], L1 Sin[θ1[t]] θ1^′[t]}

v2 = D[p2, t]

{L1 Cos[θ1[t]] θ1^′[t] + L2 Cos[θ2[t]] θ2^′[t], L1 Sin[θ1[t]] θ1^′[t] + L2 Sin[θ2[t]] θ2^′[t]}

KE = 1/2m1  v1 . v1 + 1/2m2  v2 . v2//Simplify

1/2 (L1^2 (m1 + m2) θ1^′[t]^2 + 2 L1 L2 m2 Cos[θ1[t] - θ2[t]] θ1^′[t] θ2^′[t] + L2^2 m2 θ2^′[t]^2)

(*PE = m1 g L1 Cos[θ1[t]] + m2 g (L1 Cos[θ1[t]] + L2 Cos[θ2[t]]) This is wrong, correct below*)

-g (L1 (m1 + m2) Cos[θ1[t]] + L2 m2 Cos[θ2[t]])

PE is not right

PE = m1 g y1 + m2 g y2//Simplify

-g (L1 (m1 + m2) Cos[θ1[t]] + L2 m2 Cos[θ2[t]])

c) Write the Lagrangian and construct the equations of motion.

L = KE - PE//Simplify

To find equation of motion, I can use Mathematica EulerEquations function

Needs["Calculus`VariationalMethods`"]

deps = {θ1[t], θ2[t]} ;

eqs = EulerEquations[L, deps, t] ;

Flatten[Solve[eqs, {θ1''[t], θ2''[t]}]]//Simplify//TableForm

7.  Symmetric and antisymmetric matrices can be diagonalized, i.e. if M is a symmetric matrix, there exists a matrix C such that for M= RandomSymMat[n_], C M C^T=diagonal.
C^T is the transpose of C. The rows of C are the normalized eigenvectors of M. This transformation amounts to a change in basis. The fact that diagonal matrices are easy to exponentiate provides another way to exponentiate symmetric matrices. Exp[M t]=C^T[diag] C, where the diagonal elements of diag are { e^(t λ_1),e^(t λ_2),..}, where λ_i are the eigenvalues of M.

a) Write a function RandomSymmetricMatrix[n_] which returns a random symmetric n x n matrix.

Remove["Global`*"] ;

Needs["LinearAlgebra`MatrixManipulation`"]

randomSymMat[n_] := Module[{m}, m = UpperDiagonalMatrix[Random[] &, n] ; m = m + Transpose[m] ] ;

Test the above on small matrix

randomSymMat[3]//MatrixForm

( {{0.124405, 0.920806, 0.179341}, {0.920806, 1.34487, 0.151654}, {0.179341, 0.151654, 1.06672}} )

b) Write a function Diag[mat_] that computes C and prints it in MatrixForm and returns the diagonal matrix C M C^T.  Use Dimensions if you need to know how big mat is. Remember that the Mathematica Eigensystems command does not return normalized eigenvectors. You also may want to use N to get your results in decimal form and Chop to get rid of round off trash. Test your function on several matrices generated by RandomSymmetricMatrix.

Notice: Mathematica eigenvectors returned seem to be normalized allready.

Test the function few times

Do[Print["test for n=", n] ; Print["Diag=", MatrixForm[diag[randomSymMat[n]]]], {n, 2, 4}]

test for n=2

EigenVectors C= ( {{0.552954, 0.833212}, {0.833212, -0.552954}} )

Diag= ( {{2.20064, 0}, {0, 0.424421}} )

test for n=3

EigenVectors C= ( {{-0.861761, -0.325978, -0.388722}, {-0.283218, -0.326584, 0.901738}, {0.420897, -0.887176, -0.189114}} )

Diag= ( {{2.33614, 0, 0}, {0, -0.358034, 0}, {0, 0, 0.322564}} )

test for n=4

Diag= ( {{3.0476, 0, 0, 0}, {0, 1.41189, 0, 0}, {0, 0, 1.14919, 0}, {0, 0, 0, 0.477802}} )

c) Write a function that computes Exp[M t]=C^T[diag] C. Show that it gives the same result as the MatrixExp command for a random symmetric matrix. Make sure your answer contains the variable t.

myMatrixExp[mat_, var_] := Module[{c, diag, λ}, {λ, c} = Eigensystem[mat] ; diag = DiagonalMatrix[^(λ var)] ; Transpose[c] . diag . c]

(mat = randomSymMat[3])//MatrixForm

( {{0.714617, 0.320444, 0.956927}, {0.320444, 1.55777, 0.851825}, {0.956927, 0.851825, 1.912}} )

(myMatrixExp[mat, t])//MatrixForm

(MatrixExp[mat  t])//MatrixForm

The result is the same.

8. Chemical mixing problem ( from Robertson). Tank A and tank B are connected by a series of pipes.  At time t=0, the flow rates  shown in the figure are established. The volume of each tank is as shown. Each tank also contains an initial mass of salt, mA0, and mB0. Note that the flow rates into and out of each tank are the same, so the volume of liquid in each tank is fixed, but the salt concentration varies with time.


concentration varies with time.
[Graphics:HTMLFiles/index_202.gif]

a) Denote the mass of salt in each tank as mA[t], mB[t], respectively. Write the system of differential equations which describe the rates of change of these quantities.

Remove["Global`*"]

r1 = 1 ; r2 = 2 ; r3 = 3 ; r4 = 4 ; r1 = 1 ; (*liter/min*)

Vb = 40 ;

Va = 20 ;

Null

Now for each tank, write the rate of the salt concentration in and out

(*tank A*)rateIn = r3 mB[t]/Vb ; (*notice, pure water has zero salt concentration, so not part of equation*)rateOut = r4 mA[t]/Va + r2 mA[t]/Va ; eq1 = mA '[t] == rateIn - rateOut

mA^′[t] == -(3 mA[t])/10 + (3 mB[t])/40

(*tank B*)rateIn = r4 mA[t]/Va ; rateOut = r3 mB[t]/Vb + r1 mB[t]/Vb ; eq2 = mB '[t] == rateIn - rateOut

mB^′[t] == mA[t]/5 - mB[t]/10

b)Find the general solution  using MatrixExp.

Needs["LinearAlgebra`MatrixManipulation`"]

{coeffMatrix, ddt} = LinearEquationsToMatrices[{eq1, eq2}, {mA[t], mB[t]}] ;

Print[ MatrixForm[ddt], "=", MatrixForm[coeffMatrix], MatrixForm[{mA[t], mB[t]}]]

( {{-mA^′[t]}, {-mB^′[t]}} ) = ( {{3/10, -3/40}, {-1/5, 1/10}} ) ( {{mA[t]}, {mB[t]}} )

{mA, mB} = MatrixExp[-coeffMatrix t] . {mA0, mB0} ; mA

mB

c) use the solution obtained above to find the specific solution for the case mA0=1000, mB0=0. Plot mA[t] and mB[t] on the same plot.

Needs["Graphics`Colors`"]

[Graphics:HTMLFiles/index_222.gif]

d) One of the useful features of the MatrixExp approach is that the same solution can be used for different initial conditions. Repeat the above  for mA0=0, mB0=1000

[Graphics:HTMLFiles/index_224.gif]

9. A standard representation for the delta function is
δ(x-a)=Underscript[lim, L→∞] 1/(2π)Underoverscript[∫, -L, arg3]e^(i(x - a) t)dt
Evaluate the integral for L=3, 10, 30. Convince yourself that the answer is real, and that the area under all the curves is 1.  Plot all three results on a single plot. Set a->1 and plot from {-1,2}. Conclude that the formula is in fact reasonable.

Remove["Global`*"] ;

Needs["Graphics`Colors`"]

del[x_, a_] := 1/(2 π ) Integrate[^( (x - a) t), {t, -L, L}] ;

Remove :: rmnsm : There are no symbols matching \"Global`*\". More…

sol = del[x, a]

Sin[L (a - x)]/(π (a - x))

{sol3, sol10, sol30} = sol/.L→ {3, 10, 30}

{Sin[3 (a - x)]/(π (a - x)), Sin[10 (a - x)]/(π (a - x)), Sin[30 (a - x)]/(π (a - x))}

We see that for real 'a' and 'x', the answer is real and not complex. To show that the area under the curves is 1, we integrate.

Integrate[sol3, {x, -∞, ∞}]

Integrate[sol10, {x, -∞, ∞}]

1

1

[Graphics:HTMLFiles/index_242.gif]

10. (grads)  Consider a string whose displacement y[x,t] is a function of x and t. The density of the string is ρ[x]. Let T be the tension in the string, which is constant and provided by a mass .

[Graphics:HTMLFiles/index_243.gif]

a) What is the density of kinetic energy, i.e. what function would you integrate over the length of the string to get its total kinetic energy?

b) What is the density of potential energy? A good way to think of this is to consider the constant tension supplied by a mass m hanging over a pulley. Assume that the mass of the string is negligible compared to m.

c) What is the Lagrangian density, and the equations of motion. Show that if the slope of the string is small, the equation is a wave equation. Use DForm to make it look like conventional mathematics.

Created by Mathematica  (January 30, 2007) Valid XHTML 1.1!