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

Remove["Global`*"];
maxNumberOfPeople  =  
        50;
f[s_, m_]:=s m (3 m  +  s) (E^((((-(2 s + m))/40))));
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"}];\)\)

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}}]

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

s = Round[s /. sol]

m = Round[m /. sol]

b) Using NMaximize

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

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/2) x1  y2

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

g = 
      x1 + (Sqrt[((x2^2) + (y2^2))]) + (Sqrt[(((x1 - x2))^2) + (y2^2)]) - 
        L

Aaux = A - λ g

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

y2 = y2 /. sol

areas = (1/2) x1 y2

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

areas[[4]]

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];

sphere = (x^2) + (y^2) + (z^2) - (r^2);
r = 
      5;    

p0 = 
      WireFrame[
        Graphics3D[Sphere[r, 15, 15], 
          DisplayFunction -> Identity]];
Show[{p1, p2, p0}, 
        DisplayFunction -> $DisplayFunction];

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]

MatrixForm[b]

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

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

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

minDistaceFromSphereSurface  =  % - 1

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)

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

g2 = xp + yp - zp - 11

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

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

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[%]

sol1 = 
        sol[[1, All]];

    xc = xc /. sol1;  yc = yc /. sol1; zc = zc /. sol1;

    xp = xp /. sol1;  yp = yp /. sol1; 
    zp = zp /. sol1;


    minDistance1 = 
      Sqrt[(((xc - xp))^2) + (((yc - yp))^2) + (((zc - zp))^2)] // 
        N

sol2 = 
        sol[[4, All]];

    xc = xc /. sol2;  yc = yc /. sol2; zc = zc /. sol2;

    xp = xp /. sol2;  yp = yp /. sol2; 
    zp = zp /. sol2;


    minDistance2 = 
      Sqrt[(((xc - xp))^2) + (((yc - yp))^2) + (((zc - zp))^2)] // 
        N

Min[minDistance1, minDistance2]

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

4 a) The Wronskian of a set of functions {(Subscript[y, 1])[x], (Subscript[y, 2])[x]..., (Subscript[y, n])[x] } is defined as 
det(({{(Subscript[y, 1]), (Subscript[y, 2]) ... , (Subscript[y, n])}, {((Subscript[y, 1])^′), ((Subscript[y, 2])^′) .. , ((Subscript[y, n])^′)}, {...., ... ..., ...}, {((Subscript[y, 1])^((n - 1))), ((Subscript[y, 2])^((n - 1))) .. , ((Subscript[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 {(e^x),(e^(-x)), Sin[x]} and of {(e^x),(e^(-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]

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

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 (Subscript[y, 1])and (Subscript[y, 2]) are L.d. and show that this implies that W=0
Then assume that W=0 and show that this implies that  (Subscript[y, 1])and (Subscript[y, 2]) are L.d.
Hence if W /= 0, then this must imply that (Subscript[y, 1])and (Subscript[y, 2]) are L.I.

Start of proof:
Assume (Subscript[y, 1])and (Subscript[y, 2]) are L.d then we can write
                               (Subscript[y, 1])= (Subscript[k, 1])(Subscript[y, 2])        
where (Subscript[k, 1]) is some constant

Then
     \(TraditionalForm\`\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\(\(y'\)\_1\)\)\)= (Subscript[k, 1])(Subscript[(y'), 2])  
     
Now W=(({{(Subscript[y, 1]), (Subscript[y, 2])}, {(Subscript[y, 1])', (Subscript[y, 2])'}}))

             =(({{(Subscript[k, 1])(Subscript[y, 2]) , (Subscript[y, 2])}, {(Subscript[k, 1])(Subscript[(y'), 2]), (Subscript[y, 2])'}}))

Hence the determinant of W is 

                                  (Subscript[k, 1])(Subscript[y, 2])(Subscript[y, 2])'-(Subscript[k, 1])(Subscript[y, 2])(Subscript[y, 2])' =  0
                                  
Now for the second part. Assume that det(W)=0, then this means

                                (Subscript[y, 1]) (Subscript[y, 2])'-(Subscript[y, 2])(Subscript[y, 1])' = 0                           -----------------(1)
                                
But differentiation itself is a linear operation, hence we can write 
                                     (Subscript[(y'), 1])=(Subscript[k, 1]) (Subscript[y, 1])           ------(2)
and
                                     (Subscript[(y'), 2])=(Subscript[k, 2]) (Subscript[y, 2])            -----(3)
Where (Subscript[k, 1])and (Subscript[k, 2]) are some arbitrary constanst which are unrelated.
Substitute the above 2 equations into equation (1) we obtain

                                 (Subscript[y, 1]) (Subscript[k, 2]) (Subscript[y, 2])-(Subscript[y, 2])(Subscript[k, 1]) (Subscript[y, 1]) = 0                           
Divide by (Subscript[y, 2]) we obtain
                                      (Subscript[y, 1]) (Subscript[k, 2])=(Subscript[k, 1]) (Subscript[y, 1])        
which means (Subscript[k, 1]) = (Subscript[k, 2]) (else we have (Subscript[y, 1]) = 0 everywhere)
Since (Subscript[k, 1]) = (Subscript[k, 2]), then if we divide equation (2) by (3) we obtain

                                            ((Subscript[(y'), 1])/(Subscript[(y'), 2]))=((Subscript[y, 1])/(Subscript[y, 2]))
Which implies that (Subscript[y, 1]) is L.d. to (Subscript[y, 2]). This complete the second part of the proof.

Hence if det(W)/=0, then (Subscript[y, 1]) must be L.I. to (Subscript[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 (e^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 (∇^2)f-(1/(c^2))(∂^2) (f/∂)(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 ((e^(i k r - i ω t))/r)is a solution for certain values of the scalar k.

Remove["Global`*"];
Needs["Calculus`VectorAnalysis`"]
Needs["MMHTools`DForm`"]

SetCoordinates[Spherical[r, θ, ϕ]]

waveEq[f_] := 
      Laplacian[f] - (1/(c^2)) D[f, {t, 2}]
sol = waveEq[f[r, θ, ϕ, t]]

DForm[sol]

g = ((E^(I k r  -  I ω t))/r);
sol = waveEq[g]

sol = Simplify[sol]

Reduce[sol == 0, k]

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}(e^(i(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} (e^(i(k . x - ω t))), and determine bx, by, and bz.

Remove["Global`*"];
Needs["Calculus`VectorAnalysis`"]
SetCoordinates[Cartesian[x, y, z]]

k = {0, 0, km};
xVec = {x, y, z};
E1 = {ex, 0, 
          0} (E^(I (k . xVec - ω t)))

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

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

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

b) Compute the Poynting vector S= ((\(E\^*\)×B)/((Subscript[μ, 0]))).Verify that it is proportional to the magnitude of E1 squared.

s = ComplexExpand[Cross[Conjugate[E1], B]/(Subscript[μ, 0])]

s = s /. bSolution

magS = Sqrt[Total[(s^2)]]

magS = 
      Simplify[magS, 
        Assumptions -> {Element[{km, ω, (Subscript[μ, 0]), ex}, 
              Reals]}]

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

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

propertionalAmount = (magS/magE1Squred)

7. The divergence theorem says that (∫_S) F·dS = \(TraditionalForm\`∫\+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.




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



For each integral we write \(TraditionalForm\`∫\(F\&-\) . \(n\&-\) da\) where \(TraditionalForm\`\(n\&-\)\) is a unit vector perpendiculare to the surface.

For to the top surafce, using unit vectors expressed in cylinderical coordinates,  \(TraditionalForm\`\(n\&-\)\)=\(TraditionalForm\`\(\(e\_z\)\&-\)\)={0,0,1}, while for the bottom surface  \(TraditionalForm\`\(n\&-\)\) = -\(TraditionalForm\`\(\(e\_z\)\&-\)\)={0,0,-1} and for the surrounding surface  \(n\&-\)=\(\(\(e\_r\)\&-\) = {cosθ, sinθ, 0}\)

Now before doing the dot product \(TraditionalForm\`\(F\&-\) . \(n\&-\)\) we need to convert  \(TraditionalForm\`\(\(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[θ]}

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

n = {0, 0, 1};
dot = Dot[F2, n]

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

Now do the bottom surface

n = {0, 0, -1};
dot = Dot[F2, n]

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

Now do the surrounding surface

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

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

I3 = 
      I3 /. r->1

lhs =  I1 + I2 + I3

Now do the rhs

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

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

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

lhs == rhs

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_] := {(Subscript[a, 1]) x, (Subscript[a, 2]) y, (Subscript[a, 3]) 
            z};
Bx[x_, y_, z_] := {(Subscript[b, 1]) x, (Subscript[b, 2]) y, (Subscript[b, 3]) z};

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

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

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

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

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

lhs == rhs

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]]

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

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

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

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

lhs == rhs

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