Analytic solution for damped, driven linear harmonic oscillator.
Version: 1.2 (May 15, 2006)Author: David A. Craig <http://web.lemoyne.edu/~craigda/>To do:(i) Add envelope and critical damping curves to x(t) plot.(ii) Animate phase space and 3D curves.(iii) Get initial condition solver working properlyrestart;
with(plots):physical parameters of oscillator
m := 2; # mass
k := 3; # spring constant
b := .3; # linear damping coefficient
w0 := sqrt(k/m); # natural frequency
g := b/(2*m); # damping parameter
wd := sqrt(w0^2-g^2); # damped oscillation frequency
wo := sqrt(g^2-w0^2); # overdamped oscillation 'frequency'
# comment these out if using initial condition finder [see below]
A := 1; # amplitude (underdamped)
ph := 0; # phase (underdamped)
B := A; # amplitude (critical)
C := 0; # amplitude (critical)
F := A; # amplitude (overdamped)
G := 0; # amplitude (overdamped)
transient cases
xU := t -> A*exp(-g*t)*cos(wd*t-ph); # underdamped
xC := t -> exp(-g*t)*(B + C*t); # critically damped
xO := t -> exp(-g*t)*(F*exp(wo*t) + G*exp(-wo*t) ); # overdamped
transient response
xt := t -> piecewise(g < w0, xU(t),g=w0,xC(t), g>w0,xO(t) );
physical parameters of driving force
F0 := 1; # driving force
p0 := 0; # driving phase
w := 1; # driving frequency
wr := sqrt(w0^2-2*g^2); # resonant frequency
beta := w -> (F0/m)/sqrt( (w0^2-w^2)^2+4*g^2*w^2 ); # steady-state amplitude
p := w -> arctan( 2*g*w/(w0^2-w^2) ); # steady-state phase
steady-state response
xSS := t -> beta(w)*cos(w*t-p(w)-ph);
net response of damped, driven linear oscillator
x := t -> xt(t) + xSS(t);
v := t -> diff( x(t), t);
BROKEN!
optional initial condition finder
# takes initial position and velocity, solves for unkown coefficients (A-G, above), substitutes them back into the solution
# requires that specification of the solution coefficients A-G above be commented out and the following code be uncommented
#
# broken in a number of ways:
# ic finds multiple solns in underdamped case e.g. imaginary ones; not sure what's going on yet
# x, v redefinition doesn't plug in found initial conditions to xt or xSS
# would be nice to be able to recursively redefine, but Maple complains if you try [e.g. x := t -> subs( ic, x(t)) ;]
#
# x0 := 0;
# v0 := 0;
# using subs instead of D(x) [say] because D(x) can't seem to handle piecewise
# ic := solve( {x(0)=x0,subs(t=0,v(t))=v0} );
# redefine x, v
# x := 'x';
# v := 'v';
# x := t -> subs( ic, xt(t) + xSS(t)) ;
# v := t -> diff( x(t), t);
plot responses
plot( [xt(t),xSS(t),x(t)],t=0..50, color=[green,red,blue]);
animate response
animate(plot,[x(t),t=0..T],T=0..100);amplitude and phase of steady-state response as functions of driving frequency
plot( beta(x), x=0..20);
plot( p(x), x=0..Pi);phase space
plot( [x(t),v(t),t=0..50]);
# with attractor
plot( [[x(t),v(t),t=0..50], [x(t),v(t),t=90..100]], color=[red,blue]);phase space in time
spacecurve( [x(t),v(t),t], t=0..50,axes=boxed,numpoints=1000);