ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 29 -- Application: Nonlinear Springs

Outline of Unit 29

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

> with( student ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

>

Consider the spring-mass system in which a weight is attached to one end of a massless spring and the other end is mounted to the ceiling. When the weight is pulled downwards a short distance and released (from rest), the spring-mass system oscillates up and down. The motion is completely determined by Newton's Second Law of Motion: , where is the total force on the system and is the displacement of the mass from the spring's unstretched position, i.e. , < 0 corresponds to compressing the spring and > 0 to an extension of the spring. Note that the displacement must remain with certain bounds to remain physically realistic. If the displacement becomes too large, the spring will break in two. The limit on the negative displacement occurs when the coils of the spring bunch together.

Models of this system differ in the way the forces included. For this discussion, the gravitational force, frictional force, and spring force are included. The gravitational force is . The frictional force is assumed to be proportional to velocity and acts in the direction opposite to the direction of motion. Thus, the second-order ODE for the displacement of the spring-mass system is

> F[gravity] := - m * g;

> F[friction] := - c * diff( y(t), t );

> F[spring] := S( y(t) );

> ode1 := m * diff( y(t), t$2 ) = F[spring] + F[friction] + F[gravity];

>

The corresponding first-order system is

> U := [ y(t), v(t) ]:

> f := (y,v) -> v:

> g := (y,v) -> ( S(y) - c*v ) / m - g:

> sys1 := equate( diff( U, t ), [ f( y(t), v(t) ), g( y(t), v(t) ) ] );

>

The general condition for an equilibrium solution is

> q1 := solve( eval( sys1, { y(t)=y0, v(t)=v0 } ), { y0, v0 } );

>

The most note-worthy property of this result is that the equilibrium solutions are independent of the damping force.

>

Hooke's Law assumes the spring force is linear:

> SS[linear] := y -> -k*y;

>

with > 0 (and >0, > 0, and > 0).

There is only one equilibrium solution to the linear system

> S := SS[linear]:

> equil[linear] := q1;

> for X0 in [equil[linear]] do

> b := eval( [ f(y0,v0), g(y0,v0) ], X0 );

> A := eval( evalm(jacobian( [f(y0,v0),g(y0,v0)], [y0,v0] )), X0 );

> linsys1 := equate( diff( U, t ), evalm( A &* U + b ) );

> print( `-------------------------------------------------------------` );

> print( `Linearization at [y0,v0]`=eval([y0,v0],X0) );

> print( A = evalm( A ) );

> print( lambda = [eigenvals( A )] );

> print( `-------------------------------------------------------------` );

> end do:

>

It is not difficult to see that there are three cases to consider for these eigenvalues:

Case 1: 0 < <

Case 2:

Case 3: >

Because of the complexities of the general solution to a cubic equation, this discussion is restricted to a specific example.

> param1 := { m=1, g=98/10, c=2/10, k=10, alpha=2/10 }:

> S := SS[soft]:

> sys[soft] := eval( sys1, param1 );

>

> equil[soft] := allvalues( eval( q1, param1 ) );

> evalf( equil[soft] );

>

Thus, there are three equilibria. Linearization is used to classify the equilibria

> for X0 in [equil[soft]] do

> b := eval( [ f(y0,v0), g(y0,v0) ], X0 union param1 );

> A := eval( evalm(jacobian( [f(y0,v0),g(y0,v0)], [y0,v0] )), X0 union param1 );

> linsys1 := equate( diff( U, t ), evalm( A &* U + b ) );

> print( `-------------------------------------------------------------` );

> print( `Linearization at [y0,v0]`=eval([y0,v0],X0 union param1) );

> print( A = evalm( A ) );

> print( lambda = [eigenvals( A )] );

> print( `-------------------------------------------------------------` );

> end do:

>

Based on the eigenvalues, ( -1, 0 ) is a stable spiral and ( , 0 ) and ( , 0 ) are saddle points.

To conclude the discussion of soft springs construct a plot containing the equilibrium solutions, representative solution curves, and the direction field for this example.

> Pequil := pointplot( [seq( eval( [y0,v0], E ), E=[equil[soft]] )],

> color=BLUE, symbolsize=18 ):

> icA := [ [ y(0)=-5, v(0)=-4 ], [ y(0)=-8, v(0)=8] ];

> PsolA := DEplot( sys[soft], [ y(t), v(t) ], t=-0.1..10, icA,

> stepsize=0.05, arrows=NONE ):

> icB := [ [ y(0)=-6, v(0)=-4], [ y(0)=-9, v(0)=8] ];

> PsolB := DEplot( sys[soft], [ y(t), v(t) ], t=-0.2..0.5, icB,

> stepsize=0.05, arrows=NONE ):

> icC := [ [ y(0)=-10, v(0)=25] ];

> PsolC := DEplot( sys[soft], [ y(t), v(t) ], t=-0.1..2.5, icC,

> stepsize=0.1, arrows=NONE ):

> icD := [ [ y(0)=-10, v(0)=30] ];

> PsolD := DEplot( sys[soft], [ y(t), v(t) ], t=-0.1..1.1, icD,

> stepsize=0.1, arrows=NONE ):

> icE := [ [ y(0)=10, v(0)=-10] ];

> PsolE := DEplot( sys[soft], [ y(t), v(t) ], t=-0.1..0.4, icE,

> stepsize=0.05, arrows=NONE ):

> icF := [ [ y(0)=10, v(0)=-15] ];

> PsolF := DEplot( sys[soft], [ y(t), v(t) ], t=-0.2..1.8, icF,

> stepsize=0.1, arrows=NONE ):

> Pdfield := DEplot( sys[soft], [ y(t), v(t) ], t=0..1,

> y=-15..15, v=-40..40, dirgrid=[40,40] ):

> display( Pequil, PsolA, PsolB, PsolC, PsolD, PsolE, PsolF, Pdfield, view=[-15..15, -40..40 ] );

>

Note that the solution curves that become unbounded cease to be physically realistic when they exceed the spring's ability to compress or extend.

>

Repeating the analysis for a hard spring, the equilibrium solutions are

> S := SS[hard]:

> sys[hard] := eval( sys1, param1 );

> q2 := evalf( allvalues( eval( q1, param1 ) ) );

>

Note that only one of these solutions is real; the complex-valued solutions are not physical.

> equil[hard] := op(remove( has, [q2], I ));

>

Linearization about this equilibrium yields

> for X0 in [equil[hard]] do

> b := eval( [ f(y0,v0), g(y0,v0) ], X0 union param1 );

> A := eval( evalm(jacobian( [f(y0,v0),g(y0,v0)], [y0,v0] )), X0 union param1 );

> linsys1 := equate( diff( U, t ), evalm( A &* U + b ) );

> print( `-------------------------------------------------------------` );

> print( `Linearization at [y0,v0]`=eval([y0,v0],X0 union param1) );

> print( A = evalm( A ) );

> print( lambda = [eigenvals( A )] );

> print( `-------------------------------------------------------------` );

> end do:

>

Therefore, the sole equilibrium is a stable spiral.

To conclude the discussion of hard springs construct a plot containing the equilibrium solution, representative solution curves, and the direction field for this example.

> Pequil := pointplot( [seq( eval( [y0,v0], E ), E=[equil[hard]] )],

> color=BLUE, symbolsize=18 ):

> icA := [ [ y(0)=3*j, v(0)=0 ] $ j=-3..0 ];

> PsolA := DEplot( sys[hard], [ y(t), v(t) ], t=0..10, icA,

> stepsize=0.025, arrows=NONE ):

> Pdfield := DEplot( sys[hard], [ y(t), v(t) ], t=0..1,

> y=-15..15, v=-40..40, dirgrid=[40,40] ):

> display( Pequil, PsolA, Pdfield, view=[-15..15, -40..40 ] );

>

Now, all solutions appear to remain bounded and ultimately are attracted to the equilibrium solution. This illustrates a key qualitative difference between soft and hard springs.

>