Solving Equations in Mathematica

Solve
NSolve
DSolve
NDSolve

Solving simple algebraic equations

Solve[(x^4)+(x^3)+2(x^2)+5x+12==0,x];
Solve[(x^4)+(x^3)+2(x^2)+5x+12.==0,x];
Solve[(x^5)+(x^3)+2(x^2)+5x+12==0,x];
NSolve[(x^5)+(x^3)+2(x^2)+5x+12==0,x];

Solving differential equations

(* For example: a simple harmonic oscillator *)

DSolve[y''[t]+y'[t]/10+y[t]==0,y[t],t]
DSolve[y''[t]+0.1y'[t]+y[t]==0,y[t],t]

solution1=DSolve[{y''[t]+0.1y'[t]+y[t]==0,y[0]==1,y'[0]==0},y[t],t]

Plot[y[t]/.solution1,{t,0,50},PlotRange->All,Frame->True,Axes->True]

solution2=DSolve[{y''[t]+0.1y'[t]+y[t]==a Cos[0.3 t],y[0]==1,y'[0]==0},y[t],t]

Plot[y[t]/.solution2/.a->0.01,{t,0,50},PlotRange->All,Frame->True]

DensityPlot[y[t]/.solution2,{t,0,50},{a,0,1},PlotPoints->40]

(* Introduce some extra complication *)
solution3=DSolve[{y''[t]+0.1Cos[0.5 t]y'[t]+y[t]==0.1 Cos[0.3 t],y[0]==1,y'[0]==0},y[t],t]

(* DSolve could not manage this, but we can solve it numerically using NDSolve *)
solution3=NDSolve[{y''[t]+0.1Cos[0.5 t]y'[t]+y[t]==0.1Cos[0.3 t],y[0]==1,y'[0]==0},y[t],{t,0,50}]

Plot[y[t]/.solution3,{t,0,50},PlotRange->All,Frame->True]

How to Solve Schrödinger's Equation Numerically

We can discretize the potential, and the wavefunction in space, over a suitable range.

ψ[x] becomes an array {ψ[(Subscript[x, 1])],ψ[(Subscript[x, 2])],...}
Likewise for V[x]

The step size is ((Subscript[x, i+1])-(Subscript[x, i]))=dx

Now we can make a first order approximation of the second derivative of the wavefunction:

ψ''[x] becomes ψ''[(Subscript[x, i])] = (1/(dx^2))(ψ[(Subscript[x, i+1])]-2ψ[(Subscript[x, i])]+ψ[(Subscript[x, i-1])])

Now we can write the Schrödinger equation as a matrix eigenvalue/eigenvector problem:

H ψ = E ψ

writing the hamiltonian operator -((ℏ^2)/(2m))((Subscript[∇, i])^2)+V[(Subscript[x, i])] now as a matrix.

This matrix now looks like this:

Clear[ℏ,m,dx,V];
size=5;
ExampleMatrix=(DiagonalMatrix[Table[-((ℏ^2)/(2m))((-2)/(dx^2))+V[x],{size}]]+DiagonalMatrix[Table[-((ℏ^2)/(2m))(1/(dx^2)),{size-1}],1]+DiagonalMatrix[Table[-((ℏ^2)/(2m))(1/(dx^2)),{size-1}],-1])//MatrixForm

Example 1 - finite square well potential

ℏ=1;
m=1;
a=1;
V=100;
ndivs=300;
boundary=2a;
dx=2boundary/ndivs;

(* Use 'piecewise' definition of functions *)

esc pw esc brings up the bracket to define a piecewise function
ctrl+comma moves to the right
ctrl+enter starts a new row

newfunction=(Piecewise[{{Sin[10x], x<0}, {x, x>=0}}])

Plot[newfunction,{x,-5,5},Frame->True,Axes->False,PlotRange->{Automatic,{0,3}}]

Vfunc1[x_]=(Piecewise[{{0, x<-(a/2)||(a/2)<x}, {-V, -a/2≤x≤a/2}}])

Varray1=Table[Vfunc1[x],{x,-boundary,boundary,dx}];

ListPlot[Varray1,Filling->1.1Min[Varray1],Frame->True,PlotRange->All]

MProblem1=N[DiagonalMatrix[-((ℏ^2)/(2 m))(((-2))/(dx^2))+Varray1]+DiagonalMatrix[Table[-((ℏ^2)/(2 m))(1/(dx^2)),{ndivs}],1]+DiagonalMatrix[Table[-((ℏ^2)/(2 m))(1/(dx^2)),{ndivs}],-1]];

answers1=Eigensystem[MProblem1];

{evals,evects}=Sort[Transpose[answers1],#1[[1]]<#2[[1]]&]//Transpose;

nboundstates=Select[evals,#<0&]//Length

boundstateenergies=Take[evals,nboundstates]

ListPlot[ Table[evects[[i]]+(i-1)/2,{i,1,nboundstates}],Joined->True,Filling->Table[n->(n-1)/2,{n,1,nboundstates}],PlotRange->All,ImageSize->600,Frame->True]

potplot=ListPlot[Varray1,Filling->1.1Min[Varray1],Joined->True,Frame->True,PlotRange->All,GridLines->{None,boundstateenergies},GridLinesStyle->Directive[Thickness[0.002],Red,Dashed]]

ScaleFactor=30;
stateplot=ListPlot[Table[ScaleFactor evects[[i]]+evals[[i]],{i,1,nboundstates}],Joined->True,Filling->Table[n->evals[[n]],{n,1,nboundstates}],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->All,ImageSize->600,Frame->True,Axes->False]

Show[stateplot,potplot]

Example 2: A double square well potential

Vfunc2[x_]=(Piecewise[{{0, x<-(a/2)||(a/2)<x}, {-V/2, -a/2≤x≤0}, {-V, 0<x≤a/2}}])

Varray2=Table[Vfunc2[x],{x,-boundary,boundary,dx}];

ListPlot[Varray2,Filling->1.1Min[Varray2],Frame->True,PlotRange->All]

MProblem2=N[DiagonalMatrix[-((ℏ^2)/(2 m))(((-2))/(dx^2))+Varray2]+DiagonalMatrix[Table[-((ℏ^2)/(2 m))(1/(dx^2)),{ndivs}],1]+DiagonalMatrix[Table[-((ℏ^2)/(2 m))(1/(dx^2)),{ndivs}],-1]];

answers2=Eigensystem[MProblem2];

{evals,evects}=Sort[Transpose[answers2],#1[[1]]<#2[[1]]&]//Transpose;

nboundstates=Select[evals,#<0&]//Length

boundstateenergies=Take[evals,nboundstates]

ListPlot[ Table[evects[[i]]+(i-1)/2,{i,1,nboundstates}],Joined->True,Filling->Table[n->(n-1)/2,{n,1,nboundstates}],PlotRange->All,ImageSize->600,Frame->True]

potplot=ListPlot[Varray2,Filling->1.1Min[Varray2],Joined->True,Frame->True,PlotRange->All,GridLines->{None,boundstateenergies},GridLinesStyle->Directive[Thickness[0.002],Red,Dashed]]

ScaleFactor=30;
stateplot=ListPlot[Table[ScaleFactor evects[[i]]+evals[[i]],{i,1,nboundstates}],Joined->True,Filling->Table[n->evals[[n]],{n,1,nboundstates}],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->All,ImageSize->600,Frame->True,Axes->False]

Show[stateplot,potplot]

Example 3: A harmonic oscillator

ℏ=1;
m=1;
a=1;
V=100;
ndivs=500;
boundary=5a;
dx=2boundary/ndivs;

Vfunc3[x_]=-V+30 (x^2)

Varray3=Table[Vfunc3[x],{x,-boundary,boundary,dx}];

ListPlot[Varray3,Filling->1.1Min[Varray3],Frame->True,PlotRange->All]

MProblem3=N[DiagonalMatrix[-((ℏ^2)/(2 m))(((-2))/(dx^2))+Varray3]+DiagonalMatrix[Table[-((ℏ^2)/(2 m))(1/(dx^2)),{ndivs}],1]+DiagonalMatrix[Table[-((ℏ^2)/(2 m))(1/(dx^2)),{ndivs}],-1]];

answers3=Eigensystem[MProblem3];

{evals,evects}=Sort[Transpose[answers3],#1[[1]]<#2[[1]]&]//Transpose;

(* From now on we will just arbitrarily take the lowest 10 levels *)

stateenergies=Take[evals,10]

ListPlot[ Table[evects[[i]]+(i-1)/2,{i,1,10}],Joined->True,Filling->Table[n->(n-1)/2,{n,1,10}],PlotRange->All,ImageSize->600,Frame->True]

potplot=ListPlot[Varray3,Filling->1.1Min[Varray3],Joined->True,Frame->True,PlotRange->All,GridLines->{None,stateenergies},GridLinesStyle->Directive[Thickness[0.002],Red],Axes->False]

ScaleFactor=20;
stateplot=ListPlot[Table[ScaleFactor evects[[i]]+evals[[i]],{i,1,10}],Joined->True,Filling->Table[n->evals[[n]],{n,1,10}],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->All,ImageSize->600,Frame->True,Axes->False]

Show[stateplot,potplot,PlotRange->{Automatic,{-V,0.1V}}]

(* Very nearly perfectly equally spaced energy levels *)
ListPlot[tab=Table[stateenergies[[i+1]]-stateenergies[[i]],{i,1,10-1}],PlotRange->{Automatic,{7.5,8}},PlotStyle->PointSize[0.02],GridLines->{None,{Mean[tab]}},Frame->True,Axes->False]

(* Slight deviation due to numerical errors *)