PHYS 325: Discussion 2, Problem 1

By - Sagnik Saha, 9/15/2025.
In this notebook, we will see how to use Mathematica for some basic tasks, in the context of the first problem in this week’s discussion. A couple of points to note:
​
1. Comments related to the code itself (i.e. not the physics) are provided in blue.
2. Comments related to the physics of the problem are given in black.
​
In general, this notebook will outline the syntax of the Mathematica code as clearly as possible, since it is designed to be a pedagogical introduction.

Part (0) - Skills: Defining/calling functions, rule/pattern matching.

Defining and calling functions.

(*Definethepotentialfunction*)
u[x_]:=-u0*(a^2+x^2)/(8*a^4+x^4);
(*Callthefunctionforgeneralx*)
u[x](*NotethatnowIdonotusean_or:.*)
Out[]=
-
u0(
2
a
+
2
x
)
8
4
a
+
4
x
(*Callthefunctionforaspecificvalueofx*)
u[2]
Out[]=
-
(4+
2
a
)u0
16+8
4
a

Pattern/Rule matching

(*Let'squicklylookatsomethinginMathematicacalledpattern/rulematching.*)
(*Example1*)
x/.x->5
Out[]=
5
(*Example2*)
abcde/.abcde->10
Out[]=
10
(*Seewhat'shappeninghere?Wearerequiringabunchof"pattern/rule replacements".Mathematicalooksattheexpressiongivenontheleftofthe/.,andthenusesour"rules"tomaketheproperreplacements.Thisisveryhandyaswewillseebelow.*)

Part (a) - Skills: Taking derivatives, simplifying expressions.

Taking Derivatives

The force f(x) is given as -dU/dx.
(*Let'stakethederivativetheWRONGWAYfirst.*)
In[]:=
f[x_]:=-D[u[x],{x,1}];(*WRONG*)
(*Callthefunctionwithsomegeneralx*)
f[x]
Out[]=
-
4u0
3
x
(
2
a
+
2
x
)
2
(8
4
a
+
4
x
)
+
2u0x
8
4
a
+
4
x
(*Hmm.Seemsfine-what'stheproblem?Sayyouwantf(2).Let'strythat.*)
f[2]
D
:2 is not a valid variable.
Out[]=
-
∂
{2,1}
-
(4+
2
a
)u0
16+8
4
a
(*What'shappeninghereisthatMathematicaistryingtotakeaderivativew.r.ta"variable"2.Itdoesn'tunderstandthatonemusttakethederivativeofU(x)w.r.txandTHENsetx=2.Wecanenforcethisasfollows.*)
(*Clearolddefinitionoff*)
In[]:=
Clear[f]
(*CORRECTWAY*)
In[]:=
f[x_]:=-D[u[y],{y,1}]/.y->x;(*The{y,1}isnotnecessary-canjustwritey.Higherderivativesspecifiedas{y,2/3/4...}.*)
(*Callthefunctionforgeneral'x'.*)
f[x]
Out[]=
-
4u0
3
x
(
2
a
+
2
x
)
2
(8
4
a
+
4
x
)
+
2u0x
8
4
a
+
4
x
(*Nowtryf(2)*)
f[2]
Out[]=
-
32(4+
2
a
)u0
2
(16+8
4
a
)
+
4u0
16+8
4
a
(*Works!*)

Simplifying

(*SometimeswealsoneedtoforceMathematicatocleanupbetter.UseFullSimplifyasshown.*)
(*IfonlyIhada$foreverytimeIusedFullSimplify...*)
FullSimplify[f[x]]
Out[]=
-
2u0(-8
4
a
x+2
2
a
3
x
+
5
x
)
2
(8
4
a
+
4
x
)
FullSimplify[f[2]]
Out[]=
(-2+
2
a
)(1+
2
a
)u0
2
2
(2+
4
a
)
(*Note:ThereareotherwaystosimplifyexpressionssuchasExpand,Together,TrigReduce,etc.ButFullSimplifyisthemostcommon.*)

Part (b) - Solving equations.

Solving equations

For equilibrium points, we must set U’(x) = 0 (equivalently, set f(x) = 0).
(*Findingextremumpoints.Setf(x)=0*)
expts=Solve[f[x]==0,x]
Out[]=
{x0},{x-2a},{x2a},{x-
2
a},{x
2
a}
(*Hmm.YouseethatMathematicadoesn'treallyknowthatx,aaresupposedtoberealnumbers.Let'srepeatthisstepwithsomeassumptions.*)
expts=Solve[f[x]==0,x,Assumptions->{x∈Reals,a∈Reals,a>0}]
Out[]=
{x0},{x-
2
a},{x
2
a}
(*Thislooksmuchbetter.Note,however,thattheresultisinaweirdformat.Forconvenience,let'spullouttheindividualvaluesandstorethemina``refined''listusingpatternmatching.Weuseeveryelementofouroldlist"expts"asa"rule".*)
exptsrefined={x/.expts[[1]],x/.expts[[2]],x/.expts[[3]]}
Out[]=
0,-
2
a,
2
a
(*NOTE-Inthiscase,asymbolicsolutionwaspossible.For(global)numericalsolutions,NSolveisthefunctionyouwant.YoucouldalsouseFindRoot,whichisusefulwhenyouonlycare(orknow)aboutrootsinalocalregion.*)
Ok moving on. Now, we must do the second derivative test for each of these points.
​
Define a new function for U’’(x).
In[]:=
upp[x_]:=FullSimplify[D[u[y],{y,2}]]/.y->x
(*Plugintheextremumpointsdefinedabove.*)
upp[exptsrefined[[1]]]​​upp[exptsrefined[[2]]]​​upp[exptsrefined[[3]]]
Out[]=
-
u0
4
4
a
Out[]=
u0
3
4
a
Out[]=
u0
3
4
a
Thus, assuming u0>0, x = 0 is a maxima, and x = \sqrt{2}a and -\sqrt{2}a are minima.

Part (c) - Plotting Functions.

(*Tomaketheplots,letusre-scaleouru(x)andf(x)byu0,sinceitdoesn'tchangethequalitativebehavioroftheplot.*)
In[]:=
urescaled[x_]:=u[x]/.u0->1​​frescaled[x_]:=f[x]/.u0->1
As we showed earlier, x=0 is always a maximum of the potential, while the other two are always minima.

Final Comments

In this notebook, we saw a few functionalities of Mathematica. Of course, just like any other language, there are tons of other features that I could not demonstrate here. If you decide to continue using Mathematica (highly recommended), then the best resource to learn is the Mathematica Documentation. You can look this up online:
https://reference.wolfram.com/language/
Whenever you require a certain functionality, you can look up the documentation to see if Mathematica can help you out (usually it can).
​