Unit 2: Infinite Series

Chapter 10: Fourier Series

Sections10.5: periodically driven damped oscillator

Copyright

Copyright * 2001 by Addison Wesley Longman, Inc.

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the publisher. Printed in the United States of America.

Initializations

> restart;

> with(plots):
with(inttrans):

Warning, the name changecoords has been redefined

> interface(showassumed=0);
assume(n,integer);

> alias(Y=laplace(y(t),t,s)):

> PX := proc(h, g)
local f, var, d;
if type(h, procedure) then RETURN(subs(
{'D' = rhs(g) - lhs(g), 'F' = h, 'L' = lhs(g)}, proc(
x::algebraic)
local y;
y := floor((x - L)/D); F(x - y*D)
end))
fi;
if type(g, equation) then var := lhs(g); d := rhs(g) fi;
f := unapply(h, var);
subs({'L' = lhs(d), 'D' = rhs(d) - lhs(d), 'F' = f}, proc(
x::algebraic)
local y;
y := floor((x - L)/D); F(x - y*D)
end)
end:

>

Damped Oscillator Driven Periodically

Consider the driven damped spring-mass system governed by the differential equation

> q := diff(y(t),t,t) + 2*diff(y(t),t) + 10*y(t) = f(t);

>

where the periodic forcing function is the periodic extension of

> g := piecewise(t<=Pi,t, t<=2*Pi,2*Pi-t);

>

The graph of on its domain of [ ] is given by

> plot(g,t=0..2*Pi, labels=[t,`g `], labelfont=[TIMES,ITALIC,12], xtickmarks=6, ytickmarks=3, scaling=constrained);

>

A Maple representation of is obtained with the PX command introduced in Sections10.1, so that

> f := PX(g,t=0..2*Pi):

>

with graph, Figure 10.24, given by

> plot(f(t),t=0..6*Pi, color=black, scaling=constrained, ytickmarks=2, xtickmarks=[6,12,18], labels=[t,`f `],labelfont=[TIMES,ITALIC,12]);

>

Exact Solution by Laplace Transform

A graph of the exact motion of this oscillator, started with the inert initial conditions

, y'(0) = 0

can be found via the Laplace transform. To this end, we obtain F(s), the transform of , with the formula for periodic functions from Sections6.7. Thus, with , the period of . We obtain

=

a calculation executed in Maple via

> F := simplify(int(g*exp(-s*t),t=0..2*Pi)/(1-exp(-2*Pi*s)));

>

Transforming the differential equation and solving for , the Laplace transform of , we get

calculations which we execute in Maple as follows.

The transform of the differential equation then yields

> q1 := laplace(lhs(q),t,s) = F;

>

Our notation has been simplified by aliasing the letter "Y" to the Laplace transform of . Further simplification arises from imposing the inert initial contitions. Thus,

> q2 := subs(y(0)=0, D(y)(0)=0, q1);

>

Solving for the unknown transform Y, we obtain

> q3 := factor(solve(q2,Y));

>

It is impossible to invert this transform with just a finite number of elementary functions, so we expand

in powers of , obtaining

The inverse Laplace transform of the factor

is the function

( ) +

so the solution to the differential equation is

To obtain this solution in Maple, begin by getting a series expansion of in powers of . For this, in , replace with the single letter , obtaining

> if has(q3,exp(-Pi*s)) then q4 := subs(exp(-Pi*s)=z,q3)
else q4 := subs(exp(Pi*s)=1/z,q3);fi;

>

The conditionals in the code for making this replacement are required to guarantee that each time this worksheet is executed, the result is the same. Maple can return the expression for as a fraction containing either or . The replacement of the negative exponential with is crucial for the next step, a Taylor expansion in powers of the negative exponential.

Perform a Taylor expansion about to obtain

> q5 := convert(taylor(q4,z=0,7),polynom);

>

Change back to , getting

> q6 := subs(z=exp(-Pi*s),q5);

>

then invert this truncated series representation of Y(s) to obtain what should be an approximation to . Thus,

> q7 := invlaplace(q6,s,t);

> collect( q7, [seq( Heaviside(t-k*Pi), k=1..6 ) ] );

>

A graph of this version of the solution, namely

is given in Figure 10.25, where, on the interval , the solution is exact because

is zero, since .

> p := plot(q7,t=0..20, color=black, xtickmarks=5, ytickmarks=4,labels=[t,``],labelfont=[TIMES,ITALIC,12], view=[0..20,0.. .33]):
p;

>

The first term dropped from the Taylor expansion of Y(s) was . If that term had been retained, and inverted, it would have contributed a multiple of , and this term is zero for . Hence, on the interval [0,20], we have kept enough terms for the graph to represent exactly.

>

Approximate Solution by Fourier Series

We next obtain a sequence of approximate solutions which rapidly converge to the exact solution. In the differential equation, the driving function is replaced by , a partial sum of the Fourier series of , and the corresponding solution computed. The sequence of approximate solutions converges rapidly to the exact solution.

Using the general formalism of Sections10.1, we obtain a Fourier series representation of . Hence, with the interval taken as [ ] so

> L := Pi;

>

we obtain

by writing the integrals

> q0 := Int(g,t=0..2*L)/L;
qa := Int(g*cos(n*Pi*t/L),t=0..2*L)/L;
qb := Int(g*sin(n*Pi*t/L),t=0..2*L)/L;

>

and evaluating them to get the Fourier coefficients

> a0 := value(q0);
a := simplify(value(qa));
b := value(qb);

>

We note that for all , so the sine terms disappear from the series. Also, examining the via the list

> seq(subs(n=k,a),k=1..10);

>

we see that just the cosine terms with an odd index "survive." Hence, the first three distinct partial sums of the resulting Fourier series are

> for k from 1 by 2 to 5 do
f||((k+1)/2) := a0/2 + sum(a*cos(n*Pi*t/L),n=1..k);
od;

>

Letting , be the solutions of the original IVP when is replaced by , we obtain Figure 10.26 wherein the approximations are plotted against .

To obtain the solutions in Maple, and to obtain the graphs shown in Figure 10.26, we execute the following loop which uses Maple's dsolve command to solve three IVPs whose right-hand sides are the partial sums . Within the same loop, the graphs of the solutions are generated, and displayed against the graph of the exact solution in Figure 10.25.

Also, we re-generate figure 10.25 with a thickened black line, and use it in Figure 10.26 to enhance visibility.

> p := plot(q7,t=0..20, color=black, thickness=3):
for k from 1 to 3 do
Q||k := rhs(dsolve({lhs(q) = f||k, y(0)=0, D(y)(0)=0}, y(t),method=laplace));
p||k := plot(Q||k,t=0..20, color=red, ytickmarks=2);
d||k := display([p,p||k],xtickmarks=5, ytickmarks=4, labels=[t,``], labelfont=[TIMES,ITALIC,12], view=[0..20,0.. .33]);
od:

>

Figure 10.26 is then

> display(array([d||(1..3)]),scaling=constrained);

>

The third partial sum approximates well enough for the graph of the resulting approximation to to be indistinguishable from the graph of the exact solution. In Chapter 14 we will study additional techniques for obtaining approximate analytic solutions of differential equations.

We close by separately displaying the three graphs comprising Figure 10.26. The thick black curve is the exact solution, and the approximate solutions are in red.

> d1;

> d2;

> d3;

>