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. 

 *)

Contains correction on:  
    2nd Edition:Chpt 1, Section 1.3.6.  Page 22.  rootPlot (FIO, 20 Sept. 2006)

1. Getting Started

1.1 Introduction

1.1.1  Computers as a Tool

1.1.2 A Note About Notation and Style

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

1.1.3 Notation and Symbolize 	(For experts only)

(Subscript[x, 1]) // FullForm

term  = a  (Subscript[x, 1])   +  b (Subscript[x, 2])  +  c x

term  /. {x -> 1}

<< Utilities`Notation`

    
?Symbolize

Symbolize[(Subscript[x, 1])]

(Subscript[x, 1]) //FullForm

Symbolize[(Subscript[x, 1])]
Symbolize[(Subscript[x, 2])]
term2 =a(Subscript[x, 1]) + b (Subscript[x, 2]) + c x;

term2  /. x -> 1

1.2 Arithmetic and Algebra

1.2.1 Arithmetic and Notation

2*6 (12 + 4)/(  12 - 4)

{2^3, Sqrt[2], 4!, 4!!}

N[Sqrt[2]]

Sqrt[2] // N

2  // Sqrt   // N

N[Sqrt[2],50]

Sqrt[2] // N[#,50]&

list1={I,  π , E, Infinity, GoldenRatio, EulerGamma, 
        Catalan, Degree}

list1 // N

Clear[list1]

list1

{Sqrt[-1],I,I}

(4 I + 6 I)/(2 - 4I)

{Abs[2 + I], Re[2 + I], Im[2 + I], Conjugate[2 + I] }

1.2.2 Algebraic Manipulations

eq1= Expand[(1+ x^2)^5]

eq1 // Factor

Factor[eq1,GaussianIntegers->True]

Options[Factor]

eq1= Expand[(x+1)^5 (y+1)^2 , y]

eq1 // Simplify

Sqrt[(a b c d Cos[θ])^2]  // Simplify

Simplify[ 
      Sqrt[(a b c d Cos[θ])^2], { a > 0, b ∈ Reals, 
        c ∈ 
          Complexes,  π /2 < θ  < π  }]

Simplify[{Cos[ n  π ], Sin[n π ]},  
      n ∈ Integers]

Factor[Cos[4 x]^2-2 Cos[4 x] Sin[3 x]+Sin[3 x]^2]

Factor[Cos[4 x]^2-2 Cos[4 x] Sin[3 x]+Sin[3 x]^2,Trig->True]

Cos[4 x]^2 - 2 Cos[4 x] Sin[3 x] + Sin[3 x]^2  // 
      TrigFactor

ComplexExpand[Exp[I θ]]

ComplexExpand[Exp[I θ],θ]

Options[ComplexExpand]

ComplexExpand[Exp[I θ], θ,  
      TargetFunctions -> Abs]

1.2.3 PowerExpand

{Sqrt[a^2], Log[a/b], Log[a^b]}  // PowerExpand

Solve[Sqrt[x^2] == a, x]

Solve[Sqrt[x^2] == a  // PowerExpand, x]

1.2.4 Simple Rules

{a, b} /. { a ->  b,  b -> c}

{a, b} //. { a ->  b,  b -> c}

a //.{a -> b, b -> a}

1.2.5 A Home-made Complex Conjugate with SuperStar

2 I /.{I  ->  -I}

2 I // FullForm

{I  ->  -I} // FullForm

conjugateRule= {Complex[re_,im_]->Complex[re,-im]};

2 I /. conjugateRule

Unprotect[SuperStar];

\(SuperStar\  /: \ \(exp_\^*\) := 
      exp\  /. {Complex[a_, b_] :> 
            Complex[a, \(-b\)]}\)
Protect[SuperStar];

\(\(\((a + I\ b + E\^\(I\ θ\))\)\^*\)\)

SuperStar[
      a + I b + (E^(I θ))]

Conjugate[
      a + I b + (E^(I θ))]

1.2.6 Immediate and Delayed Substitutions

{x, x,  x}  /. {x -> Random[]}

{x, x,  x}  /. {x :> Random[]}

1.2.7 Selecting Parts of Expressions

eq1 =  a1 x +  a2 x^2  +  a3 x^3 Sin[x/y];

Part[eq1, 3]

eq1[[3]]

eq1[[3, 1]]

eq1[[{3, 1}]]

position=Position[eq1,Sin[x/y]]

ReplacePart[eq1, Cos[x/y] ,position]

1.2.8 Algebraic Equations

eq1= a x^2 + b x + c == 0;

sol1= Solve[eq1,x]

x /. sol1[[2]]

eq1 /.sol1 //Simplify

Roots[eq1, x]

Reduce[eq1, x]

eq2 = {a x + b y == c, a x - b y == d};

eq3 = Reduce[eq2, {x, y}]

eq3[[4, {1, 2}]]

eq3 = 
      Reduce[{eq2, a ≠ 0, b ≠ 0} // Flatten, {x, 
          y}]

1.3 Functions and Procedures

1.3.1 Built-In Functions

list1={Sin[x], ArcCos[x], Cosh[x], ArcSinh[x], BesselJ[n, x]}

list1  /. {x -> 0.5, n -> 2}

?Sin

??Sin

list = {0, 45, 90, 180, 270 , 360}*Degree;

Sin[list]

1.3.2 User-Defined Functions

f1[x_]  =  x^2;

f1[2]

f2[x] =  x^2;
{f2[x], f2[2]}

a = 1;
f3[x_]  :=  
        If[x ≤ a, -1, 1];
f4[x_]    =  If[x ≤ a, -1, 1];

{f3[2], f4[2]}

a = 3;
{f3[2], f4[2]}

?f3

?f4

Clear[a]

Clear[fac]

fac[1] = 1;
fac[n_Integer?Positive]  :=  n fac[n - 1]

{4!,  fac[4], fac[4.0], fac[-4], fac[a]}

1.3.3 Pure Functions

f1 =  Function[{x, y},  Sin[x /y]];
f1[a, b]

f2 = (Sin[#1/#2]) &;
f2[a, b]

eq1  =  x'[t]  +  a x[t];
eq1  /. {x[t] -> y[t] + b}

eq1  /. {x -> (y[#] + b &)}

term1 =  (1 + x)^2 (1 + y)^2  // Expand

Select[term1, FreeQ[#, x] &]

Select[ term1, ! FreeQ[#, x] &]

1.3.4 Assigning Rules and Restrictions to Functions

i /: i^2  = -1;
i /: i^3  = -i;
i /: i^4  =  1;

i /: i^n_Integer := i^Mod[n,4] /; n<0 || 4<n;

Table[i^n == I^n,{n,-8,8}] /.{i->I}

Clear[step1, step2]

step1[x_]= 0 /; x<=0;
step1[x_]= 1 /; 0< x;

step2[x_ /; x<=0]= 0;
step2[x_ /; 0< x]= 1;

step1[-0.1]  // Trace

step2[-0.1]  // Trace

step3[x_]:=Which[x<=0,0,x>0,1]

step3[-0.1]  // Trace

{Table[ step1[x], {x, -10000, 10000, 1}]  // Timing  // 
        First, 
Table[ 
            step2[x], {x, -10000, 10000, 1}]  // Timing  // 
        First, 
Table[ 
            step3[x], {x, -10000, 10000, 1}]  // Timing  // 
        First}

1.3.5 Module

shift[list_] := 
      Module[{newList, a = 1, max}, 
newList = 
          list + a; 
max = 
          Max[newList]; 
Print["The max of the new list \
is: ", max]; 
Return[
          newList]
]

shift[{1, 2, 3}]

{newList, a, max}

shift2[list_, a_:  1] := 
      Module[{newList, max}, 
newList = 
          list + a; 
max = 
          Max[newList]; 
Print["The max of the new list \
is: ", max]; 
Return[
          newList]
]

shift2[{1, 2, 3}, 2]

shift2[{1, 2, 3}]

1.3.6 rootPlot   (corrected; 20 Sept. 2006)

eq1= 0 == 1 + x + 2 x^2 + 3 x^3 + 4 x^4 + 5 x^5;

NSolve[eq1,x]

(list = {Re[x], Im[x]}  /. NSolve[eq1, x]) // 
      TraditionalForm

ListPlot[list, PlotStyle -> PointSize[0.03]];

rootPlot[equation_, variable_:  x, opts___] := 
        Module[{nSol, list}, 
nSol = 
            NSolve[equation, 
              variable]; 
list =  {Re[x], Im[x]}  /. 
              nSol; 
Return[ 
            ListPlot[list, PlotStyle -> PointSize[0.03], 
              opts]];
]

rootPlot[equation_, variable_:  x, opts___]:=Module[{nSol, list},
nSol = NSolve[equation, variable];
  
list = {Re[variable], Im[variable]} /. nSol;
Return[ 
                ListPlot[list, PlotStyle -> PointSize[0.03], opts]];
]

rootPlot[(x^15) - 1 == 0, x, 
        AspectRatio -> 1];

1.4 Packages

1.4.1 Loading Packages

$Packages

Needs["Graphics`Graphics`"]

$Packages

?Graphics`Graphics`*

?PolarListPlot

Needs["Calculus`VectorAnalysis`"]

CoordinatesToCartesian[{r, θ, ϕ}, Spherical]

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

x^2 + y^2 + z^2 /. x2rRule // Simplify

1.4.2 Contexts

Contexts[]

Context[PolarListPlot]

Context[dummy]

Clear["Global`*"]

Clear["Graphics`Graphics`*"]

dummy[x_]=x;
Protect[dummy];

Clear["Global`*"]

?dummy

1.4.3 Shadowing

EasterSunday[2002]

Needs["Miscellaneous`Calendar`"]

EasterSunday[2002]

Context[EasterSunday]

Remove[EasterSunday]

Context[EasterSunday]

EasterSunday[2002]

1.5 Calculus

1.5.1 Derivatives and Integrals

f[x_] = x^n;

    
{D[x^n, {x, 1}], D[x^n, x], D[f[x], x], f'[x]}

{D[x^n, {x, 2}], D[f'[x], x], f''[x]}

term1 = x^3 y^4 + y x^5 + x/y^3;

    
D[term1, {x, 2}, y]

eq1 = 
      Integrate[(x^a), {x, 0, 1}, GenerateConditions -> True]

eq2 = Integrate[(x^(-1)), x]

eq1 = 
      Integrate[(x^a), {x, 0, 1}, GenerateConditions -> False]

Integrate[(x^a), {x, 0, 1}, 
      Assumptions -> {Re[a] > -1}]

Integrate[ 1, {x, 0, 1}, {y, 0, (Sqrt[1 - (x^2)])}]

NIntegrate[
      Sin[z Sin[x y]], {x, 0, 1}, {y, 0, x }, {z, 0, x  y}]

1.5.2 Differential Equations

eq1 = y''[t] + (ω^2) y[t] == 0;

initial = {y[0] == y0, y'[0] == v0};

eq2= Append[initial,eq1]

eq3 = DSolve[eq2, {y[t]}, t] // Flatten

1.5.3 Changing Variables and Pure Functions

eq1 = a x''[t] + b x'[t] + c x[t] == 0;

eq1 /. {x[t] -> y[t] + 2}

eq1 /. {x -> (y[#] + 2 &)}

eq1 /. {x -> Function[{t}, y[t] + 2]}

DSolve[{x[t] == y[t] + 2}, x, t]  // Flatten

1.5.4 Numerical Solutions of Differential Equations

eq1 = {x''[t] == -x[t] - x[t] y[t]^2,  y''[
              t] == -y[t] - x[t] y[t]^2};

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

ndsol = NDSolve[
          Join[eq1, initial], {x, y}, {t, 0, 20}][[1]]

Needs["Graphics`Legend`"]

\(Plot[{x[t], y[t]} /. ndsol\  // Evaluate, {t, 0, 20}, 
      PlotStyle -> {{\ }, {Dashing[{0.02}], 
            Thickness[
              0.008]}}, 
PlotLegend -> {"x[t]", "y[t]"}, LegendPosition -> {1, 0}
]\)

ParametricPlot[{x[t], y[t]} /. ndsol  // Evaluate, {t, 0, 
          20}];

Plot[{x'[t], y'[t]} /. ndsol  // Evaluate, {t, 0, 20}, 
        PlotStyle -> {{ }, {Dashing[{0.02}], 
              Thickness[0.008]}}];

1.6 Graphics

1.6.1 Using the Plot Command

Options[Plot]

\(\(pt1 = \ \ \ Plot[\ 
          1\ /\ \((1 - x^2)\), \ {x, \ \(-10\), \ 
            10}, 
\ \ \ \ \ \ 	PlotStyle -> {\ Hue[0.8], 
              Thickness[0.008]}];\)\)

colors = {red, green, blue, cyan, magenta, 
            yellow} = {RGBColor[1, 0, 0], RGBColor[0, 1, 0], 
            RGBColor[0, 0, 1], CMYKColor[1, 0, 0, 0], CMYKColor[0, 1, 0, 0], 
            CMYKColor[0, 0, 1, 0]};

Plot[Table[i x , {i, 6}]  // Evaluate , {x, 0, 1}, 
        PlotStyle -> colors ];

n = 13;
hues = Table[Hue[i/n], {i, n}];
 Plot[Table[i x , {i, n}]  // Evaluate , {x, 0, 1}, 
        PlotStyle -> hues ];

1.6.2 Animated Plots

Needs["Graphics`Animation`"]

frames=1;
Animate[
	Plot[Cos[2 π (x-t)],
	{x,-2,2},
	PlotRange->{-1,1},
	Axes-> False],
{t,1/frames,1,1/frames}]

frames = 1.1;
Animate[
      Plot[{Cos[2  π (x - t)], Cos[2  π (x + t)], 
          Cos[2  π (x - t)] + Cos[2  π (x + t)]}, {x, -2, 
          2}, 
PlotRange -> {-2, 
            2}, 
PlotStyle -> {blue, 
            green, {red, 
              Thickness[0.008] }}, 
Axes -> 
          False], 
{t, 1/frames, 1, 1/frames}]

1.6.3 Vector Field Plots

vector = {1/x, 1/y};

Needs["Graphics`PlotField`"];

PlotVectorField[
        vector, {x, -5, 5, 2}, {y, -5, 5, 
          2}, 
ScaleFunction -> (1 &), 
        PlotPoints -> 20];

Needs["Graphics`PlotField3D`"];

PlotGradientField3D[
        1/Sqrt[x^2 + y^2 + z^2], 
{x, -5, 5, 
          2}, {y, -5, 5, 2}, {z, -5, 5, 
          2}, 
ScaleFunction -> (1 &), 
VectorHeads -> True];

1.6.4 Three-dimensional Graphics using Plot3D and ParametricPlot3D

\(\(Plot3D[
        Sqrt[x + I\ y] // Im, {x, \(-1\), 1}, {y, \(-1\), 
          1}, 
PlotPoints -> 
          25, 
AspectRatio -> 
          1, 
AxesLabel -> {"Re", "Im", ""}];\)\)

Needs["Graphics`ParametricPlot3D`"]

CylindricalPlot3D[
        Sqrt[r Exp[I θ]] // Im, {r, 0, 
          10}, {θ, -π, π}, 
PlotPoints -> 25, 
AspectRatio -> 1];

SphericalPlot3D[SphericalHarmonicY[2, 0, θ, ϕ]  // 
            Abs // Evaluate, 
{θ, 
          0, π}, {ϕ, 0, 2  π}, PlotPoints -> 50, 
        Boxed -> False, Axes -> False];

ParametricPlot3D[{Cos[2t+π/4]+t/100,2 Cos[2 t]+t/200,Cos[3 t]},{t,0,100}, PlotPoints-> 3000];

1.7 Exercises

Exercise 1: Real Expressions

Exercise 2:  Series

Exercise 3: Integration

Exercise 4: Superposition of Waves

Exercise 5: Tables of Zeros for Several Functions

Exercise 6: Plot of Roots

eq1= x^7 + x^3 + x^2 + 1==0;

Exercise 7: Coupled Equations

Exercise 8: Statistical Data

r = Table[Random[Integer,{1,20}], {1000}];

Exercise 9: Differential Equation

Exercise 10: Coupled Differential Equations

Exercise 11: Trajectories and Vector Field Plots for Coupled Differential Equations

Exercise 12: Heat Equation

Exercise 13: Wave Equation

Exercise 14: Fit to Data

data = 
        Table[{x, 5 Sin[3 π  x] + Random[Real, {-0.2, 0.2}]} // 
              N // Chop[#, 10^-3] &, {x, 0, 0.5, 0.05}];

Exercise 15: Rules used to solve problems

list = Random[Integer , {0, 100}] & /@ Range[13];

Exercise 16: Procedures to select numbers

list = Table[Random[Real, {0, 20}], {i, 30}];

Exercise 17: Procedure to Select certain Pairs of Points

pairs=Table[Random[],{1000},{2}];

Exercise 18: Change

Exercise 19: Toss of a Die

Exercise 20: Three Dice

Exercise 21: Playing Darts and the Calculation of π

list = 
        Table[{Random[Real, {-1, 1}], Random[Real, {-1, 1}]}, {i, 
            10}];

Exercise 22: Dropping Needles

Exercise 23: Pacal's Wager

pairs = Table[Random[Integer, {1, 6}], {1000}, {2}]

Exercise 24: The birthday problem