Harmonic oscillator

(((d^2)x)/(dt^2))+c (dx/dt)+(k^2)x=0

sol=DSolve[x''[t]+c x'[t]+(k^2)x[t]==0,x,t]

Mathematica returns a Function.  C[1] and C[2] are the two constants.  Access the function and plug in a value for time

x[0]/.sol
x'[0]/.sol
x[t]/.sol

Solve for the constants using initial conditions:

x[0] = (Subscript[x, 0]),  x'[0] = (Subscript[v, 0])

csol=Solve[{(x[0]/.sol)==x0,(x'[0]/.sol)==v0},{C[1],C[2]}]

Plug in these values to the solution

(sol/.csol)

Then to access the function...

x[t]/.(sol/.csol)

Define a function to make this easier to access

result[t_]:=First[Flatten[x[t]/.(sol/.csol)]];

Check out the initial conditions we applied - we’ll need to use Simplify

result[0]

Simplify[result[0]]

Simplify[result'[0]]

The final form of the solution

Simplify[result[t]]

This isn’t a good form to view it in!  FullSimplify uses Exponential to Trig (ExpToTrig) and vice versa (TrigToExp) conversions

FullSimplify[result[t]]

We can look at the solution for different cases too:

Undamped:  c=0

FullSimplify[result[t]/.c->0]

Critically damped:  (c^2)=4(k^2)

FullSimplify[result[t]/.c->(Sqrt[4(k^2)])]

this gives an error... we can use Limit to evaluate the expression

Limit[result[t],c->(Sqrt[4(k^2)])]

Now lets define a function and use it to change parameters and plot results.  We can call it the same thing, and Mathematica will known which version to call based on the number of input parameters (this is called pattern matching)

result[t_,damping_,stiffness_,initialposition_,initialvelocity_]:=result[t]/.{c->damping,k->stiffness,x0->initialposition,v0->initialvelocity};

Plot[result[t,1,1,1,1],{t,0,10}]

Now define a Manipulate function to change the parameter values and view the solution result

Manipulate[
Plot[result[t,damping,stiffness,initialposition,initialvelocity],{t,0,endtime}],
{damping,0,10},{stiffness,1,10},{initialposition,0,10},{initialvelocity,0,10},{endtime,10,100}]

Note that we could have also included the initial conditions in the DSolve command, and obtained the same thing.  We’ll do this with the next example

sol2=DSolve[{x''[t]+c x'[t]+(k^2)x[t]==0,x[0]==x0,x'[0]==v0},x,t];
FullSimplify[x[t]/.sol2]

FullSimplify[result[t]]

Forced harmonic oscillator

(((d^2)x)/(dt^2))+c (dx/dt)+(k^2)x=F[t]
Different forms for F[t]:  sine wave, (sine^2), square wave, sawtooth wave, triangle wave

We’ll need to numerically solve this:  NDSolve.  First, a simple case with Sin wave forcing function, setting numerical parameter values all to 1.  We’ll include the initial conditions in the NDSolve command.

nsol=NDSolve[{x''[t]+ x'[t]+x[t]==Sin[t],x[0]==1,x'[0]==1},x,{t,0,10}]

Mathematica returns an InterpolatingFunction.  This is a function that applies over a certain range and can be used just like a normal function (e.g. you can take derivatives and integrals and plot it).  However, it doesn’ have an explicit functional form like the solution we had before.

Plot[x[t]/.nsol,{t,0,10}]

Let’s set up a Manipulate function so we can change parameter values and even forcing functions
First, define the forcing functions we want to use.  
To make the sine wave have the same period...

sinwave[t_]:=Sin[2Pi t];
sinsquared[t_]:=((Sin[2Pi t])^2);

Then SquareWave is fine.  Also others:

SquareWave[t]
TriangleWave[t]
SawtoothWave[t]

Create a manipulate to do all this.

Manipulate[
Module[{nsol},
nsol=NDSolve[{x''[t]+damping x[t]+stiffness x[t]==forcing[freq t],x[0]==initialposition,x'[0]==initialvelocity},x,{t,0,finaltime}];
Plot[{x[t]/.nsol,forcing[freq t]},{t,0,finaltime},PlotStyle->{{Blue,Thick},{Red,Thick,Dashed}}]]
,{damping,0,10},{stiffness,1,10},{initialposition,0,10},{initialvelocity,0,10},{forcing,{sinwave,sinsquared,SquareWave,TriangleWave,SawtoothWave}},{freq,.1,1},{finaltime,10,100}]

Schrodinger equation:  Finite potential well
(we didn’t get to this example during the workshop)

(non-dimensionalized version)
-(((d^2)ψ[x])/(dx^2))+(Subscript[V, 0])ψ[x]==E ψ[x], E>0

Assume that ({{E<, (Subscript[V, 0])}}), so the particle is trapped by the potential well
Boundary conditions:

(1)   (Subscript[ψ, 1])[-∞]==0
(2)   (Subscript[ψ, 1])[-a]==(Subscript[ψ, 2])[-a]
(3)   (Subscript[ψ, 1])'[-a]==(Subscript[ψ, 2])'[-a]
(4)   (Subscript[ψ, 2])[a]==(Subscript[ψ, 3])[a]
(5)   (Subscript[ψ, 2])'[a]==(Subscript[ψ, 3])'[a]
(6)   (Subscript[ψ, 3])[∞]==0

In region 1 we want

(Subscript[Lim, x->∞])(Subscript[ψ, 1])[-x]=0

First, get the general solution, then apply this boundary condition

sol1=DSolve[-psi''[x]+v0 psi[x]==e psi[x],psi,x]
Limit[psi[x]/.sol1,x->-Infinity,Assumptions->e>0&&e<v0]

For this to be true, we want C[2]=0 for region1.  Also, rename the other constant to avoid confusion

psi1[x_]:=First[(psi[x]/.sol1)/.{C[2]->0,C[1]->c3}];
psi1[x]

region 2: get solution, where (Subscript[V, 0])=0: rename the constants to prevent confusion

sol2=DSolve[{-psi''[x]==e psi[x]},psi,x]
psi2[x_]:=First[(psi[x]/.sol2)/.{C[1]->c1,C[2]->c2}];
psi2[x]

In region 3 we want

(Subscript[Lim, x->∞])(Subscript[ψ, 3])[x]=0

First, get the general solution, then apply this boundary condition

sol3=DSolve[-psi''[x]+v0 psi[x]==e psi[x],psi,x]
Limit[psi[x]/.sol3,x->Infinity,Assumptions->e>0&&e<v0]

For this to be true, we want C[1]=0 for region3.  Also, rename the constant to avoid confusion

psi3[x_]:=First[(psi[x]/.sol3)/.{C[1]->0,C[2]->c4}];
psi3[x]

Apply the first 3 boundary conditions

boundaryconditions={psi1[-a]==psi2[-a],psi1'[-a]==psi2'[-a],psi2[a]==psi3[a],psi2'[a]==psi3'[a]}

bcsol=Solve[Take[boundaryconditions,3],{c1,c2,c3}]

Apply the last boundary condition and try to solve for c4

sol4=Solve[Last[boundaryconditions]/.bcsol,c4]

This gives a trivial solution!  Need to solve for the energy level, E in order to enforce this boundary conditions

Last[boundaryconditions]/.bcsol

Simplify this assuming that c4!=0, then define a variable equal to the left hand side

eqn=FullSimplify[Last[boundaryconditions]/.bcsol,Assumptions->c4!= 0]

lhs=eqn[[1,1]]

To Solve this, we’ll need to use FindRoot.  First, let’s plot to see where the roots are for a test case

test=-2 (Sqrt[e]) (Sqrt[-e+v0]) Cos[2 a (Sqrt[e])]+(2 e-v0) Sin[2 a (Sqrt[e])];
Plot[test/.{a->1,v0->10},{e,0,10}]

FindRoot[test/.{a->1,v0->10},{e,9.9}]

Plot[lhs/.{a->1,v0->10},{e,0,10}]

Note that there is NOT a root where the asymptotes are (e.g. at ~ 1.9 and ~7).  Put in estimated values to find roots:

FindRoot[lhs/.{a->1,v0->10},{e,1}]

FindRoot[lhs/.{a->1,v0->10},{e,5}]

FindRoot[lhs/.{a->1,v0->10},{e,9.9}]

Define a Manipulate module to change the starting point, for different parameter values

Manipulate[
Module[{sol},
sol=FindRoot[lhs/.{a->size,v0->potential},{e,startingpoint}];
Show[
Plot[lhs/.{a->size,v0->potential},{e,0,potential}],
ListPlot[{{e/.sol,0}},PlotStyle->PointSize[Large]],
PlotLabel->sol]]
,{size,1,10},{potential,1,100},{startingpoint,.1,potential}]