Revitalizing Differential Equations with a Coordinated Use of Numerical, Graphical, and

Analytical Methods

Douglas B. Meade

Department of Mathematics University of South Carolina Columbia, SC 29208

Phone:  (803) 777-6183

E-mail:  meade@math.sc.edu
WWW URL: http://www.math.sc.edu/˜ meade/

18 November 1995

November 1995  1995 ICTCM – Houston  0

Revitalizing Differential Equations with a Coordinated Use of Numerical, Graphical, and Analytical Methods

Douglas B. Meade

Department of Mathematics University of South Carolina Columbia, SC 29208

Phone:  (803) 777-6183

E-mail:  meade@math.sc.edu
WWW URL: http://www.math.sc.edu/˜ meade/

18 November 1995

November 1995  1995 ICTCM – Houston  1

Outline

• Background Information

• Sampling of teaching techniques

• Examples – Direction Fields

– Picard Iterates

– Piecewise Defined Data, Systems, and Control

November 1995  1995 ICTCM – Houston  2

Background Information

Current Status • Assistant Professor of Mathematics

• Lilly Teaching Fellow (1994–1995)

• Maple Ambassador (1995–present)

• Rose-Hulman CAS Workshop (1994 and 1995)

Background Information

CAS Experience

• graduate student (Carnegie Mellon, 1984–1989)

• research (Purdue, 1989–1991; South Carolina, 1991–present)

• teaching (South Carolina, 1994–present)

– Calculus II

– Introduction to Differential Equations

– Nonlinear Optimization

– Applied Mathematics (graduate)

Including Symbolic Algebra in the Classroom

• Out-of-Class Labs

• In-Class Usage – Discuss

– Demonstrate

– Do

November 1995  1995 ICTCM – Houston  5

Discuss

There are computer packages that can assist with the calculations we are discussing. Please see the handout for a sample session.

• requires no special classroom arrangements

• fits easily into traditional lecture format

• very safe – no surprises or hardware problems

• passive participation – second-hand account of session

• minimal transfer of skill/knowledge

November 1995  1995 ICTCM – Houston  6

Demonstrate

I am now starting the computer algebra program called Maple. With this program you will be able to see how the computer can be used to help with the calculations involved in this problem. A complete transcript of this demonstration is provided on the handout.

• requires minimal, but non-trivial, hardware

• transition between lecture and demonstration takes time

• passive participation – exposure to CAS is only as a spectator • minimal transfer of skill/knowledge

• Murphy’s Law

November 1995  1995 ICTCM – Houston  7

Do

In this course we will make frequent – almost daily – use of a computer algebra system. The power of the CAS comes from its ability to manipulate mathematical expressions – not just numbers.

• special hardware needed in classroom

• restructuring of presentation is required

• active participation

• students acquire a skill for future use

• (eMurphy’s Law)^{n} where n = number of students

Differential Equations and Computer Algebra

• Objectives – Develop an understanding of what a differential equation is.

– See what it means to “solve” an initial value problem.

– First exposure to existence theory and qualitative analysis.

Differential Equations and Computer Algebra

• Techniques

– Direction Fields

(http://www.math.scarolina.edu/~meade/ictcm/dfield.ms)

– Picard Iterates

(http://www.math.scarolina.edu/~meade/ictcm/picard.ms)

– Applications – in ICTCM VIII Proceedings

(http://www.math.scarolina.edu/~meade/ictcm/parachute.ms)

November 1995  1995 ICTCM – Houston  10

Maple Example of Direction Fields

The DEtools package provides a number of useful routines for the graphical analysis of ODEs. The primary procedure is DEplot.

 > restart;
 > with( DEtools ):

Define the differential equation

 > ODE := diff( y(t), t ) = t^2 - y(t);

ODE := \frac{∂}{∂t} y( t) = t^{2} − y( t)

November 1995  1995 ICTCM – Houston  11

Maple Example of Direction Fields

The first point is that a (first-order) differential equation provides the slope of “the solution”.

 > DEplot( ODE, [t,y], t=-3..4, y=-3..10,
 >  arrows=THIN, title=`y'=t^2-y` );

November 1995  1995 ICTCM – Houston  12

Maple Example of Direction Fields

Next, it is nice to be able to talk about solutions of an ODE as a curve that follow the flow created by a direction field.

 > DEplot( ODE, [t,y], t=-3..4, { [0,1] }, y=-3..10,
 >  arrows=THIN, title=`y'=t^2-y` );

November 1995  1995 ICTCM – Houston  13

Maple Example of Direction Fields

Also, under appropriate conditions, there is exactly one solution through every point in the (t, y) plane.

> DEplot( ODE, [t,y], t=-3..4,
>  { [0,-2], [0,-1], [0,0], [0,1], [0,2] }, y=-3..10,
>  arrows=THIN, title=`y'=t^2-y` );

Maple Example of Direction Fields

• Discussion

– The third plot presents a natural discussion of qualitative behavior of solutions.

– Additional information is needed to confirm the conjectures formulated from the graphical information.

– Note that simple modifications of a single command, DEplot, provide all three views of this problem. (The seq and/or $ commands can shorten the definition of a set of initial conditions.)

– For additional information on the Maple commands, see

?DEtools and ?DEplot.

Picard Iterates and Maple

• Picard’s Method for y′(t) = f(t, y(t)), y(t_{0}) = y0 – Given ϕ_{0}

– Define, for all k = 0, 1, 2, …,

ϕ_{k+1}(t) = y_{0} + ∫_{t_{0}}^{t} f(x, ϕ_{k}(x)) dx.

• Maple implementation

> restart;
> with( DEtools ):
> with( plots ):
> picard := proc(f, init, phi)
> local t0, y0;
> t0 := init[1];
> y0 := init[2];
> unapply( y0 + int( f(x,phi(x)), x=t0..t ), t );
> end:

November 1995  1995 ICTCM – Houston  16

Picard Iterates: Example 1

Consider the initial value problem: \frac{dy}{dy} = y, y(0) = 1. A first step is to visualize the solution to this problem

 > DEplot( diff(y(t),t)=y(t), [t,y], t=0..2, { [0,1] },
 >  y=0..8, arrows=THIN, title=`y'=y, y(0)=1` );

November 1995  1995 ICTCM – Houston  17

Picard Iterates: Example 1

To use picard it is necessary to define:

• the RHS of the differential equation,

   > RHS := (t,y)->y;

RHS := ( t, y) → y

• the initial condition,

   > IC := [ 0, 1 ];

IC := [ 0, 1 ]

• the initial guess,

   > phi[0] := t -> 0;

ϕ_{0} := 0

• and the number of iterations

   > N := 8;

N := 8

Picard Iterates: Example 1

Computation of the Picard iterates can be accomplished with explicit commands, or in a loop.

 > for k from 0 to N-1 do
 > phi[k+1] := picard( RHS, IC, phi[k] );
 > od;

ϕ_{1} := 1

ϕ_{2} := t → 1 + t

ϕ_{3} := t → 1 + t + \frac{1}{2} t^{2}

⋮

ϕ_{8} := t → 1 + \frac{1}{2} t^{2} + \frac{1}{120} t^{5} + \frac{1}{5040} t^{7} + t + \frac{1}{24} t^{4} + \frac{1}{6} t^{3} + \frac{1}{720} t^{6}

November 1995  1995 ICTCM – Houston  19

Picard Iterates: Example 1

Convergence can be demonstrated graphically.

> plot( { phi[i] $ i=0..N }, 0..2,
>  title=`Picard Iterates for y'=y, y(0)=1` );

Picard Iterates: Example 1

• Discussion

– The Picard iterates should look familiar (after sorting of the powers of x).

– This provides additional justification that the solution is y(x) = e^{x}. It also can be used to motivate the seach for solutions in the form of a power series.

– Animating the iterates can also be effective

   > for k from 0 to N do
   > P[k] := plot( phi[k], 0..2 );
   > od:
   > display( [ P[i] $ i=0..N ], insequence=true );

Picard Iterates: Example 2

Consider, now, the initial value problem: \frac{dy}{dt} = y^{2}, y(0) = \frac{1}{2}.

> DEplot( diff(y(t),t)=y(t)^2, [t,y], t=0..2, { [0,1/2] },
>  y=0..20, arrows=THIN, title=`y'=y^2, y(0)=1/2` );

November 1995  1995 ICTCM – Houston  22

Picard Iterates: Example 2

Define

> RHS := (t,y) -> y^2;

RHS := ( t, y) → y^{2}

> IC := [ 0, 1/2 ];

IC := [0, \frac{1}{2}]

> phi[0] := t -> 0;

ϕ_{0} := 0

> N := 6;

N := 6

November 1995  1995 ICTCM – Houston  23

Picard Iterates: Example 2

Compute

> for k from 0 to N-1 do
> phi[k+1] := picard( RHS, IC, phi[k] );
> od;

ϕ_{1} := \frac{1}{2}

ϕ_{2} := t → \frac{1}{2} + \frac{1}{4} t

ϕ_{3} := t → \frac{1}{3} + \frac{1}{48} ( 2 + t)^{3}

Picard Iterates: Example 2

ϕ_{4} := t → \frac{1}{2} + \frac{1}{4} t+ \frac{1}{16} t^{3} + \frac{1}{48} t^{4} + \frac{1}{1152} t^{6} + \frac{1}{8} t^{2} + \frac{1}{16128} t^{7} + \frac{1}{192} t^{5}

ϕ_{5} := t → \frac{1}{2} + \frac{71}{129024} t^{8} + \frac{13}{960} t^{5} + \frac{1}{4} t + \frac{1}{1032192} t^{12} + \frac{5}{774144} t^{11}

+ \frac{11}{322560} t^{10} + \frac{29}{16128} t^{7} + \frac{43}{290304} t^{9} + \frac{1}{130056192} t^{14}

+ \frac{1}{3901685760} t^{15} + \frac{1}{192} t^{6} + \frac{1}{16} t^{3} + \frac{1}{9289728} t^{13} + \frac{1}{32} t^{4} + \frac{1}{8} t^{2}

November 1995  1995 ICTCM – Houston  25

Picard Iterates: Example 2

Plotting the iterates is only of limited use

> plot( { phi[i] $ i=0..N }, 0..2.2,
>  title=`Picard Iterates for y'=y^2, y(0)=1/2` );

Picard Iterates: Example 2

• Discussion

– Iterates converge to ∑_{n=0}^{∞} \frac{t^{n}}{2^{n+1}} , which converges for |t| < 2.

– This can serve to motivate additional discussion of the existence theory for IVPs.

– Finite time “blow-up”’ is more clearly seen when the iterates are animated, and the exact solution is included in the plots

  > EX := unapply( rhs( dsolve({diff(y(t),t)=RHS(t,y),
  >  y(IC[1])=IC[2]}, y(t)) ), t ):
  > for k from 0 to N do
  > P[k] := plot( { EX, phi[k] }, 0..2, y=0..40 );
  > od:
  > display( [ P[i] $ i=0..N ], insequence=true );

Discontinuous Data, Systems, and Control

A parachutist whose mass is 75 kg drops from a helicopter hovering 4000 m above the ground and falls toward the earth under the influence of gravity. Assume the gravitational force is constant and that the force due to air resistance is proportional to the velocity of the parachutist, with a proportionality constant k_{1} = 15 kg/sec when the chute is closed and with a constant k_{2} = 105 kg/sec when the chute is open. If the chute does not open until 1 min after the parachutist leaves the helicopter, after how many seconds will she hit the ground?

Source: Nagel and Saff (p. 94)

Discontinuous Data, Systems, and Control

• Variations – Impact Velocity

– Repeat the problem with k_{2} = 30 kg/sec. Compare the two impact velocities.

– For what values of k_{2} is the impact velocity no greater than 10 m/sec?

– Determine the shortest jump (using k_{2} = 105 kg/sec) with an impact velocity that does not exceed 7.25 m/sec.

Discontinuous Data, Systems, and Control

• More Variations – Criteria for Chute Deployment

– Analyze the jump when the chute deploys at a given altitude.

– How would the motion differ if the parachute opened when the free-fall velocity reached a prescribed velocity.