Physics 229& 100 Homework #4
Name: Nasser Abbasi

1. Here is a list of 2D vectors:

vecs = Partition[Flatten[Table[{Cos[θ]/r, Sin[θ]/r}, {θ, 0, 2π, .07}, {r, .1, 3, .5}]], 2] ;

In one line of code, construct a list containing the 5 longest vectors. Use Sort, Take and a pure function.

Take[Reverse[Sort[vecs/.{a_, b_} → {Sqrt[a^2 + b^2], a, b}]]/.{len_, a_, b_} → {a, b}, 5]//MatrixForm

( {{4.35682, 9.001}, {1.69967, 9.8545}, {-6.61179, -7.50228}, {-8.51691, 5.24044}, {-9.42222, 3.34988}} )

Take[Sort[vecs, #1 . #1>#2 . #2&], 5]

{{9.98586, -0.531602}, {9.92422, -1.22874}, {9.81398, -1.91986}, {9.65566, -2.60157}, {9.45005, -3.27055}}

2.  The Korteweg-de Vries equation is a nonlinear partial differential equation that exhibits soliton solutions (pulses that propagate without changing shape). It is ∂u/∂t+(6 u)∂u/∂x+∂^3u/∂x^3=0.
a) Write the partial differential equation as an expression in u[x,t], and as an operator using pure functions.

Remove["Global`*"] ;

pdeExpression = D[u[x, t], t] + 6u[x, t] D[u[x, t], x] + D[u[x, t], {x, 3}]

u^(0, 1)[x, t] + 6 u[x, t] u^(1, 0)[x, t] + u^(3, 0)[x, t]

pdeOperator = D[#, t] + 6# D[#, x] + D[#, {x, 3}] &

∂_t#1 + 6 #1 ∂_x#1 + ∂_ {x, 3} #1&

b) Show that u=c/2sech^2(c^(1/2)/2(x-ct)) is a solution. Compare the operator and the expression representations of the PDE for this purpose.

candidateSolution = c/2Sech[c^(1/2)/2 (x - c t)]^2 ;

exprUsingOperator = pdeOperator[candidateSolution]

If[Simplify[exprUsingOperator] == 0, Print["U is a solution"], Print["U is not solution"]]

U is a solution

expr = pdeExpression/.u→ (c/2Sech[c^(1/2)/2 (#1 - c #2)]^2&)

If[Simplify[expr] == 0, Print["U is a solution"], Print["U is not solution"]]

U is a solution

Notice the slight difference in the PDE in the two forms. In both cases ofcourse the solution is verified correctly. The expression representation is shorter, but using simplify we obtain the same final answer in both cases

3. Here is a complicated expression

φ2 = D[Underoverscript[∑, el = 0, arg3] (A[el]/r^(el + 1) + B[el] r^el) SphericalHarmonicY[el, 0, θ, φ], θ]//Expand

which is a typical form of a trial solution for many types of problems with spherical symmetry.
a) In one line of code, find all the terms of the series proportional to Cos[θ]^2

terms = Cases[φ2, any_.Cos[θ]^2]

b) in one line of code, construct a list  containing the powers of r that occur in the terms you found in part a.

terms/.any_.r^p_→p

{-4, -6, 3, 5}

4) Here is a differential equation:

Remove["Global`*"]

de1 = y ' ' '[t] + y '[t]^2 == 0

y^′[t]^2 + y^(3)[t] == 0

a)Verify that DSolve either cannot solve the equation, or produces a very complicated implicit solution, depending on which version of Mathematica you are running.

sol = DSolve[de1, y[t], t]

Solve :: ifun : Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. More…

b) Use a pure function to substitute a trial solution of the form a t^n.

poly = de1/.y→ (a #^n&)

a (-2 + n) (-1 + n) n t^(-3 + n) + a^2 n^2 t^(-2 + 2 n) == 0

c) The result should be a polynomial in t. Think carefully about the conditions for a polynomial of this type to vanish for all t. Find values of a and n which make the solution work for all t and verify that your solution works.

For the polynomial to be zero everywhere, then  f(a,n) t^-(3+n) = - g(a,n) t^(-2+2n), which happens when the powers t is raised to are equal, and when f=g. So we solve to find when this happens.

{leftTerm, firstPower, secondTerm, secondPower} = poly/.l_.t^p1_ + r_.t^p2_ == 0→ {l, p1, r, p2}

{a (-2 + n) (-1 + n) n, -3 + n, a^2 n^2, -2 + 2 n}

nSolutionFromPower = Flatten[Solve[firstPower == secondPower, n]]

{n→ -1}

aSolution = Flatten[Solve[leftTerm == -secondTerm, a]]

{a→0, a→ (-2 + 3 n - n^2)/n}

aSolution = aSolution/.nSolutionFromPower

{a→0, a→6}

Print["Values of n must be ", n/.nSolutionFromPower] ;

Print["And a has the values  ", aSolution] ;

Values of n must be -1

And a has the values   {a→0, a→6}

We see that n must be -1.  When n = -1, this means a = 6.
We also have the solution of a=0 which will work for any n.

Hence final solution is:  n=-1 and a=6, and a=0 for any n.

In fact, when a=0 the trial solution will be y=0.It is a special solution which should be ignored.

5. Explore properties of matrix Transpose: show that < Ax,y>=<x,A^Ty> for a random x, y and A.  Make a function randvec[] that makes a random 3-vector, and randmat[] that makes a random 3x3 matrix. Use these to make a function test[] which will compute the 2 inner products and test if they are equal. Verify the relation for a number of cases.

Remove["Global`*"] ;

n = 3 ;

randvec[] := Table[Random[], {n}] ;

randmat[] := Table[Random[], {n}, {n}] ;

Do[test[], {4}]

Yes, they are same inner products

Yes, they are same inner products

Yes, they are same inner products

Yes, they are same inner products

6.  Here are 2 vectors:

Remove["Global`*"] ;

v1 = {1, 2, 3}    

v2 = {2, -1, 5}

{1, 2, 3}

{2, -1, 5}

a) Find a vector perpendicular to these vectors using NullSpace, and verify that they are perpendicular.

Needs["Graphics`Colors`"]

(A = {v1, v2})//MatrixForm

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

Print["p=Vector perpendicular to these vectors is ", p = Flatten[NullSpace[A]]] ;

p=Vector perpendicular to these vectors is {-13, -1, 5}

Print["Dot product of p with v1 is ", p . v1] ;

Dot product of p with v1 is 0

Print["Dot product of p with v1 is ", p . v2] ;

Dot product of p with v1 is 0

origin = {0, 0, 0} ;

viewPoint = {3, 1.9, 1} ;

v1Line = Graphics3D[{Red,     Line[{origin, v1}]}, AspectRatio→Automatic, ViewPoint→viewPoint] ;

v2Line = Graphics3D[{Blue,   Line[{origin, v2}]}, AspectRatio→Automatic, ViewPoint→viewPoint] ;

pLine =   Graphics3D[{Black, Thickness[0.01], Line[{origin, p}]}, AspectRatio→Automatic, ViewPoint→viewPoint] ;

Show[v1Line, v2Line, pLine] ;

[Graphics:HTMLFiles/index_84.gif]

The above diagram shows geomertically that the vector p, which is the SOLID BLACK line in the diagram is perpendicular to both v1 and v2 vectors (the lines with the colors BLUE and RED). So the vector p, is perpendicular to the plane spanned by v1 and v2 as the diagram also shows.

b) Use Cross to find a vector perpendicular to v1 and v2, and verify that it is proportional to the one you found in a)

v3 = Cross[v1, v2]

{13, 1, -5}

angleBetweenV3Andp = ArcCos[v3 . p/(Norm[v3] Norm[p])]

π

If[angleBetweenV3Andp == Pi || angleBetweenV3Andp == 0, Print["Vectors are parallel as expected"], Print["Not verified. Something is wrong"]] ;

Print["Angle between v3 and p (from part a) is zero or 180, hence Cross product of v1,v2 is proportional to one found in part(a)"] ;

Vectors are parallel as expected

Angle between v3 and p (from part a) is zero or 180, hence Cross product of v1,v2 is proportional to one found in part(a)

Here are 2 vectors in 5 dimensional space:

w1 = {-3, 2, 0, 1, -1}

w2 = {1, 2, 3, -3, 0}

{-3, 2, 0, 1, -1}

{1, 2, 3, -3, 0}

c) Find a vector perpendicular to w1 and w2 using NullSpace.

(A = {w1, w2})//MatrixForm

( {{-3, 2, 0, 1, -1}, {1, 2, 3, -3, 0}} )

Print["Basis of Nullspace=\n", basisOfNullSpace = NullSpace[A]] ;

Basis of Nullspace=<br /> {{-2, 1, 0, 0, 8}, {1, 1, 0, 1, 0}, {-6, -9, 8, 0, 0}}

Now generate some scales to stretch the basis vectors with. These scales are what we consider to be the coordinates of the position vector itself. Multiply each basis vector by some integer number (I could also have used Real number, but it is more clear if I use Integer multiplier.

scales = Table[Random[Integer, 10], {Length[basisOfNullSpace]}] ;

Print["The following is a random scalling of the basis \n", s = scales * basisOfNullSpace] ;

Print["The following is a vector s perpendicular to w1 and w2  \n", Total[s]] ;

The following is a random scalling of the basis <br /> {{0, 0, 0, 0, 0}, {3, 3, 0, 3, 0}, {-30, -45, 40, 0, 0}}

The following is a vector s perpendicular to w1 and w2  <br /> {-27, -42, 40, 3, 0}

To verify, dot multiply 's' with w1 and w2 at a time, the result should be zero in both cases

s . w1

s . w2

{0, 0, 0}

{0, 0, 0}

7.  Here is a matrix

A = ({{1, 3, -2, 0, 0}, {-1, 2, 2, -1, 2}, {3, 1, 4, 0, 3}, {0, 5, 2, 1, 4}, {6, -5, -4, 1, -5}}) ;

a) What is the rank of A?

MatrixRank[A]

4

b) Find a value of rhs so that A.x==rhs is a consistent equation.

{nRow, nCol} = Dimensions[A]

{5, 5}

Det[A]

0

basisOfRangeOfA = NullSpace[NullSpace[Transpose[A]]]

{{1, 0, 0, 0, 1}, {1, 0, 0, 1, 0}, {-1, 0, 1, 0, 0}, {2, 1, 0, 0, 0}}

scales = Table[Random[Integer, 10], {Length[basisOfRangeOfA]}] ;

Print["rhs = ", rhs = Total[scales * basisOfRangeOfA]]

rhs = {23, 8, 3, 0, 10}

c) Using your choice of rhs, find the general solution using LinearSolve and NullSpace.

nullSpaceBasis = NullSpace[A]

{{3, -27, -39, -19, 58}}

x = LinearSolve[A, rhs]

{259/58, 221/58, -103/29, -693/58, 0}

Clear[t] ;

generalSolution = x + t * Random[Integer, 10] nullSpaceBasis[[1]]

{259/58 + 30 t, 221/58 - 270 t, -103/29 - 390 t, -693/58 - 190 t, 580 t}

Where 't' above is any scalar

8.  Here is a table of data in {x,y} pairs:

Remove["Global`*"] ;

a) Fit the data to the form y=a x^2+b x +c. Do this by requiring that this form be true for every data point. This will produce many equations in the three unknowns a, b, c. Use PseudoInverse to solve this overdetermined set of equations to get the best least squared estimates of a,b,c. Compare the data and the fit on the same plot.

x = data[[All, 1]] ;

y = data[[All, 2]] ;

A = Table[{x[[i]]^2, x[[i]], 1}, {i, 1, Length[x]}] ;

{a, b, c} = PseudoInverse[A] . y ;

2 Print["This is the fit: ", a , x +, b, " x +", c] ;

2 This is the fit: 0.985304 x +2.24634 x +1.69617

p1 = ListPlot[data, PlotStyle->PointSize[.02], DisplayFunction→Identity] ;

p2 = Plot[a x^2 + b x + c, {x, -3, 3}, DisplayFunction→Identity] ;

Show[{p1, p2}, AxesLabel→ {"x", "y"}, DisplayFunction→$DisplayFunction] ;

[Graphics:HTMLFiles/index_142.gif]

b) Now use Fit to find the values of a,b,c. The method you used above is the numerical method that Fit uses. Show that the answers obtained are the same both ways.

Clear[x] ;

eq = Fit[data, {1, x, x^2}, x]

1.69617 + 2.24634 x + 0.985304 x^2

We see that the Fit command gives the same answer as part(a)

9. (grads) Residues can be used to compute integrals that would otherwise be difficult or impossible. The same is true of sums. We have previously seen that Mathematica can do sums such as

Underoverscript[∑, n = 1, arg3] 1/n^2

π^2/6

The trick for doing sums of the form  Underoverscript[∑, n = 1, arg3]g[n]  using complex function theory is to find a contour integral that is reasonably easy to do and whose residues are the values g[1], g[2] .., g[n]. The first step is to find a function that has poles at integer values along the real line. Such a function is π Cot[π z].

a) Expand  π Cot[π z] in a series centered at z=n, with n an integer for several values of n to convince yourself that the poles are there.

b) Construct a proof (or at least a plausibility argument ) that these are the only poles the function has.

c) Now consider the integral  ∮  (π Cot[π z])/z^2dz over a contour consisting of a circle centered at the origin of radius R as R->∞. The value of the integral is zero. Show why this is true by examining the behavior of the real and imaginary parts of π Cot[π z] in polar coordinates for fixed θ and large r. Find an upper bound for the integral and show that the upper bound tends to zero for large r.

d) Finally, use the residue theorem to express the integral (which=0) in terms of the sum of the residues. Construct a table using Table of the residues of (π Cot[π z])/z^2 for n=-10,-9...9,10. From the structure of the table, explain why Underoverscript[∑, n = 1, arg3]1/n^2=π^2/6.

e) Calculate the sum Underoverscript[∑, n = 1, arg3]1/(n^2 + a^2) using the above technique. (note that this time there are poles not on the real line) Calculate the appropriate residues, and check your answer using Sum.

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