mth3060402m9.html

Forced Harmonic Oscillators, I

November 11, 1998

Example: Sinusoidal Forcing

We consider the following underdamped, forced oscillator.

> eqn1:=diff(y(t),t,t)+6*diff(y(t),t)+13*y(t)=sin(3*t);

dsolve will find the solution, but MAPLE seems to need a little help cleaning up the mess.

> sol1:=dsolve(eqn1,y(t));

> sol11:=combine(sol1,trig);

Subexample I

If we put in some values for the constants of integration, we get something we can plot.

> y1:=subs({_C1=1,_C2=1},subs(sol11,y(t)));

> y1prime:=diff(y1,t);

Here are the graphs of the solution (red) and its derivative (blue).

> plot([y1,y1prime],t=0..10,color=[red,blue]);

We get a phase plane picture by plotting the solution and its derivative parametrically.

> plot([y1,y1prime,t=0..10]);

Subexample 2

Here's the same thing with different values.

> y2:=subs({_C1=0,_C2=-1},subs(sol11,y(t)));

> y2prime:=diff(y2,t);

> plot([y2,y2prime],t=0..10,color=[red,blue]);

> plot([y2,y2prime,t=0..10]);

Using phaseportrait

We don't get the same sort of output because we no longer have an autonomous system. Phaseportrait just plots the solutions (numerically).

> with(DEtools):

> phaseportrait(diff(y(t),t$2)+6*diff(y(t),t)+13*y(t)=sin(3*t),
y(t),t=0..10,
[[y(0)=-1,D(y)(0)=3],[y(0)=1,D(y)(0)=-1]],
stepsize=0.1,
linecolor=[red,blue]);

To get a phase plane picture, it seems that we need to write the equation as a system. Notice that since the system is not autonomous, we don't get a direction field.

> phaseportrait([diff(y(t),t)=v(t),
D(v)(t)=-6*v(t)-13*y(t)+sin(3*t)],
[y(t),v(t)],t=0..10,
[[y(0)=-1,v(0)=3],[y(0)=1,v(0)=-1]],
stepsize=0.05,
linecolor=[red,blue]);

Example: Undamped Forcing

Here we consider an undamped oscillator forced at a frequency .

> eqn2:=diff(y(t),t,t)+8*y(t)=cos(omega*t);

dsolve will find the general solution, but Maple seems to need a little help in cleaning up the answer.

> sol2:=dsolve(eqn2,y(t));

> sol21 := combine(sol2,trig);

Let's look at the particular solution for which and .

> sol2part := dsolve({eqn2,y(0)=0,D(y)(0)=0},y(t));

> sol22part := combine(sol2part,trig);

> y3 := subs(sol22part,y(t));

Note that the particular solution is the some of a free response and a steady-state response . The addition of two such sinusoidal functions with different frequencies results in "beats".

Here are the beats when . The solution is graphed in red, and the beats in blue.

> plot([subs(omega=2,y3),
(1/2)*sin((1-sqrt(2))*t),
-(1/2)*sin((1-sqrt(2))*t)],
t=0..30,
color=[red,blue,blue]);

When we also have beats, but they have a much longer period.

> plot([subs(omega=3,y3),
2*sin((3-2*sqrt(2))*t/2),
-2*sin((3-2*sqrt(2))*t/2)],
t=0..80,
color=[red,blue,blue]);