Legal Notice, Disclaimer, and Book Ordering Information:

(* Legal Notice, Disclaimer, and Book Ordering Information:

Last Modified: FIO June 8, 2002

Mathematica For Physics: 2nd Edition

Robert L. Zimmerman and Fredrick I. Olness

WebSite:  
  www.physics.smu.edu/~olness
  darkwing.uoregon.edu/~phys600/

MathSource Number: 0206-862

ISBN 0-201-53796-6
For ordering information, call 1-800-282-0693
Addison-Wesley Publishing Company

Communication with the authors:
Fredrick I. Olness:  olness@mail.physics.smu.edu
Robert L. Zimmerman: bob@zim.uoregon.edu

Copyright 2002, Addison-Wesley Publishing Company, Inc. The material in
this file may be distributed freely  so long as the content remains
unchanged, and the copyright and  reference notices are included. The
publisher grants permission for the noncommercial use of these programs
and program segments.  All other uses require the prior written consent
of the publisher.    

All rights are reserved with regard to the  material in the
"Mathematica For Physics" text, and may not be reproduced, stored in a
retrieval system, or transmitted,  in any form or by any means,
electronic, mechanical, photocopying,  recording, or otherwise, without
the prior written permission  of the publisher. 

The programs and applications presented in this book have been included
for  their instructional value. They have been tested with care but are
not  guaranteed for any particular purpose. The publisher does not
offer any  warranties or representations, nor does it accept any
liabilities with  respect to the programs or applications. 

 *)

2. GENERAL PHYSICS

Introduction

Off[General::spell1];   
Off[General::spell];

2.1 Newtonian Mechanics in Inertial Frames

Overview of Newtonian Mechanics in Inertial Frames

Problem 1: Escape Velocity

Remarks and Outline

Solution

Clear["Global`*"];

T= (1/2) m v^2;     
U= -G M m/r;
energy=T+U

sol = Solve[energy == 0, v][[2]]

Needs["Miscellaneous`PhysicalConstants`"];

?Miscellaneous`PhysicalConstants`*

values = {G -> 
            GravitationalConstant, 
r -> 
            EarthRadius, 
M -> 
            EarthMass, 
Newton -> 
            Kilogram Meter/Second^2};

v0 = v /. sol //. values  // PowerExpand

Needs["Miscellaneous`Units`"];

?Convert

Convert[v0,Kilo Meter/Second]

Convert[v0,Mile/Second]

Problem 2: Projectile in a Uniform Gravitional Field

Remarks and Outline

Solution

Clear["Global`*"];

Part a

eq1 ={m x''[t]==0, m z''[t]==-m g };

initial1={x[0]==0, x[t1]==x1,z[0]==0, z[t1]==z1};

dsol1=DSolve[Join[eq1,initial1], {x[t],z[t]},t] //Flatten

Map[Collect[#, t] &, dsol1, {2}]

Map[collect[#, t] &, dsol1, {2}]

values1 = {x1 -> 300, z1 -> 0, g -> 9.8};

coord1[t_,t1_]= {x[t],z[t]} /.dsol1 /.values1

points[t1_]  := Table[coord1[t, t1], {t, 0, t1, 0.5}]

plot1[t1_] := 
      ListPlot[ 
        points[t1] // 
          Evaluate, 
PlotStyle -> {Hue[0.8], 
            PointSize[0.02]}, 
GridLines -> 
          Automatic, 
DisplayFunction -> 
          Identity]

\(\(Show[{plot1[5], plot1[10], 
          plot1[15]}, 
Epilog -> {
Text["t1= 5", coord1[\ \ 5/2, \ \ 5], 
              TextStyle -> {FontWeight -> "Bold"}]\ , 
Text["t1=10", coord1[10/2, 10], 
              TextStyle -> {FontWeight -> "Bold"}], 
Text["t1=15", coord1[15/2, 15], 
              TextStyle -> {FontWeight -> "Bold"}]}\ , 
DisplayFunction -> $DisplayFunction\ \ , 
PlotRange -> All];\)\)

Part b

initial2={x[0]==0 ,x'[0]==v0 Cos[θ], z[0]==0, z'[0]==v0 Sin[θ]};

dsol2=DSolve[Join[eq1,initial2],{x,z},t] //Flatten

eq2 =   z'[t] == 0 /. dsol2

z'[t] == 0 /. dsol1

tsol= Solve[eq2,t] //Flatten

{x[t],z[t]} /.dsol2 /.tsol //Simplify

values2 = {g->9.8, v0->100};

coord2[t_,θ_]= {x[t],z[t]} /.dsol2 /.values2

plot2[θ_] := 
      ListPlot[ 
        Table[coord2[t, θ], {t, 0, 20, 1/2}] // 
          Evaluate, 
PlotStyle -> {Hue[ .8], 
            PointSize[0.02]}, 
GridLines -> 
          Automatic, 
DisplayFunction -> 
          Identity]

\(\(Show[\ \ \(plot2[#] &\)\  /@ \ {π/8, \ 2\ π/8, \ 
            3\ π/
                8}, 
Epilog -> \ {
Text["θ= π/8", \ coord2[5, \ \ \ π/8]\ , 
              TextStyle -> {FontWeight -> "Bold"}], 
            Text["θ=2π/8", \ coord2[5, 2  π/8]\ , 
              TextStyle -> {FontWeight -> "Bold"}], 
Text["θ=3π/8", \ coord2[5, 3  π/8]\ , 
              TextStyle -> {FontWeight -> "Bold"}]}, 
        PlotRange -> {{0, 1000}, {0, 
              500}}, 
DisplayFunction -> $DisplayFunction, 
FrameLabel -> {"x[t]", "z[t]"}
];\)\)

eq3 =  coord2[t, θ] == {xMax, 0}  // Thread

sol3 =  Solve[eq3, {t, xMax}]

xMax = xMax  /. sol3  // Last

(θ/π)  /. 
        Solve[ D[xMax, θ] == 0 , θ]  // 
      Rationalize

Problem 3: Reflecting Trajectories

Solution

Clear["Global`*"];

{xEq, yEq} = {x[t] == x0 +  vx0 t  +  (1/2) ax (t^2), 
          y[t] == y0 +  vy0 t  +  
              (1/2) ay (t^2)};
rules = {vy0  ->  0,  ax -> 0,  
          ay -> -g};

    
{xEq, yEq} = {xEq, yEq}  /. rules

ySol =  DSolve[yEq, y, t]  // Flatten

xSol =  DSolve[xEq, x, t]  // Flatten

position[
        t_] := If[ t < tWall, {+x[t], y[t]}, {-x[t], y[t]}] /. 
          ySol  /. xSol

tWall =  
      t  /. Solve[x[t] == 0  /. xSol, t][[1]]

values = {x0 -> -100, g -> 9.8,  
        y0 -> 100}

plot[v_] := 
        ParametricPlot[position[t]  //. values  /. {vx0 -> v}  // 
            Evaluate, {t, 0, 5}, 
          DisplayFunction -> Identity];

Show[ 
        Table[plot[v], {v, 5, 50, 
            5}], 
DisplayFunction -> $DisplayFunction, 
PlotRange -> {{-125, 0}, {0, 
              100}}, 
Epilog -> {Thickness[0.02], 
            Hue[1.0], 
            Line[{{0, 0}, {0, 100}}]}
];

Problem 4: Falling Projectile with Linear Air Resistance Drag

Remarks and Outline

Solution

Clear["Global`*"];

Part a

eq1 =  {m z''[t] == -m g - 
              b z'[t]};

initial = {z[0] == z0, z'[0] == 0};

z[t_] = z[t] /. DSolve[ Join[eq1, initial], {z}, t]  // 
          First//Simplify

v[t_] = z'[t]  // Simplify

vTerminal = 
      Limit[v[t]  /. {t -> τ m/
                  b}, τ -> ∞]

{zff,  vff} = Limit[{z[t], v[t]}, b -> 0]

Part b

values = {g -> 9.8, m -> 1, z0 -> 0};

pointZ[b_, t_] = {t, z[t]} /. values

pointZ[0, t_] = {t, zff} /. values

plotZ[b_, hue_] := 
        ParametricPlot[ 
          pointZ[b, t] , {t, 0, 
            10}, 
PlotStyle -> {Hue[hue], 
              Thickness[0.008]}, 
DisplayFunction -> 
            Identity
];

part1 = 
        Show[ Table[ plotZ[b, b/2], {b, 0, 1, 1/2}], 
          DisplayFunction -> $DisplayFunction
];

textZ = {
Text[
              FontForm["b=0  ", {"Times-Bold", 12}], 
              pointZ[0, 7] ], 
Text[
              FontForm["b=1/2", {"Times-Bold", 12}], 
              pointZ[0.5, 7]], 
Text[
              FontForm["b=1  ", {"Times-Bold", 12}], 
              pointZ[1, 7]]
}  // Graphics;

\(\(Show[{part1, textZ}, 
Frame -> 
          True, 
GridLines -> 
          Automatic, 
FrameLabel -> {"time", "height", "Effect of Air Resistance", " "}, 
RotateLabel -> 
          False, 
DisplayFunction -> $DisplayFunction
];\)\)

Part c

pointV[b_, t_] = {t, v[t]} /. values

pointV[0, t_] = {t, vff} /. values

plotV[b_, hue_] := 
        ParametricPlot[ 
          pointV[b, t] , {t, 0, 
            10}, 
PlotStyle -> {Hue[hue], 
              Thickness[0.008]}, 
DisplayFunction -> 
            Identity
];

part2 = 
        Show[ Table[ plotV[b, b/2], {b, 0, 1, 1/2}], 
          DisplayFunction -> $DisplayFunction
];

textV = {
Text[
              FontForm["b=0  ", {"Times-Bold", 12}], 
              pointV[0, 7] ], 
Text[
              FontForm["b=1/2", {"Times-Bold", 12}], 
              pointV[0.5, 7]], 
Text[
              FontForm["b=1  ", {"Times-Bold", 12}], 
              pointV[1, 7]]
}  // Graphics;

\(\(Show[{part2, textV}, 
Frame -> 
          True, 
GridLines -> 
          Automatic, 
FrameLabel -> {"time", "speed", "Effect of Air Resistance", " "}, 
RotateLabel -> 
          False, 
DisplayFunction -> $DisplayFunction
];\)\)

Part d

<<Statistics`NormalDistribution`

gaussian :=  
      Random[NormalDistribution[1, 0.1]]

    
Table[gaussian, {i, 10}]

data =  Table[{t, z[t] gaussian}, {t, 0, 40, 1}]  /. 
            values  /. {b -> 1/2};

    
listplot = 
        ListPlot[data, PlotStyle -> {PointSize[0.03]}];

<< Statistics`NonlinearFit`

?NonlinearFit

fit = 
      NonlinearFit[data, 
        a1 + a2 t  + a3 Exp[-a4 t], {t}, {a1, a2, a3, 
          a4} ]

fitPlot =  
        Plot[{fit, z[t]} /. values  /. {b -> 1/2}  // 
            Evaluate, {t, 0, 40}, 
          PlotStyle -> {{RGBColor[0, 0, 1]}, {RGBColor[1, 0, 0], 
                Thickness[0.008]}}];

Plot[fit/z[t] /. values  /. {b -> 1/2}  // Evaluate, {t, 
          0, 40}];

\(\(Show[listplot, 
        fitPlot, 
GridLines -> 
          Automatic, 
Frame -> 
          True, 
FrameLabel -> {"time", "height", "Experimental Data", ""}, 
RotateLabel -> 
          False
];\)\)

Problem 5: Falling Projectile with Quadratic Air Resistance Drag

Remarks and Outline

Solution

Clear["Global`*"];

Part a

eq1[
          sign_] =  {m z''[t] == -m g -  
              sign b z'[t]^2};

dsolUp =  DSolve[ Join[eq1[+1], {z'[0] == vz0}], {z}, 
            t]  // First   // Simplify

c2sol= Solve[z[0] == z0  /. dsolUp, (Subscript[c, 2])] //First  //Simplify

dsolUp = dsolUp  /. c2sol;

{zUp, 
        vUp} =  {z[t], z'[t]}  /. 
          dsolUp  //FullSimplify

{zUp0, vUp0} = Limit[{zUp, vUp}, b -> 0]

Part b

dsolDown =  DSolve[ 
            Join[eq1[-1], {z'[timeTop] == 0}], {z}, t]  // 
          First   // Simplify

c2sol= Solve[z[timeTop] == zMax  /. dsolDown, (Subscript[c, 2])] //First  //Simplify

dsolDown = dsolDown  /. c2sol;

{zDown, 
        vDown} =  {z[t], z'[t]}  /. 
          dsolDown //FullSimplify

{zDown0, vDown0} =  Limit[{zDown, vDown}, b -> 0]

vTerminal =  
      Limit[ vDown /. {t -> (((Sqrt[m]) τ)/((Sqrt[b]) (Sqrt[g])))}, τ -> ∞]

Part c

values = {g -> 9.8, m -> 1, z0 -> 0,  
          b -> 0.1,  vz0 -> +20};

timeTop = 
      t  /.  Solve[vUp == 0  /. values, t]  // First

zMax = zUp /. values  /. {t -> timeTop}

zFull[
        tIn_] := If[tIn < timeTop,  zUp, zDown]  /. 
          values  /. {t -> tIn}
vFull[
        tIn_] := If[tIn < timeTop,  vUp, vDown]  /. 
          values  /. {t -> tIn}

zFull0[
        tIn_] := zUp0 /. 
          values  /. {t -> tIn}
vFull0[tIn_] := vUp0 /. values  /. {t -> tIn}

Plot[{zFull[t], zFull0[t]} , {t, 0, 
          7}, 
PlotStyle -> {{Thickness[
                0.008]}, {}}];

Plot[{vFull[t], vFull0[t]} , {t, 0, 
          7}, 
PlotStyle -> {{Thickness[
                0.008]}, {}}];

{vFull[100], vTerminal } /. values

Problem 6: Rocket with varying Mass

Remarks and Outline

Solution

Clear["Global`*"];

Part a

p[t   ]=  m[t] v;
p[t+δt]= (m[t]-δm)(v+δv) + (+δm)(v-ve);

δv = x''[t] δt;
δm = -m'[t]   δt;

δp =(p[t+δt] - p[t])//ExpandAll

eq1 = -m[t] g == δp/δt  // 
        ExpandAll

eq1 = eq1 /.  δt->0

Part b

dsol1 = 
      DSolve[{m'[t] == -Δm, m[0] == m0}, m, 
          t] // Flatten

eq2 = eq1 /. dsol1 // Simplify

eq3 = (Solve[eq2, x''[t]] // First  // 
            First) /. Rule->Equal

position = x[t] /. 
            DSolve[{eq3, x[0] == 0, x'[0] == 0}, x, t] // 
          First // FullSimplify

speed =  D[position, t] // Simplify

height= Limit[position,t->m0/Δm]

Part c

values = {m0 -> 100, Δm -> 10, 
          g -> 9.8};

tf = (m0/Δm) /. values

plot1[
        veIn_] := 
ParametricPlot[{position, speed} /. 
              values /. {ve -> veIn}  // 
          Evaluate, 
{t, 0, 
          tf}, 
PlotStyle -> 
          Thickness[ .010], 
DisplayFunction -> 
          Identity
]

elevation[
          ve_] = 
{Thickness[0.004], Hue[0.8], 
              Line[{{height, 0}, {height, 10000}}]} /. values // 
          Graphics;

\(\(Show[
        Table[{elevation[ve], plot1[ve]}, {ve, 100, 400, 
            100}], 
Frame -> 
          True, 
FrameLabel -> {"position", "speed", "", " "}, 
RotateLabel -> 
          False, 
PlotRange -> {{0, 4000}, {0, 
              4000}}, 
DisplayFunction -> $DisplayFunction
];\)\)

Part d

initial = {x[0] -> 0, x'[0] -> 
            0};
order = 6;

eq4 =  Series[eq3[[1]] - eq3[[2]], {t, 0, order}] /. 
          initial//Simplify

variables = 
      Table[D[x[t], {t, n}], {n, 0, order + 1}] /. {t -> 0}

sol1 = Solve[eq4 == 0, variables][[1]]

series1 = (x[t] + O[t]^order)  /. sol1 /. initial // 
        ExpandAll

term[
          n_] = (((t^(n + 1)) ve (((Δ(m/m0)))^n))/(n (n + 1)));

series2= -(1/2) g t^2 + Sum[term[n],{n,1,order-2}];

series1==series2 //Normal //ExpandAll

exact2 =  -(1/2) g t^2  +  
          Sum[term[n], {n, 1, ∞}]  //  FullSimplify

((position == exact2 )
 // 
            MapAll[Together, #] & 
 // 
          PowerExpand
 //  
        FullSimplify
)

Problem 7: Keplerian Orbits

Remarks and Outline

Solution

Clear["Global`*"];

Part a

Mercury  ={3.300*10^23,    0.5791*10^11,     0.241,  
              0.206};
Venus       = {4.870*10^24,    1.0821*10^11,     0.615,  0.007};
Earth       = {5.974*10^24,    1.4960*10^11,     1.000,  0.017};
Mars         = {6.419*10^23,    2.2794*10^11,     1.881,  0.093};
Jupiter  = {18.99*10^26,    7.7840*10^11,   11.860,  
            0.048};
Saturn    = {5.680*10^26,  14.2700*10^11,   29.420,  
            0.054};
Uranus   ={8.680*10^25,  28.7100*10^11,   83.750,  0.047};
Neptune  = {10.20*10^25,  44.9800*10^11, 163.700,  
            0.009};
Pluto      = {1.270*10^22,   59.5060*10^11, 248.000,  
            0.249};

data = {Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, 
          Neptune, Pluto};

\(\(names = {"Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto"};\)\)

\(TableForm[data, 
      TableHeadings -> {names, {"Mass", "Axes", "Period", "ϵ"}}]\)

{mass, axes, period,  ϵ} = Transpose[data];

period

<<Miscellaneous`Units`
<<Miscellaneous`PhysicalConstants`

values = {
year2sec -> 
            Convert[Year, Second],  
G -> 
            GravitationalConstant,  
Newton  ->  
            Kilogram Meter/(Second^2), 
AstronomicalUnit -> 
            Convert[AstronomicalUnit, Meter]};
values  // TableForm

kepler = ((4 (a^3) (π^2))/((T^2) G))  /. {a -> 
              axes * Meter, T ->  period * year2sec}   //. 
        values

solarMass = 1.99 * (10^30) Kilogram; 
    TableForm[(kepler/solarMass), TableHeadings -> {names}] // 
      PaddedForm[#, {4, 3}] &

Part b

r = ((a (1 - (e^2)))/((1 + e Cos[ θ])))  /. {a -> axes * Meter, T ->  period * year2sec,  
            e -> ϵ}   //. values

points = 
      Outer[#1 * #2 &, 
          (r/AstronomicalUnit) , {Cos[θ], Sin[θ]}]  /. 
        values

plot1 = 
        ParametricPlot[ 
          points[[{1, 2, 3, 4}]]  // Evaluate, {θ, 0, 2 π}, 
          PlotRange -> All, AspectRatio -> 1];

plot2 = 
        ParametricPlot[ 
          points[[{5, 6, 7, 8, 9}]]  // Evaluate, {θ, 0, 
            2 π}, PlotRange -> All, 
          AspectRatio -> 1];

Show[GraphicsArray[{plot1, plot2}]];

2.2 Newtonian Mechanics in Rotating Frames

Overview of Newtonian Mechanics in Rotating Frames

Problem 1:  Projectile Motion as Measured by an Observer on Earth

Remarks and Outline

Solution

Clear["Global`*"];

Part a

r[t_] = {x[t], y[t], z[t]};
ω[t_] = {-ω0 Sin[θ], 
          0, ω0 Cos[θ]};

fTransverse =  -    m Cross[ω'[t], 
          r[t]]

fCoriolis =    -   2 m Cross[ω[t],  r'[
            t]]

fCentrifugal =  -m Cross[ ω[t],  
            Cross[ω[t], r[t]] ]  // Simplify

fInertial =  {0, 0, -m g};

eq1 = 
      m  r''[t] - fInertial  - fTransverse - fCoriolis - 
        fCentrifugal

initialRule =  
      Thread  /@  {r[0] ->  0, r'[0] ->  {1, 0, 0}}  // 
        Flatten

nOrder = 0;

eq2 =  ( Series[eq1/m, {t, 0, nOrder}] == 
                  0 
 /. initialRule // Normal  // 
            Thread  // Simplify )

vars =  
Table[
            D[{x[t], y[t], z[t]}, {t, i}], {i, 2, nOrder + 2}] /. 
          t -> 0 //Flatten

sol =  Solve[eq2 , vars]  // First

point[
        t_, θ_] = 
Series[ 
              r[t], {t, 0, nOrder + 2}]  /. sol  /. initialRule // 
        Normal

plot[θ_]:=ParametricPlot[ point[
                    t, θ][[{1, 2}]] /. {ω0 -> 1, 
                  g -> 9.8}  // Evaluate,
{t, 0, 20},
PlotStyle -> { Thickness[0.016]},
\(AxesLabel -> {"x,South", "y,East"}\),
PlotRange -> {{0, 
                20}, {+200, -200}},
Ticks -> False ,
 DisplayFunction -> Identity]

{plot[π/4], plot[3  π/4]}  // GraphicsArray  // 
        Show;

Part b

initial =  
      Thread  /@  {r[0] == 0, r'[0] == {vx0, 0, 0}}  // 
        Flatten

eqNorth = 
      Join[Thread[eq1 == 0], 
          initial]  /. θ->0

solNorth =  DSolve[eqNorth, {x[t], y[t], z[t]}, t]  // First // 
        FullSimplify

values = {g -> 9.8, z0 -> 500, vx0 -> 1};

pointNorth[t_, ω0_] = r[t] /. solNorth  /. 
        values

plot1 = 
        ParametricPlot3D[
          pointNorth[t, 0] // Evaluate, 
{t, 0, 
            10}, 
BoxRatios -> {1, 1, 
              1}];

plot2 = 
        ParametricPlot3D[
          Join[pointNorth[t, 0.01], {Thickness[ .02]} ] // Evaluate, {t, 0, 
            10}, 
BoxRatios -> {1, 1, 
              1}];

\(\(Show[\ {plot1, 
          plot2}, 
FaceGrids -> {{0, 1, 0}, {\(-1\), 0, 
              0}}, 
BoxRatios -> {1, 1, 
            1}, 
ViewPoint -> {1.5, \(-1.0\), 
            0.7}, 
AxesLabel -> {"x-south", "west", "z-vertical"}
];\)\)

Part c

eqEquator = 
      Join[Thread[eq1 == 0], 
          initial]  /. θ -> π/2

solEquator =  DSolve[eqEquator, {x[t], y[t], z[t]}, t]  // First // 
        FullSimplify

pointEquator[t_, ω0_] = r[t] /. solEquator /. 
        values

pointEquator[t_, 0] =  
      Limit[pointEquator[t, ω0], ω0 -> 0]

plotEquator[ω0_, 
          style_] := 
ParametricPlot3D[
Join[ pointEquator[t, ω0], {style}]  /. values  // 
            Evaluate, 
{t, 0, 
            10}, 
BoxRatios -> {1, 1, 
              1}, 
DisplayFunction -> 
            Identity];

\(\(Show[\ 
{\ 
          plotEquator[
            0\ \ \ \ \ , {Thickness[0.02], 
              Hue[0.6]}], 
plotEquator[
            0.01, {Thickness[0.01], 
              Hue[0.4]}]\ }, 
FaceGrids -> {{0, \(-1\), 
              0}, \ {\(-1\), 0, 0}}, 
BoxRatios -> {1, 1, 
            1}, 
ViewPoint -> {1, 0.7, 
            0.7}, 
AxesLabel -> {"x-South", "y-East", "z-Up"}, 
DisplayFunction -> $DisplayFunction\ \ ];\)\)

Problem 2:  Foucault Pendulum

Remarks and Outline

Solution

Clear["Global`*"];

Part a

r[t_] = {x[t], y[t], z[t]};
ω[t_] = {-ω0 Sin[θ], 
          0, ω0 Cos[θ]};

fTransverse =  -    m Cross[ω'[t], 
          r[t]]

fCoriolis =    -   2 m Cross[ω[t],  r'[
            t]]

fCentrifugal =  -m Cross[ ω[t],  
            Cross[ω[t], r[t]] ]  // Simplify

fGravity = {0, 0, -m g};
fTension = 
        tension ( {0, 0, 1} - r[t]/L );
fInertial = fGravity + fTension

eq1 = 
      0 == 
          m  r''[t] - fInertial  - fTransverse - fCoriolis - 
            fCentrifugal  // Thread

initialRule =  
      Thread  /@  {r[0] -> 0, r'[0] -> {1, 0, 0}}  // 
        Flatten

zSol = DSolve[{z[t] == 0}, z, t]  // First

eq2 =  eq1 /. {ω(0^2) -> 0,  Sin[θ] -> 0}  /. 
        zSol

eq3 =  
      Solve[eq2, {x''[t], y''[t]}, 
            tension]  /. {Rule  ->  Equal}  // First

eq4 =  eq3 /. { g ->  (ω^2) L}  // Simplify

DSolve[eq4, {x, y}, t];

Part b

imEq =   I eq4[[2]] // Thread[#, Equal] &

reEq= eq4[[1]]

eq5 =  reEq  +  imEq // Thread[#, Equal] &  // 
        Simplify

eq6 =  
      eq5  /. DSolve[u[t] == x[t] + I y[t], x, t][[1]]  // 
        Simplify

initial = {u[0] == I y0, u'[0] == 0};

dsol = DSolve[Join[{eq6}, initial], {u[t]}, t]  // First

uSol = ( u[t] 
 /. 
                  dsol 
 /.  ω0^2 -> 
                  0
 // PowerExpand
 // 
            ExpToTrig  
 // ExpandAll )

x =  reSol =  
        uSol /. {Complex[a_, b_] :> a}

y =  imSol =  (uSol -  reSol)/I  // Simplify

point = {x, y};

point /. {ω(0/ω) -> 0}  // FullSimplify // 
      TableForm

Part c

values = {ω -> Sqrt[g/L],  L -> 1,  
          g -> 9.8,  ω0 -> 1/10,  
          y0 -> 1};

\(\(plot[θin_, 
          options_: {}] := 
ParametricPlot[\ \(point\ \  //. 
                values\) /. {θ -> θin} // 
            Evaluate, 
{t, 0, 
            10}, 
Ticks -> 
            False, 
PlotStyle -> {options}, 
PlotRange -> {{\(-1\), 1}, {\(-1\), 1}}, 
          AxesLabel -> {"x-axis, south", "y-axis, east"}];\)\)

plot[0];

plot[π];

plot[π/2, Thickness[0.02]];

2.3 Electricity and Magnetism

Overview of Electricity and Magnetism

Problem 1: Charged Disk

Remarks and Outline

Solution

Clear["Global`*"];

Part a

dV = k  (dq/(Sqrt[(z^2) + (r^2)]));

dq = 2 π  r dr σ;

(dV/dr)

diskPotential = 
      Integrate[(dV/dr), {r, 0, R}]   // PowerExpand

Needs["Calculus`VectorAnalysis`"];

Edisk = -Grad[diskPotential, Cartesian[x, y, z]]

Part b

Eplane = Limit[Edisk, R -> ∞]

Eplane  /. {k -> (1/(4  π (Subscript[ϵ, 0])))}

Problem 2: Uniformly Charged Sphere

Remarks and Outline

Solution

Clear["Global`*"];

Part a

eq1 = {
4 π (r^2) Eoutside == 
            qTotal        , 
4 π (r^2) Einside == 
            qEnclosed   };

qTotal = ((4  π (R^3))/(3 )) ρ;
qEnclosed = ((4  π (r^3))/3) ρ;

sol = Solve[eq1, {Einside, Eoutside}]  // First

eField[rr_, 
        RR_:  1, ρρ_:  
          1] :=  If[ rr ≤ RR,  Einside, Eoutside] /. 
          sol /. {r -> rr, 
          R -> RR, ρ -> ρρ}

p1 = Plot[eField[r], {r, (10^(-6)), 3}];

medium1 = {Hue[ .2], 
          Polygon[{{1,  .51}, {1, 0}, {3, 
                0}, {3,  .51}}]};
medium2 = {Hue[ .5], 
          Polygon[{{0,  .51}, {0, 0}, {1, 0}, {1,  .51}}]};

\(\(Show[
{Graphics[{
medium1, 
medium2, 
Text["Inside ", {0.5, 
                  0.4}], 
Text["Outside", {2.0, 
                  0.4}]}], 
p1
}, 
AxesLabel -> {"Distance ", "Electric Field"}, 
Axes -> True];\)\)

Part b

vOut[r_, 
        R_] = -Integrate[Eoutside /. sol, {r, ∞, r}]

vIn[r_, R_] = (vOut[R, R] - Integrate[Einside /. sol, {r, R, r}]) // 
        Simplify

v[r_,R_]:= If[r<R,vIn[r,R],vOut[r,R]];

p2 = 
        Plot[v[r, 1] /. ρ -> 1, {r, (10^(-6)), 3}];

\(\(Show[{
Graphics[{
medium1, 
medium2, 
Text["Inside ", {0.5, 
                  0.3}], 
Text["Outside", {2.0, 
                  0.3}]}], 
p2
}, 
AxesLabel -> {"Distance ", "Potential Field"}, 
Axes -> 
          True
];\)\)

Part c

vv[x_, y_, z_:  0, 
        R_:  1, ρρ_:  1] := 
If[r < R, 
          vIn[r, R] // Evaluate, 
          vOut[r, R] // 
            Evaluate] /.{r -> (Sqrt[(x^2) + (y^2) + (z^2)]), ρ -> ρρ}

vPlot = 
        ContourPlot[
          vv[x, y] , {x, -2, 2}, {y, -2, 
            2}, 
PlotPoints -> 
            30, 
ColorFunction -> 
            Hue, 
ContourSmoothing -> 
            True, 
DisplayFunction -> 
            Identity];

Needs["Graphics`PlotField`"];

ePlot = 
        PlotGradientField[-vv[x, y], {x, -2, 2}, {y, -2, 
            2}, 
ScaleFunction -> (1 &), 
DisplayFunction -> Identity];

Show[{vPlot, ePlot}, 
        DisplayFunction -> $DisplayFunction];

Part d

Needs["Graphics`PlotField3D`"];

ePlot3D = 
        PlotGradientField3D[
-vv[x, y, z], {x, 0, 
            2}, {y, 0, 2}, {z, 0, 2}, 
VectorHeads -> 
            True,  
ScaleFunction -> (1 &), 
DisplayFunction -> Identity];

point[r_, R_:  1] := {r Sin[θ] Cos[φ], 
        r Sin[θ] Sin[φ], r Cos[θ]}

potPlot3D[
        r_] := 
ParametricPlot3D[
point[
            r] // Evaluate, {θ, 0, π/2}, {φ, 
          0, π/2}, 
PlotPoints -> 
          12, 
DisplayFunction -> Identity]

Show[
ePlot3D, 
Table[
          potPlot3D[r], {r, 1/2, 2, 
            1/2}], 
PlotRange -> {0, 
            1}, 
DisplayFunction -> $DisplayFunction];

Problem 3: Electric Dipole

Remarks and Outline

Solution:

Clear["Global`*"];

Part a

monopole[q_, p0_: {0, 0, 0}, p1_: {x, y, z}] := 
      (q/(Sqrt[(p0 - p1) . (p0 - p1)]))

dipole = (monopole[-q, {0, 0, -d/2}] + 
          monopole[+q, {0, 0, +d/2}])

ContourPlot[
        dipole /. {q -> 1, d -> 1, y -> 0}  // 
          Evaluate, 
{x, -2, 2}, {z, -2, 
          2}, 
PlotPoints -> 
          45, 
Contours -> 
          21, 
ColorFunction -> 
          Hue, 
Axes -> 
          True, 
AxesLabel -> {"x-axis", "z-axis\
"}];

Needs["Graphics`PlotField3D`"];

plot = 
PlotGradientField3D[-dipole /. {q -> 1, d -> 1}  // Evaluate, 
{x, -0.5, 
            0.5}, {y, -0.5, 0.5}, {z, -1.0, 
            1.0}, 
PlotPoints -> 
            6, 
ColorFunction -> 
            Hue, 
VectorHeads -> 
            True,  
ScaleFunction -> (1 &), 
BoxRatios -> {1, 1, 
              2}
];

Show[plot, ViewPoint -> {0, 3, 0} ];

Part b

Needs["Calculus`VectorAnalysis`"];

x2rRule = 
      Thread[{x, y, z} -> 
          CoordinatesToCartesian[{r, θ, φ}, 
            Spherical]]

potential = dipole /. x2rRule // Simplify

dipoleR = Series[potential, {r, ∞, 2}] // Normal // 
        Simplify

eField = -Grad[dipoleR, 
          Spherical[r, θ, φ]]

Part c

rotDipole[θ_] = (  monopole[+1, 
              1/2 {+Sin[θ], 0, +Cos[θ]}] + 
            monopole[-1, 
              1/2 {-Sin[θ], 
                  0, -Cos[θ]}] );

frames = 1;
Do[
PlotGradientField[
-rotDipole[θ] /. {y -> 0}, 
{x, -2, 2}, {z, -2, 
          2}, 
PlotPoints -> 
          14, 
VectorHeads -> 
          True,  
ScaleFunction -> (1 &)]
, {θ, 2 π/frames, 2 π, 
        2 π/frames}]

Problem 4: Magnetic Vector Potential for a Long Straight Wire

Remarks and Outline

Solution

Clear["Global`*"];

Part a

Az[r_]=Integrate[
        (I0/(Sqrt[(z^2) + (r^2)])), {z, -L, L}]

A[x_, y_, z_] = {0, 0, Az[(Sqrt[(x^2) + (y^2)])]}

Part b

Needs["Calculus`VectorAnalysis`"];

    
B[x_, y_, z_] = 
      Curl[A[x, y, z], Cartesian[x, y, z]] // Simplify

B∞[x_, y_, z_] = 
      Limit[B[x, y, z], L -> ∞]

Part c

Needs["Graphics`PlotField3D`"];

p1 = 
        PlotVectorField3D[
B∞[x, y, z] /. 
            I0 -> 1, 
{x, -1, 1}, {y, -1, 
            1}, {z, -1, 1}, 
PlotPoints -> 
            6, 
VectorHeads -> 
            True, 
ViewPoint -> {1, 1, 
              3}, 
ColorFunction -> 
            Hue,  
ScaleFunction -> (1 &), 
ViewPoint -> {2, 2, 
              16}
];

Needs["Graphics`PlotField`"];

    
PlotVectorField[
B∞[x, y, 
              z][[{1, 2}]] /. 
          I0 -> 1, 
{x, -3, 3}, {y, -3, 
          3}, 
PlotPoints -> 
          20, 
VectorHeads -> 
          True, 
ColorFunction -> 
          Hue,  
ScaleFunction -> (1 &)
];

Part d

Acylind[r_]={0, 0, Az[r]}

Bcylind[r_] =  
      Curl[Acylind[r], Cylindrical[r, θ, z]] // 
        Limit[#, L -> ∞] &

Problem 5: Motion of a Charged Particle in a Uniform B Field

Solution

Clear["Global`*"];

Part a

acceleration = {x''[t], y''[t], z''[
            t]};
velocity = {x'[t], y'[t], z'[
            t]};
magneticField = {0, 0, B};

eq1 = acceleration == 
          q/m  Cross[velocity, 
              magneticField]//Thread

eq1 = eq1 /. {B -> m/q  ω}

initial = {
x[0] == 0, x'[0] == 
            v0x, 
y[0] == 0, y'[0] == 
            v0y, 
z[0] == 0, z'[0] == 
            v0z};

dsol =  DSolve[Join[{eq1, initial}], {x[t], y[t], z[t]}, t]  // 
          Flatten // FullSimplify

Part b

values = {B -> m/q  ω,  ω -> 1, 
          v0x -> 1, v0y -> 1, v0z -> 1/10 };

plot1 =  
        ParametricPlot3D[{x[t], y[t], z[t], {Thickness[0.02]}} /. dsol /. 
              values  // Evaluate, 
{t, 0, 
            8  π}, 
PlotPoints -> 
            100 ];

Needs["Graphics`PlotField3D`"]

    
 plot2 = 
        PlotVectorField3D[ {0, 0, 1} // 
            Evaluate, 
{x, -2, 2}, {y, -2, 2}, {z, 
            0, 2}, 
VectorHeads -> 
            True , 
PlotPoints -> 
            4
];

Show[plot1, plot2];

Problem 6: Motion of a Charged Particle in a Uniform B Field and Time Varying E Field

Solution

Clear["Global`*"];

F = {Fx, Fy, Fz}; 
r = {x[t], y[t], z[t]}; 
EF = {Ex, Ey, Ez}; 
BF = {Bx, By, Bz};  v = D[r, t]

F = q EF + q Cross[v, BF]

eq1 =  
        F - m D[r, {t, 2}] == 0  // 
          Thread;
eq1  // ColumnForm

initial = {
x[0] == 0, x'[0] == 
            0, 
y[0] == 0, y'[0] == 
            0, 
z[0] == 0, z'[0] == 
            0};

fields =  
      Thread  /@  {BF -> {0, 0, m ω0/q}, 
            EF -> {E0  Cos[ω t ], 0, 0}}  // 
        Flatten

eqs = Join[eq1, initial] /. fields

dsol = DSolve[eqs, {x[t], y[t], z[t]}, t]  // Flatten // 
        FullSimplify

values = {m -> 1, E0 -> 1, ω0 -> 1,  
          q ->  1 };

point[t_, ω_] = {x[t], y[t], z[t]}  /. dsol  /. 
        values

point[t_, 1] = 
      Limit[ point[t, ω], ω -> 1]

point[t_, 0] = 
      Limit[ point[t, ω], ω -> 0]

plot[ω_ ] := 
        ParametricPlot[ point[t, ω][[{1, 2}]]  // 
            Evaluate, {t, 0, 
            20 π   }, 
PlotLabel -> ω, 
Ticks -> 
            None, 
DisplayFunction -> 
            Identity];

table = Table[ plot[ω], {ω, 0, 4, 1/3}];

Show[GraphicsArray[Partition[table, 3]]];

2.4 Modern Physics

Problem 1:  Carbon Dating

Solution

Clear["Global`*"];

n[t_] = n0 Exp[-a t]

eq1 =  1/2 n[0] ==  n[ 5750]

aSol =  Solve[eq1, a]  // Flatten

n[t_] = n[t]  /. aSol

eq2 = Eliminate[{n1 == n[t],  ratio == n1/n0}, n0]

ageRule = Solve[eq2, t]  // First

time[ratio_] =  t  /. ageRule

data =  
        Table[ {ratio,  time[ratio]}, {ratio, 0, 1, 1/10}]  // 
          N;

\(TableForm[\ data\ , 
\ 
        TableHeadings -> {{}, {"Fraction Left"\ , "Age(yrs)"}}, 
TableAlignments -> {Right, 
            Bottom}
]\  // 
      PaddedForm[#, {6, 1}] &\)

t0 = 5750;

    
Plot[ 
        n[t]/n0, {t, 0, 5 t0}, 
AspectRatio -> 
          1, 
Background -> 
          Hue[0.2, 0.5, 
            1], 
Ticks -> {Range[0, 5 t0, t0], 
            Table[1/2^i, {i, 0, 
                3}]}, 
AxesLabel -> {"Age(yrs)", "\
ratio=n1/n0"}, 
PlotStyle -> {Thickness[0.02], 
            Hue[0.9]}, 
GridLines -> {Range[0, 5 t0, 
              t0], Table[
              1/2^i, {i, 0, 5}]}
];

Problem 2:  Stable Isotopes

Solution

Clear["Global`*"];

Part a

Needs["Miscellaneous`ChemicalElements`"]

Remove[Global`Mercury]

?Miscellaneous`ChemicalElements`*

{Elements[[26]],  ElementAbbreviations[[26]]}

StableIsotopes[Fe]

Length[Elements]

StableIsotopes[Fe] - 26

StableIsotopes[Fe]*0 + 26

neutrons[n_] :=  
      StableIsotopes[ElementAbbreviations[[n]] ]       - 
        n
protons[  n_] :=  
      StableIsotopes[ElementAbbreviations[[n]] ]*0  + n

pairs[n_] := 
      Transpose[{protons[n], neutrons[n]}]

    
pairs[26]

\(\(ListPlot[
        pairs[26]\  // Evaluate, 
PlotStyle -> 
          PointSize[
            0.03], 
AxesLabel -> {"proton", "neutron"}, 
PlotLabel -> "Iron (Fe)"
];\)\)

data =  
        Table[ pairs[i], {i, Length[Elements]}]  // 
          Flatten[#, 1] &;

\(\(plot1 = 
        ListPlot[
          data, 
PlotStyle -> 
            PointSize[
              0.015], 
AxesLabel -> {"proton", "neutron"}, 
PlotLabel -> "All Elements", 
AspectRatio -> 
            1, 
DisplayFunction -> $DisplayFunction, 
Epilog -> 
            Line[{{0, 0}, {100, 100}}]
];\)\)

Show[plot1, PlotRange -> {{20, 30}, {20, 40}}];

Part b

Select[
      data , (OddQ[#[[1]]]  &&  OddQ[#[[2]]]) &]

Select[data , (OddQ[#[[1]]]  &&  OddQ[#[[2]]]) &]  // 
      Length

Select[data , (EvenQ[#[[1]]]  &&  EvenQ[#[[2]]]) &]  // 
      Length

Length[
      Cases[data, {x_, y_ } /; OddQ[y]  &&   EvenQ[x]] ]

Length[
      Cases[data, {x_, y_ } /; EvenQ[y]  &&   OddQ[x]] ]

Problem 3: The Bohr Atom

Remarks and Outline

Solution

Clear["Global`*"];

Part a

coulomb = k  ((Z (e^2))/(r^2));
force = m a;
\(\(\(a = v\^2\/r;\)\(
\)
    \)\)
eq1 = coulomb == force

angularMomentum = m v r;

eq2=  angularMomentum == n ℏ

rule1=Solve[{eq1,eq2},{v,r}] //First

kineticEnergy = (1/2) m (v^2);
potentialEnergy = - k (( Z (e^2))/r);
totalEnergy = kineticEnergy + potentialEnergy

energy[n_]= totalEnergy /.rule1 //Simplify

Needs["Miscellaneous`PhysicalConstants`"];

values = {
e -> 
            ElectronCharge, 
m -> 
            ElectronMass, 
k -> 
            1/(4 π ϵ0), 
ϵ0 -> VacuumPermittivity, 
ℏ -> 
            PlanckConstantReduced, 
Z -> 
            1};

Needs["Miscellaneous`Units`"];

energyEV[n_]:= Convert[energy[n] //.values,ElectronVolt];

energyEV[1]  // NumberForm[#, 3] &

Plot[
        Table[ energyEV[n]/ElectronVolt, {n, 10}]  // 
          Evaluate, 
{x, 0, 
          1}, 
PlotRange -> 
          All, 
Axes -> {False, True}];

Part b

eq3 = ℏ ω == energy[(Subscript[n, i])] - energy[(Subscript[n, f])]

Solve[eq3, ω]  // First  // Simplify

Problem 4:  Relativistic Collision

Remarks and Outline

Solution:

Clear["Global`*"];

p1 = m1 γ1 {1, 0, 0, β1};
p2 = m2 {1, 0, 0, 0};

p3 = m3 γ3 {1, 0, 0, β3};

eq1 = p1 + p2 == p3 // Thread

eq2= {γ(3^2) == (1/(1 - β(3^2))), γ(1^2) == 
              (1/(1 - β(1^2)))};

eq3 = Join[{eq1, eq2  }]  // Flatten

sol = Solve[ 
            eq3, {β3, 
              m3 }, {γ3 } ][[2]]//Simplify

2.5 Exercises

Exercise 1: Exploding Projectiles

Exercise 2: Intersecting Trajectories (Shoot the Monkey)

Exercise 3: Principle of Galilean Invariance

Exercise 4: Elastic Collision

Exercise 5: Projectile and Earth's Rotation

Exercise 6: Projectile Thrown up an Incline Plane

Exercise 7: Falling Object with Air Resistance as Observed  from the Rotating Earth

Exercise 8: Projectile Thrown up with Air Resistance

Exercise 9: Quadrupole

Exercise 10: Motion of a Charged Particle

Exercise 11:  (Overscript[B, ⇀])  Field of a Spinning Disk