Theory and Applications of Boundary Value Problems for ODEs

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/

17 October 1997

October 1997  AMS Southeastern Section Meeting (Atlanta)  0

Theory and Applications of Boundary Value Problems for ODEs

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/

17 October 1997

Outline

• Introduction

• BVPs in Introductory ODE Courses – Past and Present

– Computer Algebra Tools for Solving BVPs – Future

Play at the Plate

Bottom of the ninth; one out; tying run on third. A fly ball is hit to an outfielder positioned 300 feet from the plate. The catch is made on the run; 0.6 seconds later the ball starts a perfect flight to the catcher (i.e., on target with no bounces). The speedy baserunner runs from third to home in 4.0 seconds. How hard must the ball be thrown to prevent the runner from scoring?

— based on “The Physics of Baseball”, Robert K. Adair, 1994.

October 1997  AMS Southeastern Section Meeting (Atlanta)  3

Play at the Plate

• Assumptions – throw arrives at home plate without bouncing

– catcher needs 0.3 seconds to apply the tag

– no wind

– quadratic air resistance (coefficient k)

– the ball has mass 5 oz.

– the ball is released from a height of 8 feet above the ground – g = 9.8j m/sec^{2}

– the umpire is in a good position to make the call

October 1997  AMS Southeastern Section Meeting (Atlanta)  4

Play at the Plate

• Problem Formulation:

determine x′(0) and ∥x′(0)∥ where x is a solution to

mx′′ = −k∥x′∥x′ − g
x(0) = 300i + 8j
x(3.1) = 0

Play at the Plate

• Further Discussion

– Why does the fielder try to make the catch while running towards home plate?

– Are there advantages of bouncing the ball to the catcher?

– Given the maximum initial velocity for a player, where should the fielder be positioned?

Past and Present

Boundary Value Problems. When the values of a solution to a differential equation are specified at two different points, these conditions are called boundary conditions. (In contrast, initial conditions specify the values of a function and its derivative at the same point.) Given that every solution to y′′ + y = 0 is of the form

y(x) = c_{1} cos x + c_{2} sin x,

where c_{1} and c_{2} are arbitrary constants, show that

(a) there is a unique solution that satisfies the boundary conditions y(0) = 2 and y(π/2) = 0;

(b) there is no solution that satisfies y(0) = 2 and y(π) = 0;

(c) there are infinitely many solutions that satisfy y(0) = 2 and y(π) = −2.

— Fundamentals of Differential Equations, Nagle & Saff (# 30, p. 159)

Computer Algebra Tools for Solving BVPs

> ODE := diff( y(x), x$2 ) + y(x) = 0:
> BCa := y(0)=2, y(Pi/2)=0:
> SOLa := dsolve( { ODE, BCa }, y(x) );

SOLa := y(x) = 2 cos(x)

> BCb := y(0)=2, y(Pi)=0:
> SOLb := dsolve( { ODE, BCb }, y(x) );

SOLb :=

> BCc := y(0)=2, y(Pi)=-2:
> SOLc := dsolve( { ODE, BCc }, y(x) );

SOLc := y(x) = C1 sin(x) + 2 cos(x)

Customized symbolic BVP solver: bvpsolve

> bvpsolve := proc( BVP::set(equation), FN )
> local ODE, BC, gSOL, Y, DVAR, IVAR, CONST, cSOL;
> ODE := select( has, BVP, diff );
> BC := BVP minus ODE;
> DVAR := op(0,FN); IVAR := op(1,FN);
> gSOL := dsolve( ODE, FN );
> if gSOL=NULL then
>  RETURN( `No general solution found` )
> fi;
> if nops([gSOL])>1 then
>  RETURN( `Nonunique general solution found` )
> fi;
> CONST := indets(gSOL,name) minus {IVAR};
> Y := unapply( rhs( gSOL ), IVAR );
> cSOL := solve( eval( subs( DVAR=Y, BC ) ), CONST );
> if cSOL=NULL then
>  RETURN( `No particular solution satisfies all BCs` )
> fi;
> subs( cSOL, gSOL );
> end:

Simple Examples

> bvpsolve( { ODE, BCa }, y(x) );

y(x) = 2 cos(x)

> bvpsolve( { ODE, BCb }, y(x) );

No particular solution satisfies all BCs

> bvpsolve( { ODE, BCc }, y(x) );

y(x) = 2 cos(x) + C2 sin(x)

October 1997  AMS Southeastern Section Meeting (Atlanta)  10

More Challenging Problems

• multi-point boundary conditions y′′ + y = 0, y(π/4) + y(π/2) = 1, y′(π/4) + y′(π/2) = −1

• nonlinear boundary conditions y′′ + y = 0, y(0) = 2, y(π/2)^{2} + y′(π) = 1

• nonlinear differential equation y′′ + y^{2} = 0, y(0) = 2, y(π/2) = 0

Multi-point Boundary Conditions

> ODE := diff( y(x), x$2 ) + y(x) = 0:
> BC2p := D(y)(Pi/2)+D(y)(Pi/4)=-1, y(Pi/4)+y(Pi/2)=1;

BC2p := D(y)(\frac{1}{2} π) + D(y)(\frac{1}{4} π) = −1, y(\frac{1}{4} π) + y(\frac{1}{2} π) = 1

> dsolve( { ODE, BC2p }, y(x) );
> bvpsolve( { ODE, BC2p }, y(x) );

y(x) = \frac{1}{2} \sqrt{2} cos(x) + (1 − \frac{1}{2} \sqrt{2}) sin(x)

Nonlinear Boundary Conditions

> ODE := diff( y(x), x$2 ) + y(x) = 0:
> BCnl := y(0)=2, y(Pi/2)^2 + D(y)(Pi)=1;

BCnl := y(0) = 2, y(\frac{1}{2} π)^{2} + D(y)(π) = 1

> dsolve( { ODE, BCnl }, y(x) );
> bvpsolve( { ODE, BCnl }, y(x) );

y(x) = 2 cos(x) + RootOf( Z^{2} − Z − 1) sin(x)

 > SS := allvalues( " );

SS := y(x) = 2 cos(x) + ( \frac{1}{2} \sqrt{5} + \frac{1}{2} ) sin(x), y(x) = 2 cos(x) + ( \frac{1}{2} − \frac{1}{2} \sqrt{5}) sin(x)

 > plot( map( rhs, [SS] ), x=0..Pi );

Nonlinear ODE I

> ODE2 := diff( y(x), x$2 ) + y(x)^2 = 0;

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

> dsolve( { ODE2, BCa }, y(x) );
> bvpsolve( { ODE2, BCa }, y(x) );

Nonunique general solution found

> dsolve( ODE2, y(x) );

x = ∫_{0}^{y(x)}− \frac{3}{\sqrt{−6 y1^{3} + 9 C1}} dy1 − C2 , x = ∫_{0}^{y(x)} \frac{3}{\sqrt{−6 y2^{3} + 9 C1}} dy2 − C2

October 1997  AMS Southeastern Section Meeting (Atlanta)  15

Nonlinear ODE II

> ODE3 := diff( y(x),x$2 ) + y(x)^2 = 7*x^(-4);

ODE3 := (\frac{∂^{2}}{∂x^{2}} y(x)) + y(x)^{2} = \frac{7}{x^{4}}

> IC := y(1)=1; BC := y(3)=1/3;

IC := y(1) = 1

BC := y(3) = \frac{1}{3}

> dsolve( { ODE3, IC, BC }, y(x) );
> bvpsolve( { ODE3, IC, BC }, y(x) );

No general solution found

Numerical Solutions to Second-Order BVPs

> bvpsolve2 := proc( alpha )
> local SYS, Y, dY, bcerr;
> global n, SOL;
> n := n+1;

>

> SYS := { ODE, IC, D(y)(1)=alpha };
> SOL := dsolve( SYS, y(x), type=numeric,
>  output=listprocedure );

>

> Y  := subs( SOL, y(x) );
> dY := subs( SOL, diff(y(x),x) );
> bcerr := eval( subs( D(y)=dY, y=Y, lhs(BC) - rhs(BC) ) );

>

> print(printf( `iteration %d: alpha=%10.6f, bcerr=%10.6f\n`,
>  n, alpha, bcerr ));
> RETURN( bcerr )
> end:

Example 1: Linear Two-Point BVP

> ODE := diff( y(x),x$2 ) + y(x) = cos(x);
> IC := y(1)=3;
> BC := y(3)=1;
> Digits := 5:
> n:=0:
> fsolve( 'bvpsolve2(alpha)'=0, alpha, 0..5);

iteration 1: alpha=  0.000000, bcerr= -2.489894
iteration 2: alpha=  5.000000, bcerr=  2.056592
iteration 3: alpha=  2.500000, bcerr=  -.216650
iteration 4: alpha=  2.738262, bcerr=   .000000
iteration 5: alpha=  2.738253, bcerr=  -.000007

2.7383

> P1 := plots[odeplot]( SOL, [x,y(x)], 0..Pi):
> P2 := plot( [[1,3],[3,1]], style=POINT,
>  color=BLUE, symbol=BOX ):
> plots[display]( [P1,P2] );

October 1997  AMS Southeastern Section Meeting (Atlanta)  19

> dsolve( { ODE, IC, BC }, y(x) );

y(x) = \frac{1}{2} sin(x) x

+ \frac{1}{2} \frac{(−3 sin(1) − 8 + 4 sin(1)^{3} + 12 sin(1)^{2}) cos(1) cos(x)}{sin(1)^{2} − 1}

+ \frac{(2 sin(1)^{3} + 6 sin(1)^{2} − 2 sin(1) − 1) sin(x)}{sin(1)}

> bvpsolve( { ODE, IC, BC }, y(x));

y(x) = \frac{1}{2} sin(x) x

+ \frac{1}{2} \frac{(−3 sin(1) − 8 + 4 sin(1)^{3} + 12 sin(1)^{2}) cos(1) cos(x)}{sin(1)^{2} − 1}

+ \frac{(2 sin(1)^{3} + 6 sin(1)^{2} − 2 sin(1) − 1) sin(x)}{sin(1)}

> evalf( subs( x=1, diff( rhs("), x ) ));

2.738261577

Example 2: Linear BVP with Robin BC

> ODE := diff( y(x),x$2 ) + y(x) = cos(x);
> IC := y(1)=3;
> BC := y(3)+D(y)(3)=1;
> Digits := 5:
> n:=0:
> fsolve( 'bvpsolve2(alpha)'=0, alpha, 0..15);

iteration 1: alpha=  0.000000, bcerr= -5.962131
iteration 2: alpha= 15.000000, bcerr=  1.435127
iteration 3: alpha=  7.500000, bcerr= -2.263501
iteration 4: alpha= 12.089880, bcerr=   .000000
iteration 5: alpha= 12.089840, bcerr=  -.000019

12.090

> plots[odeplot]( SOL,
>  [ [x,y(x)], [x,y(x)+diff(y(x),x)-1] ],
>  0..5 );

> dsolve( { ODE, IC, BC }, y(x) );
> bvpsolve( { ODE, IC, BC }, y(x));

y(x) = \frac{1}{2} sin(x) x + \frac{1}{2} (12 sin(1)^{4} + 24 sin(1)^{3} + 8 cos(1) sin(1)^{3} + 24 cos(1) sin(1)^{2} − 2 cos(1) sin(1) − 6 cos(1) − 9 sin(1)^{2} − 16 sin(1)) cos(x) /(2 sin(1)^{2} − 1 − 2 cos(1) sin(1)) + \frac{1}{2} (−12 sin(1)^{2} + 3 + 11 cos(1) sin(1) − 12 cos(1) sin(1)^{3} − 24 cos(1) sin(1)^{2} + 8 sin(1)^{4} − 18 sin(1) + 24 sin(1)^{3} + 4 cos(1))sin(x) /( 2 sin(1)^{2} − 1 − 2 cos(1) sin(1))

> evalf( subs( x=1, diff(rhs(exSOL),x) ) );

12.08987956

Example 3: Nonlinear Two-Point BVP

> ODE := diff( y(x),x$2 ) + y(x)^2 = 3*x^(-2);
> IC := y(1)=1;
> BC := y(3)=1/3;
> Digits := 5:
> n:=0:
> fsolve( 'bvpsolve2(alpha)'=0, alpha, 0..5);

iteration 1: alpha=  0.000000, bcerr=   .606250
iteration 2: alpha=  5.000000, bcerr= -2.090706
iteration 3: alpha=  2.500000, bcerr=  -.318299
iteration 4: alpha=  1.858074, bcerr=   .031241
iteration 5: alpha=  1.918843, bcerr=   .000459
iteration 6: alpha=  1.919744, bcerr=  -.000000
iteration 7: alpha=  1.919743, bcerr=   .000000

1.9197

      > plots[odeplot]( SOL, [x,y(x)], 0..4 );

Ongoing Efforts

• General Goals – BVPs deserve more attention than they currently receive – should some material be omitted from traditional courses to

allow discussion of BVPs?

– how do BVPs fit into newer courses that emphasize

dynamical systems?

• Implementation

– continued development of specialized and general CAS tools for symbolic, numeric, and graphic analysis of BVPs

– collect realistic and interesting examples of BVPs for use as homework, in-class examples, and projects

– generalize and extend bvpsolve and bvpsolve2 (and shoot)

October 1997  AMS Southeastern Section Meeting (Atlanta)  26