Discontinuous Coefficients and Forcing Functions

discont.mws

Douglas B. Meade (meade@math.sc.edu)

1 Overview

Differential equations with piecewise-defined coefficients (including the forcing function) arise in a number of natural circumstances. For example, the air resistance encountered by a skydiver is (dramatically) different between the free-fall and full deployment stages of the jump (see Example 3, pp. 112–113)

We will use graphical solution methods to first see that the solution satisfies the fundamental principles for all differential equations: solutions must follow the direction field. In particular, even though the direction field may not be continuous, the solution can still be continuous (as would be expected in the skydiving example). The issue of finding an explicit formula for the solution is a little different. The key is to break the problem into pieces where the solution can be determined, then to piece these solutions together with appropriate transition conditions.

The discussion follows the general steps found in Problems 31 and 32 from Section 2.3 Nagle and Saff.

2 Objectives

• understand how ODEs with discontinuous coefficients might arise in a physical application

• develop conceptual and mechanical skills for the analysis of ODEs with piecewise-defined coefficients

• continue to develop Maple skills for the manipulation of differential equations

• continue to built a geometric intuition about ODEs and IVPs

3 New Maple Commands

• piecewise creates a piecewise-defined function, i.e., a function defined by different rules for different parts of the domain

4 Problem 31 (Section 2.3, p. 53) — Discontinuous Coefficients

There are times when the coefficient P(x) in a linear equation may not be continuous but may have a jump discontinuity. Fortunately, we may still obtain a “reasonable” solution. For example, consider the initial value problem: y’( x) + P(x) y = x, y(0) = 1 where P(x) = 1 if 0<= x<=2 and P(x) = 3 if x>2.

4.1 Solution

   > restart;
   > with( plots ):
   > with( DEtools ):

Prior to following the book’s outline for the computation of an explicit solution, let’s spend a little time considering exactly what the solution should look like. For 0 <= x <= 2, the problem is simply a standard IVP for a first-order linear equation with constant coefficients. For x>2 a different ODE is involved, but it is still a first-order linear equation with constant coefficients. The general solution to either problem is simple to determine.

The question is, exactly which of these solutions should be selected to form the solution to the IVP?

The initial condition y(0) = 1 identifies the specific solution curve for 0 <= x <= 2. The solution curve followed for x>2 is the one that meets the solution at the end of the first stage. That is, the solution from the first stage is used to determine the initial condition for the second stage. (The precise definition of this will be clearer later.)

But first, let’s develop some intuition by considering the direction field for this equation. The ODE can be represented as

   > ODE := diff( y(x), x ) + P*y(x) = x;

ODE := (\frac{∂}{∂x} y(x)) + P y(x) = x

The coefficient of y(x) can be defined, as a function, using the piecewise command as follows:

> P := piecewise( x<=2, 1, 3 );

P := \{ \begin{matrix}1 & x ≤ 2 \\ 3 & otherwise \end{matrix}

We can check that this definition is correct by evaluating P for a variety of x values:

> seq( eval(subs( x=X, P )), X=[ 0, 1, 2, 3, 4 ] );

1, 1, 1, 3, 3

Another check is to explicitly view the differential equation:

> ODE;

(\frac{∂}{∂x} y(x)) + (\{ \begin{matrix}1 & x ≤ 2 \\ 3 & otherwise \end{matrix} ) y(x) = x

4.1.1 Alternate Approach (using functions)

Note that if you wish to implement P as a function, it would be necessary to replace P with P(x) in ODE and to define the coefficient using P := x -> piecewise( x<=2, 1, 3 );.

The initial condition can be represented as.

> IC := y(0)=1;

IC := y(0) = 1

The direction field for this problem is then created exactly as all previous direction fields have been created:

> DEplot( ODE, y, x=0..5, { [0,1] }, arrows=NONE,
> title=`Problem 31 (Section 2.3) -- Direction Field and Solution`
);

4.1.2 Bug in Release 4

Note that the arrows in the direction field are not displayed in the above plot. This is because DEplot appears to be mishandling the piecewise defined coefficient. (See if you see anything strange when the arrows=NONE option is removed.) This worked flawlessly in Release 3. This problem has been reported to Maple’s Technical Support.

The solution obviously experiences a change at x = 2. This is a direct result of the discontinuity in the direction field at x = 2. Note also how the solution, which follows the direction field at all points, is continuous.

The remainder of this problem is concerned with finding an explicit solution to this problem. This will be done in a series of steps that require the application of our prior knowledge about solving first-order linear ODEs:

4.1.3 (a) Find the general solution for 0 <= x <= 2

For these values of x, P(x) = 1 and the integrating factor is μ(x) = e^{x}. The general solution is:

   > ODE1 := diff( y(x), x ) + y(x) = x;

ODE1 := (\frac{∂}{∂x} y(x)) + y(x) = x

   > GSOLN1 := dsolve( ODE1, y(x) );

GSOLN1 := y(x) = x − 1 + e^{(−x)} C1 4.1.4 (b) Choose the constant in the solution of part (a) so that the

initial condition is satisfied The initial condition is y(0) = 1. The appropriate value of C1 can be found as follows

> CEQN1 := subs( x=0, IC, GSOLN1 );

CEQN1 := 1 = −1 + e^{0} C1

> CSOLN1 := solve( CEQN1, { _C1 } );

CSOLN1 := { C1 = 2}

> SOLN1 := subs( CSOLN1, GSOLN1 );

SOLN1 := y(x) = x − 1 + 2 e^{(−x)}

or, dsolve could be used to find the solution in one clean step

> dsolve( { ODE1, IC }, y(x) );

y(x) = \frac{e^{x} x − e^{x} + 2}{e^{x}}

Well, a second step is required to put the solution in a more familiar form

> expand( " );

y(x) = x − 1 + \frac{2}{e^{x}}

4.1.5 (c) Find the general solution for x>2 The differential equation for these values of x is

> ODE2 := diff( y(x), x ) + 3*y(x) = x;

ODE2 := (\frac{∂}{∂x} y(x)) + 3 y(x) = x

The integrating factor is μ(x) = e^{(3 x)} and the general solution is

> GSOLN2 := dsolve( ODE2, y(x) );

GSOLN2 := y(x) = \frac{1}{3} x − \frac{1}{9} + e^{(−3 x)} C1

4.1.6 (d) choose the constant in the general solution from part (c) so that the solution from part (b) and the solution from part (c) agree at x=2

In other words, find C1 so that

> Limit( rhs(SOLN1), x=2, left ) = Limit( rhs(GSOLN2), x=2, right
);

lim_{x→2−} x − 1 + 2 e^{(−x)} = lim_{x→2+} \frac{1}{3} x − \frac{1}{9} + e^{(−3 x)} C1

Note that this condition forces the continuity of the solution for all x>=0.

Another way of thinking about this is to say that the initial condition for the second stage ( x>2) is obtained from the value of the solution at the end of the first stage. That is,

> IC2 := subs( x=2, SOLN1 );

IC2 := y(2) = 1 + 2 e^{(−2)}

This initial condition can then be used to find the value of C1 that yields the solution which agrees with the solution during the first stage. The compatibility condition at x = 2 is:

> CEQN2 := subs( x=2, IC2, GSOLN2 );

CEQN2 := 1 + 2 e^{(−2)} = \frac{5}{9} + e^{(−6)} C1

Solving for the unknown constant and inserting this value into the general solution proceeds as follows:

> CSOLN2 := solve( CEQN2, { _C1 } );

CSOLN2 := \{ C1 = \frac{2}{9} \frac{2 + 9 e^{(−2)}}{e^{(−6)}}\}

> SOLN2 := simplify( subs( CSOLN2, GSOLN2 ) );

SOLN2 := y(x) = \frac{1}{3} x − \frac{1}{9} + \frac{4}{9} e^{(−3 x+6)} + 2 e^{(−3 x+4)}

A more direct solution to this problem can be obtained by directly solving the appropriate IVP for x>2:

> dsolve( { ODE2, IC2 }, y(x) );

y(x) = \frac{1}{3} x − \frac{1}{9} + \frac{2}{9} \frac{e^{(−3 x)} (2 + 9 e^{(−2)})}{e^{(−6)}}

   > simplify( " );

y(x) = \frac{1}{3} x − \frac{1}{9} + \frac{4}{9} e^{(−3 x+6)} + 2 e^{(−3 x+4)}

Thus, the continuous solution to this problem is the solution formed by piecing together the solutions for 0 <= x <= 2 and for x>2. This is most conveniently accomplished using the piecewise command:

   > SOLN := piecewise( x<=2, rhs(SOLN1), rhs(SOLN2) );

SOLN := \{ \begin{matrix}x − 1 + 2 e^{(−x)} & x ≤ 2 \\ \frac{1}{3} x − \frac{1}{9} + \frac{4}{9} e^{(−3 x+6)} + 2 e^{(−3 x+4)} & otherwise \end{matrix}

4.1.7 (e) Sketch the graph of the solution from x=0 to x=5

There are a number of ways in which the graph can be constructed. One of the most efficient, and straightforward, is to use the previously defined piecewisedefined solution, SOLN. Then the requested plot is obtained with a single plot command

   > plot( SOLN, x=0..5, title=`Problem 31 (Section 2.3) -- Explicit
   Solution` );

Note that this solution agrees with the graphical solution obtained at the beginning of the problem.

4.2 Update for Release 4

One of the new features in Release 4 is the ability to solve ODEs and IVPs with piecewise-defined coefficients directly with dsolve. We illustrate for this example; see ?dsolve,piecewise for more details.

   > IVP := { ODE, IC };

IVP := \{(\frac{∂}{∂x} y(x)) + (\{ \begin{matrix}1 & x ≤ 2 \\ 3 & otherwise \end{matrix} ) y(x) = x, y(0) = 1\}

   > SOLNp := dsolve( IVP, y(x) );

SOLNp := y(x) = \{ \begin{matrix}−1 + 2 e^{(−x)} & x ≤ \frac{4}{3} \\ −1 + 2 e^{(−x)} & x ≤ 2 \\ −\frac{2}{3} x + \frac{4}{9} e^{(−3 x+6)} − \frac{1}{9} + 2 e^{(−3 x+4)} & 2 < x \end{matrix}

4.2.1 Bug In Release 4

It appears as though there are problems with this solution; the first piece is redundant. The cause of the false transition at t = \frac{4}{3} has not been fully investigated but there are other problems with this solution. The output from dsolve is not always the same: sometimes the solution is correct, sometimes the solution is incorrect and sometimes dsolve returns no value. The following plot allows a simple comparison of the solutions. This problem was reported to Maple’s Technical Support on 12 November 1996.

   > plot( rhs(SOLNp), x=0..5 );

5 Problem 32 (Section 2.3, p. 54) — Discontinuous Forcing Terms

There are times when the forcing term Q(x) in a linear equation may not be continuous but may have a jump discontinuity. Fortunately, we may still obtain a reasonable solution imitating the procedure discussed in Problem 31. Use this procedure to find the continuous solution to the initial value problem y’( x) + 2 y = Q(x), y(0) = 0 where Q(x) = 2 if 0 <= x <= 3 and Q(x) = −2 if x > 3.

Sketch the graph of the solution from x = 0 to x = 7.

5.1 Solution

> restart;
> with( plots ):
> with( DEtools ):

One Maple implementation of this problem is

> ODE := diff( y(x), x ) + 2*y(x) = Q(x);

ODE := (\frac{∂}{∂x} y(x)) + 2 y(x) = Q(x)

> IC := y(0) = 0;

IC := y(0) = 0

where the piecewise-defined right-hand side is defined as

> Q := x -> piecewise( x<=3, 2, -2 );

Q := x → piecewise(x ≤ 3, 2, −2)

That is,

> ODE;

(\frac{∂}{∂x} y(x)) + 2 y(x) = \{ \begin{matrix}2 & x ≤ 3 \\ −2 & otherwise \end{matrix}

Note that the RHS has been defined as a function (compare this with the approach used in the previous example).

The direction field and graphical solution to this IVP are produced by one call to DEplot:

> DEplot( ODE, y, x=0..7, { [0,0] }, arrows=NONE,
> title=`Problem 32 (Section 2.3) -- Direction Field & Graphical
Solution` );

The direction field has a similar problem as was noted in the first example.

The low resolution of this graph could be improved by specifying the stepsize= optional argument. (The dirgrid= option is often useful when trying to get a decent picture of the direction field.) Regardless, it seems as though the solution to the first solution wants to be asymptotic to y = 1 initially, but the jump in Q pushes the solutions towards y = −1 for large values of t. The explicit solution should confirm these observations.

The explicit solution is found by following the outline introduced in Problem 31. Both the “long” and the “short” methods will be demonstrated.

5.1.2 Solution by “long approach”

The “long approach” begins by finding the general solution on 0 <= x <= 3:

   > ODE1 := diff( y(x), x ) + 2*y(x) = 2;

ODE1 := (\frac{∂}{∂x} y(x)) + 2 y(x) = 2

Then, find the value of the constant, C1, that satisfies the initial condition

   > GSOLN1 := dsolve( ODE1, y(x) );

GSOLN1 := y(x) = 1 + e^{(−2 x)} C1

   > CEQN1 := subs( x=0, IC, GSOLN1 );

CEQN1 := 0 = 1 + e^{0} C1

> CSOLN1 := solve( CEQN1, { _C1 } );

CSOLN1 := { C1 = −1}

> SOLN1 := subs( CSOLN1, GSOLN1 );

SOLN1 := y(x) = 1 − e^{(−2 x)}

The general solution to the differential equation for x>3 is found as follows:

> ODE2 := diff( y(x), x ) + 2*y(x) = -2;

ODE2 := (\frac{∂}{∂x} y(x)) + 2 y(x) = −2

> GSOLN2 := dsolve( ODE2, y(x) );

GSOLN2 := y(x) = −1 + e^{(−2 x)} C1

To ensure continuity of the final solution, the constant C1 must be chosen so the solution for x>3 agrees with the solution for 0 <= x <= 3 at the boundary point x = 3. Thus,

> IC2 := subs( x=3, SOLN1 );

IC2 := y(3) = 1 − e^{(−6)}

> CEQN2 := subs( x=3, IC2, GSOLN2 );

CEQN2 := 1 − e^{(−6)} = −1 + e^{(−6)} C1

> CSOLN2 := solve( CEQN2, { _C1 } );

CSOLN2 := \{ C1 = −\frac{−2 + e^{(−6)}}{e^{(−6)}}\}

> SOLN2 := simplify( subs( CSOLN2, GSOLN2 ) );

SOLN2 := y(x) = −1 + 2 e^{(−2 x+6)} − e^{(−2 x)}

5.1.3 Solution by “short approach”

The “short approach” is characterized by working with the solution to two initial value problems, rather than the general solutions to an ODE. The first step is to solve the initial value problem for 0 <= x <= 3:

> SOLN1a := dsolve( { ODE1, IC }, y(x) );

SOLN1a := y(x) = 1 − e^{(−2 x)} The initial condition for the second stage is found as before:

> IC2a := subs( x=3, SOLN1a );

IC2a := y(3) = 1 − e^{(−6)}

> SOLN2a := dsolve( { ODE2, IC2 }, y(x) );

SOLN2a := y(x) = −1 − \frac{e^{(−2 x)} (−2 + e^{(−6)})}{e^{(−6)}}

> SOLN2a := simplify( SOLN2a );

SOLN2a := y(x) = −1 + 2 e^{(−2 x+6)} − e^{(−2 x)}

In either case, the complete solution can be expressed using piecewise:

> SOLN := piecewise( x<=3, rhs(SOLN1), rhs(SOLN2) );

SOLN := \{ \begin{matrix}1 − e^{(−2 x)} & x ≤ 3 \\ −1 + 2 e^{(−2 x+6)} − e^{(−2 x)} & otherwise \end{matrix}

> plot( SOLN, x=0..7, title=`Problem 32 (Section 2.3) -- Explicit
Solution` );

Note that this plot has better resolution, and better exhibits the previously identified asymptotic behaviors, than the plot generated with DEplot (the DEplot could be improved by specifying a sufficiently small stepsize).

5.1.4 Update for Release 4

The direct solution from dsolve appears to be correct for this problem

> SOLNp := dsolve( { ODE, IC }, y(x) );

SOLNp := y(x) = \{ \begin{matrix}1 − e^{(−2 x)} & x ≤ 3 \\ −1 + 2 e^{(−2 x+6)} − e^{(−2 x)} & 3 < x \end{matrix}

> plot( rhs(SOLNp), x=0..7 );

6 Copyright

discont.mws – a Maple worksheet illustrating Maple solution techniques for IVPs with discontinuous coefficients and forcing terms (for use with section 2.3 of Nagle and Saff) Prepared by: Douglas B. Meade (7 March 1996) as part of the Maple Supplement for Nagle/Saff (Addison-Wesley) Updated to Maple V, Release 4, by Douglas B. Meade (27 December, 1996) For more information, please contact the author at one of following addresses: Department of Mathematics, Univ. of South Carolina, Columbia, SC 29208 (803) 777-6183 meade@math.sc.edu (copyright, Douglas B. Meade, 1996)