Unit 3: Ordinary Differential Equations - Part 2

Chapter 14: Series Solutions

Sections14.5: the nonlinear spring and Lindstedt's method

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(DEtools):

Warning, the name changecoords has been redefined

> interface(showassumed=0);

>

A Nonlinear Oscillator

The differential equation

that is,

> q := diff(x(t),t,t) + x(t) + a*x(t)^3 = 0;

>

governs an undamped oscillator whose spring obeys a nonlinear restoration law, . For , the spring is called "hard," while for , it is called "soft." The magnitude of the restorative force for the hard spring is greater than that of the linear spring (whose restorative force is just ). Likewise, the magnitude of the restorative force for the soft spring is less than that of the linear spring for small enough displacements.

When the system is just a simple harmonic oscillator which clearly supports periodic solutions. We wish to show that for the system also exhibits periodic solutions. We show this by obtaining a first integral (a function constant on the trajectories in the phase plane) and showing that its level curves, the trajectories, are closed curves. To this end, multiply the differential equation by x'( t ), and integrate with respect to , obtaining

> q2 := map(Int,expand(q*diff(x(t),t)),t);

>

The antiderivative on the right is a constant we will call E , and each term on the left is an exact differential. Hence, we obtain

> q3 := lhs(value(q2)) = E;

>

For each value of , the level curves of define solutions of the differential equation. Indeed, differentiating E with respect to yields

> q4 := collect(diff(q3,t),diff);

>

which is the original differential equation multiplied by .

Hence, a solution of the differential equation is a level curve of and a level curve of defines a solution of the nonlinear differential equation.

It is easier to plot the level curves of if we change variables according to

u = x( t )

v = x'( t )

so that

u ' = x' = v

v ' = x'' =

which then gives in the form

+

that is,

> E := subs(diff(x(t),t)=v,x(t)=u,lhs(q3));

>

For each value of , the level curves of E( u , v ) = constant define solutions of the differential equation. Some notion of how the level curves (which are phase plane trajectories) vary with is obtained from the following animation in which three contours of E are shown while varies from 0.2 to 2.

> f := z -> contourplot(subs(a=z,E),u=-1..1,v=-1..1, contours=[.1,.2,.3], color=black):
display([seq(f(k/5),k=1..10)],insequence=true, scaling=constrained);

>

We rely on the graphs in the animation to support the claim that the trajectories for this differential equation are closed curves. Hence, for the undamped oscillator with the nonlinear spring supports periodic solutions.

>

Dependence of Frequency on Initial Displacement

In the nonlinear case, the (angular) frequency depends on the initial displacement. If we let that intial displacement be first 1, then 2, (along with an initial velocity x'(0) = 0), and obtain graphs of the numeric solutions when , we obtain Figure 14.10, the following graph, which shows two solutions corresponding to .

> p1:=DEplot(subs(a=1,q),x(t),t=0..10,[[x(0)=1,D(x)(0)=0], [x(0)=2,D(x)(0)=0]], stepsize=.1,linecolor=black):
p2:=subs(THICKNESS(3)=THICKNESS(1),p1):
display(p2, scaling=constrained, xtickmarks=5, ytickmarks=5, labels=[`t `,x],labelfont=[TIMES,ITALIC,12]);

>

The graph clearly indicates different periods, and hence frequencies, for the two solutions which start with different initial displacements.

This is in distinction to the behavior of the linear oscillator where, for , the angular frequency is

= 1

a constant depending on just the spring constant and the mass . Indeed, for we have

that is,

> dsolve(subs(a=0,q),x(t));

>

and the angular frequency is the multiplier of in the trig functions in the solution.

In the nonlinear case, the dependence of the angular frequency on the initial displacement can be computed from the first integral . First, evaluate the constant E from the initial conditions = , obtaining

> E0 := subs(v=0,u=A,E);

>

Next, solve the equation for v = , obtaining the two solutions

= +

that is,

> V := solve(E=E0,v);

>

If we set = and separate variables, we get

On a closed trajectory in the phase plane, motion is clockwise. At , the differential equations

u ' = x' = v

v ' = x'' =

give

v ' =

so motion at this point is downward, and hence, clockwise around the trajectory. Thus, we pick the negative square root for , and integrate to , one-quarter of the way around the trajectory. Thus, the total time for one orbit is

> T := Int(4/V[2],u=A..0);

>

To third-order, a series expansion of in powers of is then

obtained in Maple via

> assume(A>0):
omega = series(2*Pi/T,a,4);

>

Thus, for the nonlinear oscillator, the angular frequency depends on the initial displacement, and we have obtained a series expansion of this dependence. Shortly, we will compute this same expression by Lindstedt's perturbation technique.

>

Regular Perturbation

As we saw in Sections14.4, a regular perturbation will lead to an expansion containing secular terms which are not periodic. Hence, the regular perturbation series of Poincare cannot approximate the periodic solutions of the nonlinear oscillator. We therefore fail to learn about the depencence of the frequency on the initial displacement.

However, for the sake of completeness, we will develop the regular perturbation. The purpose is to justify the claim that it contains secular terms which prevent it from approximating the periodic solutions of the nonlinear oscillator.

Assume a regular perturbation series of the form

> X := add(x[k](t)*a^k,k=0..2);

>

and substitute this into the differential equation, obtaining ,

that is,

> q5 := collect(simplify(subs(x(t)=X,q)),a);

>

where the coefficients are given by

> for k from 0 to 2 do
A[k] = map(coeff,q5,a,k);
od;

>

The identical vanishing of the gives differential equations for the . The first equation

depends only on , so can be solved immediately. For simplicity, we choose the initial point as , thereby obtaining

Each successive differential equation contains terms depending on the solution of the previous equations. We solve this sequence of equations, taking the initial point as for each successive approximation.

Clearly, this is a job for Maple. To facilitate our work in Maple, we replace the symbols with via the loop

> for k from 0 to 2 do
Q||k := map(coeff,q5,a,k);
od;

>

We then obtain , the solution of the first equation and initial condition , as

> s0 := dsolve({Q0,x[0](0)=1,D(x[0])(0)=0},x[0](t));

>

Each successive differential equation contains terms depending on the solutions of previous equations. We solve this sequence of equations, taking the initial point as (0,0) for each successive approximation. We thereby obtain

> for k from 1 to 2 do
s||k := combine(dsolve({subs(seq(s||n,n=0..k-1),Q||k),x[k](0)=0,D(x[k])(0)=0}, x[k](t)),trig);
od;

>

The solutions for both and contain secular terms, as predicted. The regular perturbation series is not useful for studying the dependence of the frequency on the initial displacement. We therefore consider Lindstedt's method for generating a more useful asymptotic expansion.

>

Lindstedt's Method

In Lindstedt's method, the angular frequency is made to depend on by assuming

= 1 +

for which a truncated version is

> W := 1 + add(omega[k]*a^k,k=1..3);

>

We then carry out a perturbation solution of the differential equation by writing

Since the angular frequence is related to the period by

the two expansions in the Lindstedt method are merged by the change of variables . In terms of , the solution assumes the form

=

and the differential equation becomes

where now, primes denote differentiation with respect to .

The change of variables rests on the chain rule which gives

A second application of the chain rule then gives the multiplying U''( ).

In fact, with the differential equation given by

> q;

>

the change of variables is effected in Maple by

> q6 := Dchangevar({t=tau/omega,x(t)=U(tau)},q,t,tau);

>

We next write the first few terms in the expansion for U( ).

> UU := add(u[k](tau)*a^k,k=0..3);

>

Then, substitute into the differential equation

> q6;

>

the expansions for both and , collect coefficients of like powers of , and thus obtain

> q7 := collect(simplify(subs(omega=W,U(tau)=UU,q6)),a);

>

Putting this into the form , we get for the expresions

> for k from 0 to 3 do
sigma[k] = map(coeff,q7,a,k);
od;

>

The identical vanishing of the give differential equations for the . To facilitate this calculation in Maple, we rewrite the equations as via the following loop.

> for k from 0 to 3 do
Q||k := map(coeff,q7,a,k);
od;

>

The first equation contains only , so this equation, and the initial conditions , '(0) = 0, lead to the solution

> s0 := dsolve({Q0,u[0](0)=A,D(u[0])(0)=0},u[0](tau));

>

The equation for contains terms in . Incorporating this solution and the initial conditions , '(0) = 0, we obtain

> s1 := combine(dsolve({simplify(subs(s0,Q1)),u[1](0)=0,D(u[1])(0)=0}, u[1](tau)),trig);

>

As with the regular perturbation, we have secular terms appearing. However, the coefficient is as yet undetermined, so we choose a value for it which will guarantee the vanishing of the secular term. Thus, collect all secular terms with

> q8 := collect(s1,[tau,sin]);

>

and set the coefficient of the secular term equal to zero.

> q9 := map(coeff,q8,sin(tau));

>

Solving for the value of which eliminates the secular term from , we find

> w[1] := solve(q9,omega[1]);

>

so that becomes

> subs(omega[1]=w[1],s1);

>

The equation for contains terms involving and . Thus, making the appropriate substitutions, we solve for the function satisfying the initial conditions '(0) = 0, obtaining

> s2 := combine(dsolve({simplify(subs(s0,s1,omega[1]=w[1],Q2)), u[2](0)=0, D(u[2])(0)=0}, u[2](tau)),trig);

>

Again, there are secular terms which we collect via

> q10 := collect(s2,[tau,sin]);

>

We then determine by the requirement that there be no secular terms in . Hence, form the equation

> q11 := map(coeff,q10,sin(tau));

>

from which we determine

> w[2] := solve(q11,omega[2]);

>

The function thus becomes

> subs(omega[2]=w[2],s2);

>

The equation for contains terms in and . Making the appropriate substitutions, and solving for the which satisfies the initial conditions '(0) = 0, we obtain

> s3 := combine(dsolve({simplify(subs(s0,s1,s2,omega[1]=w[1], omega[2]=w[2], Q3)), u[3](0)=0,D(u[3])(0)=0}, u[3](tau)),trig);

>

Once again, there are secular terms which we gather together with

> q12 := collect(subs(omega[1]=w[1],omega[2]=w[2],s3),[tau,sin]);

>

To choose the which eliminates the secular term from , form the equation

> q13 := map(coeff,q12,sin(tau));

>

and solve, obtaining

> w[3] := solve(q13,omega[3]);

>

Consequently, the solution for is

> subs(omega[3]=w[3],s3);

>

and, the first few terms of the expansion for are

> omega = 1 + add(w[k]*a^k,k=1..3);

>

This compares favorably to the expansion

> series(2*Pi/T,a,4);

>

found in an earlier section.

The Lindstedt solution is then given by

> U(tau) = subs(s||(0..3),seq(omega[k]=w[k],k=1..3),UU);

>

where = , and is as above. To judge the viability of the Lindstedt expansion, write the approximation as a function of , obtaining, for

> q14 := subs(s||(0..3),seq(omega[k]=w[k],k=1..3),
tau=t*(1+add(w[k]*a^k,k=1..3)),UU);

>

Since our test of the utility of the expansion will be a graph, assign the values and to the initial displacement and spring parameter, respectively. Thus, the approximate solution for is

> q15 := subs(A=2,a=1/10,q14);

>

To establish as a numeric reference solution against which to test Lindstedt's solution, execute

> F := dsolve({subs(a=1/10,q),x(0)=2,D(x)(0)=0},x(t),numeric):
xx := z -> subs(F(z),x(t)):

>

and then draw Figure 14.11, the graph of the Lindstedt solution (the dotted red curve), and a numeric solution (the solid black curve), for .

> plot(['xx'(t),q15],t=100..110,color=[black,red], linestyle=[1,2], xtickmarks=6, ytickmarks=5, labels=[t,`x `], labelfont=[TIMES,ITALIC,12]);

>

In the neighborhood of , the two solutions are still in agreement. Figure 14.12, the next plot, shows these two solutions for . In the neighborhood of , the two solutions start separating.

> plot(['xx'(t),q15],t=150..160,color=[black,red], linestyle=[1,2], xtickmarks=6, ytickmarks=5, labels=[t,`x `], labelfont=[TIMES,ITALIC,12]);

>

However, with and , the Lindstedt expansion gives remarkable accuracy.

>