Part_06.html

Theory of Oscillations
Part 6: Velocity and Acceleration of Oscillations and Phase Curves

Alexei V. Tikhonenko

Email: avtikhon@okclub.org

General Physics Department, Institute of Nuclear Power Engineering, Russia

This worksheet demonstrates the use of Maple for analyzing velocity and acceleration oscillations and illustrates properties of oscillations represented by phase curves.

Introduction

Velocity and acceleration of oscillations

We will receive analytical expressions for velocity and acceleration of oscillations:

, ,

and visualize plots of x(t), v(t) and w(t).

Phase curves of oscillations

Phase curves of oscillations are plots of velocity v(t) depending on oscillation variable x(t).

DEtools Package and plots Package

We will give demonstration of analytical and numerical calculations of velocity and acceleration and phase curves of oscillations with " DEtools " and " plots " Packages.

> restart; with(DEtools): with(plots):

Warning, the name changecoords has been redefined

Section 6.1. FREE DAMPED OSCILLATIONS

Consider free damped oscillations with two different values of damping factors.

MAPLE executable expression for

- proper frequency w0 ,

- two values of damping factor g0 , G0 ,

- initial conditions x0 , v0 ,

- differential equation for free damped oscillation,

- resulting analytical solutions for two values of damping factor,

- plots of solutions for three values of damping factor: a) x_01a(t) , b) x_01b(t),

are:

> omega0:=2; gamma0:=0.05; Gamma0:=0.15; x0:=4; v0:=0;

Eqn_01a:=diff(x(t),t,t)+2*gamma0*diff(x(t),t)+omega0^2*x(t)=0;
Eqn_01b:=diff(x(t),t,t)+2*Gamma0*diff(x(t),t)+omega0^2*x(t)=0;

sol1a:=dsolve({Eqn_01a,x(0)=x0,D(x)(0)=v0},x(t));
sol1b:=dsolve({Eqn_01b,x(0)=x0,D(x)(0)=v0},x(t));

plot(rhs(sol1a),t=0..60);
plot(rhs(sol1b),t=0..60);

Now we will calculate velocity and acceleration of oscillations:

> x_01a(t):= rhs(sol1a);
v_01a(t):=diff(x_01a(t),t);
w_01a(t):=diff(v_01a(t),t);

x_01b(t):=rhs(sol1b);
v_01b(t):=diff(x_01b(t),t);
w_01b(t):=diff(v_01b(t),t);

plot([x_01a(t),v_01a(t),w_01a(t)],t=0..30,color=[red,blue,gold],linestyle=[1,4,3],thickness=[4,2,2]);
plot([x_01b(t),v_01b(t),w_01b(t)],t=0..30,color=[red,blue,gold],linestyle=[1,4,3],thickness=[4,2,2]);

Next plots are phase curve of oscillations ( gold color ), variable of oscillations x(t) ( red color ) and velocity v(t) ( blue color ):

> G_01a:=array(1..2):
G_01a[1]:=plot([v_01a(t),x_01a(t),t=0..60],color=magenta, axes=box):
G_01a[2]:=plot([t,x_01a(t),t=0..60],color=red, axes=box):
display(G_01a, linestyle=1, thickness=2, axes=box);
G_011a:=array(1..2):
G_011a[1]:=plot([v_01a(t),-t,t=0..60],color=blue, thickness=2, axes=box):
G_011a[2]:=plot([t,-t,t=0..60],color=green, thickness=4, axes=box):
display(G_011a, linestyle=1, axes=box);

The last plot does not have any information for free oscillations).

> G_01b:=array(1..2):
G_01b[1]:=plot([v_01b(t),x_01b(t),t=0..60],color=magenta, axes=box):
G_01b[2]:=plot([t,x_01b(t),t=0..60],color=red, axes=box):
display(G_01b, linestyle=1, thickness=2, axes=box);
G_011b:=array(1..2):
G_011b[1]:=plot([v_01b(t),-t,t=0..60],color=blue, thickness=2, axes=box):
G_011b[2]:=plot([t, -t,t=0..60],color=green, thickness=4, axes=box):
display(G_011b, linestyle=1, axes=box);

Next plots are a comparison of the esults for a) and b):

> G2d_011ab := array(1..2):
G2d_011ab[1]:=plot([[v_01a(t),x_01a(t),t=0..30],[v_01b(t),x_01b(t),t=0..30]],color=[magenta,pink], axes=box, linestyle=[1,3]):
G2d_011ab[2]:=plot([[t,x_01a(t),t=0..30], [t,x_01b(t),t=0..30]],color=[red,pink], axes=box,linestyle=[1,3]):
display(G2d_011ab, linestyle=1, thickness=2, axes=box);
G2d_012ab := array(1..2):
G2d_012ab[1]:=plot([[v_01a(t),-t,t=0..30], [v_01b(t),-t,t=0..30]],color=[blue,pink], thickness=2, axes=box,linestyle=[1,3]):
G2d_012ab[2]:=plot([t, -t,t=0..30],color=green, thickness=4, axes=box):
display(G2d_012ab, linestyle=1, axes=box);

At last we represent phase curves of oscillations with time:

> G3d_01:= array(1..2):
G3d_01[1]:=spacecurve([v_01a(t),x_01a(t),t],t=0..30,numpoints=1000,color=red):
G3d_01[2]:=spacecurve([v_01b(t),x_01b(t),t],t=0..30,numpoints=1000,color=blue):
display(G3d_01,orientation=[0,80],thickness=2,axes=none);

Geometric properties of phase curves (curvature and torsion) correspond to damping characteristics of oscillations. The blue line corresponds to oscillations with greater damping factor.

Section 6.2. FORCED DAMPED OSCILLATIONS 1 (HARMONIC EXTERNAL FORCE)

Consider forced damped oscillations with two external forces.

MAPLE executable expressions for

- frequency of external force W ,

- amplitude of external force f0 ,

- proper frequency w0 ,

- damping factor g0 ,

- initial conditions: x0 , v0 ,

- external forces Force_02a and Force_02b ,

- differential equation of oscillation with external forces Force_02a and Force_02b ,

- resulting analytical solutions for two external forces: a) x_03a(t) , b) x_03b(t) ,

- plots of solutions for or two external forces: a) x_03a(t) , b) x_03b(t) ,

- plots of forces Force_02a and Force_02b

are:

> Omega:=3; f_02:=10;
Force_02a:=f_02*cos(Omega*t);
Force_02b:=f_02*sin(Omega*t);

omega0:=2; gamma0:=0.1; x0:=0; v0:=0;

Eqn_02a:=diff(x(t),t,t)+2*gamma0*diff(x(t),t)+omega0^2*x(t)=Force_02a;
Eqn_02b:=diff(x(t),t,t)+2*gamma0*diff(x(t),t)+omega0^2*x(t)=Force_02b;

sol2a:=dsolve({Eqn_02a,x(0)=x0,D(x)(0)=v0},x(t));
sol2b:=dsolve({Eqn_02b,x(0)=x0,D(x)(0)=v0},x(t));

plot(rhs(sol2a),t=0..60);
plot(Force_02a,t=0..60, color=green, linestyle=1, thickness=3);
plot(rhs(sol2b),t=0..60);
plot(Force_02b,t=0..60, color=green, linestyle=1, thickness=3);

Now we will calculate velocity and acceleration of oscillations and extra functions (amplitudes of the oscillations and damped function):

> x_02a(t):=rhs(sol2a);
v_02a(t):=diff(x_02a(t),t);
w_02a(t):=diff(v_02a(t),t);
A_02a:=sqrt((1625/126483*sqrt(399))^2+(625/317)^2):
x_02expa:=A_02a*exp(-1/10*t):

x_02b(t):=rhs(sol2b);
v_02b(t):=diff(x_02b(t),t);
w_02b(t):=diff(v_02b(t),t);
A_02b:=sqrt((6275/42161*sqrt(399))^2+(75/317)^2):
x_02expb:=A_02b*exp(-1/10*t):

plot([x_02a(t),v_02a(t),x_02expa,-x_02expa,w_02a(t)],t=0..60,color=[red,blue,violet,violet,magenta],linestyle=[1,3,3,3,2],thickness=[4,1,1,1,2]);
plot(Force_02a,t=0..60, color=green, linestyle=1, thickness=3);

plot([x_02b(t),v_02b(t),x_02expb,-x_02expb,w_02b(t)],t=0..60,color=[red,blue,violet,violet,magenta],linestyle=[1,3,3,3,2],thickness=[4,1,1,1,2]);
plot(Force_02b,t=0..60, color=green, linestyle=1, thickness=3);

The next plots are phase curve of oscillations ( magenta color ), variable of oscillations x(t) ( red color ), velocity v(t) ( blue color ) and external force ( green color ):

> G2d_021a := array(1..2):
G2d_021a[1]:=plot([v_02a(t),x_02a(t),t=0..30],color=magenta, axes=box):
G2d_021a[2]:=plot([t,x_02a(t),t=0..30],color=red, axes=box):
display(G2d_021a, linestyle=1, thickness=2, axes=box);
G2d_022a := array(1..2):
G2d_022a[1]:=plot([v_02a(t),-t,t=0..30],color=blue, thickness=2, axes=box):
G2d_022a[2]:=plot([t, Force_02a,t=0..30],color=green, thickness=3, axes=box):
display(G2d_022a, linestyle=1, axes=box);

> G2d_021b := array(1..2):
G2d_021b[1]:=plot([v_02b(t),x_02b(t),t=0..30],color=magenta, axes=box):
G2d_021b[2]:=plot([t,x_02b(t),t=0..30],color=red, axes=box):
display(G2d_021b, linestyle=1, thickness=2, axes=box);
G2d_022b := array(1..2):
G2d_022b[1]:=plot([v_02b(t),-t,t=0..30],color=blue, thickness=2, axes=box):
G2d_022b[2]:=plot([t, Force_02b,t=0..30],color=green, thickness=3, axes=box):
display(G2d_022b, linestyle=1, axes=box);

Next plots are a omparison of the results for a) and b) :

> G2d_021ab := array(1..2):
G2d_021ab[1]:=plot([[v_02a(t),x_02a(t),t=0..20],[v_02b(t),x_02b(t),t=0..20]],color=[magenta,pink], axes=box, linestyle=[1,3]):
G2d_021ab[2]:=plot([[t,x_02a(t),t=0..20], [t,x_02b(t),t=0..20]],color=[red,pink], axes=box,linestyle=[1,3]):
display(G2d_021ab, linestyle=1, thickness=2, axes=box);
G2d_022ab := array(1..2):
G2d_022ab[1]:=plot([[v_02a(t),-t,t=0..20], [v_02b(t),-t,t=0..20]],color=[blue,pink], thickness=2, axes=box,linestyle=[1,3]):
G2d_022ab[2]:=plot([[t, Force_02a,t=0..20],[t, Force_02b,t=0..20]],color=[green,pink], thickness=2, axes=box):
display(G2d_022ab, linestyle=1, axes=box);

At last we represent phase curves of oscillations:

> G3d_02:= array(1..2):
G3d_02[1]:=spacecurve([v_02a(t),x_02a(t),t],t=0..60,numpoints=1000,color=red):
G3d_02[2]:=spacecurve([v_02b(t),x_02b(t),t],t=0..60,numpoints=1000,color=blue):
display(G3d_02,orientation=[0,80],thickness=2,axes=none);

Geometric properties of phase curves (curvature and torsion) correspond to damping characteristics of oscillations. The blue line corresponds to oscillations with Force_02b.

Section 6.3. FORCED DAMPED OSCILLATIONS 2 (PIECEWISE EXTERNAL FORCE)

Consider forced damped oscillations with two different values of initial conditions.

MAPLE executable expressions for

- external force Force_03 ,

- plot of external force Force_03

- proper frequency w0 ,

- damping factor g0 ,

- differential equation of oscillation (in form of the system of two differential equation of the first order),

- two types of initial conditions: a) x0a , v0a , b) x0b , v0b ,

- resulting analytical solutions for two types of initial conditions : a) x_03a(t) , b) x_03b(t) ,

- plots of solutions for two types of initial conditions: a) x_03a(t) , b) x_03b(t) ,

are:

> restart;

> Force_03:=PIECEWISE([15, t <= 7.5],[0, t< 12.5],[15, 12.5 < t]);
plot(Force_03,t=0..60, color=green, linestyle=1, thickness=4);
omega0:=2; gamma0:=0.1;x0a:=0; v0a:=0; x0b:=0; v0b:=12;
System_03:= {v(t)=diff(x(t),t),diff(v(t),t)+2*gamma0*v(t)+omega0^2*x(t)=Force_03};
IC_03a := {x(0)= x0a,v(0)=v0a};
IC_03b := {x(0)= x0b,v(0)=v0b};
sol_03a := combine(dsolve(System_03 union IC_03a,{x(t),v(t)}),trig):
sol_03b := combine(dsolve(System_03 union IC_03b,{x(t),v(t)}),trig);

$System_03 := {diff(v(t),t)+.2*v(t)+4*x(t) = PIECEWISE([15, t <= 7.5],[0, t < 12.5],[15, 12.5 < t]), v(t) = diff(x(t),t)}$

$IC_03a := {x(0) = 0, v(0) = 0}$

$IC_03b := {x(0) = 0, v(0) = 12}$

$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$
$sol_03b := {x(t) = 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 1...$

Now we will calculate velocity and acceleration of oscillations and extra functions (amplitude of proper oscillations and damped function):

> x_03a(t):= -15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))-5/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 15/2],[150*sin(3/4*399^(1/2))*exp(3/4)*399^(1/2), t <= 25/2],[150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2))+150*sin(3/4*399^(1/2))*exp(3/4)*399^(1/2)-150*sin(5/4*399^(1/2))*exp(5/4)*399^(1/2), 25/2 < t])*399^(1/2)*sin(1/10*t*399^(1/2))-10/399*exp(-1/10*t)*399^(1/2)*PIECEWISE([-3/8*399^(1/2)*exp(1/10*t)*cos(1/10*t*399^(1/2))+3/8*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 15/2],[-3/8*399^(1/2)*exp(3/4)*cos(3/4*399^(1/2))+3/8*exp(3/4)*sin(3/4*399^(1/2)), t <= 25/2],[-3/8*399^(1/2)*exp(1/10*t)*cos(1/10*t*399^(1/2))+3/8*exp(1/10*t)*sin(1/10*t*399^(1/2))-3/8*399^(1/2)*exp(3/4)*cos(3/4*399^(1/2))+3/8*exp(3/4)*sin(3/4*399^(1/2))+3/8*exp(5/4)*399^(1/2)*cos(5/4*399^(1/2))-3/8*exp(5/4)*sin(5/4*399^(1/2)), 25/2 < t])*cos(1/10*t*399^(1/2))-10/399*exp(-1/10*t)*sin(1/10*t*399^(1/2))*PIECEWISE([-3/8*399^(1/2)*exp(1/10*t)*cos(1/10*t*399^(1/2))+3/8*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 15/2],[-3/8*399^(1/2)*exp(3/4)*cos(3/4*399^(1/2))+3/8*exp(3/4)*sin(3/4*399^(1/2)), t <= 25/2],[-3/8*399^(1/2)*exp(1/10*t)*cos(1/10*t*399^(1/2))+3/8*exp(1/10*t)*sin(1/10*t*399^(1/2))-3/8*399^(1/2)*exp(3/4)*cos(3/4*399^(1/2))+3/8*exp(3/4)*sin(3/4*399^(1/2))+3/8*exp(5/4)*399^(1/2)*cos(5/4*399^(1/2))-3/8*exp(5/4)*sin(5/4*399^(1/2)), 25/2 < t]):
v_03a(t):=diff(x_03a(t),t):
w_03a(t):=diff(v_03a(t),t):

x_03b(t):= 155/532*exp(-1/10*t)*399^(1/2)*sin(1/10*t*399^(1/2))-15/4*exp(-1/10*t)*cos(1/10*t*399^(1/2))+10/159201*exp(-1/10*t)*PIECEWISE([150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 15/2],[150*sin(3/4*399^(1/2))*exp(3/4)*399^(1/2), t <= 25/2],[150*399^(1/2)*exp(1/10*t)*sin(1/10*t*399^(1/2))+150*sin(3/4*399^(1/2))*exp(3/4)*399^(1/2)-150*sin(5/4*399^(1/2))*exp(5/4)*399^(1/2), 25/2 < t])*399^(1/2)*sin(1/10*t*399^(1/2))-10/399*exp(-1/10*t)*399^(1/2)*PIECEWISE([-3/8*399^(1/2)*exp(1/10*t)*cos(1/10*t*399^(1/2))+3/8*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 15/2],[-3/8*399^(1/2)*exp(3/4)*cos(3/4*399^(1/2))+3/8*exp(3/4)*sin(3/4*399^(1/2)), t <= 25/2],[-3/8*399^(1/2)*exp(1/10*t)*cos(1/10*t*399^(1/2))+3/8*exp(1/10*t)*sin(1/10*t*399^(1/2))-3/8*399^(1/2)*exp(3/4)*cos(3/4*399^(1/2))+3/8*exp(3/4)*sin(3/4*399^(1/2))+3/8*exp(5/4)*399^(1/2)*cos(5/4*399^(1/2))-3/8*exp(5/4)*sin(5/4*399^(1/2)), 25/2 < t])*cos(1/10*t*399^(1/2))-10/399*exp(-1/10*t)*sin(1/10*t*399^(1/2))*PIECEWISE([-3/8*399^(1/2)*exp(1/10*t)*cos(1/10*t*399^(1/2))+3/8*exp(1/10*t)*sin(1/10*t*399^(1/2)), t <= 15/2],[-3/8*399^(1/2)*exp(3/4)*cos(3/4*399^(1/2))+3/8*exp(3/4)*sin(3/4*399^(1/2)), t <= 25/2],[-3/8*399^(1/2)*exp(1/10*t)*cos(1/10*t*399^(1/2))+3/8*exp(1/10*t)*sin(1/10*t*399^(1/2))-3/8*399^(1/2)*exp(3/4)*cos(3/4*399^(1/2))+3/8*exp(3/4)*sin(3/4*399^(1/2))+3/8*exp(5/4)*399^(1/2)*cos(5/4*399^(1/2))-3/8*exp(5/4)*sin(5/4*399^(1/2)), 25/2 < t]):
v_03b(t):=diff(x_03b(t),t):
w_03b(t):=diff(v_03b(t),t):

plot([x_03a(t),v_03a(t),w_03a(t)],t=0..60,color=[red,blue,magenta],linestyle=[1,3,3],thickness=[4,1,1]);
plot(Force_03,t=0..60, color=green, linestyle=1, thickness=3);

plot([x_03b(t),v_03b(t),w_03b(t)],t=0..60,color=[red,blue,magenta],linestyle=[1,3,3],thickness=[4,1,1]);

Next plots are phase curve of oscillations ( magenta color ), variable of oscillations x(t) ( red color ), velocity v(t) ( blue color ) and external force ( green color ):

> G2d_031a := array(1..2):
G2d_031a[1]:=plot([v_03a(t),x_03a(t),t=0..40],color=magenta, axes=box):
G2d_031a[2]:=plot([t,x_03a(t),t=0..40],color=red, axes=box):
display(G2d_031a, linestyle=1, thickness=2, axes=box);
G2d_032a := array(1..2):
G2d_032a[1]:=plot([v_03a(t),-t,t=0..40],color=blue, thickness=2, axes=box):
G2d_032a[2]:=plot([t, Force_03,t=0..40],color=green, thickness=4, axes=box):
display(G2d_032a, linestyle=1, axes=box);

> G2d_031b := array(1..2):
G2d_031b[1]:=plot([v_03b(t),x_03b(t),t=0..40],color=magenta, axes=box):
G2d_031b[2]:=plot([t,x_03b(t),t=0..40],color=red, axes=box):
display(G2d_031b, linestyle=1, thickness=2, axes=box);
G2d_032b := array(1..2):
G2d_032b[1]:=plot([v_03b(t),-t,t=0..40],color=blue, thickness=2, axes=box):
G2d_032b[2]:=plot([t, Force_03,t=0..40],color=green, thickness=4, axes=box):
display(G2d_032b, linestyle=1, axes=box);

Next plots are comparison of results for a) and b) situations:

> G2d_031ab := array(1..2):
G2d_031ab[1]:=plot([[v_03a(t),x_03a(t),t=0..30],[v_03b(t),x_03b(t),t=0..30]],color=[magenta,pink], axes=box, linestyle=[1,3]):
G2d_031ab[2]:=plot([[t,x_03a(t),t=0..30], [t,x_03b(t),t=0..30]],color=[red,pink], axes=box,linestyle=[1,3]):
display(G2d_031ab, linestyle=1, thickness=2, axes=box);
G2d_032ab := array(1..2):
G2d_032ab[1]:=plot([[v_03a(t),-t,t=0..30], [v_03b(t),-t,t=0..30]],color=[blue,pink], thickness=2, axes=box,linestyle=[1,3]):
G2d_032ab[2]:=plot([t, Force_03,t=0..30],color=green, thickness=4, axes=box):
display(G2d_032ab, linestyle=1, axes=box);

At last we represent phase curves of oscillations:

> G3d_03:= array(1..2):
G3d_03[1]:=spacecurve([v_03a(t),x_03a(t),t],t=0..60,numpoints=1000,color=red):
G3d_03[2]:=spacecurve([v_03b(t),x_03b(t),t],t=0..60,numpoints=1000,color=blue):
display(G3d_03,orientation=[0,80],thickness=2,axes=none);

The blue curve corresponds to oscillations with initial conditions b).

Section 6.4. FORCED DAMPED OSCILLATIONS 3 (PIECEWISE EXTERNAL FORCE)

Consider forced damped oscillations with two different values of initial conditions.

MAPLE executable expressions for

- external force Force_0 4,

- plot of external force Force_04

- proper frequency w0 ,

- damping factor g0 ,

- differential equation of oscillation (in form of the system of two differential equation of the first order),

- initial conditions: x0 , v0 ,

- resulting analytical solution,

are:

> Force_04:=PIECEWISE([0, t <= 5],[10, t <= 10],[-10, t <= 15],[10, t <= 20],[-10, 20 < t]);
plot(Force_04,t=0..60, color=green, linestyle=1, thickness=4);
x0:=0;v0:=0; omega0:=2; gamma0:=0.1;
System_04:= {v(t)=diff(x(t),t),diff(v(t),t)+2*gamma0*v(t)+omega0^2*x(t)=Force_04};
IC_04 := {x(0)= x0,v(0)=v0};
Solution_04 := combine(dsolve(System_04 union IC_04,{x(t),v(t)}),trig):

$System_04 := {diff(v(t),t)+.2*v(t)+4*x(t) = PIECEWISE([0, t <= 5],[10, t <= 10],[-10, t <= 15],[10, t <= 20],[-10, 20 < t]), v(t) = diff(x(t),t)}$

$IC_04 := {x(0) = 0, v(0) = 0}$

> x_04(t):=-10/399*exp(-1/10*t)*sqrt(399)*PIECEWISE([0, t <= 5],[-1/4*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/4*exp(1/10*t)*sin(1/10*t*sqrt(399))+1/4*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399))-1/4*exp(1/2)*sin(1/2*sqrt(399)), t <= 10],[1/4*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))-1/4*exp(1/10*t)*sin(1/10*t*sqrt(399))-1/2*sqrt(399)*exp(1)*cos(sqrt(399))-1/4*exp(1/2)*sin(1/2*sqrt(399))+1/2*exp(1)*sin(sqrt(399))+1/4*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399)), t <= 15],[-1/4*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/4*exp(1/10*t)*sin(1/10*t*sqrt(399))+1/2*sqrt(399)*exp(3/2)*cos(3/2*sqrt(399))-1/4*exp(1/2)*sin(1/2*sqrt(399))-1/2*exp(3/2)*sin(3/2*sqrt(399))+1/4*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399))-1/2*sqrt(399)*exp(1)*cos(sqrt(399))+1/2*exp(1)*sin(sqrt(399)), t < 20],[-1/4*sqrt(399)*exp(2)*cos(2*sqrt(399))-1/4*exp(1/2)*sin(1/2*sqrt(399))+1/4*exp(2)*sin(2*sqrt(399))+1/4*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399))+1/2*sqrt(399)*exp(3/2)*cos(3/2*sqrt(399))+1/2*exp(1)*sin(sqrt(399))-1/2*exp(3/2)*sin(3/2*sqrt(399))-1/2*sqrt(399)*exp(1)*cos(sqrt(399)), t = 20],[1/4*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))-1/4*exp(1/10*t)*sin(1/10*t*sqrt(399))-1/2*sqrt(399)*exp(2)*cos(2*sqrt(399))-1/4*exp(1/2)*sin(1/2*sqrt(399))+1/2*exp(2)*sin(2*sqrt(399))+1/4*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399))+1/2*sqrt(399)*exp(3/2)*cos(3/2*sqrt(399))+1/2*exp(1)*sin(sqrt(399))-1/2*exp(3/2)*sin(3/2*sqrt(399))-1/2*sqrt(399)*exp(1)*cos(sqrt(399)), 20 < t])*cos(1/10*t*sqrt(399))-10/399*exp(-1/10*t)*sin(1/10*t*sqrt(399))*PIECEWISE([0, t <= 5],[-1/4*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/4*exp(1/10*t)*sin(1/10*t*sqrt(399))+1/4*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399))-1/4*exp(1/2)*sin(1/2*sqrt(399)), t <= 10],[1/4*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))-1/4*exp(1/10*t)*sin(1/10*t*sqrt(399))-1/2*sqrt(399)*exp(1)*cos(sqrt(399))-1/4*exp(1/2)*sin(1/2*sqrt(399))+1/2*exp(1)*sin(sqrt(399))+1/4*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399)), t <= 15],[-1/4*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/4*exp(1/10*t)*sin(1/10*t*sqrt(399))+1/2*sqrt(399)*exp(3/2)*cos(3/2*sqrt(399))-1/4*exp(1/2)*sin(1/2*sqrt(399))-1/2*exp(3/2)*sin(3/2*sqrt(399))+1/4*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399))-1/2*sqrt(399)*exp(1)*cos(sqrt(399))+1/2*exp(1)*sin(sqrt(399)), t < 20],[-1/4*sqrt(399)*exp(2)*cos(2*sqrt(399))-1/4*exp(1/2)*sin(1/2*sqrt(399))+1/4*exp(2)*sin(2*sqrt(399))+1/4*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399))+1/2*sqrt(399)*exp(3/2)*cos(3/2*sqrt(399))+1/2*exp(1)*sin(sqrt(399))-1/2*exp(3/2)*sin(3/2*sqrt(399))-1/2*sqrt(399)*exp(1)*cos(sqrt(399)), t = 20],[1/4*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))-1/4*exp(1/10*t)*sin(1/10*t*sqrt(399))-1/2*sqrt(399)*exp(2)*cos(2*sqrt(399))-1/4*exp(1/2)*sin(1/2*sqrt(399))+1/2*exp(2)*sin(2*sqrt(399))+1/4*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399))+1/2*sqrt(399)*exp(3/2)*cos(3/2*sqrt(399))+1/2*exp(1)*sin(sqrt(399))-1/2*exp(3/2)*sin(3/2*sqrt(399))-1/2*sqrt(399)*exp(1)*cos(sqrt(399)), 20 < t])+10/159201*exp(-1/10*t)*sqrt(399)*sin(1/10*t*sqrt(399))*PIECEWISE([0, t <= 5],[100*sqrt(399)*exp(1/10*t)*sin(1/10*t*sqrt(399))-100*exp(1/2)*sqrt(399)*sin(1/2*sqrt(399)), t <= 10],[-100*sqrt(399)*exp(1/10*t)*sin(1/10*t*sqrt(399))+200*exp(1)*sqrt(399)*sin(sqrt(399))-100*exp(1/2)*sqrt(399)*sin(1/2*sqrt(399)), t <= 15],[100*sqrt(399)*exp(1/10*t)*sin(1/10*t*sqrt(399))-200*exp(3/2)*sqrt(399)*sin(3/2*sqrt(399))-100*exp(1/2)*sqrt(399)*sin(1/2*sqrt(399))+200*exp(1)*sqrt(399)*sin(sqrt(399)), t < 20],[100*sqrt(399)*exp(2)*sin(2*sqrt(399))-100*exp(1/2)*sqrt(399)*sin(1/2*sqrt(399))-200*exp(3/2)*sqrt(399)*sin(3/2*sqrt(399))+200*exp(1)*sqrt(399)*sin(sqrt(399)), t = 20],[-100*sqrt(399)*exp(1/10*t)*sin(1/10*t*sqrt(399))+200*sqrt(399)*exp(2)*sin(2*sqrt(399))-100*exp(1/2)*sqrt(399)*sin(1/2*sqrt(399))-200*exp(3/2)*sqrt(399)*sin(3/2*sqrt(399))+200*exp(1)*sqrt(399)*sin(sqrt(399)), 20 < t]):
v_04(t):=diff(x_04(t),t):
w_04(t):=diff(v_04(t),t):
plot([x_04(t),v_04(t),w_04(t)],t=0..60,color=[red,blue,magenta],linestyle=[1,3,3],thickness=[4,1,1]);
plot(Force_04,t=0..60, color=green, linestyle=1, thickness=3);

Next plots are phase curve of oscillations ( magenta color ), variable of oscillations x(t) ( red color ), velocity v(t) ( blue color ) and external force ( green color ):

> with(plots):
G2d_041 := array(1..2):
G2d_041[1]:=plot([v_04(t),x_04(t),t=0..60],color=magenta, axes=box):
G2d_041[2]:=plot([t,x_04(t),t=0..60],color=red, axes=box):
display(G2d_041, linestyle=1, thickness=2, axes=box);
G2d_042 := array(1..2):
G2d_042[1]:=plot([v_04(t),-t,t=0..60],color=blue, thickness=2, axes=box):
G2d_042[2]:=plot([t, Force_04,t=0..60],color=green, thickness=4, axes=box):
display(G2d_042, linestyle=1, axes=box);

Warning, the name changecoords has been redefined

Section 6.5. FORCED DAMPED OSCILLATIONS 4 (PIECEWISE EXTERNAL FORCE)

Consider forced damped oscillations with two different values of initial conditions.

MAPLE executable expression for

- external force Force_05 ,

- plot of external force Force_05 ,

- proper frequency w0 ,

- damping factor g0 ,

- differential equation of oscillation (in form of the system of two differential equation of the first order),

- two types of initial conditions: a) x0a , v0a , b) x0b , v0b ,

- resulting analytical solutions for two types of initial conditions : a) x_05a(t) , b) x_05b(t) ,

are:

> Force_05:=PIECEWISE([t, t <= 10],[20-t, t <= 20],[0, 20 < t]);
plot(Force_05,t=0..40, color=green, linestyle=1, thickness=4);
omega0:=2; gamma0:=0.1; x0a:=0; v0a:=0; x0b:=0; v0b:=2;
System_05:= {v(t)=diff(x(t),t),diff(v(t),t)+2*gamma0*v(t)+omega0^2*x(t)=Force_05};
IC_05a := {x(0)= x0a,v(0)=v0a};
IC_05b := {x(0)= x0b,v(0)=v0b};
Solution_05a := combine(dsolve(System_05 union IC_05a,{x(t),v(t)}),trig):
Solution_05b := combine(dsolve(System_05 union IC_05b,{x(t),v(t)}),trig):

$System_05 := {v(t) = diff(x(t),t), diff(v(t),t)+.2*v(t)+4*x(t) = PIECEWISE([t, t <= 10],[20-t, t <= 20],[0, 20 < t])}$

$IC_05a := {x(0) = 0, v(0) = 0}$

$IC_05b := {x(0) = 0, v(0) = 2}$

> x_05a(t):=1/80*exp(-1/10*t)*cos(1/10*t*sqrt(399))-199/31920*exp(-1/10*t)*sqrt(399)*sin(1/10*t*sqrt(399))-10/399*exp(-1/10*t)*sqrt(399)*PIECEWISE([-1/40*exp(1/10*t)*cos(1/10*t*sqrt(399))*t*sqrt(399)+1/800*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/40*exp(1/10*t)*sin(1/10*t*sqrt(399))*t+199/800*exp(1/10*t)*sin(1/10*t*sqrt(399)), t <= 10],[-401/800*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/40*exp(1/10*t)*cos(1/10*t*sqrt(399))*t*sqrt(399)+201/800*exp(1/10*t)*sin(1/10*t*sqrt(399))-1/40*exp(1/10*t)*sin(1/10*t*sqrt(399))*t+1/400*sqrt(399)*exp(1)*cos(sqrt(399))+199/400*exp(1)*sin(sqrt(399)), t < 20],[-1/800*sqrt(399)*exp(2)*cos(2*sqrt(399))+199/400*exp(1)*sin(sqrt(399))-199/800*exp(2)*sin(2*sqrt(399))+1/400*sqrt(399)*exp(1)*cos(sqrt(399)), 20 <= t])*cos(1/10*t*sqrt(399))-10/399*exp(-1/10*t)*PIECEWISE([-1/40*exp(1/10*t)*cos(1/10*t*sqrt(399))*t*sqrt(399)+1/800*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/40*exp(1/10*t)*sin(1/10*t*sqrt(399))*t+199/800*exp(1/10*t)*sin(1/10*t*sqrt(399)), t <= 10],[-401/800*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/40*exp(1/10*t)*cos(1/10*t*sqrt(399))*t*sqrt(399)+201/800*exp(1/10*t)*sin(1/10*t*sqrt(399))-1/40*exp(1/10*t)*sin(1/10*t*sqrt(399))*t+1/400*sqrt(399)*exp(1)*cos(sqrt(399))+199/400*exp(1)*sin(sqrt(399)), t < 20],[-1/800*sqrt(399)*exp(2)*cos(2*sqrt(399))+199/400*exp(1)*sin(sqrt(399))-199/800*exp(2)*sin(2*sqrt(399))+1/400*sqrt(399)*exp(1)*cos(sqrt(399)), 20 <= t])*sin(1/10*t*sqrt(399))+10/159201*exp(-1/10*t)*sqrt(399)*sin(1/10*t*sqrt(399))*PIECEWISE([399/4*exp(1/10*t)*cos(1/10*t*sqrt(399))+10*exp(1/10*t)*sin(1/10*t*sqrt(399))*t*sqrt(399)-1/4*sqrt(399)*exp(1/10*t)*sin(1/10*t*sqrt(399)), t <= 10],[-399/4*exp(1/10*t)*cos(1/10*t*sqrt(399))+801/4*sqrt(399)*exp(1/10*t)*sin(1/10*t*sqrt(399))-10*exp(1/10*t)*sin(1/10*t*sqrt(399))*t*sqrt(399)+399/2*exp(1)*cos(sqrt(399))-1/2*exp(1)*sqrt(399)*sin(sqrt(399)), t < 20],[-399/4*exp(2)*cos(2*sqrt(399))-1/2*exp(1)*sqrt(399)*sin(sqrt(399))+1/4*sqrt(399)*exp(2)*sin(2*sqrt(399))+399/2*exp(1)*cos(sqrt(399)), 20 <= t]):
v_05a(t):=diff(x_05a(t),t):
w_05a(t):=diff(v_05a(t),t):

x_05b(t):=1/80*exp(-1/10*t)*cos(1/10*t*sqrt(399))+467/10640*exp(-1/10*t)*sqrt(399)*sin(1/10*t*sqrt(399))-10/399*exp(-1/10*t)*sqrt(399)*PIECEWISE([-1/40*exp(1/10*t)*cos(1/10*t*sqrt(399))*t*sqrt(399)+1/800*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/40*exp(1/10*t)*sin(1/10*t*sqrt(399))*t+199/800*exp(1/10*t)*sin(1/10*t*sqrt(399)), t <= 10],[-401/800*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/40*exp(1/10*t)*cos(1/10*t*sqrt(399))*t*sqrt(399)+201/800*exp(1/10*t)*sin(1/10*t*sqrt(399))-1/40*exp(1/10*t)*sin(1/10*t*sqrt(399))*t+1/400*sqrt(399)*exp(1)*cos(sqrt(399))+199/400*exp(1)*sin(sqrt(399)), t < 20],[-1/800*sqrt(399)*exp(2)*cos(2*sqrt(399))+199/400*exp(1)*sin(sqrt(399))-199/800*exp(2)*sin(2*sqrt(399))+1/400*sqrt(399)*exp(1)*cos(sqrt(399)), 20 <= t])*cos(1/10*t*sqrt(399))-10/399*exp(-1/10*t)*PIECEWISE([-1/40*exp(1/10*t)*cos(1/10*t*sqrt(399))*t*sqrt(399)+1/800*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/40*exp(1/10*t)*sin(1/10*t*sqrt(399))*t+199/800*exp(1/10*t)*sin(1/10*t*sqrt(399)), t <= 10],[-401/800*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/40*exp(1/10*t)*cos(1/10*t*sqrt(399))*t*sqrt(399)+201/800*exp(1/10*t)*sin(1/10*t*sqrt(399))-1/40*exp(1/10*t)*sin(1/10*t*sqrt(399))*t+1/400*sqrt(399)*exp(1)*cos(sqrt(399))+199/400*exp(1)*sin(sqrt(399)), t < 20],[-1/800*sqrt(399)*exp(2)*cos(2*sqrt(399))+199/400*exp(1)*sin(sqrt(399))-199/800*exp(2)*sin(2*sqrt(399))+1/400*sqrt(399)*exp(1)*cos(sqrt(399)), 20 <= t])*sin(1/10*t*sqrt(399))+10/159201*exp(-1/10*t)*sqrt(399)*sin(1/10*t*sqrt(399))*PIECEWISE([399/4*exp(1/10*t)*cos(1/10*t*sqrt(399))+10*exp(1/10*t)*sin(1/10*t*sqrt(399))*t*sqrt(399)-1/4*sqrt(399)*exp(1/10*t)*sin(1/10*t*sqrt(399)), t <= 10],[-399/4*exp(1/10*t)*cos(1/10*t*sqrt(399))+801/4*sqrt(399)*exp(1/10*t)*sin(1/10*t*sqrt(399))-10*exp(1/10*t)*sin(1/10*t*sqrt(399))*t*sqrt(399)+399/2*exp(1)*cos(sqrt(399))-1/2*exp(1)*sqrt(399)*sin(sqrt(399)), t < 20],[-399/4*exp(2)*cos(2*sqrt(399))-1/2*exp(1)*sqrt(399)*sin(sqrt(399))+1/4*sqrt(399)*exp(2)*sin(2*sqrt(399))+399/2*exp(1)*cos(sqrt(399)), 20 <= t]):
v_05b(t):=diff(x_05b(t),t):
w_05b(t):=diff(v_05b(t),t):

plot([x_05a(t),v_05a(t),w_05a(t)],t=0..40,color=[red,blue,magenta],linestyle=[1,3,1],thickness=[4,1,1]);
plot(Force_05,t=0..40, color=green, linestyle=1, thickness=3);
plot([x_05b(t),v_05b(t),w_05b(t)],t=0..40,color=[red,blue,magenta],linestyle=[1,3,1],thickness=[4,1,1]);

Next plots are phase curve of oscillations ( magenta color ), variable of oscillations x(t) ( red color ), velocity v(t) ( blue color ) and external force ( green color ):

> G2d_051a := array(1..2):
G2d_051a[1]:=plot([v_05a(t),x_05a(t),t=0..40],color=magenta, axes=box):
G2d_051a[2]:=plot([t,x_05a(t),t=0..40],color=red, axes=box):
display(G2d_051a, linestyle=1, thickness=2, axes=box);
G2d_052a := array(1..2):
G2d_052a[1]:=plot([v_05a(t),-t,t=0..40],color=blue, thickness=2, axes=box):
G2d_052a[2]:=plot([t, Force_05,t=0..40],color=green, thickness=4, axes=box):
display(G2d_052a, linestyle=1, axes=box);

> G2d_051b := array(1..2):
G2d_051b[1]:=plot([v_05b(t),x_05b(t),t=0..40],color=magenta, axes=box):
G2d_051b[2]:=plot([t,x_05b(t),t=0..40],color=red, axes=box):
display(G2d_051b, linestyle=1, thickness=2, axes=box);
G2d_052b := array(1..2):
G2d_052b[1]:=plot([v_05b(t),-t,t=0..40],color=blue, thickness=2, axes=box):
G2d_052b[2]:=plot([t, Force_05,t=0..40],color=green, thickness=4, axes=box):
display(G2d_052b, linestyle=1, axes=box);

Next plots are comparison of results for a) and b) situations:

> G2d_051ab := array(1..2):
G2d_051ab[1]:=plot([[v_05a(t),x_05a(t),t=0..30],[v_05b(t),x_05b(t),t=0..30]],color=[magenta,pink], axes=box, linestyle=[1,3]):
G2d_051ab[2]:=plot([[t,x_05a(t),t=0..30], [t,x_05b(t),t=0..30]],color=[red,pink], axes=box,linestyle=[1,3]):
display(G2d_051ab, linestyle=1, thickness=2, axes=box);
G2d_052ab := array(1..2):
G2d_052ab[1]:=plot([[v_05a(t),-t,t=0..30], [v_05b(t),-t,t=0..30]],color=[blue,pink], thickness=2, axes=box,linestyle=[1,3]):
G2d_052ab[2]:=plot([t, Force_05,t=0..30],color=green, thickness=4, axes=box):
display(G2d_052ab, linestyle=1, axes=box);

At last we represent phase curves of oscillations with time evolvent:

> G3d_05:= array(1..2):
G3d_05[1]:=spacecurve([v_05a(t),x_05a(t),t],t=0..40,numpoints=1000,color=red):
G3d_05[2]:=spacecurve([v_05b(t),x_05b(t),t],t=0..40,numpoints=1000,color=blue):
display(G3d_05,orientation=[10,80],thickness=2,axes=none);

Section 6.6. FORCED DAMPED OSCILLATIONS 5 (PIECEWISE EXTERNAL FORCE)

Consider forced damped oscillations with two external forces.

MAPLE executable expression for

- external forces Force_06a and Force_06b ,

- plots of external forces Force_06a and Force_06b,

- proper frequency w0 ,

- damping factor g0 ,

- differential equations of oscillations (in form of the system of two differential equation of the first order),

- initial conditions: x0 , v0 ,

- resulting analytical solutions for two external forces: a) x_06a(t) , b) x_06b(t)

are:

> Force_06a:=PIECEWISE([0, t <= 5],[t^2, t <= 15],[0, 15 < t]);
Force_06b:=PIECEWISE([0, t <= 5],[-0.525*t^2, t <= 25],[0, 25 < t]);
plot([Force_06a,Force_06b],t=0..60, color=green, linestyle=[1,3], thickness=4);
x0:=0; v0:=5; omega0:=2; gamma0:=0.1;
System_06a:= {v(t)=diff(x(t),t),diff(v(t),t)+2*gamma0*v(t)+omega0^2*x(t)=Force_06a};
System_06b:= {v(t)=diff(x(t),t),diff(v(t),t)+2*gamma0*v(t)+omega0^2*x(t)=Force_06b};
IC_06 := {x(0)= x0,v(0)=v0};
Solution_06a := dsolve(System_06a union IC_06,{x(t),v(t)}):
Solution_06b := dsolve(System_06b union IC_06,{x(t),v(t)}):

$System_06a := {v(t) = diff(x(t),t), diff(v(t),t)+.2*v(t)+4*x(t) = PIECEWISE([0, t <= 5],[t^2, t <= 15],[0, 15 < t])}$

$System_06b := {diff(v(t),t)+.2*v(t)+4*x(t) = PIECEWISE([0, t <= 5],[-.525*t^2, t <= 25],[0, 25 < t]), v(t) = diff(x(t),t)}$

$IC_06 := {v(0) = 5, x(0) = 0}$

> x_06a(t):=50/399*exp(-1/10*t)*sqrt(399)*sin(1/10*t*sqrt(399))-10/399*exp(-1/10*t)*sqrt(399)*PIECEWISE([0, t <= 5],[-1/40*exp(1/10*t)*cos(1/10*t*sqrt(399))*sqrt(399)*t^2+1/400*exp(1/10*t)*cos(1/10*t*sqrt(399))*t*sqrt(399)+99/8000*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/40*exp(1/10*t)*sin(1/10*t*sqrt(399))*t^2+199/400*exp(1/10*t)*sin(1/10*t*sqrt(399))*t-299/8000*exp(1/10*t)*sin(1/10*t*sqrt(399))+4801/8000*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399))-24601/8000*exp(1/2)*sin(1/2*sqrt(399)), t <= 15],[-44601/8000*sqrt(399)*exp(3/2)*cos(3/2*sqrt(399))-24601/8000*exp(1/2)*sin(1/2*sqrt(399))+104401/8000*exp(3/2)*sin(3/2*sqrt(399))+4801/8000*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399)), 15 < t])*cos(1/10*t*sqrt(399))-10/399*exp(-1/10*t)*PIECEWISE([0, t <= 5],[-1/40*exp(1/10*t)*cos(1/10*t*sqrt(399))*sqrt(399)*t^2+1/400*exp(1/10*t)*cos(1/10*t*sqrt(399))*t*sqrt(399)+99/8000*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))+1/40*exp(1/10*t)*sin(1/10*t*sqrt(399))*t^2+199/400*exp(1/10*t)*sin(1/10*t*sqrt(399))*t-299/8000*exp(1/10*t)*sin(1/10*t*sqrt(399))+4801/8000*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399))-24601/8000*exp(1/2)*sin(1/2*sqrt(399)), t <= 15],[-44601/8000*sqrt(399)*exp(3/2)*cos(3/2*sqrt(399))-24601/8000*exp(1/2)*sin(1/2*sqrt(399))+104401/8000*exp(3/2)*sin(3/2*sqrt(399))+4801/8000*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399)), 15 < t])*sin(1/10*t*sqrt(399))+10/159201*exp(-1/10*t)*sqrt(399)*sin(1/10*t*sqrt(399))*PIECEWISE([0, t <= 5],[399/2*exp(1/10*t)*cos(1/10*t*sqrt(399))*t-399/40*exp(1/10*t)*cos(1/10*t*sqrt(399))+10*exp(1/10*t)*sin(1/10*t*sqrt(399))*sqrt(399)*t^2-1/2*exp(1/10*t)*sin(1/10*t*sqrt(399))*t*sqrt(399)-199/40*sqrt(399)*exp(1/10*t)*sin(1/10*t*sqrt(399))-39501/40*cos(1/2*sqrt(399))*exp(1/2)-9701/40*sqrt(399)*sin(1/2*sqrt(399))*exp(1/2), t <= 15],[119301/40*cos(3/2*sqrt(399))*exp(3/2)+89501/40*sqrt(399)*sin(3/2*sqrt(399))*exp(3/2)-9701/40*sqrt(399)*sin(1/2*sqrt(399))*exp(1/2)-39501/40*cos(1/2*sqrt(399))*exp(1/2), 15 < t]):
v_06a(t):=diff(x_06a(t),t):
w_06a(t):=diff(v_06a(t),t):

x_06b(t):=50/399*exp(-1/10*t)*sqrt(399)*sin(1/10*t*sqrt(399))-10/399*exp(-1/10*t)*sqrt(399)*PIECEWISE([0, t <= 5],[21/1600*exp(1/10*t)*cos(1/10*t*sqrt(399))*sqrt(399)*t^2-21/16000*exp(1/10*t)*cos(1/10*t*sqrt(399))*t*sqrt(399)-2079/320000*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))-21/1600*exp(1/10*t)*sin(1/10*t*sqrt(399))*t^2-4179/16000*exp(1/10*t)*sin(1/10*t*sqrt(399))*t+6279/320000*exp(1/10*t)*sin(1/10*t*sqrt(399))-100821/320000*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399))+516621/320000*exp(1/2)*sin(1/2*sqrt(399)), t <= 25],[2612421/320000*sqrt(399)*cos(5/2*sqrt(399))*exp(5/2)+516621/320000*exp(1/2)*sin(1/2*sqrt(399))-4708221/320000*sin(5/2*sqrt(399))*exp(5/2)-100821/320000*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399)), 25 < t])*cos(1/10*t*sqrt(399))-10/399*PIECEWISE([0, t <= 5],[21/1600*exp(1/10*t)*cos(1/10*t*sqrt(399))*sqrt(399)*t^2-21/16000*exp(1/10*t)*cos(1/10*t*sqrt(399))*t*sqrt(399)-2079/320000*sqrt(399)*exp(1/10*t)*cos(1/10*t*sqrt(399))-21/1600*exp(1/10*t)*sin(1/10*t*sqrt(399))*t^2-4179/16000*exp(1/10*t)*sin(1/10*t*sqrt(399))*t+6279/320000*exp(1/10*t)*sin(1/10*t*sqrt(399))-100821/320000*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399))+516621/320000*exp(1/2)*sin(1/2*sqrt(399)), t <= 25],[2612421/320000*sqrt(399)*cos(5/2*sqrt(399))*exp(5/2)+516621/320000*exp(1/2)*sin(1/2*sqrt(399))-4708221/320000*sin(5/2*sqrt(399))*exp(5/2)-100821/320000*sqrt(399)*exp(1/2)*cos(1/2*sqrt(399)), 25 < t])*exp(-1/10*t)*sin(1/10*t*sqrt(399))+10/159201*exp(-1/10*t)*sqrt(399)*sin(1/10*t*sqrt(399))*PIECEWISE([0, t <= 5],[-8379/80*exp(1/10*t)*cos(1/10*t*sqrt(399))*t+8379/1600*exp(1/10*t)*cos(1/10*t*sqrt(399))-21/4*exp(1/10*t)*sin(1/10*t*sqrt(399))*sqrt(399)*t^2+21/80*exp(1/10*t)*sin(1/10*t*sqrt(399))*t*sqrt(399)+4179/1600*sqrt(399)*exp(1/10*t)*sin(1/10*t*sqrt(399))+829521/1600*cos(1/2*sqrt(399))*exp(1/2)+203721/1600*sqrt(399)*sin(1/2*sqrt(399))*exp(1/2), t <= 25],[-4181121/1600*cos(5/2*sqrt(399))*exp(5/2)-5235321/1600*sqrt(399)*sin(5/2*sqrt(399))*exp(5/2)+203721/1600*sqrt(399)*sin(1/2*sqrt(399))*exp(1/2)+829521/1600*cos(1/2*sqrt(399))*exp(1/2), 25 < t]):
v_06b(t):=diff(x_06b(t),t):
w_06b(t):=diff(v_06b(t),t):

plot([x_06a(t),v_06a(t),w_06b(t)],t=0..60,color=[red,blue,magenta],linestyle=[1,3,1],thickness=[4,1,1]);
plot([Force_06a,Force_06b],t=0..60, color=green, linestyle=[1,3], thickness=3);

plot([x_06b(t),v_06b(t),w_06b(t)],t=0..60,color=[red,blue,magenta],linestyle=[1,3,1],thickness=[4,1,1]);

Next plots are phase curve of oscillations ( magenta color ), variable of oscillations x(t) ( red color ), velocity v(t) ( blue color ) and external force ( green color ):

> G2d_061a := array(1..2):
G2d_061a[1]:=plot([v_06a(t),x_06a(t),t=0..60],color=magenta, axes=box):
G2d_061a[2]:=plot([t,x_06a(t),t=0..60],color=red, axes=box):
display(G2d_061a, linestyle=1, thickness=2, axes=box);
G2d_062a := array(1..2):
G2d_062a[1]:=plot([v_06a(t),-t,t=0..60],color=blue, thickness=2, axes=box):
G2d_062a[2]:=plot([t, Force_06a,t=0..60],color=green, thickness=4, axes=box):
display(G2d_062a, linestyle=1, axes=box);

> G2d_061b := array(1..2):
G2d_061b[1]:=plot([v_06b(t),x_06b(t),t=0..60],color=magenta, axes=box):
G2d_061b[2]:=plot([t,x_06b(t),t=0..60],color=red, axes=box):
display(G2d_061b, linestyle=1, thickness=2, axes=box);
G2d_062b := array(1..2):
G2d_062b[1]:=plot([v_06b(t),-t,t=0..60],color=blue, thickness=2, axes=box):
G2d_062b[2]:=plot([t, Force_06b,t=0..60],color=green, thickness=4, axes=box):
display(G2d_062b, linestyle=1, axes=box);

Next plots are comparison of results for a) and b) situations:

> G2d_061ab := array(1..2):
G2d_061ab[1]:=plot([[v_06a(t),x_06a(t),t=0..60],[v_06b(t),x_06b(t),t=0..60]],color=[magenta,pink], axes=box, linestyle=[1,3]):
G2d_061ab[2]:=plot([[t,x_06a(t),t=0..60], [t,x_06b(t),t=0..60]],color=[red,pink], axes=box,linestyle=[1,3]):
display(G2d_061ab, linestyle=1, thickness=2, axes=box);
G2d_062ab := array(1..2):
G2d_062ab[1]:=plot([[v_06a(t),-t,t=0..60], [v_06b(t),-t,t=0..60]],color=[violet,pink], thickness=2, axes=box,linestyle=[1,3]):
G2d_062ab[2]:=plot([[t, Force_06a,t=0..60],[t, Force_06b,t=0..60]],color=green, thickness=4, axes=box,linestyle=[1,3]):
display(G2d_062ab, linestyle=1, axes=box);

At last we represent phase curves of oscillations with time:

> G3d_06:= array(1..2):
G3d_06[1]:=spacecurve([v_06a(t),x_06a(t),t],t=0..60,numpoints=1000,color=red):
G3d_06[2]:=spacecurve([v_06b(t),x_06b(t),t],t=0..60,numpoints=1000,color=blue):
display(G3d_06,orientation=[10,80],thickness=2,axes=none);

>

Geometrical properties of phase curves (curvature and torsion) correspond to damping characteristics of oscillations. Blue line corresponds to oscillations with Force_06b.

� 2002. Alexei V. Tikhonenko