Introduction to Mathematica 
Luca Amendola, SS2025

Thanks to Prof. Miguel Quartin (CNBPF-Brazil) for several contributions to this course

Course web page:
https : // www . thphys . uni - heidelberg . de/~amendola/intromath - ss2025 . html

Lecture 1

First Law: Put the cursor in the cell and enter Shift - Return to execute it!

Second Law: always start with:

ClearAll["Global`*"];SetDirectory[NotebookDirectory[]]

Quick facts about Mathematica

First release by Stephen Wolfram in 1988. Current version 14.2
Composed of a Notebook (this graphical interface) and a Kernel (which does the computation and is automatically launched at the first operation)
More than 6000 built-in functions, from Sin[x] to InverseContinuousWaveletTransform[cwd], from AASTriangle to ZTransform
The first manual of Mathematica weighed 3kg!
Many more community-provided functions at Wolfram Function Repository
Support for automatic Parallel Computing
Main advantage: numerical and symbolic manipulations, formatted text, and plots, all in the same interface; intuitive
Main disadvantage: it’s not open source and it’s not free

Start!

Everything is a function

Assignment

a=2

Now a is always 2 until you clear or change assignment

a

Clear[a];a

Use ";" if you don' t want to see the output

a=2;
a
a=3;
Clear[a];

Functions

f[x_]=x^2

or also (use Palettes -> Basic Math Assistant) or ctrl 6

f[x_]=(x^2)

The argument x_ is a "slot", i.e. it applies to everything you insert there

f[y]
f[a*b]
f[{a,b}]

But if you forget the "_" then it is only applied to the named argument

Clear[f];f[x]=x^2;
f[x]
f[y]

Three or four different "="

With “=” the definition is executed immediately (immediate assignment) :

Clear[f,x];x=3;f[x_]=x^2;
f[x]
f[y]
f[4]

With “:=” the definition is executed when the function is called (delayed assignment)

Clear[f,x];x=3;f[x_]:=x^2;
f[x]
f[y]
f[4]

another example

Clear[f,x];f[x_]:=Integrate[x^3,x]

f[y]

However, be careful with delayed assignments:

f[1]

If there' s no analytical form, it just gives back the output

f[x_]:=Integrate[Exp[Sin[x]],x];f[x]

but it' s still doing its job...

D[f[x],x]

The expression /. x -> y is a rule: means replace every x with y (local assignment)

Clear[f,x];f[y_]:=Integrate[x^3,x]/.x->y
f[z]

Notice also

f[1]
f[1.]

So typically : 
use "=" for a general/immediate assignment; 
use “x->y” for a local assignment;
use ":=" for functions

Some more functions :

f[x_]:=Cos[(x^2)];
f[Pi]
g[x_]:=Integrate[y^3,{y,0,x}];
g[Pi]

Yet another "=" : Boolean test

Clear[f];f[x_,a_,b_]:=If[a==b,x^2,x^3]
f[z,a,b]
f[z,3,4]

Boolean operators : && (AND), || (OR) , != (Not Equal To) , >= , < , <= , == ...

Clear[f];f[x_,a_,b_]:=If[a≥ b && (b>0 || b==2),x^2,x^3]
f[z,5,4]
f[z,5,-1]

The output of a test is just True or False

2==3

If and Which

f[x_]=If[x>3,x,0]
g[x_]=Which[x>3,x,x>0,1,x>-2,0]

g[2.5]

Numeric vs symbolic

f1[x_]=Integrate[Cos[x],x]

f2=NIntegrate[Cos[x],{x,0,1}]

... but this doesn' t work:

f3[x_]=NIntegrate[Cos[y],{y,0,x}]

... because everything must be numeric now:

f3=NIntegrate[Cos[x],{x,0,3}]

...but it is sufficient that is numerical when you call it:

f4[y_]:=NIntegrate[Cos[x],{x,0,y}];
f4[Pi/4]
f4[Pi]
Chop[f4[Pi]]

Clearly, sometimes there is no analytical solution :

f5[y_]=Integrate[Exp[Cos[Log[x]]],{x,1,y}]

Derivatives :

f1[x_]=Integrate[Cos[x],x];
f6[y_]=D[f1[y],y]

f6[y_]=D[f5[y],y]

sometimes Mathematica needs a little help :

NIntegrate[a*Cos[x],{x,0,1}]
a*NIntegrate[Cos[x],{x,0,1}]

Multidimensional integrals and partial derivatives

f[x_,y_]:=Integrate[Cos[x]*Sin[y],x,y];
f[x,y]

f[a_,b_]:=NIntegrate[Cos[x]*Sin[y],{x,0,a},{y,0,b}];
f[2,3]

f[3,5]

f[x_,y_]:=Integrate[Cos[x]*Sin[y],x,y];
D[f[x,y],x,y]
D[f[x,y],x,x,x]

Generic expression for partial derivatives

D[g[x,y],x]

D[g[x,y,z,w],x,y]

Cell types (this is a subsection)

Input, Output, Text, Section, sub - section, etc.

This cell is in Text Format, and cannot be executed

This is an Input cell

"This is an Input cell"

Mathematica assumes each word is a variable, and by default orders them alphabetically

((%)^2)

(* this is a commment that is not printed *)
Print["this is a commment that is printed"]
"this is a commment that is printed"
1+1

Brackets and Capital Letters

Third Law: ALL Mathematica built - in functions start with Capital Letters
If you define yours with smallcase first letter, you are guaranteed not to have conflicts of definitions

sin=3;

Sin=3;

Palettes

Palettes are helpful for beginners
You can write equations using a graphical interface (similar to, say, Office programs)

(Sqrt[(Subscript[a, x]^2)])

(Sqrt[π E 8.1])

Pi

1*Pi

N[Pi,40]

Clear[x]
(∂_x)(x^2)

ψ λ

Space = multiplication

a*b
a b
3*2
3 2

Exercise :
write a function with delayed assignment containing Sin, Exp, and Sqrt; differentiate it; integrate it numerically in some range

A function for every need

3*I

Some functions

Ceiling[7.4]
Floor[8.1]
Conjugate[I]
Sqrt[16]
Mod[17,4]
Random[]
RandomInteger[{1,100}]
RandomVariate[NormalDistribution[0,1],3]

very useful to build lists of random numbers from various distributions :

Table[RandomVariate[NormalDistribution[0,1]],{i,10}]

RandomVariate[NormalDistribution[0,1],100]

More functions

Re[2+3*I]
Im[2+3*I]
Im[a+I*b]
Log10[100]
LegendreP[2,x]
PrimeQ[13]

PrimeQ[20]

Table[{6*n-1,6*n+1},{n,1,10}]

PrimeQ[Table[{6*n+1,6*n-1},{n,1,100}]]

Do[Print[6*n+1,":",PrimeQ[6*n+1],"; ",6*n-1,":",PrimeQ[6*n-1]],{n,1,20}]

Simplification: one of the most important commands

D[Integrate[1/(x^3+1),x],x]

Simplify[%]

Mathematica does not assume symbols to represent real variables : they might also be complex numbers. So sometimes an obvious simplification in the real domain is not carried out:

Simplify[Sqrt[x^2]]

It helps then to specify what we know about the symbols we use:

assume={Im[x]==0};
Simplify[Sqrt[x^2],Assumptions->assume]

assume={Im[x]==0,x>0};
Simplify[Sqrt[x^2],Assumptions->assume]

Hotkeys (shortcuts)

a*b/.b-> β

Press "ESC-key expression ESC-key" to get a compact form of the expression
E.g.:  
ESC beta ESC -> β ;   
ESC x ESC (or ESC xi ESC) -> ξ ;   
ESC X ESC -> Ξ ;   
ESC Gamma ESC (or ESC G ESC) -> Γ;   
ESC m ESC -> μ
ESC pd ESC -> ∂
ESC inf ESC -> ∞

"CTRL 6" -> enter exponentiation power 
"CTRL 2" -> sqrt
“CTRL /“ -> start fraction
"CTRL -" -> subindex (typesetting only, no computational value)

Strings

text1="This";
text2=" is";
text3=" text";
sentence=text1<>text2<>text3

"1"==1
ToString[1]
ToExpression["1"]
ToString[1]=="1"

Many operations on strings : StringTake, StringDrop, StringSplit, StringCount, etc

StringReplace[sentence,"t"->"p"]

Clear[j];Do[{Print["i=",i," , ",j],j=i},{i,5}]

Strings can be used as arguments as well:

plot[sx_,sy_]:=Plot[x^2,{x,0,2},Frame->True,FrameLabel->{sx,sy}]

plot["distance","redshift"]
plot["x","y"]

Simple plots

Mathematica can do many different plots : Plot, LogPlot, LogLogPlot, ListPlot, Plot3D, ContourPlot, Histogram etc .

?Plot
?LogLinearPlot

Plot[x,{x,-1,3}]

Plot[x,{x,-1,3},Frame->True]

even plots are functions!

plot[f_]:=Plot[f,{x,-1,3},Frame->True];
p1=plot[x^2]
p2=plot[Sin[x]]

(* combine plots *)
Show[p1,p2]

LogLogPlot[x^5,{x,-5,5},Frame->True]

p1=Plot[{x^2,x^3+1,Sin[x]},{x,-1,2},Frame->True,FrameLabel->{"time","space"},PlotLabel->"Trajectory",PlotStyle->{{Thick,Red},Dotted,{Thick,Dashed}},GridLines->Automatic,PlotLegends->Automatic]

(* you can directly export your plots, or select them and go to File->Save Selection As...*)
Export["newplot.pdf",p1]

Plot3D[Exp[-x^2-y^2],{x,-3,3},{y,-3,3},PlotRange->All]

?Plot3D

ContourPlot[x^2+3*y^2+2*x*y,{x,-3,3},{y,-3,3},Contours->{.5,2,10},PlotPoints->50]

data1=Table[RandomVariate[NormalDistribution[0,1]],{i,5000}];
data2=Table[RandomVariate[NormalDistribution[3,1]],{i,5000}];
Histogram[{data1,data2},{0.25},Frame->True]

hist[n_]:=Histogram[Table[Random[],{i,n}],{0.1}];
hist[30]
hist[10000]

Animations

Animate[Plot3D[Exp[-(x-Cos[a])^2-(y-Sin[a])^2],{x,-3,3},{y,-3,3},PlotRange->All,PlotPoints->30],{a,-Pi,Pi},AnimationRunning->False]

The Gaussian race

Animate[Plot3D[Exp[-(x-Cos[a])^2-(y-Sin[a])^2]+.5*Exp[-(x-2*Cos[-a])^2-(y-2*Sin[-a])^2],{x,-3,3},{y,-3,3},PlotRange->All,PlotPoints->30],{a,-Pi,Pi},AnimationRunning->False]

Main commands this lecture :

Clear,Integrate, NIntegrate, D, If, Which,Table, Random, 
RandomVariate,Simplify,Plot, Plot3D, Show,Animate

Lecture 2

Lists

Main commands :

Length, Head, Total, Sum, Drop, TableForm, TeXForm

a=2;
ClearAll["Global`*"]; (* this clears all variables! *)

a

Lists

Introduction

Arguments of functions are always in  [... ]; list elements are in {...}

mylist={1,1,2,3,5,8,13,21};
Length[mylist]

Specific elemensts of lists are in [[ ...]]

mylist[[4]]
mylist[[-1]]
mylist[[Length[mylist]]]
mylist[[-3]]

(* most functions, but not all, can be applied directly to lists *)

Sin[mylist]
Sin[N[mylist]]

General use of index brackets

(x+1)

(x+1)[[1]]

(* "Head" of an expression *)
((a+b)^3)[[0]]
(x+1)[[0]]

{1,2}[[0]]
{1,2}[[1]]

Head[{1,2,3}]
Head[40]
Head[Sin[x]]

Parts of lists

mylist

mylist[[{2,7}]]
mylist[[2;;7]]

"Table" is the main way to produce a list

rlist=Table[RandomReal[{0,1}],{i,1,5}]

Forms of iterations

Table[i,{i,10}]
Table[i,{i,1,10}]
Table[i,{i,1,10,1}]
Table[i,{i,1,10,2}]
Table[i,{10}]
list={1,2,4,8,9,16};
Table[i,{i,list}]

Operation on lists

mylist={1,1,2,3,5,8,13,21};

Total[mylist]
Sum[mylist[[i]],{i,1,8}]

Timing (important lesson: always try to use built-in functions!)

listsize=(10^6);
mylist3=Table[RandomReal[],{listsize}];

Timing[Total[mylist3]]
Timing[Sum[mylist3[[i]],{i,1,listsize}]]

Matrices

mymatrix1={{1,2},{3,4}};
mymatrix2[dim_]:=Table[Random[],{i,dim},{j,dim}];
mymatrix1
mymatrix1//TableForm (* only for visualization! *)
mymatrix2[5]//TableForm (* only for visualization! *)

mymatrix2[5]

This command can save you LOTS of work:

(x^2+α*3*Sqrt[x])

(x^2+α*3*Sqrt[x])//TeXForm

mymatrix2[5]//TableForm//TeXForm

Parts of matrices

mymatrix1//TableForm

mymatrix1[[1]]
mymatrix1[[1,2]]
mymatrix1[[1,3]]

(* select a row or column *)
mymatrix1//TableForm
mymatrix1[[2,All]]
mymatrix1[[All,1]]

Combining functions and list elements

mymatrix3[factor_]=factor*Table[Random[],{i,5},{j,5}];
mymatrix3[5][[1]]
mymatrix3[5][[1,2]]

?Drop

list={a,b,c,d,e};
Drop[list,3]
Drop[list,{3}]

ClearAll["Global`*"];

Vector and matrix operations

Main commands :

Outer, Cross, Norm, Join, Dimensions, Inverse, Det, 
Transpose, Eigenvalues, Eigenvectors, IdentityMatrix,DiagonalMatrix, ListPlot

?Outer

vector={a,b,c}

vector*vector

vector.vector

mymatrix1={{1,2.},{3,4}};

T_ik = M_ij M_jk

mymatrix1.mymatrix1

3*vector
Sin[vector]
Total[vector]
vector*vector
vector.vector (* scalar product *)
Cross[vector,vector] (* cross product *)
Outer[Times,vector,vector]//TableForm (* outer product T_ij = v_i v_j *)
Norm[vector]
Join[vector,{d,e,f}]

mymatrix=Table[Random[],{i,3},{j,3}];
mymatrix

Dimensions[mymatrix]
Length[mymatrix] 
invmymatrix=Inverse[mymatrix];
ident=(invmymatrix.mymatrix)//TableForm  
Chop[ident]

vector

Transpose[mymatrix]//TableForm
Det[mymatrix]
Eigenvectors[mymatrix]
Chop[Eigenvectors[mymatrix]]
Eigenvalues[mymatrix]
(* v_i M_ij v_j *)
Simplify[vector.mymatrix.vector] (* v_i M_ij v_j*)

Also symbolically :

sm={{a,b},{c,d}};
Det[sm]
Eigenvectors[sm]
Eigenvalues[sm]
Eigenvalues[sm]//TeXForm

Exercise : given an arbitrary 3 x3 random matrix , construct a symmetric matrix by using Transpose ...

Special matrices

IdentityMatrix[4]//TableForm
vector={a,b,c}
DiagonalMatrix[vector]//TableForm

ListPlot

ListPlot[{7,2,5,9}]

list=Table[{i,i^2},{i,5,50,2.5}]
p1=ListPlot[list,Frame->True,PlotStyle->PointSize[0.02]]
p2=ListPlot[list,Frame->True,Joined->True]
Show[p1,p2]

Testing, searching, ordering Lists

Main commands :

Position, Count, Cases, Min, Max, Sort, Ordering, Partition, Flatten, DeleteDuplicates

mytab=Table[i^2-6*i-2,{i,0,14}]

pos=Position[mytab,25]

pos=Position[mytab,25][[1,1]]

mytab[[pos]]

?Count

Count[mytab,-7]

Patterns and tests

f[x_]=x^2

(* "_?Test" this applies the function only if the test gives True*)

p[x_?NumberQ]:=x^2

p[4.5]+p[3/2]+p[u]

mytab

Count[mytab,_?OddQ]

Count[mytab,_?PrimeQ]
Cases[mytab,_?PrimeQ]

Min[mytab]
Max[mytab]
MinMax[mytab]

mytab

Sort[mytab]
(* ordering gives the position according to Sort *)
Ordering[mytab]

(* changing the sorting *)
Sort[mytab,Greater]

?Partition

mytab

Length[mytab]
Partition[mytab,3]
mytabpart=Partition[mytab,5]
Flatten[mytabpart]

mytab
DeleteDuplicates[mytab]

Loops

Main commands :

AppendTo, AddTo, Throw, Catch, Monitor, Solve, NSolve

“Do” loops can perform many operations in each iteration

Do[
{var1=4i+1;
If[PrimeQ[var1], var2=var1,var2=ToString[var1]<>" it's not a prime number"];
Print[var2]},
{i,1,10}]

Table[i,{i,1,10}]

AppendTo: create tables when you do not know how big the table will be:

table={};Do[If[PrimeQ[i],AppendTo[table,i]],{i,100}];table

AddTo : adding elements

sumOfPrimes=0;Do[If[PrimeQ[i],AddTo[sumOfPrimes,i]],{i,100}];sumOfPrimes

Exercise : plot a unit circular region filled with random points

More about loops

Throw and Catch an exception

?Throw

(* factorial *)

3!

f[x_]:=If[x>10,Throw["overflow"],x!]

Catch[Print[f[3]]]
Catch[Print[f[11]]]

f[11]

Monitoring progress

?NSolve

Solve[x^2==a,x]

Clear[x]
tab=Table[NSolve[(x^2)+j x+(j^4)==0],{j,1,10}]

Monitor[
mytab=Table[NSolve[(x^2)+j x+(j^4)==0,x],{j,1,5000}];
,j]

... and make it more readable

imax=10000;
Monitor[mytab=Table[NSolve[(x^2)+i x+(i^4)==0,x],{i,1,imax}];,1.i/imax]

imax=10000;
Monitor[mytab=Table[NSolve[(x^2)+i x+(i^4)==0,x],{i,1,imax}];,100.i/imax"%"]

imax=15000;
Monitor[mytab=Table[NSolve[(x^2)+i x+(i^4)==0,x],{i,1,imax}];,Round[100.i/imax,1]"%"]

Nested Loops

imax=30;jmax=50;
mytab=Table[NSolve[(x^(i+2))+(x^i)+(x^j)-3==0,x]
,{i,1,imax},{j,1,jmax}];
Dimensions[mytab]
mytab[[1,1]]

TableForm and TeXForm again

imax=6;jmax=6;
mytab=Table[i*j,{i,1,imax},{j,1,jmax}];
Dimensions[mytab]

mytab//TableForm

mytab//TeXForm

Arbitrary Precision Arithmetic

(2.^(1/10))
(2^(1/10))

1.*π

N[π,500]

N[Exp[-100]]
N[Exp[-(10^3)]]-N[Exp[-(10^3)-1]]
N[Exp[-(10.^16)]]

Numeric vs symbolic

Clear[x];intsymbolic=Integrate[BesselJ[0,10000x],{x,0,10}]
N[intsymbolic]

Integrate[BesselJ[0,10000.x],{x,0,10}]

Integration parameters:
MaxRecursion, WorkingPrecision, AccuracyGoal, PrecisionGoal...

MaxRecursion is important for fast oscillating functions

Plot[BesselJ[0,10000*x]*BesselJ[1,10001x],{x,0,10}]

nint=NIntegrate[BesselJ[0,10000x]BesselJ[1,10001x],{x,0,10},MaxRecursion->5]
nint=NIntegrate[BesselJ[0,10000x]BesselJ[1,10001x],{x,0,10},MaxRecursion->10]
nint=NIntegrate[BesselJ[0,10000x]BesselJ[1,10001x],{x,0,10},MaxRecursion->15]
nint=NIntegrate[BesselJ[0,10000x]BesselJ[1,10001x],{x,0,10},MaxRecursion->20]

nint=NIntegrate[BesselJ[0,10000x]BesselJ[1,10001x],{x,0,10},MaxRecursion->20,WorkingPrecision->10,AccuracyGoal->10,PrecisionGoal->10]

nint=NIntegrate[BesselJ[0,10000x]BesselJ[1,10001x],{x,0,10},MaxRecursion->20,WorkingPrecision->20,AccuracyGoal->15,PrecisionGoal->15]

With AccuracyGoal->a and PrecisionGoal->p 
Mathematica attempts to make the numerical error in a result of size x be less than 
((10)^(-a))+|x|*((10)^(-p))

Main commands this lecture :

Length, Head, Total, Sum, Drop, TableForm, TeXForm

Outer, Cross, Norm, Join, Dimensions, Inverse, Det, 
Transpose, Eigenvalues, Eigenvectors, IdentityMatrix,DiagonalMatrix, ListPlot

Position, Count, Cases, Min, Max, Sort, Ordering, Partition, Flatten, DeleteDuplicates

AppendTo, AddTo, Throw, Catch, Monitor, Solve, NSolve

Lecture 3

Anything can be a Function!

Simple Functions

f[n_]:=Table[i^2,{i,n}];
f[5]
f[15]

f[a_,b_]:=NIntegrate[x^2,{x,a,b}]
f[1,2]
f[-3,3]

Clear[f];f[a_,b_]:=Plot[x^2,{x,-1,b},Frame->True,PlotLabel->a]
f["what a nice plot!",3]
f["but this is even better!",4]

Pure Functions

g1[x_]:=x+3
(* anonymous functions *)
g2=Function[3+#]
{g1[2],g2[2]}

g2=Function[3+#]

Shorthand notation for pure functions:

g3=(3+#)&

(3+#)&[6]

{g1[2],g2[2],g3[2]}

Pure functions of many variables:

g4=(#1^2+#2^4)&

Applying pure functions to values:

Clear[a];
(#1^2+#2^4)&[a,b]
g4[x,y]

Pure functions are very useful when using built-in functions like Map:

?Map

Map[Sin,{0,1,2,3,π/2,π,I}]

Sin[{0,1,2,3,π/2,π,I}]

Map[(((#)^2)+1)&,{0,1,2,3,π/2,π,x,I}]

Map[Re[#]&,{0,1+4*I,2,π,x,I}]

g5=NIntegrate[#1,{x,0,2}]&;
Map[g5,{x,x^2,Sin[x]}]

g5=NIntegrate[#1,{x,0,#2}]&;
g5[x,2]

g5=NIntegrate[#[[1]],{x,0,#[[2]]}]&;
Map[g5,{{x,1},{x^2,2},{Sin[x],2}}]

Warning! Map applies to elements!

Clear[f];Map[f,a+b+c+d]
Map[f,{a+b+c+d,e+f,11}]

(* order matters *)

g4=(#1^2+#2^4)&;
Map[g4[a,#]&,{0,1,2,3,π/2,π,I}]
Map[g4[#,a]&,{0,1,2,3,π/2,π,I}]

Recursive functions

factorial[x_]:=If[x==0,1,x*factorial[x-1]];
factorial[5]

factorial[100]

Interpolating functions

This is one of the most  time-saving trick!

Clear[x];
prec=12;f[a_]:=NIntegrate[Exp[x]*y,{x,-1,a},{y,-1,a},MaxRecursion->prec,AccuracyGoal->prec,PrecisionGoal->prec,WorkingPrecision->prec]
Timing[f[3]]

Timing[Plot[f[a],{a,1,10}]]

Timing[Plot[f[a],{a,1,10},PlotPoints->5]]

Table[{a,f[a]},{a,1.,10.,1}]

?InterpolationOrder

fint=Interpolation[Table[{a,f[a]},{a,1.,10.,1}],InterpolationOrder->3]

fint[9.4]

(* see what happen by increasing "InterpolationOrder" *)

Timing[Plot[fint[a],{a,1,10}]]

An interpolating function is a full function :

fdint[x_]:=D[fint[a],a]/.a->x

fdint[3]

Choosing the interpolation scheme

Range[0,10]/10

{Range[0,10]/10,RandomReal[1,11]}//TableForm

Transpose[{Range[0,10]/10,RandomReal[1,11]}]//TableForm

data=Transpose[{Range[0,10]/10,RandomReal[1,11]}];
sp=Interpolation[data,Method->"Spline"]
he=Interpolation[data,Method->"Hermite"]

GraphicsRow[{Show[Plot[{sp[x],he[x]},{x,0,1},PlotRange->All],ListPlot[data,PlotStyle->{Red,PointSize[.02]}]],Plot[{sp'[x],he'[x]},{x,0,1},PlotRange->All]}]

data=Transpose[{Range[0,10]/10,RandomReal[1,11]}];
sp=Interpolation[data,InterpolationOrder->1]
he=Interpolation[data,InterpolationOrder->3]

GraphicsRow[{Show[Plot[{sp[x],he[x]},{x,0,1},PlotRange->All],ListPlot[data,PlotStyle->{Red,PointSize[.02]}]],Plot[{sp'[x],he'[x]},{x,0,1},PlotRange->All]}]

Exercises :

 1) write down a polynomial of order 4 with random real coefficients between 0, 1. 
 Select 10 random points between 0 and 10 from this polynomial. Interpolate the points with increasing 
InterpolationOrder and plot the results .
  
2) Solve x^3 + a*x - 2 = 0 with NSolve, for a integer from 0 to 10. Use Flatten to create a one-dimensional list. Create a table of the real roots, interpolate them and plot them

Solution

tab=Flatten[Table[NSolve[x^3+a*x-2==0,x],{a,0,10,1}]]
rtab={};Do[If[Im[tab[[i,2]]]==0,AppendTo[rtab,tab[[i,2]]]],{i,Length[tab]}];
rtab
ir=Interpolation[rtab];Plot[ir[x],{x,1,11}]

Modules

The Module command defines a function with parameters that remain local:

?Module

?While

Clear[f];f[x0_]:=Module[{x=x0},While[x>0,x=Log[x]];
x]

f[3]
f[3.]

Clear[f,np];f[y0_,np0_]:=Module[{np=np0,a=Pi},x=NIntegrate[Exp[a*Sin[y]^np],{y,0,y0}];
x]

f[3,3]

np

Numerical and Symbolic Computation

Integrals with Singularities

NIntegrate[Log[x-1],{x,0,2}]
NIntegrate[(1/(Sqrt[x-1])),{x,0,2}]

Clear[x];Integrate[Log[x-1],{x,0,2}]
Integrate[(1/(Sqrt[x-1])),{x,0,2}]

Helping Mathematica by indicating singularities

NIntegrate[Log[x-1],{x,0,2},Exclusions->1]
NIntegrate[(1/(Sqrt[x-1])),{x,0,2},Exclusions->1]
NIntegrate[(1/(Sqrt[x-1])),{x,0,2},Exclusions->Automatic]

NIntegrate[(1/(Sqrt[x]))(1/(Sqrt[x-1])),{x,-2,2},Exclusions->{0,1}]

Integrals with singular integrands cannot be computed directly
But we can compute the Cauchy Principal Values

Integrate[(1/(x-1)),{x,0,2}]
Integrate[(1/(x-1)),{x,0,2},PrincipalValue->True]
Integrate[(1/(x-1)),{x,0,3}]
Integrate[(1/(x-1)),{x,0,3},PrincipalValue->True]

Solve, NSolve, DSolve & NDSolve

Solve, NSolve, DSolve and NDSolve have transformation rules as outputs

sol=Solve[x^2+a x+1==0,x]

sol[[2]]

xsols=x/.sol
xsols[[1]]
x/.sol[[1]]

sol=Solve[x^3+a x+1==0,x]
sol=Simplify[Solve[x^3+a x+1==0,x],TimeConstraint->10]
sol=FullSimplify[Solve[x^3+a x+1==0,x]]

sol/.a->1.37
sol/.a->1.
sol/.a->-2.
sol/.a->-2
FullSimplify[sol/.a->-2]

Chop again :

x/.sol/.a->-2.
Chop[x/.sol/.a->-2.]

sol/.a->1.
x/.sol/.a->1.

5^th order polynomial: Solve and NSolve

Clear[poly];poly[x_,a_]=x^5-x^4-a*x+1
Plot[poly[x,2],{x,-2,2},Frame->True]

Solve[-x^5-a*x+1==0,x]

NSolve[poly[x,2]==0,x]

FindRoot versus NSolve

NSolve[Exp[x]==x^2+x,x]

FindRoot[Exp[x]==x^2+x,{x,2}]

newpoly[x_]=3*x^5-x^4-4*x^3+1;Plot[newpoly[x],{x,-1.5,2}]

Try with various initial values

FindRoot[newpoly[x]==0,{x,-0.3}]

FindRoot[newpoly[x]==0,{x,-10}]

FindRoot[newpoly[x]==0,{x,10}]

NSolve[newpoly[x]==0,x]

Differential equations:
DSolve, NDSolve

sol=DSolve[{y'[x]+y[x]==a Sin[x],y[0]==2},y[x],x]

Plot[y[x]/.sol/.a->2,{x,0,12},Frame->True]

(* or create a new function *)
Clear[z]
a=2;
z[x_]=(2*y[x]/.sol[[1]])
z[1]
z[1.]

Plot[{z[x],z'[x]},{x,0,12},Frame->True]

Clear[sol];sol=NDSolve[{y'[x]+y[x]== Sin[x],y[3]==2},y,{x,0,12}]

Clear[z,y];z:=y/.sol[[1]]

Plot[z[x],{x,0,13},Frame->True]

System of differential equations

Clear[x,y];s=NDSolve[{x'[t]==-y[t]-x[t]^2,y'[t]==2 x[t]-y[t]^3,x[0]==y[0]==1},{x,y},{t,50}]

ParametricPlot[{x[t],y[t]}/. s,{t,0,50},Frame->True,FrameLabel->{x,y}]

Exercise

Exercise : Solve numerically the growth equation in FRW, pure CDM, from N = log a = -5 to N = 0
delta'' + 1/2*delta' - 3/2*(om*Exp[-3*n]/(om*Exp[-3*n]+1-om))*delta = 0
with om=0.3 and
delta[-5]=1 and delta’[-5]=1
and plot the result

Solution

om=0.3;sol=NDSolve[{delta''[n]+1/2*delta'[n]-3/2*om*Exp[-3*n]/(om*Exp[-3*n]+1-om)==0,delta[-5]==1,delta'[-5]==1},delta,{n,-5,0}]

delta/.sol[[1]]

LogPlot[(delta/.sol[[1]])[n],{n,-5,2},Frame->True]

Reduce

Reduce[expr, vars, dom] -> reduces the statement expr by solving equations or inequalities for vars in domain dom. 
Gives full Boolean solution instead of rules.

Reduce[x^2-7y^2==1,{x,y},Reals]

FullSimplify[%]

Reduce[x^2-7y^2==1&&x>0&&y>0,{x,y},Reals]
Reduce[{x^2-7y^2==1,x>0,y>0},{x,y},Reals]

Reduce[(x^4-1)(x^4-4)==0,x,Complexes]

Reduce[(x^4-1)(x^4-4)==0,x,Reals]

Reduce[(x^4-1)(x^4-4)==0,x,Integers]

Reduce[Sin[x]==1/2,x]

Reduce[{Sin[x]==1/2,0≤x≤2π},x]

Reduce[{Sin[x y]==1/2,0≤x≤π,0≤y≤π},{x,y}]

Reduce[{((Sin[x])^y)==1/2,0≤x≤π,0≤y≤4},{x,y}]

Reduce[{Sin[x]==1/2,0≤x≤π,4≤x≤10},{x,y}]

Limits

Limit[ ]  ->  usual mathematical lim
(by default the limit is taken from the right)

Limit[((Sin[x])/x),x->0]

Limit[(((1+(x/n)))^n),n->Infinity]

Clear[lim,x,a]
lim:=Limit[(x^a),x->Infinity]
FullSimplify[lim,a<0]

Limit[(x^a),x->Infinity,Assumptions->a>0]

Limit[Exp[x](Exp[(1/x)-Exp[-x]]-Exp[(1/x)]),x->∞]

Limit[1/x,x->0]
Limit[1/x,x->0,Direction->"FromAbove"]
Limit[1/x,x->0,Direction->-1]
Limit[1/x,x->0,Direction->"FromBelow"]
Limit[1/x,x->0,Direction->1]

Limit[Sin[x],x->Infinity]

The golden ratio ϕ

Clear[k]
Limit[((Fibonacci[k+1])/(Fibonacci[k])),k->∞]
FullSimplify[%-GoldenRatio]
1.GoldenRatio

(* Golden Spiral r= (ϕ^(2ϑ/π)) *)Graphics[Point[Table[{Exp[theta*Log[GoldenRatio]/(Pi/2)]*Cos[ theta ],Exp[theta*Log[GoldenRatio]/(Pi/2)]*Sin[ theta ]},{theta,.3,30,.1}]]]

Simplifications and Assumptions

We often make implicit assumptions about many variables
Examples:  x ∈ Reals, x ∈ Integers,  x > 0, x ≠ 0, etc.
Mathematica by default make no such assumptions

Important Functions: ‘Simplify’ and ‘FullSimplify’
Other functions: ‘Expand’ ‘PowerExpand’, ‘TrigExpand’ and ‘TrigReduce’

Exact Simplifications: Simplify, Expand,  TrigExpand etc

Integrate[1/x,{x,2,y}]

?NotElement

NotElement[Pi,Integers]

Integrate[1/x,{x,2,y},Assumptions->y∉Reals]

Integrate[1/x,{x,2,y},Assumptions->y>2]

Integrate[1/x,{x,2,y},Assumptions->1>y>0]

Integrate[1/x,{x,2,y},Assumptions->y<0]

Integrate[1/x,{x,2,I-1}]

Simplify[Sqrt[x^2]]

Simplify[Sqrt[x^2],x∈Reals]

Simplify[Sqrt[x^2],x<0]

Simplify[(1/(3 (1+x)))-((-1+2 x)/(6 (1-x+(x^2))))+(2/(3 (1+(1/3) (((-1+2 x))^2))))]

Expand[(((1+x))^5)]

TrigExpand[Sin[2x]]

TrigReduce[2 Cos[x] Sin[x]]

texp=TrigExpand[Expand[(((1+x))^5)Sin[2x]]]

Simplify[((Sin[x])^2)+((Cos[x])^2)]
Simplify[2Tan[x]/(1+((Tan[x])^2))]

FullSimplify[...]  knows about more special functions than Simplify[...], but is slower

Simplify[Cosh[x]-Sinh[x]]
FullSimplify[Cosh[x]-Sinh[x]]

Main commands this lecture :

Interpolation, Map, Module, While, FullSimplify, Expand, TrigExpand,
Solve, NSolve, FindRoot, Reduce, DSolve, NDSolve

Lecture 4

Series and  Transforms

Taylor Series

Series[ ]  ->  Taylor expands a function around any point

s0=Series[(((1-x))^n),{x,0,3}]

Normal[s0]+x^6

Normal[ ]  ->  Changes the header of expressions into more “normal” headers

s0[[0]]

Normal[s0][[0]]

s0+(x^6)

Series[(((1-x))^n),{x,5,2}]

Series[Exp[x],{x,0,5}]

Series[Sin[x],{x,0,4}]

ser=Series[ArcSin[x],{x,0,4}]

InverseSeries[ser]

Series around different points

Use of Evaluate to force evaluation even when normally the evaluation should be delayed:

Plot[{Sin[x],Table[Normal[Series[Sin[x],{x,2,n}]],{n,3}]},{x,0,2Pi},PlotRange->{-5,5},Frame->True]

Plot[{Sin[x],Evaluate[Table[Normal[Series[Sin[x],{x,2,n}]],{n,4}]]},{x,0,2Pi},PlotRange->{-5,5},Frame->True,PlotStyle->{Thick,Dashed,Dashed,Dashed,Dashed}]

Expand around infinity

Clear[n,s]
s=Series[(((1+(1/n)))^n),{n,Infinity,5}]

SeriesCoefficient[series,n] -> gives the coefficient of the n-th power
CoefficientList[ ]  ->  gives the coefficients of all the powers of a SeriesData

SeriesCoefficient[s,4]

s=Series[Sin[x],{x,0,7}]
CoefficientList[s,x]

Legendre Polynomials from a generating function:

s2=Series[(1/(Sqrt[1-2 t x +(t^2)])),{t,0,5}]
CoefficientList[s2,t]

Table[LegendreP[n,x],{n,0,5}]

Multi-variate Series:

s1=Series[Sin[x+y],{x,0,3},{y,0,3}]
s2=Series[Sin[x+y],{y,0,3},{x,0,3}]

FullSimplify[s1-s2]

Simplify[Normal[s1]-Normal[s2]]

Fourier Series

?FourierSeries

?Piecewise

Plot[Piecewise[{{x/2+Pi,-3*Pi<x<-Pi},{x/2,-Pi<x<Pi},{x/2-Pi,Pi<x<3*Pi}}],{x,-3*Pi,3*Pi},PlotRange->All]

?FourierSeries

fs=FourierSeries[(t/2),t,2] (* by default is between -Pi,+Pi*)
Plot[fs,{t,-3Pi,3Pi}]

fs=FourierSeries[t/2,t,7]
Plot[fs,{t,-3Pi,3Pi}]

fs=FourierSeries[t/2,t,20]
Plot[fs,{t,-3Pi,3Pi}]

fff[n_]:=Plot[Evaluate[{FourierSeries[t/2,t,n],Mod[t/2,Pi,-Pi/2]},{t,-3Pi,3Pi}],PlotRange->1.9];
Manipulate[fff[n],{n,1,13,1}]

Discuss the Gibbs phenomenon

From Wikipedia :

squarewave[x_]:=Piecewise[{{-1,x<-π},{1,-π<x<0},{-1,π>x>0},{1,x>π}}]
fourier=FourierSeries[squarewave[x],x,10];
Plot[{fourier,squarewave[x]},{x,-2π,2π}]
fourier=FourierSeries[squarewave[x],x,20];
Plot[{fourier,squarewave[x]},{x,-2π,2π}]

fouriern[n_]:=FourierSeries[squarewave[x],x,n];
Plot[{Evaluate[fouriern[25]],squarewave[x]},{x,-.5π,.5π}]

Find numerically the overshoot :

FindMaximum[Evaluate[fouriern[15]],{x,-.3}]

FindMaxValue[Evaluate[fouriern[15]],{x,-.3}]

Clear[diff]
diff[n_]:=FindMaxValue[Evaluate[fouriern[n]],{x,-.2}]-1;
diff[12]
diff[20]
diff[25]

difftab=Table[{n,diff[n]},{n,12,28,1}];
ListPlot[difftab]

Fourier Transform

FourierTransform[Cos[t],t, ω]

FT of a Gaussian

g[ω_]:=FourierTransform[Exp[-(t^2)/2] ,t,ω]
g[x]
Plot[g[x],{x,-8,8}]

FourierTransform[Cos[Ω t],t,ω]

Inverse FT :

InverseFourierTransform[FourierTransform[Cos[Ω t],t,ω],ω,t]
ExpToTrig[InverseFourierTransform[FourierTransform[Cos[Ω t],t,ω],ω,t]]
Simplify[ExpToTrig[InverseFourierTransform[FourierTransform[Cos[Ω t],t,ω],ω,t]]]

g[ω_]=FourierTransform[1,t,ω]
Integrate[g[x],{x,-1,∞}]

Integrate[DiracDelta[ω],ω]

Power spectrum of a damped sinusoid:
(power spectrum is defined as g[ω]*Conjugate[g[ω]])

?UnitStep

Plot[Sin[(10^3) t]Exp[-1t]UnitStep[t],{t,-3,3},PlotRange->All,Frame->True]

(* Here using := instead of = leads to very slow plotting below  *)
Clear[g]
g[ω_]=FourierTransform[Sin[(10^3) t]Exp[-1t] UnitStep[t],t,ω]

We expect a peak at a frequency of 1000, a constant spectrum for smaller frequencies, and a decay for larger frequencies

LogLogPlot[Evaluate[g[ω]*Conjugate[g[ω]]],{ω,(10^(-1)),(10^5)},PlotRange->All,Frame->True,FrameLabel->{ω,A}]

FourierTransform[(1/((t^2)+1)),t,ω]
InverseFourierTransform[%, ω,t]

Numerical FT of a random signal

?Needs

"download a package"
Needs["FourierSeries`"]

signal=Table[{t,RandomReal[{-1,1}]},{t,-50,50}];

f=Interpolation[signal];
fext[t_]=If[t<-50 || t>50,0,f[t]];
Plot[fext[t],{t,-100,100},Frame->True]

g[ω_]=NFourierTransform[fext[t],t,ω]

g[1.]

Plot[g[ω]*Conjugate[g[ω]],{ω,1,5},PlotPoints->10,PlotRange->All,Frame->True,FrameLabel->{ω,A}]

gex[1.]

gnum[ω_]=NFourierTransform[Exp[-t^2],t,ω];
gex[ω_]=FourierTransform[Exp[-t^2],t,ω];
Plot[{1.01*gnum[ω]*Conjugate[gnum[ω]],Evaluate[gex[ω]*Conjugate[gex[ω]]]},{ω,-5,5},PlotPoints->10,PlotRange->All,Frame->True,FrameLabel->{ω,A}]

Laplace transform

Clear[s];LaplaceTransform[(t^4)Sin[t],t,s]

LaplaceTransform[HeavisideTheta[t],t,s]

LaplaceTransform[DiracDelta[t],t,s]

Importing and Exporting Data

SetDirectory[NotebookDirectory[]]

Image Formats

Image files can be Vectorial or Rasterized
Vectorial images have “infinite resolution” -> they are described by analytical functions (interpolated polynomials + fonts + ...)
	Best vectorial format: “PDF”
Rasterized images are defined pixel-by-pixel. They have finite resolution -> pixel size
	Best rasterized format: “PNG” (you can also consider “JPEG” for photos)	
		
PNG is a lossless format -> great for plots and text, but can lead to large filesizes for photos  
JPEG uses lossy compression -> good for photos with small filesizes, but can lead to artifacts in plots

p1=Plot[{2(x^2)-(x^3),2(x^2)-(x^3)+.3,2(x^2)-(x^3)-.3},{x,0,2},Frame->True,FrameStyle->20,FrameLabel->{"Temperature","ζ[x]"},PlotRange->{-.5,1.6},Filling->{2->{3}},ImageSize->500]

Export["myplot.png",p1]
Export["myplot.jpg",p1]
Export["myplot.pdf",p1]

ImageResolution changes the resolution of exported images
ImageResolution -> r   sets the resolution to r dpi (default = 72 dpi = 72 dots per inch)

Export["myplot-hires.png",p1,ImageResolution->300]

op=0.3;
colorlist={Black,Red,Orange,Yellow,White};
shadelist=Table[{colorlist[[i]],Opacity[op]},{i,1,5}];
cplot=ContourPlot[BesselJ[n,y],{n,0,13},{y,0,20},Contours->{-.2,0,.2,.6},FrameStyle->20,PlotPoints->40,ContourStyle->{Red,White,Green,Yellow},ContourShading->shadelist,GridLines->{{2},{10}},GridLinesStyle->Thick]

Export["myplot2.png",cplot]
Export["myplot2.pdf",cplot]

p2=Plot3D[Sin[x y],{x,0,3},{y,0,3},Mesh->None]

p2=Plot3D[Sin[x y],{x,0,3},{y,0,3},ColorFunction->"Rainbow",Mesh->None]

Export["myplot4.png",p2,ImageResolution->200]
Export["myplot4.jpg",p2,ImageResolution->200]
Export["myplot4.pdf",p2]

Numeric Formats

Numeric formats can be exported in ASCII strings (common text) or in binary formats
ASCII can be opened easily in word processors, spreadsheet programs and in fact in most software
Binary formats can only be opened in specific software, but have much smaller file sizes

mydata2=Table[RandomReal[{0,i}](j^2),{i,1,200},{j,1,5}];
mydata2[[1]]
Export["myfile.txt",mydata2]
Export["myfile.dat",mydata2,"Table"]
Export["myfile.r32",mydata2,"Real32"]

Dimensions[mydata2]

mydata3=Import["myfile.dat","Table"];
Dimensions[mydata3]

mydata3[[1]]

To write a Mathematica expression on a file as it is, use Put or >>; then retrieve it with Get or << (but this only returns the last expression)

?Put

?Get

table={{2,a},{3,c}};table>>"mytable.txt"
table1=<<"mytable.txt"

Import["mytable.txt"]

Put[Sin[x],Sqrt[y],"myexpr.txt"]
<<"myexpr.txt"

Put[{Sin[x],Sqrt[y]},"myexpr.txt"]
<<"myexpr.txt"

Put[{a=3,b=5},"mynewexpr.txt"]
<<"mynewexpr.txt"
a
b

Main commands this lecture :

Series, SeriesCoefficient, CoefficientList,
FourierSeries, FourierTransform,NFourierTransform,LaplaceTransform,
Import, Export, Put,Get

Lecture 5

?ImplicitRegion

gx[x_,y_]:=x+y;
reg=ImplicitRegion[gx[x,y]<1,{{x,-4,4},{y,-4,4}}];

NIntegrate[x^2+y^2,{x,y}∈reg]

The funniest Mathematica command...

Drawing Professional Plots

2D Plots

Types of Plots: 
Plot, LogPlot, LogLogPlot, LogLinearPlot, ContourPlot, ParametricPlot, PolarPlot, ...
ListPlot, ListLogPlot, ListContourPlot, ListParametricPlot, ...
Graphics

Options to be discussed:  PlotStyle, Frame, FrameLabel, FrameStyle, PlotLegends, Filling, FillingStyle

Common Plot Options

Plot[{BesselY[0,x],BesselY[1,x],BesselY[2,x]},{x,0,12}]

Plotstyle adjusts Style of the curves (color, thickness, dashed lines, etc.)

Plot[{BesselY[0,x],BesselY[1,x],BesselY[2,x]},{x,0,12},PlotStyle->{Red,{Blue,Thick},{Orange,Dashed}}]

Frame puts the plot on a frame

Plot[{BesselY[0,x],BesselY[1,x],BesselY[2,x]},{x,0,12},Frame->True]

FrameStyle adjusts Style of the Frame (size, color, etc.)

Plot[{BesselY[0,x],BesselY[1,x],BesselY[2,x]},{x,0,12},Frame->True,FrameStyle->15]

FrameStyle -> {{left, right}, {bottom, top}}    [specify a different style for each border]

Plot[{BesselY[0,x],BesselY[1,x],BesselY[2,x]},{x,0,12},Frame->True,FrameStyle->{{20,Red},{Blue,{Thick,Orange}}}]

We might need  to repeat 4 times the Style...

Plot[{BesselY[0,x],BesselY[1,x],BesselY[2,x]},{x,0,12},Frame->True,FrameStyle->{{{20,Blue,Thick},{20,Blue,Thick}},{{20,Blue,Thick},{20,Blue,Thick}}}]

Use Directive[] to pass on various styles simultaneously in a simpler way

Plot[{BesselY[0,x],BesselY[1,x],BesselY[2,x]},{x,0,12},Frame->True,FrameStyle->Directive[20,Blue,Thick]]

Use PlotLegends to distinguish different curves

Plot[{BesselY[0,x],BesselY[1,x],BesselY[2,x]},{x,0,12},PlotLegends->Automatic]

Plot[{BesselY[0,x],BesselY[1,x],BesselY[2,x]},{x,0,12},PlotLegends->Automatic,Frame->True,FrameStyle->19]

Plot[{BesselY[0,x],BesselY[1,x],BesselY[2,x]},{x,0,12},PlotLegends->Placed[{"This","B","C"},{Right,Bottom}],Frame->True,FrameStyle->19]

FrameLabel -> {bottom, left}
FrameLabel -> {{left, right}, {bottom, top}}

myplot=Plot[{BesselJ[0,x],BesselJ[1,x],BesselJ[2,x]},{x,0,16},Frame->True,FrameStyle->15,FrameLabel->{"Temperature","Time"}]

Dashed, Dotted, DotDashed, Thick, Thin
Dashing & Thickness -> sets specific values

Plot[{BesselJ[0,x],BesselJ[1,x]},{x,0,16},Frame->True,FrameStyle->19,FrameLabel->{"Temperature","Time"},PlotStyle->{{Red,Dashed,Thin},{Purple,Thick}}]

Plot[{BesselJ[0,x],BesselJ[1,x]},{x,0,16},Frame->True,FrameStyle->19,FrameLabel->{"Temperature","Time"},PlotStyle->{{Red,Dashing[.03],Thickness[.01]},{LightBlue,Thickness[.02]}}]

Filling and FillingStyle -> paints region below / above / between curves

Plot[2(x^2)-(x^3)+.3,{x,0,2},Frame->True,FrameStyle->20,FrameLabel->{"Temperature","ζ[x]"},PlotRange->{-.5,1.6},Filling->1.1]

Plot[{2(x^2)-(x^3),2(x^2)-(x^3)+.3,2(x^2)-(x^3)-.3},{x,0,2},Frame->True,FrameStyle->20,FrameLabel->{"Temperature","ζ[x]"},PlotRange->{-.5,1.6},Filling->{2->{3}}]

Plot[{2(x^2)-(x^3),2(x^2)-(x^3)+.3,2(x^2)-(x^3)-.3},{x,0,2},Frame->True,FrameStyle->16,FrameLabel->{"Temperature","ζ[x]"},PlotStyle->{Gray,Black,Red},PlotRange->{-.5,1.6},Filling->{3->{2}}]

myplot=Plot[{2(x^2)-(x^3),2(x^2)-(x^3)+.3,2(x^2)-(x^3)-.3},{x,0,2},Frame->True,FrameStyle->20,FrameLabel->{"Temperature","ζ[x]"},Filling->{2->{3}},PlotRange->{-.5,1.6},FillingStyle->Opacity[.2,Blue],PlotStyle->{Red,{Red,Dashing[.1]},{Red,Dotted}}]

Table[.3i,{i,-10,10}]

myticks={{Table[.3i,{i,-10,10}],None},{{0,1,2},None}}
myplot=Plot[{2(x^2)-(x^3),2(x^2)-(x^3)+.3,2(x^2)-(x^3)-.3},{x,0,2},Frame->True,FrameStyle->16,FrameLabel->{"Temperature","ζ[x]"},Filling->{2->{3}},PlotRange->{-.5,1.6},FillingStyle->Opacity[.1,Blue],PlotStyle->{Red,{Red,Dashing[.1]},{Red,Dotted}},Axes->False,FrameTicks->myticks,GridLines->{{.7,.9873,1.3354,1.49873},{.5}}]

Error Plots

?Around

Around[x,dx]

(* Error propagation with Around*)
 Sin[Around[x,dx]]
PowerExpand[Sin[Around[x,dx]]]

err=Table[{x,Around[x^2,.1*x^2]},{x,1,10,.5}]

ListPlot[err,Frame->True]

err=Table[{x,Around[x^2,{.05*x^2,.2*x^2}]},{x,1,10,.5}]
ListPlot[err,Frame->True]

Exercise:
generate a list of points distributed as a Normal around x with a relative error of 10% for x=1 to 10, in steps of 0.5. 
Plot the points with errorbars given by the standard deviation around each point

relerr=.1;rv=Table[{x,Around[RandomVariate[NormalDistribution[x,.1*x]],.1*x]},{x,1,10,.5}]
ListPlot[rv,Frame->True]

Beyond the simple Plot[] function

Log Plots: LogLogPlot, LogPlot, LogLinearPlot

LogLogPlot[(1/x),{x,.01,10},Frame->True,FrameStyle->18,FrameTicks->{{{.1,1,10,100},Automatic},{{.01,.03,.1,.3,1,3,10},Automatic}}]

LogPlot[(1/x),{x,.01,10},Frame->True,FrameStyle->18]

LogLinearPlot[(1/x),{x,.01,10},Frame->True,FrameStyle->18]

Parametric Plots

ParametricPlot[{Sin[u],Sin[3u]},{u,0,2Pi},Frame->True]

Manipulate[ParametricPlot[{Sin[u],Sin[n u]},{u,0,2Pi}],{n,.5,10},Frame->True]

ParametricPlot[{{2Cos[t],2Sin[t]},{2Cos[t],Sin[t]},{Cos[t],2Sin[t]},{Cos[t],Sin[t]}},{t,0,2Pi},PlotLegends->"Expressions"]

ParametricPlot[{{2Cos[t],2Sin[t]},{2Cos[t],Sin[t]},{Cos[t],2Sin[t]},{Cos[t],Sin[t]}},{t,0,2Pi},PlotLegends->{"circle","ellipsis","ellipsis","circle"}]

Contours -> N
Contours -> List

ContourPlot[Sin[x]Cos[y],{x,0,3π},{y,0,3π},FrameStyle->18,Contours->11]

ContourPlot[Sin[x]Cos[y],{x,0,3π},{y,0,3π},Contours->{-.5,-.2,0,.2,.5,.9,1.1}]

Manipulate[ContourPlot[Sin[x]Cos[y],{x,0,xmax},{y,0,ymax},Contours->{-.5,-.2,0,.2,.5}],{xmax,1,50},{ymax,1,50}]

Plot[{BesselJ[0,y],BesselJ[1,y],BesselJ[7,y]},{y,0,20}]

colorlist={Black,Red,Orange,Yellow,White}

myblue=RGBColor[0.18,0.5,0.38]
Darker[myblue,0]
Darker[myblue,.5]
Lighter[myblue,.3]
Darker[Yellow,.4]

op=.6;
shadelist={{Black,Opacity[op]},{Red,Opacity[op]},{Orange,Opacity[op]},{Yellow,Opacity[op]},{White,Opacity[op]}};
colorlist={Black,myblue,Orange,Yellow,White}
shadelist=Table[{colorlist[[i]],Opacity[op]},{i,1,5}];
cplot=ContourPlot[BesselJ[n,y],{y,0,20},{n,0,13},Contours->{-.2,0,.2,.6},FrameStyle->18,PlotPoints->50,ContourStyle->{Red,White,Green,Yellow},ContourShading->shadelist,GridLines->{{{2,{Red,Thick}},4,6,8,10},{3,6,9,12}},GridLinesStyle->Black]

?ColorFunctionScaling

GraphicsGrid[{{ContourPlot[((BesselJ[n,y]+1)/2),{y,0,10},{n,0,6},Contours->80,FrameStyle->20,PlotPoints->40,ContourStyle->None,ColorFunction->"ThermometerColors",ColorFunctionScaling->False],ContourPlot[((BesselJ[n,y]+1)/2),{y,0,10},{n,0,6},Contours->80,FrameStyle->20,PlotPoints->40,ContourStyle->None,ColorFunction->"ThermometerColors",ColorFunctionScaling->True]}}]

Combining Plots

g[x_]:=(x^(1/3));


mylist=Table[{1.Log[i],1.(i^2)},{i,1,6}];

myplot1=ListPlot[mylist,PlotStyle->{Red},PlotMarkers->{{Graphics[Circle[{0,0},.01]],.04}}]
myplot2=Plot[20g[x],{x,0,2.3},Frame->True,FrameStyle->Red]

Show[myplot2,myplot1,Frame->True,PlotRange->All,FrameStyle->Blue]

myplot3=ContourPlot[Tanh[x]Cos[((y^2)/30)],{x,-1,3},{y,-2,20},ContourShading->None]

Show[myplot1,myplot2,myplot3,PlotRange->{0,19},Frame->True,AspectRatio->1]

Fitting

Linear fits

f (x) = A*x +B*x^2 + C*Exp[x]

Random[]

RandomInteger[{0,100}]

RandomVariate[NormalDistribution[0,1]]

RandomVariate[NormalDistribution[0,1],10]

(* create random data *)

ndata=50;
xtab=Table[i,{i,ndata}];
ytab=3*xtab^2;
relerror=0.1;

data=Table[{xtab[[i]],RandomVariate[NormalDistribution[ytab[[i]],relerror*ytab[[i]]]]},{i,ndata}]
lp=ListPlot[data,Frame->True]

?Fit

Linear in the parameters:
a*1 + b*x
a + b*x+c*x^2
a+b*Sin[x]
Sin[b*x]

Clear[f1,f2];f1[x_]=Fit[data,{1,x},x]

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

f2[x_]=Fit[data,{x^2},x]

pl=Plot[{f1[x],f2[x]},{x,xtab[[1]],xtab[[-1]]}];
Show[lp,pl]

?LinearModelFit

lm1=LinearModelFit[data,x,x]

lm2=LinearModelFit[data,{x,x^2},x]

f3[x_]=Normal[lm2]

lm2["BasisFunctions"]

lm2["ParameterConfidenceIntervals"]  (* default is 95% c.l. *)

lm2["CovarianceMatrix"]//TableForm

lm2["CorrelationMatrix"]//TableForm

lm2["BIC"]
lm1["BIC"]

Fit in several dimensions

?ListPlot3D

(* z=x^2+y^2; *)

ndata=100;
xtab=Table[RandomVariate[UniformDistribution[{-3,3}]],{i,ndata}];
ytab=Table[RandomVariate[UniformDistribution[{-3,3}]],{i,ndata}];
ztab=xtab^2+ytab^2;
relerror=0.1;
data=Table[{xtab[[i]],ytab[[i]],RandomVariate[NormalDistribution[ztab[[i]],relerror*ztab[[i]]]]},{i,ndata}];
ListPlot3D[data]
ListPointPlot3D[data]

z = a + b*x + c*x^2 + d*y + e*y^2 + f*x*y

Fit[data,{x^2,y^2},{x,y}]
Fit[data,{1,x,x^2,y,y^2,x*y},{x,y}]

Exercise:
generate a list of points distributed as a Normal around x with a relative error of 10% for x=1 to 10, in steps of 0.5. 
Plot the points with errorbars given by the standard deviation around each point. Fit with a straight line.

relerr=.1;rvshow=Table[{x,Around[RandomVariate[NormalDistribution[x,.1*x]],.1*x]},{x,1,10,.5}]
rv=Table[{x,RandomVariate[NormalDistribution[x,.1*x]]},{x,1,10,.5}]
fit=Fit[rv,{1,x},x]
Show[ListPlot[rvshow,Frame->True],Plot[fit,{x,1,10}]]

Non linear fits

(* create random data *)

ndata=50;
xtab=Table[i/10,{i,ndata}];
ytab=xtab*Exp[-(xtab-1)^2];

data=Table[{xtab[[i]],RandomVariate[NormalDistribution[ytab[[i]],.1*ytab[[i]]]]},{i,ndata}];

lp=ListPlot[data,Frame->True]

?NonlinearModelFit

y = x*Exp[-(x - 1)^2]

Clear[a,b,c,d];lm=NonlinearModelFit[data,x^a*Exp[c*(x-b)^2],{a,b,c},x]
f[x_]=Normal[lm]

f[x_]=Normal[lm];
pl=Plot[f[x],{x,xtab[[1]],xtab[[-1]]}];
Show[lp,pl]

lm["CovarianceMatrix"]//TableForm

Fit with errorbars

ndata=100;
xtab=Table[i,{i,0,10,10./ndata}];
ytab=10+3*xtab^2;
relerror=0.1;
errors=relerror*ytab;

data=Table[{xtab[[i]],RandomVariate[NormalDistribution[ytab[[i]],relerror*ytab[[i]]]]},{i,Length[xtab]}];
pl=ListPlot[data]
lm=NonlinearModelFit[data,a+b*x^2,{a,b},x,Weights->1/errors^2]
fit=Normal[lm]
pp=Plot[fit,{x,xtab[[1]],xtab[[-1]]}];
Show[pp,pl]

Main commands this lecture :

Plot,LogPlot,LogLogPlot,LogLinearPlot,ContourPlot,ParametricPlot,
ListPlot,ListLogPlot,
(Options: PlotStyle,Frame,FrameLabel,FrameStyle,PlotLegends,Filling,FillingStyle...)
Around,Fit, Normal,LinearModelFit,NonlinearModelFit

Lecture 6

Wolfram Data Base

Wolfram hosts a huge  database at https : // datarepository . wolframcloud . com/
that can be accessed directly from Mathematica (usually you need a data connection). Data about almost everything, from economy to astronomy, from movies to weather.

Natural language: ctrl = opens a “natural language” cell
e.g. ctr=elephant

(Entity["African savanna elephant", Entity["TaxonomicSpecies","LoxodontaAfricana::qtsb3"], "Entity[\"TaxonomicSpecies\", \"LoxodontaAfricana::qtsb3\"]", "species specification"])

(Entity["African savanna elephant", Entity["TaxonomicSpecies","LoxodontaAfricana::qtsb3"], "Entity[\"TaxonomicSpecies\", \"LoxodontaAfricana::qtsb3\"]", "species specification"])

or ctrl = italy

(Entity["Italy", Entity["Country","Italy"], "Entity[\"Country\", \"Italy\"]", "country"])

(Entity["Italy", Entity["Country","Italy"], "Entity[\"Country\", \"Italy\"]", "country"])

(Entity["Germany", Entity["Country","Germany"], "Entity[\"Country\", \"Germany\"]", "country"])

(Entity["Germany", Entity["Country","Germany"], "Entity[\"Country\", \"Germany\"]", "country"])["HighestElevation"]

(Entity["Germany", Entity["Country","Germany"], "Entity[\"Country\", \"Germany\"]", "country"])[{"HighestElevation","Population","PopulationDensity"}]

CountryData[ "Germany","GDP"]
CountryData["Germany","HighestElevation"]
pop=CountryData["Germany",{"Population",{1914,2020}}]

DateListLogPlot[pop,FrameStyle->16]

your location :

mygeo=FindGeoLocation[]

FindGeoLocation["Heidelberg"]

Weather

?WeatherData

WeatherData[mygeo,"Temperature"]
WeatherData[mygeo,"WindSpeed"]
WeatherData[mygeo,"WindDirection"]

min=WeatherData[mygeo, "MinTemperature", {{2000, 1, 1}, {2000, 12, 31}, "Day"}];
max=WeatherData[mygeo, "MaxTemperature", {{2000, 1, 1}, {2000, 12, 31}, "Day"}];
p2000=DateListPlot[{min,max}, Joined->True,Filling->{1->{2}},FillingStyle->LightGray,PlotRange->{-10,40}]

min=WeatherData[mygeo, "MinTemperature", {{2022, 1, 1}, {2022, 12, 31}, "Day"}];
max=WeatherData[mygeo, "MaxTemperature", {{2022, 1, 1}, {2022, 12, 31}, "Day"}];
p2022=DateListPlot[{min,max}, Joined->True,Filling->{1->{2}},FillingStyle->LightGray,PlotRange->{-10,40}]

GraphicsRow[{p2000,p2022}]

Planets and meteorites

MinorPlanetData["Pluto","Radius"]
MinorPlanetData["Eris","Radius"]
MinorPlanetData["Pluto","Color"]

(* download external databases *)

ResourceObject["Meteorite Landings"]

ResourceData["Meteorite Landings"]

CountryData["Italy"]

?GeoGraphics

GeoGraphics[(Entity["Italy", Entity["Country","Italy"], "Entity[\"Country\", \"Italy\"]", "country"]),ImageSize->200]
GeoGraphics[Polygon[(Entity["Italy", Entity["Country","Italy"], "Entity[\"Country\", \"Italy\"]", "country"])],ImageSize->200]
GeoGraphics[Polygon[CountryData["Italy"]],ImageSize->200]

GeoGraphics[Entity["SolarSystemFeature","OlympusMonsMars"],GeoRange->Quantity[400,"Miles"]]

(* use of @ to apply functions *)
Sin@BesselJ[1.,2]
Sin[BesselJ[1.,2]]

GeoGraphics[{Red,PointSize[0.015],Opacity[0.8],Point@DeleteMissing[ResourceData["Meteorite Landings"][All,"Coordinates"]]},GeoRange->CountryData["Germany"],ImageSize->300]

?Gather

Clear[d];Gather[{a,b,a,d,b}]

?EntityClass

?EntityValue

?EntityList

EntityClass["Country","Europe"]
EntityList[EntityClass["Country","Europe"]]

types of entities :
   reference . wolfram . com/language/guide/EntityTypes . html

?DeleteMissing

SNdata=Gather[DeleteMissing[EntityValue[EntityClass["Supernova",All],"DiscoveryYear"]]]

SNdata[[30]]

SNdata[[30,1]]["Year"]

(* number of SN discoveries per year *)

SNlist=Table[{SNdata[[i,1]]["Year"],Length[SNdata[[i]]]},{i,Length[SNdata]}]
Total[SNlist[[All,2]]]

snshortlist=Drop[SNlist,50]

nlm=NonlinearModelFit[snshortlist,b Exp[a*x/1000],{a,b},x,MaxIterations->10000];
fun[x_]=Normal[nlm]

plot1=ListPlot[snshortlist];
plot2=Plot[fun[x],{x,1970,2025},PlotStyle->Red];
Show[plot1,plot2,FrameStyle->20,Frame->True]

BarChart[snshortlist[[All,2]],ChartLabels->Placed[snshortlist[[All,1]],Axis,Rotate[#,90Degree]&],PlotLabel->"number of Supernovae discoveries per year"]

Exercise : use natural language to find the properties of an animal (eg  dog, cat, elephant, etc)
and find their maximum speed and weight

Exercise : use the command WeatherData to find the mean (maximum) summer temperature from 1990 to 2020, and plot the result, including the rms error (square root of the variance) in the plot

?Mean

?StandardDeviation

Solutions

(Entity["Cape mountain zebra", Entity["Species","Infraspecies:EquusZebraZebra"], "Entity[\"Species\", \"Infraspecies:EquusZebraZebra\"]", "species specification"])

s=(Entity["Cape mountain zebra", Entity["Species","Infraspecies:EquusZebraZebra"], "Entity[\"Species\", \"Infraspecies:EquusZebraZebra\"]", "species specification"])["MaximumSpeed"]
w=(Entity["Cape mountain zebra", Entity["Species","Infraspecies:EquusZebraZebra"], "Entity[\"Species\", \"Infraspecies:EquusZebraZebra\"]", "species specification"])["Weight"]

Clear[max];
mygeo=FindGeoLocation[];
max[year_]:=DeleteMissing[WeatherData[mygeo, "MaxTemperature", {{year, 1, 1}, {year, 12, 31}, "Day"}]];
Mean[max[2000]]
StandardDeviation[max[2000]]

(* this produces a "time series" *)
max[2000]

tabmax=Table[{y,Mean[max[y]]},{y,1990,2020,1}]

ListPlot[tabmax,Frame->True]

Clear[max];summermax[year_]:=DeleteMissing[WeatherData[mygeo, "MaxTemperature", {{year, 6,21}, {year, 9, 21}, "Day"}]];
ListPlot[Table[{y,Around[Mean[summermax[y]],Sqrt[Variance[summermax[y]]]]},{y,1990,2020,1}],PlotRange->{10,40},Frame->True]

"TimeSeriesInsert" joins two time series

Clear[max,year];wintermax[year_]:=TimeSeriesInsert[DeleteMissing[WeatherData[mygeo, "MaxTemperature", {{year, 1,1}, {year, 3, 31}, "Day"}]],DeleteMissing[WeatherData[mygeo, "MaxTemperature", {{year, 10,1}, {year, 12, 31}, "Day"}]]];
ListPlot[Table[{y,Around[Mean[wintermax[y]],Sqrt[Variance[wintermax[y]]]]},{y,1990,2020,1}],PlotRange->{-10,20},Frame->True]

Main commands this lecture :

CountryData,WeatherData,Entity, EntityClass,GeoGraphics,TimeSeries,DeleteMissing

Lecture 7

Thanks to Prof . D . Rodrigues for a first version of this lecture

ClearAll["Global`*"];SetDirectory[NotebookDirectory[]]

Tips for faster codes

Use built - in functions; use Interpolation...

Memorization (or memoization)

Functions that remember the computed values

Clear[fib];
fib[1] = 1;
fib[2] = 2;
fib[n_]:=fib[n-1]+fib[n-2];
fib[40] // EchoTiming

Clear[fib];
fib[1] = 1;
fib[2] = 2;
fib[n_]:=fib[n]=fib[n-1]+fib[n-2];
fib[40] // EchoTiming

Use real numbers

* For numerical computations in which  MachinePrecision is sufficient, use float point numbers as soon as possible. That is, write “10.” instead of “10”, or 0.5 instead of 1/2. It is  faster to compute with approximate (float point) numbers than with exact numbers.

N[Det[Table[1/(1+Abs[i-j]),{i,1,120},{j,1,120}]]]//AbsoluteTiming

Det[Table[1./(1.+Abs[i-j]),{i,1.,120.},{j,1.,120.}]]//AbsoluteTiming

Parallelization: ParallelTable, ParallelDo

tab=Table[i,{i,10}]//AbsoluteTiming
(* "First" picks the first element *)
First[tab]

t[n_]:=First[(Table[i, {i,n}]; // AbsoluteTiming)]
d[n_]:=First[(Do[i, {i,n}]; // AbsoluteTiming)]

tab = Table[t[n],{n,1000,100000, 2000}];
dtab = Table[d[n],{n,1000,100000, 2000}];

tab[[10]]

Clear[pt, pd]
pt[n_]:=First[(ParallelTable[i, {i,n}]; // AbsoluteTiming)]
pd[n_]:=First[(ParallelDo[i, {i,n}]; // AbsoluteTiming)]

ptab = Table[pt[n],{n,1000,100000, 2000}];
pdtab = Table[pd[n],{n,1000,100000, 2000}];

(* doesn't seem so good...*)

ListPlot[{ dtab,  pdtab,ptab}, PlotLegends->{ "Do",  "ParallelDo","ParallelTable"},Frame->True]

?$KernelCount

$KernelCount

list=10^Table[i,{i,3,8,.25}]

list=10^Table[i,{i,3,8,.25}]
dtabExt = Table[{list[[i]],d[list[[i]]]},{i,Length[list]}];
pdtabExt = Table[{list[[i]],pd[list[[i]]]},{i,Length[list]}];

ListLogLogPlot[{ dtabExt, pdtabExt}, PlotLegends->{ "Do",  "ParallelDo"},Frame->True]

External variables/functions are automatically distributed to all subkernels; variables/function defined internally are only valid for the subkernels where have been defined

ClearAll["Global`*"];a=Table[0,{n,10}];

ParallelDo[{sol=NDSolve[{d''[x]==-n*d[x],d[0]==1,d'[0]==0},d,{x,-100000.,100000.}];f=d/.sol[[1]];Print[f[1]];a[[$KernelID]]=f[1];b=f[1]},{n,10}]

a
b

Share variables

ClearAll["Global`*"];a=Table[0,{n,10}];
SetSharedVariable[a,b];
ParallelDo[{sol=NDSolve[{d''[x]==-n*d[x],d[0]==1,d'[0]==0},d,{x,-100000.,100000.}];f=d/.sol[[1]];a[[$KernelID]]={$KernelID,n,f[1]};
Print[{a[[$KernelID]]}];b=f[1]},{n,10}]

a
b

Tips for ParallelDo : write results that you need to use outside the Do in a file, and read it afterwards

Main commands this lecture :

ParallelTable,ParallelDo

Lecture 8

Linking Mathematica to code in other languages

Compile[] -> compiles simple functions on the fly in order to make your code faster.

Compiling simple functions

?Compile

Essentially, "Compile" rewrites your code internally as a single optimized code, rather than running command after command. The problem is that it is hard to predict when it is really more efficient.

"WVM" is Wolfram Virtual Machine

Clear[fx];
fraw[x_]=Sin[x] + x^2 - 1/(1+x);
fx=Compile[{{x, _Real}}, Sin[x] + x^2 - 1/(1+x),CompilationTarget->"WVM"]

?CompilationTarget

xi=RandomReal[{0,1},1000000];

Table[fraw[xi[[i]]],{i,1000000}];//AbsoluteTiming
Table[fx[xi[[i]]],{i,1000000}];//AbsoluteTiming

Listable functions

ClearAll["Global`*"];

list={a,b};
(* non listable *)
fun[list]
(* listable *)
Sin[list]
Sqrt[list]
list^2
Attributes[Sin]

SetAttributes[fun,Listable];
fun[list]

?Listable

RuntimeAttributes

?RuntimeAttributes

?Parallelization

fx=Compile[{{x, _Real}}, Sin[x] + x^2 - 1/(1+x),CompilationTarget->"WVM",RuntimeAttributes->{Listable},Parallelization->True]

fxc=Compile[{{x, _Real}}, Sin[x] + x^2 - 1/(1+x),CompilationTarget->"C",RuntimeAttributes->{Listable},Parallelization->True]

xi=RandomReal[{0,1},10000000];

Total[fx[xi]]//AbsoluteTiming

Total[fxc[xi]]//AbsoluteTiming

CForm & FortranForm Example

CForm[PDF[NormalDistribution[mu,sigma],x]]

FortranForm[PDF[NormalDistribution[mu,sigma],x]]

Probability Distributions

The Wolfram Language is very versatile when it comes to statistical distributions
A distribution is a special kind of function. It is nonnegative-definite and its integral over the whole parameter space is unity.

Many functions operate directly on distributions, like PDF, CDF, Probability, Mean, Variance etc.

{μ=2,σ=1};
dist1=NormalDistribution[μ,σ]
PDF[dist1]
GraphicsRow[{Plot[PDF[dist1][y],{y,-2,6},Frame->True,FrameStyle->15],
Plot[CDF[dist1][y],{y,-2,6},Frame->True,FrameStyle->15]},ImageSize->600]

{μ=2,σ=1};
dist2=LogNormalDistribution[μ,σ]
{μ=2,σ=2};
dist3=LogNormalDistribution[μ,σ]
Clear[μ,σ];PDF[LogNormalDistribution[μ,σ]];
GraphicsRow[{Plot[{PDF[dist2][y],PDF[dist3][y]},{y,-2,26},Frame->True,FrameStyle->15,PlotRange->All],
Plot[{CDF[dist2][y],CDF[dist3][y]},{y,-2,26},Frame->True,FrameStyle->15]},ImageSize->600]

{Mean[dist2],Variance[dist2]}
{Mean[NormalDistribution[μ,σ]],Variance[NormalDistribution[μ,σ]]}
{Mean[LogNormalDistribution[μ,σ]],Variance[LogNormalDistribution[μ,σ]]}
{Mean[PoissonDistribution[μ]],Variance[PoissonDistribution[μ]]}
StandardDeviation[PoissonDistribution[μ]]

E.g.: The Maxwell-Boltzmann velocity distribution

?MaxwellDistribution

dist0=MaxwellDistribution[sqrtT];
PDF[dist0]
(* find mean and mode *)
{Mean[dist0],Simplify[Solve[{PDF[dist0]'[x]==0,x>0},x],sqrtT>0]}

Manipulate[{Plot[PDF[dist0][x]/.sqrtT->(Sqrt[T]),{x,0,10},PlotRange->{0,.6},Frame->True,FrameStyle->15]},{T,1,10}]

Two important operators for random variables:
x \[Distributed] dist  ->  the random variable “x” is distributed according to the distribution dist   (hotkey dist or esc dist esc)
expr \[Conditioned] cond  ->  the expression expr is conditioned to the hypothesis cond   (hotkey cond or esc cond esc)

\[Distributed]
\[Conditioned]

\[Distributed]
\[Conditioned]

\[Distributed]
\[Conditioned]

\[Distributed]

αβ

Probability[x<5,x\[Distributed]dist1]

Probability[x<5.,x\[Distributed]dist1]
Probability[x Tanh[x]<1.,x\[Distributed]dist1]

Probability[(x^2)<16.,x\[Distributed]dist1]
Probability[(x^2)<16.\[Conditioned]x>3,x\[Distributed]dist1]
Probability[(x^2)<16.&& x>3,x\[Distributed]dist1]

Show[Plot[PDF[dist1][x],{x,-6,6},Frame->True],RegionPlot[x>3 && x<4,{x,3,6},{y,0,2},BoundaryStyle->None]]

(* conditional probability *)

Show[Plot[If[x>3,PDF[dist1][x],0],{x,-6,6},Frame->True,PlotRange->All],RegionPlot[x>3 && x<4,{x,3,6},{y,0,.4},BoundaryStyle->None]]

Transformed Distributions

We often need to carry out change of variables in probability
TransformedDistribution is an built-in function for that
If the resulting distribution is not tabulated, further computation will be slow

TransformedDistribution[3x+1,x\[Distributed]NormalDistribution[μ,σ]]
TransformedDistribution[Exp[x],x\[Distributed]NormalDistribution[μ,σ]]

dist4=TransformedDistribution[(x^2),x\[Distributed]NormalDistribution[0,1]]
Plot[PDF[dist4][x],{x,-.1,.5},PlotRange->{0,9},Frame->True,FrameStyle->15]

Manipulate[Plot[PDF[ChiSquareDistribution[n]][x],{x,0,30},PlotRange->{{0,30},{0,1}}],{n,1,30}]

dist5=TransformedDistribution[(x^2)+(y^2)+(z^2),{x\[Distributed]NormalDistribution[0,1],y\[Distributed]NormalDistribution[0,1],z\[Distributed]NormalDistribution[0,1]}]
Plot[PDF[dist5][x],{x,-1,5},PlotRange->All,Frame->True,FrameStyle->15]

dist7=TransformedDistribution[(x^a),x\[Distributed]UniformDistribution[]]
Manipulate[Plot[PDF[dist7][x]/.a->aa,{x,-.1,1.1},PlotRange->All,Frame->True,FrameStyle->15],{aa,.2,2}]

TransformedDistribution[x/y,{x\[Distributed]NormalDistribution[0,σ1],y\[Distributed]NormalDistribution[0,σ2]}]

TransformedDistribution[x/Sqrt[y],{x\[Distributed]NormalDistribution[0,σ1],y\[Distributed]ChiSquareDistribution[1]}]

Exercise 1. Use Manipulate to visualize how the PDF and the CDF of the chi - squared distribution change with the number of degrees of freedom

Creating your own distribution from data

?Histogram

data=RandomVariate[NormalDistribution[],10^4];Histogram[data,{0.25}]

Export["data"<>ToString[.25]<>".dat",data]

ToExpression["0.25"]+1

newdist=HistogramDistribution[data];
{Mean[newdist],StandardDeviation[newdist]}
Table[Plot[f[newdist,x],{x,-4,4},PlotLabel->f,Frame->True,FrameStyle->14],{f,{PDF,CDF}}]

(* smoothed version *)
newdist=SmoothKernelDistribution[data];
{Mean[newdist],StandardDeviation[newdist]}
Table[Plot[f[newdist,x],{x,-4,4},PlotLabel->f,Frame->True,FrameStyle->14],{f,{PDF,CDF}}]

Lecture 9

SetDirectory[NotebookDirectory[]]

Resource Functions

The Wolfram Function Repository is a curated cloud repository of functions set up to be instantly usable in the Wolfram Language . The Function Repository includes a growing number of functions designed for a wide range of application areas . The Wolfram Language supports creation of private resource functions, which can then be submitted for inclusion in the Wolfram Function Repository .

https://resources.wolframcloud.com/FunctionRepository

PlotGrid

ResourceFunction["PlotGrid"]

[{{Plot[x,{x,0,1},Frame->True],Plot[x^2,{x,0,1},Frame->True]}}]

ResourceFunction["PlotGrid"][{{Plot[x,{x,0,1},Frame->True],Plot[x^2,{x,0,1},Frame->True]}}]

[{{Plot[x,{x,0,1},Frame->True]},{Plot[x^2,{x,0,1},Frame->True]}},ImageSize->200]

Metropolis Hastings

ResourceFunction["MetropolisHastingsSequence"]

(* generate a MCMC sampling of a likelihood *)

mc=[Exp[-1./2.(x^2.+y^2.)],{{x,0},{y,0}},{50000,0,1}];

mc[[3]]

Histogram3D[mc]

SmoothHistogram3D[mc]

(* convert into a prob distribution *)

dist=SmoothKernelDistribution[mc]

Plot3D[PDF[dist,{x,y}],{x,-3,3},{y,-3,3},PlotRange->All,PlotPoints->50]

ContourPlot[PDF[dist,{x,y}],{x,-3,3},{y,-3,3},PlotRange->All]

Covariance[dist]//MatrixForm

(* generate a MCMC sampling of a likelihood *)

cov={{1,.5},{.5,1}};invcov=Inverse[cov];var={x,y};mc=[Exp[-1./2.*var.invcov.var],{{x,0},{y,0}},{50000,0,1}];
(* convert into a prob distribution *)
dist=SmoothKernelDistribution[mc];
ContourPlot[PDF[dist,{x,y}],{x,-3,3},{y,-3,3},PlotRange->All]
Covariance[dist]//MatrixForm

max=NMaximize[PDF[dist,{x,y}],{x,y}]
maxl=max[[1]]

(* integrate over a region high than level *)

level=0.15;r=ImplicitRegion[PDF[dist,{x,y}]>level && -2<x<2 && -2<y<2,{x,y}];
NIntegrate[PDF[dist,{x,y}],{x,y}∈r]

(* do that for several levels *)

res={};tablev=Table[level,{level,0,.95*maxl,maxl/5}];
Do[{r=ImplicitRegion[PDF[dist,{x,y}]>tablev[[i]] && -2<x<2 && -2<y<2,{x,y}];
AppendTo[res,{NIntegrate[PDF[dist,{x,y}],{x,y}∈r],tablev[[i]]}]},{i,Length[tablev]}]

res

hei=Interpolation[res]

ContourPlot[{PDF[dist,{x,y}]==hei[.68],PDF[dist,{x,y}]==hei[.9]},{x,-3,3},{y,-3,3},PlotRange->All]

RandomMaze

ResourceFunction["RandomMaze"][40]

Magic Squares

ms=ResourceFunction["MagicSquare"][10];ms//TableForm

ms[[All,1]]
ms[[1,All]]

Table[Total[ms[[All,i]]],{i,Length[ms]}]
Table[Total[ms[[i,All]]],{i,Length[ms]}]
Total[Diagonal[ms]]

QR code

BarcodeImage["https://www.thphys.uni-heidelberg.de/~amendola/teaching.html","QR",200]

BarcodeImage["You can store any random text in a QR code, up to 4000 characters","QR",200]

Parsing text in internet

text=Import["https://edition.cnn.com/"];

?StringCount

StringCount[text,"Biden"]
StringCount[text,"Trump"]
StringCount[text,"Musk"]

An association is a list in which key a is associated with value x, key b with value y, etc . :

(* associations *)

assoc=<|a->x,b->y,c->z|>

<|"a"->x,"b"->{y,z}|>[["b"]]

(* convert lists into associations*)

list2Assoc[in_]:=Which[VectorQ[in,Head[#]===Rule&],Apply[Association,Map[list2Assoc,in]],ListQ[in],Map[list2Assoc,in],Head[in]===Rule&&ListQ[in[[2]]],in[[1]]->list2Assoc[in[[2]]],True,in];

list2Assoc[{"name"->"ab","surname"->"cd"}]

(* get a response from the inspirehep API (Application Programming Interface); the response is a list in JSON format *)

authorsn="Witten";authorn="E.";
authorsn="Einstein";authorn="A.";
url="https://inspirehep.net/api/literature?sort=mostrecent&size=500&page=1&q=find%20au%20"<>authorsn<>"%2C%20"<>authorn;

text=Import[url,"JSON"];
list=list2Assoc[text];
data=list[["hits","hits"]];
cit=Table[data[[i]][["metadata","citation_count"]],{i,Length[data]}];
Print["name: ",authorsn<>", "<>authorn," , papers: ",Length[cit],"  , citations: ",Total[cit]]

text
list

URLEncode["<"]
URLEncode[">"]
URLEncode[" "]
URLEncode[","]

(* get a response from the inpirehep API (Application Programming Interface); the response is a list in JSON format *)

authorsn="Einstein";authorn="A.";
authorsn="Witten";authorn="E.";
datein="1950";dateout="2025";
url="https://inspirehep.net/api/literature?sort=mostrecent&size=500&page=1&q=find%20au%20"<>authorsn<>"%2C%20"<>authorn<>"%20and%20date%3E"<>datein<>"%20and%20date%3C"<>dateout;
text=Import[url,"JSON"];
list=list2Assoc[text];
data=list[["hits","hits"]];
cit=Table[data[[i]][["metadata","citation_count"]],{i,Length[data]}];
Print["name: ",authorsn<>", "<>authorn," , papers: ",Length[cit],"  , citations: ",Total[cit]]

(* get a response from the inpirehep API (Application Programming Interface); the response is a list in JSON format *)

authorsn="Einstein";authorn="A.";
datein="1950";dateout="2025";
url="https://inspirehep.net/api/literature?sort=mostrecent&size=500&page=1&q=find%20au%20"<>authorsn<>"%2C%20"<>authorn<>"%20and%20date%3E"<>datein<>"%20and%20date%3C"<>dateout;
text=Import[url,"JSON"];
list=list2Assoc[text];
data=list[["hits","hits"]];
cit=Table[data[[i]][["metadata","citation_count"]],{i,Length[data]}];
Print["name: ",authorsn<>", "<>authorn," , papers: ",Length[cit],"  , citations: ",Total[cit]]

Machine Learning

?FindClusters

cov={{1,0},{0,1}};np=20;c1={0,0};c2={5,5};

points=Join[RandomVariate[MultinormalDistribution[c1,cov],np],RandomVariate[MultinormalDistribution[c2,cov],np]];
ListPlot[points,Frame->True]

cl=FindClusters[points];
ListPlot[cl,Frame->True]

cov={{1,0},{0,1}};np=20;c1={0,0};c2={5,5};c3={-5,-5};c4={2,-2};
set[c_]:=RandomVariate[MultinormalDistribution[c,cov],np];
points=Join[set[c1],set[c2],set[c3],set[c4]];
ListPlot[points,Frame->True,PlotRange->All,PlotLabel->"Original"]

(* free number of clusters *)
cl=FindClusters[points];
ListPlot[cl,Frame->True,PlotLabel->"Free number of clusters"]

(* require exactly 4 clusters *)
cl=FindClusters[points,4];
ListPlot[cl,Frame->True,PlotLabel->"Four clusters"]

ImageIdentify[]

figure=Import["galaxy.jpg"];
figure

ImageIdentify[figure]

Strange symbols: @ and &

@ : short for applying a function to an argument

Sin[Pi/4.]
Sin@(Pi/4.)

Plot[x^2,{x,0,3},PlotStyle->Thickness@.03]

Pure (or anonymous) functions

f=#^2&

f@4
#^2&@4

Clear[f];list={a,b,c};
list2={d,e,f};

Map[#^2&,list]

#^2&/@list

Apply[Plus,list]
Apply[Times,list]

Plus@@list
Times@@list
Times@@@{list,list2}

Music

We discuss three simple musical commands : Play, SoundNote, and Sound

?Play

?SoundNote

(* SoundNote creates a note, and Sound reproduces it *)

Sound[{SoundNote["C",1,"Violin"],SoundNote["G",0.5,"Piano",SoundVolume->0.8]}]

You can choose many instruments (where did they find them??) :

p=0.25;Sound[{SoundNote["C",{0,1},"Violin"],SoundNote["G",{0,.5},"Piano"],SoundNote["C",{1,1.25},"Piano"],SoundNote["D",{1.25,1.75},"Piano"],SoundNote["A",{1.25,2},"Violin"],SoundNote["D",{1.75,2},"Piano"],
SoundNote["C",2.+p+{0,1},"Violin"],SoundNote["G",2.+p+{0,.5},"Piano"],SoundNote["C",2.+p+{1,1.25},"Piano"],SoundNote["D",2.+p+{1.25,1.75},"Piano"],SoundNote["B",2.+p+{1.25,2},"Violin"],SoundNote["D",2.+p+{1.75,2},"Piano"],SoundNote["C",2+p+{2.25,3},"Violin"]}]

p=0.25;instrument="Percussion";
instrument="Harp";
instrument="Percussion";Sound[{SoundNote["C",{0,1},instrument],SoundNote["G",{0,.5},instrument],SoundNote["C",{1,1.25},instrument],SoundNote["D",{1.25,1.75},instrument],SoundNote["A",{1.25,2},instrument],SoundNote["D",{1.75,2},instrument],
SoundNote["C",2.+p+{0,1},instrument],SoundNote["G",2.+p+{0,.5},instrument],SoundNote["C",2.+p+{1,1.25},instrument],SoundNote["D",2.+p+{1.25,1.75},instrument],SoundNote["B",2.+p+{1.25,2},instrument],SoundNote["D",2.+p+{1.75,2},instrument],SoundNote["C",2+p+{2.25,3},instrument]}]

?Play

Play[Sin[1500 t]+Sin[1520 t],{t,0,2}]

There is also a function that converts mathematical expressions into music

ResourceFunction["MathToMusicTransform"][x^2,{x,-7,7},.5]

ResourceFunction["MathToMusicTransform"][Exp[x/2]+Exp[-x/3],{x,-7,7},.3]

Composing Random music

?RandomSample

(* create a random melody *)

instruments={"Harp","Piano","Xylophone"};
Clear[dt0];
dt={0.05,0.1,0.2};
dt={0.1,0.2,0.4};

score={};t0new=0;
nn=30;
notes={"C4","D4","E4","F4","G4","A4","B4"};

Do[{dt0=RandomSample[dt,1][[1]];t0=t0new;AppendTo[score,SoundNote[RandomSample[notes,1][[1]],{t0,t0+dt0},RandomSample[instruments,1][[1]]]],t0new=t0+dt0},{i,nn}];
tt=t0
Sound[score]

(* add a random base with lower volume *)

instruments={"BlownBottle"};
Clear[dt0];
dt={1};
basescore={};t0new=0;
nn=Floor[tt];
notes={"C3","E3"};
Do[{dt0=RandomSample[dt,1][[1]];t0=t0new;AppendTo[basescore,SoundNote[RandomSample[notes,1][[1]],{t0,t0+dt0},RandomSample[instruments,1][[1]],SoundVolume->.8]],t0new=t0+dt0},{i,nn}];
basescore
Sound[basescore]

(* add instead a not-so-random base  *)

instruments={"BlownBottle"};
Clear[dt0];
dt={.5};
basescore={};t0new=0;
nn=Floor[tt]/dt[[1]];
notes={"C3","E3","G3"};
j=0;
Do[{dt0=dt[[1]];t0=t0new;
note=notes[[1+Mod[i,3]]];
AppendTo[basescore,SoundNote[note,{t0,t0+dt0},RandomSample[instruments,1][[1]],SoundVolume->1.]],t0new=t0+dt0},{i,nn}];
basescore
Sound[basescore]

(* combine *)

Sound[Join[score,basescore]]

Sound[{SoundNote["C",{0,1},instrument],SoundNote["G",{0,.5},instrument],SoundNote["C",{1,1.25},instrument],SoundNote["D",{1.25,1.75},instrument],SoundNote["A",{1.25,2},instrument],SoundNote["D",{1.75,2},instrument],
SoundNote["C",2.+p+{0,1},instrument],SoundNote["G",2.+p+{0,.5},instrument],SoundNote["C",2.+p+{1,1.25},instrument],SoundNote["D",2.+p+{1.25,1.75},instrument],SoundNote["B",2.+p+{1.25,2},instrument],SoundNote["D",2.+p+{1.75,2},instrument],SoundNote["C",2+p+{2.25,3},instrument]}]