Physics 229& 100 Homework #5
Name: Nasser Abbasi

1. Suppose  that the profit your company makes is given by P=s m (3m+s)E^(-(2s +m)/40) , where s and m are the number of sales and manufacturing people, respectively. If the total number of people is equal to 50, how do you distribute the human resources to maximize your profit.  By how much would profit increase if you hire an additional person?  Do the problem two ways:
a) using Lagrange multipliers and FindRoot

The problem formulation is the following
Maximize the objective function P (which represents profit) subject to contraint  g(x,m)=50, where in this problem g(x,m) is defined as x+m

If we consider a 3D space, where the z axis represents the profit, the x axis is s, and y axis is m, then we want to locate the point above the line s=50-m (marked as RED in the diagram, which is the constraint) where this point has the largest distance from the x-y plane, as shown in this diagram
[Graphics:HTMLFiles/index_1.gif]

Remove["Global`*"] ;

maxNumberOfPeople = 50 ;

RowBox[{RowBox[{f[s_, m_], :=, RowBox[{s,  , m,  , (3 m + s),  , ^(-(2 s + m)/40), Cell[]}]}], ;}]

g[s_, m_] := s + m - maxNumberOfPeople ;

faux[s_, m_, λ_] := f[s, m] - λ g[s, m] ;

eq1 = {D[faux[s, m, λ], m] == 0, D[faux[s, m, λ], s] == 0, D[faux[s, m, λ], λ] == 0} ;

MatrixForm[%]

So we have 3 equations to solve for 3 unknowns,(s,m,λ).Since this is not a polynomials,we can't use Solve or NSolve,we need to find FindRoot[],but for that we need a guess.So lets make a plot of the profit to get a guess for s and m

Plot3D[f[s, m], {s, 0, 50}, {m, 0, 50}, PlotLabel->"Profit vs m and s", AxesLabel→ {"s", "m", "profit"}] ;

[Graphics:HTMLFiles/index_11.gif]

We see from the above that guess for s=15and m=35 looks promissing.Lets use that (what about λ?,I tried many different values,and it did not make a difference)

sol = FindRoot[eq1, {{s, 15}, {m, 35}, {λ, 1}}]

{s→14.375, m→35.625, λ→345.415}

Hence, since the solution must be an integer, we round the numbers

s = Round[s/.sol]

14

m = Round[m/.sol]

36

b) Using NMaximize

Clear[s, m] ;

NMaximize[{f[s, m], g[s, m] == 0}, {s, m}]

{12419.4, {m→35.625, s→14.375}}

We see that the answer agrees with the method of the lagrange multiplier

2. Find the triangle of perimeter L with maximum area.  A triangle is specified by three points in a plane. Put the first point at the origin, the second point along the x axis at {x1,0} , and the third point at {x2,y2}.
a) Write an expression for the area of the triangle in terms of x1,x2,y2 , neglecting the perimeter constraint.

Remove["Global`*"] ;

A = 1/2x1  y2

(x1 y2)/2

b) What is the maximal area in terms of L?

g = x1 + (x2^2 + y2^2)^(1/2) + ((x1 - x2)^2 + y2^2)^(1/2) - L

-L + x1 + ((x1 - x2)^2 + y2^2)^(1/2) + (x2^2 + y2^2)^(1/2)

Aaux = A - λ g

(x1 y2)/2 - (-L + x1 + ((x1 - x2)^2 + y2^2)^(1/2) + (x2^2 + y2^2)^(1/2)) λ

eq1 = {D[Aaux, x1] == 0, D[Aaux, x2] == 0, D[Aaux, y2] == 0, D[Aaux, λ] == 0} ;

MatrixForm[eq1]

sol = Solve[eq1, {x1, x2, y2, λ}] ;

MatrixForm[sol]

x1 = x1/.sol

{0, 0, L/3, L/3}

y2 = y2/.sol

{0, 0, -L/(2 3^(1/2)), L/(2 3^(1/2))}

areas = 1/2x1 y2

{0, 0, -L^2/(12 3^(1/2)), L^2/(12 3^(1/2))}

Hence maximum area in terms of L is the last case, which is

areas[[4]]

L^2/(12 3^(1/2))

3.  Find the minimum distance between points  in the intersection of
3x-4y+5z==8       and      x+y-z==11  
and points on the unit sphere. Do it 2 ways, and show that you get the same answer both ways:
a) using PseudoInverse and geometry ( utilizing the fact that the PseudoInverse solution to the simultaneous equations will have minimum length)

The constraint is the line of the interection. The objective function is the distance between this line and points on the unit sphere.

Remove["Global`*"] ;

Needs["Graphics`Shapes`"] ;

Needs["Graphics`Colors`"] ;

plan1[x_, y_] := 1/5 (8 + 4 y - 3 x) ;

plan2[x_, y_] := 11 - y - x ;

p1 = Plot3D[plan1[x, y], {x, -4, 16}, {y, -4, 16}, PlotPoints→50, ViewPoint→ {-2, -2.4, 2.}, DisplayFunction→Identity] ;

p2 = Plot3D[plan2[x, y], {x, -4, 16}, {y, -4, 16}, PlotPoints→30, ViewPoint→ {-2, -2.4, 2.}, DisplayFunction→Identity] ;

Show[{p1, p2}, DisplayFunction→$DisplayFunction] ;

[Graphics:HTMLFiles/index_50.gif]

sphere = x^2 + y^2 + z^2 - r^2 ;

r = 5 ; (*I enlarge the sphere just to be able to see it*)p0 = WireFrame[Graphics3D[Sphere[r, 15, 15], DisplayFunction→Identity]] ;

Show[{p1, p2, p0}, DisplayFunction→$DisplayFunction] ;

[Graphics:HTMLFiles/index_54.gif]

Now use PsudoInverse to find the point on the line of the intersection of the 2 planes where the point is closes to the origin

Clear[x, y, z] ;

Needs["LinearAlgebra`MatrixManipulation`"] ;

eq1 = 3 x - 4 y + 5 z == 8 ;

eq2 = x + y - z == 11 ;

 {A, b} = LinearEquationsToMatrices[{eq1, eq2}, {x, y, z}] ; MatrixForm[A]

( {{3, -4, 5}, {1, 1, -1}} )

MatrixForm[b]

( {{8}, {11}} )

{x, y, z} = PseudoInverse[A] . b//N

{7.61404, 2.08772, -1.29825}

distaceFromOrigin = Sqrt[x^2 + y^2 + z^2]

8.0011

The above is the distance to the origin. Hence the distance to the surface of the sphere is

minDistaceFromSphereSurface = % - 1

7.0011

b) using Lagrange multipliers.

Let xc,yc,zc represent coordinates of a point on the sphere. Let xp,yp,zp represent coordinates of a point on the surafces. Now setup the object function and the constrains and the auxilary function.

Remove["Global`*"] ;

r = 1 ;

f = (xc - xp)^2 + (yc - yp)^2 + (zc - zp)^2

(xc - xp)^2 + (yc - yp)^2 + (zc - zp)^2

g1 = 3 xp - 4 yp + 5 zp - 8

-8 + 3 xp - 4 yp + 5 zp

g2 = xp + yp - zp - 11

-11 + xp + yp - zp

g3 = xc^2 + yc^2 + zc^2 - r^2

-1 + xc^2 + yc^2 + zc^2

fAux = f - (λ1 g1 + λ2 g2 + λ3  g3)

(xc - xp)^2 + (yc - yp)^2 + (zc - zp)^2 - (-8 + 3 xp - 4 yp + 5 zp) λ1 - (-11 + xp + yp - zp) λ2 - (-1 + xc^2 + yc^2 + zc^2) λ3

eq = Thread[{D[fAux, xc], D[fAux, yc], D[fAux, zc], D[fAux, xp], D[fAux, yp], D[fAux, zp]} == 0]

MatrixForm[eq]

maxeqs = Join[eq, {g1 == 0, g2 == 0, g3 == 0}] ;

MatrixForm[maxeqs]

sol = NSolve[maxeqs, {xc, yc, zc, xp, yp, zp, λ1, λ2, λ3}] ;

MatrixForm[%]

9.0011

7.0011

Min[minDistance1, minDistance2]

7.0011

We see that the solution agrees with part(a).

4 a) The Wronskian of a set of functions {y_1[x], y_2[x]..., y_n[x] } is defined as
det(

y_1 y_2 ... y_n
y_1^′ y_2^′ .. y_n^′
.... ...... ...
y_1^(n - 1) y_2^(n - 1) .. y_n^(n - 1)
). Write a function wronski[funcs_List, var_] that will compute the Wronskian of any set of functions. Use it to compute the Wronskian of {^x,^(-x), Sin[x]} and of {^x,^(-x), Sinh[x]}

Remove["Global`*"] ;

wronski[funcs_List, var_] := Det[NestList[D[#, var] &, funcs, Length[funcs] - 1]] ;

wronski[{E^x, E^-x, Sin[x]}, x]

4 Sin[x]

wronski[{E^x, E^-x, Sinh[x]}, x]

0

b) One of the uses of the Wronskian is to determine if a set of functions is linearly independent. Recall that two vectors v1 and v2 are independent if A v1+B v2=0   implies that A and B are zero. Similarly,  if the only solution to A y1[x]+B y2[x]=0   ( for all x) is the trivial one {A->0,B->0}, y1 and y2 are independent functions. Note that if   A y1[x]+B y2[x]=0, then by taking a derivative, we conclude that  A y1' [x]+B y2' [x]  must also vanish.  Show that the condition for independence is that the Wronskian of y1 and y2 does not vanish  for all x.

sketch of proof:
Assume that y_1and y_2 are L.d. and show that this implies that W=0
Then assume that W=0 and show that this implies that  y_1and y_2 are L.d.
Hence if W /= 0, then this must imply that y_1and y_2 are L.I.

Start of proof:
Assume y_1and y_2 are L.d then we can write
                               y_1= k_1y_2        
where k_1 is some constant

Then
                                       y ' _1= k_1y ' _2  
     
Now W=(

y_1 y_2
y_1 ' y_2 '
)

             =(
k_1 y_2 y_2
k_1 y ' _2 y_2 '
)


Hence the determinant of W is

                                  k_1y_2y_2'-k_1y_2y_2' =  0
                                  
Now for the second part. Assume that det(W)=0, then this means

                                y_1y_2'-y_2y_1' = 0                           -----------------(1)
                                
But differentiation itself is a linear operation, hence we can write
                                     y ' _1=k_1 y_1           ------(2)
and
                                     y ' _2=k_2 y_2            -----(3)
Where k_1and k_2 are some arbitrary constanst which are unrelated.
Substitute the above 2 equations into equation (1) we obtain

                                 y_1k_2y_2-y_2k_1 y_1 = 0                           
Divide by y_2 we obtain
                                      y_1k_2=k_1 y_1        
which means k_1=k_2 (else we have y_1=0 everywhere)
Since k_1=k_2, then if we divide equation (2) by (3) we obtain

                                            y ' _1/y ' _2=y_1/y_2
Which implies that y_1 is L.d. to y_2. This complete the second part of the proof.

Hence if det(W)/=0, then y_1 must be L.I. to y_2

This completes the proof.                                            

c) Consider a second order DE with coefficients that are not necessarily constant:
y''[x]+p1[x] y'[x]+p0[x] y[x]=0.  Let y1[x] and y2[x] be two solutions to this equation. Let W be the Wronskian of the two solutions. Show that dW/dx=-p1[x] W[x]. If p1[x] is constant, show that if W is non zero for any x, it will be non zero for all x.

I will just write everything without [x] next to it for ease of reading.

W=det (

y1 y2
y1' y2'
)=y1 y2'-y2 y1'

Hence dW/dx=y1'  y2'  + y1  y2''- y2'  y1'- y2  y1''   ------(1)

Since y1[x] is a solution, then
                           y1''=-p1 y1' -p0 y1
Similarly, since y2[x] is a solution then
                           y2''=-p1 y2' -p0 y2
Substitute the above 2 equations into (1) we obtain

dW/dx=y1'  y2'  + y1  (-p1 y2' -p0 y2)- y2'  y1'- y2  (-p1 y1' -p0 y1)

dW/dx=y1'  y2'  - y1  p1 y2' -y1 p0 y2- y2'  y1'+ y2  p1 y1' + y2 p0 y1

dW/dx=y1'  y2'  - y1  p1 y2' - y2'  y1'+ y2  p1 y1'

dW/dx=- p1 (y1 y2' - y2 y1' )

But the last term above is W, hence

dW/dx=- p1[x]  W[x]

Which is what we are asked to show.          

Last part. If p1[x] is constant, say m then the above is dW/dx=mW[x] which has the solution of W(x)=n ^mx where n is some other constant.
Hence we see that if we can find one x where W(x) is not zero, then this implies that the constant n is not zero (since the exponential term is not zero unless x=-∞). So, if n≠0 and since the exponential term is not zero then W(x) can not be zero for any other x.

5. Write a function waveEq[f_] which computes the wave operator ∇^2f-1/c^2∂^2f/∂t^2 in spherical coordinates f[r,θ,φ,t]  ( use the VectorAnalysis Add-on package). The wave operator is identically zero if f satisfies the wave equation. Use it to show that ^( k r -  ω t)/ris a solution for certain values of the scalar k.

Remove["Global`*"] ;

Needs["Calculus`VectorAnalysis`"]

Needs["MMHTools`DForm`"]

SetCoordinates[Spherical[r, θ, φ]]

Spherical[r, θ, φ]

waveEq[f_] := Laplacian[f] - 1/c^2D[f, {t, 2}]

sol = waveEq[f[r, θ, φ, t]]

DForm[sol]

g = ^( k r -  ω t)/r ;

sol = waveEq[g]

sol = Simplify[sol]

(^( (k r - t ω)) (-c^2 k^2 + ω^2))/(c^2 r)

Reduce[sol == 0, k]

c≠0&& ((k == ω/c&&r≠0) || (k == -ω/c&&r≠0))

We see that the solution is zero form specific values of k, which are +- frequency/speed of light

6.  A plane wave polarized in the x direction and propagating in the z direction has an electric field E1={ex, 0, 0}^(k . x - ω t), where ex is the magnitude of the field and  k is a constant vector k={0, 0, km}. According to Maxwell's equations, this wave has a magnetic field associated with it which obeys the equation ∇×E=-∂B/∂t.
a) Assume B={bx, by, bz} ^(k . x - ω t), and determine bx, by, and bz.

Remove["Global`*"] ;

Needs["Calculus`VectorAnalysis`"]

SetCoordinates[Cartesian[x, y, z]]

Cartesian[x, y, z]

k = {0, 0, km} ;

xVec = {x, y, z} ;

E1 = {ex, 0, 0} ^( (k . xVec - ω t))

{^( (km z - t ω)) ex, 0, 0}

B = {bx, by, bz} ^( (k . xVec - ω t))

{bx ^( (km z - t ω)), by ^( (km z - t ω)), bz ^( (km z - t ω))}

eq = Curl[E1] == - D[B, t]

bSolution = Flatten[Solve[ComplexExpand[eq], {bx, by, bz}]]

{bx→0, by→ (ex km)/ω, bz→0}

b) Compute the Poynting vector S= (E^* ×B)/μ_0.Verify that it is proportional to the magnitude of E1 squared.

s = ComplexExpand[Cross[Conjugate[E1], B]/μ_0]

{0, -(bz ex)/μ_0, (by ex)/μ_0}

s = s/.bSolution

{0, 0, (ex^2 km)/(ω μ_0)}

magS = Sqrt[Total[s^2]]

(ex^4 km^2)/(ω^2 μ_0^2)^(1/2)

magS = Simplify[magS, Assumptions→ {Element[{km, ω, μ_0, ex}, Reals]}]

(ex^2 Abs[km])/(Abs[ω] Abs[μ_0])

(*Is it ok to assume km≥0 ?, this is just so not to have Abs[] pop up*)magS = Simplify[magS, Assumptions→ {km≥ 0, ω≥ 0, μ_0≥ 0}]

(ex^2 km)/(ω μ_0)

magE1 = Sqrt[Total[ComplexExpand[Re[E1]]^2 + ComplexExpand[Im[E1]]^2]] ;

magE1Squred = FullSimplify[magE1]^2

ex^2

propertionalAmount = magS/magE1Squred

km/(ω μ_0)

7. The divergence theorem says that Underscript[∫, S] F·dS = Underscript[∫, V] ∇·F dV . Verify the theorem by computing both sides of the identity for the cartesian vector field  F= {z x y^2, x^2 y, y x z^2}  , where V is the interior of a cylinder centered at the origin of radius 1 and extending from z=-1 to z=1. Calculate ∇·F  in Cartesian coordinates, and then perform a change of variables and do the integral in  cylindrical coordinates. Express dS in cartesian coordinates, and then compute the dot product.

Remove["Global`*"] ;

Needs["Calculus`VectorAnalysis`"]

Needs["MMHTools`DForm`"]


Working on the lhs. We need to integrate F over the suraface area of the cylinders. There are 3 integrals. One for the top surface, one for the bottom surface, and one integral for the surroundaing surface.


[Graphics:HTMLFiles/index_253.gif]

To help me understand how to use Mathematica to perform coordinates transformation, I made the following diagram

[Graphics:HTMLFiles/index_254.gif]

For each integral we write Overscript[F, -].Overscript[n, -]da where Overscript[n, -] is a unit vector perpendiculare to the surface.

For to the top surafce, using unit vectors expressed in cylinderical coordinates,  Overscript[n, -]=Overscript[e_z, -]={0,0,1}, while for the bottom surface  Overscript[n, -] = -Overscript[e_z, -]={0,0,-1} and for the surrounding surface  Overscript[n, -]=Overscript[e_r, -]={cosθ,sinθ,0}

Now before doing the dot product Overscript[F, -].Overscript[n, -] we need to convert  Overscript[F, -]  to cylindrical coordinates, which we can do using the transformation x=rcosθ,y=rsinθ,z=z

Now we need to find the scalar quantity da. For the top and botton surfaces, da=r dr dθ, and for the surrounding surface, da=dθ dz

Start by converting vector F coordinates from cartesian to cylinderical

F = {z x y^2, x^2 y, y x z^2} ;

F2 = F/.{x→r Cos[θ], y→ r Sin[θ]}

{r^3 z Cos[θ] Sin[θ]^2, r^3 Cos[θ]^2 Sin[θ], r^2 z^2 Cos[θ] Sin[θ]}

Now get the cylinderical coordinates of the unit normal vector to the top surface

n = {0, 0, 1} ;

dot = Dot[F2, n]

r^2 z^2 Cos[θ] Sin[θ]

I1 = Integrate[dot  r, {r, 0, 1}, {θ, 0, 2π}]

0

Now do the bottom surface

n = {0, 0, -1} ;

dot = Dot[F2, n]

-r^2 z^2 Cos[θ] Sin[θ]

I2 = Integrate[dot  r, {r, 0, 1}, {θ, 0, 2π}]

0

Now do the surrounding surface

n = { Cos[θ], Sin[θ], 0} ;

dot = Dot[F2, n]

r^3 Cos[θ]^2 Sin[θ]^2 + r^3 z Cos[θ]^2 Sin[θ]^2

I3 = Integrate[(dot ), {z, -1, 1}, {θ, 0, 2π}]

(π r^3)/2

I3 = I3/.r→1 (*since radius is 1*)

π/2

lhs = I1 + I2 + I3

π/2

Now do the rhs

div = Div[F, Cartesian[x, y, z]]

x^2 + 2 x y z + y^2 z

div = div/.{x→r Cos[θ], y→ r Sin[θ]}

r^2 Cos[θ]^2 + 2 r^2 z Cos[θ] Sin[θ] + r^2 z Sin[θ]^2

rhs = Integrate[(div r), {r, 0, 1}, {z, -1, 1}, {θ, 0, 2π}]

π/2

lhs == rhs

True

8. Verify the theorem from vector calculus that ∇·(A×B)=B·(∇×A)-A·(∇×B)  for vector functions A and B by computing both sides and showing that they are equal. Use the VectorAnalysis package.
a) Do the calculation using Cartesian coordinates, i.e. A &B are functions of  x,y,z.

Remove["Global`*"] ;

Needs["Calculus`VectorAnalysis`"]

Needs["MMHTools`DForm`"]

A = {Ax[x, y, z], Ay[x, y, z], Az[x, y, z]}

B = {Bx[x, y, z], By[x, y, z], Bz[x, y, z]}

Make 2 general vector functions of coordinates x,y,z to work with.

Ax[x_, y_, z_] := {a_1 x, a_2y, a_3z} ;

Bx[x_, y_, z_] := {b_1 x, b_2y, b_3z} ;

Cross[Ax[x, y, z], Bx[x, y, z]]

{-y z a_3 b_2 + y z a_2 b_3, x z a_3 b_1 - x z a_1 b_3, -x y a_2 b_1 + x y a_1 b_2}

lhs = Div[%, Cartesian[x, y, z]]

0

curlA = Curl[Ax[x, y, z], Cartesian[x, y, z]]

{0, 0, 0}

curlB = Curl[Bx[x, y, z], Cartesian[x, y, z]]

{0, 0, 0}

rhs = Dot[Bx[x, y, z], curlA] - Dot[Ax[x, y, z], curlB]

0

lhs == rhs

True

This is not a general case

b) Do the calculation using cylindrical coordinates, i.e. A&B are functions of r, θ,z.

Ar[r_, θ_, z_] := Ax[x, y, z]/.{x→r Cos[θ], y→r Sin[θ], z→z}

Br[r_, θ_, z_] := Bx[x, y, z]/.{x→r Cos[θ], y→r Sin[θ], z→z}

Cross[Ar[r, θ, z], Br[r, θ, z]]

{-r z Sin[θ] a_3 b_2 + r z Sin[θ] a_2 b_3, r z Cos[θ] a_3 b_1 - r z Cos[θ] a_1 b_3, -r^2 Cos[θ] Sin[θ] a_2 b_1 + r^2 Cos[θ] Sin[θ] a_1 b_2}

lhs = Div[%, Cylindrical[r, θ, z]]//FullSimplify

z Sin[θ] (-a_3 (b_1 + 2 b_2) + (a_1 + 2 a_2) b_3)

curlA = Curl[Ar[r, θ, z], Cylindrical[r, θ, z]]

{0, 0, (r Sin[θ] a_1 + 2 r Sin[θ] a_2)/r}

curlB = Curl[Br[r, θ, z], Cylindrical[r, θ, z]]

{0, 0, (r Sin[θ] b_1 + 2 r Sin[θ] b_2)/r}

rhs = Dot[Br[r, θ, z], curlA] - Dot[Ar[r, θ, z], curlB]//FullSimplify

z Sin[θ] (-a_3 (b_1 + 2 b_2) + (a_1 + 2 a_2) b_3)

lhs == rhs

True

Suggestion: try to use somtthing like  Ax[x,y,z] or Ar[r,φ,z]


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