(*
This notebook demonstrates some of the capabilities of Mathematica. Cells are grouped into sections. Each section is started off by a purple coloured cell which contains a brief description of the contents.

Very likely, when you open this notebook in Mathematica, only this cell will be showing. You can open the notebook by clicking on the outermost grouping bar (thin blue line) on the right of this cell,then pulling down the menu item "Cell" from the action menu bar on top, selecting the submenu "Grouping", and then clicking on "Open all subgroups".

This notebook was created by Sourendu Gupta (sgupta@tifr.res.in). I would appreciate hearing from you if you use this for teaching or have suggestions for improvements.
*)

(* 
This section is an overview of the kind of things you can do with Mathematica: do calculations, control the precision, solve equations, define and use functions, use lists, make plots, get help, etc.

Work your way through these cells, evaluating each cell in order.

General instructions: to evaluate, click on a cell and press [shift][enter]
*)

(* A simple calculation *)
2+2

(* Several calculations can go into a single cell *)
2+2(1+1)
2+2(1+1)^2
(2+2(1+1)^2)/5

(* Solve a simple equation: note the difference between '=' and '==' *)
fahrenheit = 32 + 9 * celsius /5
Solve[ fahrenheit == celsius, celsius]

(* Defining a function:
Note the difference between '=' and ':='
Note also that Mathematica is object-oriented, the function "Celsius" applies to integer, real and complex numbers, as well as lists.
*)
Celsius[fahrenheit_] := 5( fahrenheit -32)/9 (* Defining a function *)
Celsius[-40] (* Function applied to an integer *)
Celsius[98.4] (* Function applied to areal number *)
Celcius[33.0+I] (* Function applied to a complex number *)
fahrenheitvalues = {17, 19, 23, 31, 37} (* Defining a list *)
celsiusvalues=Celsius[fahrenheitvalues] (* Function applied to a list *)

(* Numerical evaluation: controlling the precision *)
N[celsiusvalues]      (* The default numerical evaluation *)
N[celsiusvalues,2] (* Evaluation with 2 digit precision *)
N[celsiusvalues,30] (* More precision *)


?N  (* Help on the built in function N *)

(* Making tables *)
conversiontable = Table[{fahrenheitvalues[[i]], celsiusvalues[[i]]}, {i,Length[fahrenheitvalues]}]
MatrixForm[conversiontable]


?Table (* Help on the built in function Table *)

(* Plotting *)
ListPlot[conversiontable, PlotStyle -> {Red}, Joined -> True]

?ListPlot  (* Help on the built in function ListPlot *)

(*
This section demonstrates how Mathematica can be used to manipulate algebraic expressions, evaluate derivatives, and perform logical operations based on evaluation of expressions.
*)

(* Simple algebra *)
a=(x+y)^4 (* Define an expression *)
Expand[a /. x -> 1] (* A substitution and an expansion *)
Factor[a - y^4] (* Factorization works *)
Plot3D[Tanh[a],{x,-1,1},{y,-1,1}] (* 3-dimensional plots *) 


?Expand
?Factor
?Plot3D

(* Derivatives, extrema and conditional statements *)
expr = x^4 + 2 x^3 - x^2 -2 x 
derv = D[expr, x]
extrema = Solve[ derv == 0, x]
derv2 = D[ expr, {x,2}]

?D

(* Iteration and conditional statements *)
Do[ If [ N[derv2 /. extrema[[i,1]]] < 0 , Print[extrema[[i,1]], "  Maximum"], Print[extrema[[i,1]],"  Minimum"]], {i,3}]

?Do
?If

(* Checking the result of the above manipulations *)
Plot[expr,{x,-5/2,3/2}]

(* 
This section shows how a computation can be developed into into a program, using the example given above of finding the extrema of a polynomial. The main new Mathematica structure used is the "Module".

Here is the program logic. 
First check whether the input expression is indeed a polynomial in the given variable.
Next check that it is not a constant.
If both tests are passed then repeat the manipulations of the previous section.
*)

(* The function definition *)
Extremaof[poly_,var_] := 
If[ PolynomialQ[poly,var],

(* This is executed when "poly" is a polynomial in "var" *)
If [Exponent[poly,var] ≤ 1,

(* This is executed when "poly" is constant in "var" *)
Print[poly, " has no finite extrema in ", var],

(* This is executed when "poly" is a non-trivial polynomial in "var" *)
Module[{derv, derv2, xtrma, nxtrm, i}, (* These variables are internal to the module *)
Print[poly, " is a non-trivial polynomial in ", var];
derv=D[poly,var]; (* The derivative *)
derv2=D[derv,var]; (* The second derivative *)
xtrma=Solve[derv == 0, var]; (* Find the extrema *)
nxtrm=Length[xtrma];(* Number of extrema found *)
Do[ 
If [ N[derv2 /. xtrma[[i,1]]] < 0 , 
(* This is executed when the 2nd derivative is negative *)
Print["\t Maximum:\t", xtrma[[i,1]]], 
(* This is executed when the 2nd derivative is positive *)
Print["\t Minimum:\t", xtrma[[i,1]]]
] (* End of innermost If *), 
{i,nxtrm}] (* End of Do *)
] (* End of Module *)
] (* End of second If *),
(* This is executed when "poly" is not a polynomial in "var" *)
Print[poly, " is not a polynomial in ", var]
] (* End of outermost If *)

(* Now use the function *)
Extremaof[x^2+w Tanh[y], w]
Extremaof[x^2+w Tanh[y], x]
Extremaof[x^2+w Tanh[y], y]
Extremaof[x^2+w Tanh[y], z]
Extremaof[x^5+2x^4 -x-5,x]

(* 
Exercise: there is still one problem with the function, as you saw above. Fix it. 
*)

(*
The purpose of this section is to demonstrate some matrix techniques built into Mathematica
*)

(* Define a matrix and a function on the matrix *)

HermiteanConjugate[m_] :=Conjugate[Transpose[m]]
m = {{0,1,0,0},{0,0,1,0},{0,0,0,1},{1,0,0,0}}  (* A 4-by-4 matrix *)
Det[m] (* The determinant *)
Tr[m] (* The trace *)

(* Some matrix manipulation and display *)

MatrixForm[m] (* Writing the matrix in the usual form *)
md=HermiteanConjugate[m]; (* The Hermitean conjugate: the semicolon stops display of results *)
MatrixForm[md] (* Writing the matrix in the usual form *)
MatrixForm[m.md-md.m] (* The matrix commutes with its Hermitean conjugate *)
MatrixForm[m.md] (* In fact, the matrix m is unitary *)

(* A vector and matrices acting on vectors: m permutes the components *)

v={va,vb,vc,vd} (* A vector *)
m.v (* The matrix acting on the vector *)
md.v (* The conjugate matrix acting on the vector *)
v.m.v (* A quadratic form *)

(* Eigensystems: the eigenvalues of an unitary matrix have modulus 1 *)

l=Eigenvalues[m] (* The eigenvalues of the matrix *)
u=Eigenvectors[m] (* The eigenvectors of the matrix: the semicolon stops display of results *)
MatrixForm[u] (* The eigenvectors of the matrix *)
MatrixForm[u.HermiteanConjugate[u]] (* The eigenvectors of the matrix are not normalized *)
eig=Eigensystem[m] (* The eigenvalues and eigenvectors of the matrix *)

(* Help on the functions used *)

?Det
?Tr
?MatrixForm
?Transpose
?Conjugate
?Eigenvalues
?Eigenvectors
?Eigensystem

(* Every function is overloaded: sometimes that is unfortunate *)

msq=m.m; (* The square of m: both unitary and Hermitean *)
MatrixForm[msq]
MatrixForm[Exp[msq]] (* This is NOT the exponential of a matrix *)

(* This is the definition of a matrix exponential: rotate into diagonal form, exponentiate the diagonal entries, rotate back. There is a built-in Mathematica function which does the same thing; it is called MatrixExp. *)
Matrixexp[m_] := Module[{u,l,n=Length[m]},
 l=Eigenvalues[m];
u=Eigenvectors[m];
ud=HermiteanConjugate[u];
Table[
Sum[ud[[i,k]] Exp[l[[k]]] u[[k,j]],{k,n}],
{i,n},{j,n}]
]

(* Check *)
MatrixForm[Matrixexp[msq]]

(*
This section demonstrates how to solve differential equations and to examine the solutions.
*)

(* The equations for a simple harmonic oscillator *)
soln = DSolve[
 { x'[t] ==  p[t], 
p'[t] ==  -x[t], (* Two differential equation *)
x[0] == 0,
p[0] == 1  (* Two initial conditions *)
},
{x[t], p[t]}, (* The functions to be solved for *)
t (* The independent variable *)
]

(* Plot the solutions as functions of t *)
Plot[{x[t]/. soln ,p[t] /. soln}, {t,0,2 Pi}, Filling -> Axis, AxesLabel -> Automatic]

(* Plot the trajectory in phase space *)
ParametricPlot[{x[t],p[t]} /. soln, {t,0,2 Pi}, PlotStyle -> Thick, AxesLabel -> {x,p}]


(* Is the energy conserved? *)
Plot[(x[t]^2 + p[t]^2)/2 /. soln, {t,0,2 Pi}, AxesLabel-> {time, energy}, PlotStyle -> Thick]

(* Help on every new function used here *)

?DSolve
?ParametricPlot