ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 34 -- Piecewise-Defined and Impulse Functions

Outline of Unit 34

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

> with( inttrans ):

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

Warning, the name hilbert has been redefined

> #assume( a>0 ):

> Alist := [ 0, 1, 2, 4, 5, 10, 15, a ]:

> alias( seq( u[A](t) = Heaviside(t-A), A=Alist ) );

>

A piecewise-defined function is a function that cannot be defined by one elementary formula for all points in the domain. For example, the piecewise linear function that linearly interpolates the origin, ( 1, 1 ), and ( 2, 0 ) and is zero for all t > 2 can be defined in Maple using the piecewise command :

> f := piecewise( t<1, t, t<2, 2-t, t>=2, 0 );

>

To verify the correctness of this definition, consider the plot of the function :

> plot( f, t=0..3, thickness=3, scaling=constrained );

>

If the function is discontinuous, e.g. ,

> g := piecewise( t<1, t, t<2, t-1, 0 );

> plot( g, t=0..3, thickness=3, scaling=constrained );

>

The vertical line segments in this plot are not part of the graph of this function; they arise when Maple connects the sample of points on the graph of the function.

> plot( g, t=0..3, thickness=3, scaling=constrained, style=point );

>

The discont=true option instructs Maple to attempt to identify discontinuities in the plotting interval and to not connect points on either side of the discontinuity.

> plot( g, t=0..3, thickness=3, scaling=constrained, discont=true );

>

The Heaviside function is the piecewise-defined function with definition

> convert( Heaviside(t), piecewise, t );

>

More generally, the Heaviside function at is

> `u`[a] = convert( Heaviside(t-a), piecewise, t );

>

A piecewise-defined function can be written in terms of Heaviside functions with appropriate shifts. For example,

> f = collect( convert( f, Heaviside ), Heaviside );

> g = collect( convert( g, Heaviside ), Heaviside );

>

Let > 0. The Laplace transform of the Heaviside function at can be obtained by a simple change of variables ( ) in the integral definition.

= = = = = .

The result can also be obtained directly from Maple

> assume( a>0 );

> L( 'u'[a](t) ) = laplace( u[a](t), t, s );

>

Consider the IVP

> ode1 := 2*diff( x(t), t ) + x(t) = h(a,t);

> ic1 := x(0) = 5;

>

where the forcing function is initially zero, jumps to the value 1 at , remains at 1 until , decreases linearly to zero between and , then stays at zero for > 15 :

> h := unapply( u[a](t) + (10-t)/5*u[10](t) + (t-15)/5*u[15](t), (a,t) );

>

For example,

> convert( h(5,t), piecewise, t );

>

In advance of finding an explicit solution to this IVP, take a moment to think about the direction field for this ODE. When , the direction field is

> dir_field1 := DEplot( eval(ode1,a=5), x(t), t=0..20, x=0..6, arrows=THIN ):

> dir_field1;

>

The solution curve satisfying the initial condition must follow the direction field. The fact that there are drastic changes in the directions at several times does not change this basic fact. In addition, note that the solution will be continuous even when the directions are discontinuous.

The solution to the general IVP is

> sol1 := dsolve( { ode1, ic1 }, x(t), method=laplace );

>

Or, as a piecewise-defined function when :

> convert( rhs(eval(sol1,a=5)), piecewise, t );

>

A graph of the solution, with , superimposed on the direction field shows how this solution follows the direction field

> sol_plot1 := plot( eval(rhs(sol1),a=5), t=0..20, color=BLUE ):

> display( [dir_field1,sol_plot1] );

>

An animation of the solution and direction fields for integer values of between 0 and 10 shows how the solutions vary with the parameter :

> P := A -> display( [ DEplot( eval(ode1,a=A), x(t), t=0..20,

> x=0..6, arrows=THIN ),

> plot( eval(rhs(sol1),a=A), t=0..20,

> color=BLUE ) ] ):

> display( [ 'P'(i) $ i=0..10 ], insequence=true );

>

The Dirac function is defined by the seemingly inconsistent pair of properties : for all and . Because the integral of any function that is zero almost everywhere must be zero, is not really a function (it is a distribution). Because of this fact, it can be a little difficult to work with the Dirac function. Most properties of the Dirac function are obtained using the fact that

where

> h := Delta -> 1/(2*Delta) * ( Heaviside(t+Delta) - Heaviside(t-Delta) );

>

For all > 0 observe that has the property that :

> Int( h[Delta], t=-infinity..infinity )

> = int( h(Delta), t=-infinity..infinity );

>

and, for all :

> Limit( h[Delta], Delta=0, right ) = limit( h(Delta), Delta=0, right );

>

The following animation provides additional information about this sequence of functions and its limit.

> display( [seq( plot(h(Delta), t=-1.1..1.1), Delta=[2^(-i) $ i=0..8] )], insequence=true );

>

Superimposing each frame of this animation in a single plot is also useful.

> plot( [seq( h(Delta), Delta=[2^(-i) $ i=0..8] )], t=-1.1..1.1 );

>

Note that it is not possible to plot the Dirac function.

More generally, the Dirac function at time is

= =

where

> h := (a,Delta) -> 1/(2*Delta) * ( Heaviside(t-a+Delta) - Heaviside(t-a-Delta) );

>

The Laplace transform of the Dirac function at is derived from the Laplace transforms of the sequence of approximating functions:

That is,

> assume( Delta>0, Delta<a );

> q1 := laplace( h(a,Delta), t, s );

> L( delta[a] ) = Limit( q1, Delta=0, right );

> `` = limit( q1, Delta=0, right );

>

The Maple name for the Dirac function (at time 0) is Dirac . Thus, the previous result could have been obtained directly with

> L( delta[a] ) = laplace( Dirac(t-a), t, s );

>

A common physical interpretation of the Dirac function at time is as an instantaneous impulse when . To illustrate, consider an undamped harmonic oscillator with unit mass ( ) and spring constant ( ) that is subject to an instantaneous impulse at and a second impulse, in the opposite direction and with three times the force, at . Assuming the system starts at rest the IVP is

> ode2 := diff( x(t), t$2 ) + 3*x(t) = Dirac(t-1) - 3*Dirac(t-4);

> ic2 := x(0)=0, D(x)(0)=0;

>

The solution to this problem is

> infolevel[dsolve] := 3:

> sol2 := dsolve( {ode2,ic2}, x(t), method=laplace );

> infolevel[dsolve] := 0:

dsolve/inttrans/solveit:   Transform of eqns is   {s^2*laplace(x(t),t,s)+3*laplace(x(t),t,s)-exp(-s)+3*exp(-4*s)}

dsolve/inttrans/solveit:   Algebraic Solution is   {laplace(x(t),t,s) = -(-exp(-s)+3*exp(-4*s))/(s^2+3)}

>

Note the appearance of Heaviside functions at and in the solution. The graph of the solution is

> plot( rhs(sol2), t=0..10 );

>

The solution is continuous for all time and that the amplitude of the oscillation depends on the magnitude of the impulse.

>