4.3 and 4.4: Forcing and Resonance

November 13, 1998

The undamped osillator forced at resonant frequency

The solutions here are linear combinations of the real and imaginary parts of the general solution to

diff(y, `$`(t, 2))+omega[0]*y = exp(i*omega*t) ,

that is,

t*exp(i*omega[0]*t)/2*i*omega[0]+k[1]*exp(i*omega[0]*t)+k[2]*exp(-i*omega[0]*t) .

For example, with initial conditions y = 0 and dy/dt = 0 at t = 0 , we get

> ivalues := solve({subs(t=0,t*exp(I*omega[0]*t)/(2*I*omega[0])+k[1]*exp(I*omega[0]*t)+k[2]*exp(-I*omega[0]*t))=0,subs(t=0,diff(t*exp(I*omega[0]*t)/(2*I*omega[0])+k[1]*exp(I*omega[0]*t)+k[2]*exp(-I*omega[0]*t),t))=0},{k[1],k[2]});

ivalues := {k[2] = -1/4/omega[0]^2, k[1] = 1/4/omega[0]^2}

and the solution is

> ycomplex:=subs(ivalues,t*exp(I*omega[0]*t)/(2*I*omega[0])+k[1]*exp(I*omega[0]*t)+k[2]*exp(-I*omega[0]*t));

ycomplex := -1/2*I*t*exp(I*omega[0]*t)/omega[0]+1/4*exp(I*omega[0]*t)/omega[0]^2-1/4*exp(-I*omega[0]*t)/omega[0]^2

We get two real solutions by taking real and imaginary parts.  The real part corresponds to cos(omega[0]*t) forcing, the imaginary part to sin(omega[0]*t) forcing.

> yreal := evalc(Re(ycomplex));

yreal := 1/2*t*sin(omega[0]*t)/omega[0]

> yimaginary := evalc(Im(ycomplex));

yimaginary := -1/2*t*cos(omega[0]*t)/omega[0]+1/2*sin(omega[0]*t)/omega[0]^2

Here is a picture of the two functions.

> plot([subs(omega[0]=1,yreal),subs(omega[0]=1,yimaginary)],
t=0..10*Pi,color=[red,blue]);

[Plot]

Amplitude and Phase Angle

Amplitude

> A := 1/sqrt((q-omega^2)^2+p^2*omega^2);

A := 1/sqrt(q^2-2*q*omega^2+omega^4+p^2*omega^2)

Here's a picture of the graph of A when omega = 1 .

> plot3d(subs(omega=1,A),q=0..2,p=0..1,axes=boxed,view=[0..2,0..1,0..10]);

[Plot]

Here's a picture of how A varies when p = 1 .

> plot3d(subs(p=1,A),q=0..2,omega=0..2,axes=boxed,view=[0..2,0..2,0..2]);

[Plot]

And when q = 1 .

> plot3d(subs(q=1,A),p=0..2,omega=0..2,axes=boxed,view=[0..2,0..2,0..5]);

[Plot]

Phase Angle

> phi := arctan(-p*omega,(q-omega^2));

phi := arctan(-p*omega, q-omega^2)

Fixed omega = 1 :

> plot3d(subs(omega=1,phi),p=0..1,q=0..2,axes=boxed,view=-Pi..0,
orientation=[45,-45]);

[Plot]

Fixed q = 1 :

> plot3d(subs(q=1,phi),p=0..1,omega=0..2,axes=boxed,view=-Pi..0);

[Plot]

Fixed p = 1 :

> plot3d(subs(p=1,phi),omega=0..2,q=0..4,axes=boxed,view=-Pi..0);

[Plot]