PHYS 2063 Wave Physics:
Solving the Equation of Motion

Differential equations and Newton’s Second Law

Remove["Global`*"]

Solving Differential Equations

The syntax is very similar to solving algebraic equations, but the result is a function

eq=y'[x]==a x
sol=DSolve[eq,y[x],x][[1]]

ysol=y[x]/.sol

You can also write the equation using the “Differentiate” function

eq=D[y[x],x]==a y[x]
sol=DSolve[eq,y[x],x][[1]]

ysol=y[x]/.sol

eq=y''[x]==a^2 y[x]
sol=DSolve[eq,y[x],x][[1]]

ysol=y[x]/.sol

Adding Boundary (or Initial) Conditions

The syntax is almost obvious. Boundary conditions are just additional equations

eq=y'[x]==a x;
bc=y[0]==b;
sol=DSolve[{eq,bc},y[x],x][[1]];

ysol=y[x]/.sol

eq=y'[x]==a y[x];
bc=y[0]==b;
sol=DSolve[{eq,bc},y[x],x][[1]];

ysol=y[x]/.sol

eqs={y''[x]==a^2 y[x],y[0]==A,y'[0]==B};
sol=DSolve[eqs,y[x],x][[1]];

ysol=y[x]/.sol

ysol/.x->0
D[ysol,x]/.x->0

Simplify[%,a>0]

Newton’s Second Law as a Differential Equation

Just write “F=ma” with F=F(x) and a=D[D[x[t],t],t]

The simplest possible example is F(x)=f=ma, x[0]=x0, and v[0]=v0

f=m a;
xsol=x[t]/.DSolve[{f==m D[D[x[t],t],t],x[0]==x0,x'[0]==v0},x[t],t][[1]]

This is just the equation for position as a function of time for constant acceleration

Example: A Mass Falling with Drag

Remove["Global`*"]

Let y be the vertical position, and the mass m falls with drag -bv from a height h, starting from rest

eq=m y''[t]==-m g-b y'[t];
ic1=y[0]==h;
ic2=y'[0]==0;

y=y[t]/.DSolve[{eq,ic1,ic2},y[t],t][[1]]

Does this result make sense for small and large values of b?

Limit[y,b->0]

Limit[y,b->∞,Assumptions->{m>0&&t>0&&g>0}]

Now plot the result, allowing us to manipulate the value of b

vals={m->0.1,h->10,g->9.8};

yplot=y/.vals;
tMaxP=tMax/.vals;
Manipulate[
Plot[yplot/.b->bSlide,{t,0,5},PlotRange->{0,10}],
{bSlide,0.001,0.5}]

Example: The Simple Harmonic Oscillator

Remove["Global`*"]

Let it start out at x(0)=A (the “amplitude”) and with zero velocity

f=-k x[t];
sol=DSolve[{m D[D[x[t],t],t]==f,x[0]==A,x'[0]==0},x[t],t];
x1=x[t]/.sol[[1]]

Write it in terms of the angular frequency ω

x2=x1/.k->ω^2m
x=Simplify[x2,m>0&&ω>0]

Plot the motion for a few cycles

xplot=x/.{A->5,t->2 Pi τ/ω};
Plot[xplot,{τ,0,4}]

Example: The Damped Harmonic Oscillator

Remove["Global`*"]

Follow the same procedure, just add a damping term to the force. Convert to frequency parameters.

f=-k x[t]-b x'[t];
eq=m D[D[x[t],t],t]==f/.{k->ω^2m,b->2r m}

Set the “natural” (undamped) period to unity

ω=2 Pi;

Do the underdamped, overdamped, and critically damped cases separately, with example values for β

sol=DSolve[{eq/.r->ω/10,x[0]==A,x'[0]==0},x[t],t];
xU=x[t]/.sol[[1]]

sol=DSolve[{eq/.r->3ω,x[0]==A,x'[0]==0},x[t],t];
xO=x[t]/.sol[[1]]

sol=DSolve[{eq/.r->ω,x[0]==A,x'[0]==0},x[t],t];
xC=x[t]/.sol[[1]]

Make a plot of the solutions

A=5;
Plot[{xU,xO,xC},{t,0,5},
PlotRange->{-A,A},PlotLegends->{Under,Over,Critical}]

Energy for the Damped Oscillator

Remove["Global`*"]

Set up and solve the under damped oscillator, and keep the undamped solution

f=-k x[t]-b x'[t];subs={k->ω^2m,b->2r m};
eq=m D[D[x[t],t],t]==f/.subs;
x=x[t]/.DSolve[{eq/.r->α ω,x[0]==A,x'[0]==0},x[t],t][[1]];

ωA={A->5,ω->2Pi};α0=α->0;α1=α->0.002;α2=α->0.01;
x0=x/.ωA/.α0;
x1=x/.ωA/.α1;
x2=x/.ωA/.α2;
tMax=25;
Plot[{x0,x1,x2},{t,0,tMax},PlotLegends->{α/.α0,α/.α1,α/.α2}]

Now plot the energy as a function of time

e0=(1/2)m D[x0,t]^2+(1/2)k x0^2/.subs/.ωA/.m->1;
e1=(1/2)m D[x1,t]^2+(1/2)k x1^2/.subs/.ωA/.m->1;
e2=(1/2)m D[x2,t]^2+(1/2)k x2^2/.subs/.ωA/.m->1;
Plot[{e0,e1,e2},{t,0,tMax},PlotLegends->{α/.α0,α/.α1,α/.α2}]