Mathematical Models with Boundary Value

Problems

Douglas B. Meade  Allan A. Struthers
Department of Mathematics  Department of Mathematical Sciences
University of South Carolina  Michigan Technical University
Columbia, SC 29208  Houghton, MI 49931
(803) 777-6183  (906) 487-3541

     meade@math.sc.edu                 struther@math.mtu.edu
http://www.math.sc.edu/˜ meade/    http://www.math.mtu.edu/˜ struther/

8 November 1997

November 1997  ICTCM X (Chicago)  0

Mathematical Models with Boundary Value Problems

Douglas B. Meade  Allan A. Struthers
Department of Mathematics  Department of Mathematical Sciences
University of South Carolina  Michigan Technical University
Columbia, SC 29208  Houghton, MI 49931
(803) 777-6183  (906) 487-3541

     meade@math.sc.edu                 struther@math.mtu.edu
http://www.math.sc.edu/˜ meade/    http://www.math.mtu.edu/˜ struther/

8 November 1997

November 1997  ICTCM X (Chicago)  1

Outline

• Past and Present

• Computer Algebra Tools for Solving BVPs

• Application — Infinite Rotating Disk (Chemical Engineering)

• Future

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::function )
> local ODE, BC, gSOL, Y, DVAR, IVAR, CONST, cSOL;
> ODE := select( has, BVP, diff );
> BC := BVP minus ODE;
> DVAR, IVAR := op(0,FN), 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)

November 1997  ICTCM X (Chicago)  6

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

November 1997  ICTCM X (Chicago)  11

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

> fBVPsolve := proc( BVP::set(equation),
>  FN::function,
>  RANGE::function=range(numeric) )
> local i, auxIC, ODE, BC, IC, IVP, IVAR, RANGE2,
>  RANGE3, PARAM, POINTS, RESULT, X0;
> global n, BCERR;

>

> IVAR   := op(1,FN);
> ODE    := select( has, BVP, IVAR );
> BC     := BVP minus ODE;
> X0     := op(1,lhs(RANGE));

>

> IC     := select( (l,r) -> evalb(op(1,lhs(l))=r), BC, X0 );
> BC     := convert(map(e->lhs(e)-rhs(e), BC minus IC ),list);

>

> POINTS := convert(map2( op, 1, indets(BC,function) ),list);

>

> RANGE2 := [RANGE];
> auxIC  := {seq(lhs(RANGE2[i])=_alpha[i],i=1..nops(RANGE2))};
> RANGE3 := {seq(_alpha[i]=rhs(RANGE2[i]),i=1..nops(RANGE2))};

>

> IVP    := ODE union IC union auxIC;

Numerical Solutions to Second-Order BVPs

> BCERR  := subs( _IVP=IVP, _FN=FN, _POINTS=POINTS, _BC=BC,
>  proc( )
>  local i, j, ARGS, FF, PTS, SOL, SYS, VALS, bcerr;
>  global n;
>  n := n+1;

>

>  ARGS := seq( _alpha[i]=args[i], i=1..nargs );
>  SYS := subs( ARGS, _IVP );
>  PTS := convert(map2( op, 1, indets(_BC,function) ),list);

>

>  SOL := dsolve( SYS, _FN, type=numeric,
>  value=array(evalf(PTS)) );

>

>  FF := convert(map( convert, SOL[1,1], D ),list);
>  VALS := [ seq( seq( op(0,FF[j])(PTS[i]) = SOL[2,1][i,j],
>  i=1..nops(PTS) ), j=2..nops(FF) ) ];

>

>  bcerr := evalf( op( subs( VALS, _BC ) ) );
>  userinfo( 3, BCERR,
>  sprintf(`iteration %2d: alpha=%10.6f bcerr=%10.6f`,
>  n, rhs(ARGS), bcerr) );
>  RETURN( bcerr );
>  end );

Numerical Solutions to Second-Order BVPs

> PARAM  := op( map( lhs, RANGE3 ) );

>

> Digits := 5:
> n      := 0:
> RESULT := fsolve( 'BCERR'(PARAM)=0, {PARAM}, op(RANGE3) );
> Digits := 10:

>

> RETURN( [RESULT,dsolve(subs(RESULT,IVP),FN,type=numeric)] );
> end:

>

> infolevel[BCERR]:=3:

Example 1: Linear Two-Point BVP

> ODE := diff( y(x),x$2 ) + y(x) = cos(x);
> IC := y(1)=3;
> BC := y(3)=1;
> PARAM, SOL := op(fBVPsolve( {ODE,IC,BC}, y(x), D(y)(1)=0..5)):

BCERR:   "iteration  1: alpha=  2.500000 bcerr=  -.213300"
BCERR:   "iteration  2: alpha=  0.000000 bcerr= -2.489100"
BCERR:   "iteration  3: alpha=   .100000 bcerr= -2.398200"
BCERR:   "iteration  4: alpha=  -.010000 bcerr= -2.498800"
BCERR:   "iteration  5: alpha=  2.580837 bcerr=  -.139900"
BCERR:   "iteration  6: alpha=  2.734313 bcerr=  -.000500"
BCERR:   "iteration  7: alpha=  2.734863 bcerr=  -.001400"
BCERR:   "iteration  8: alpha=  2.734012 bcerr=  -.001800"
BCERR:   "iteration  9: alpha=  2.734292 bcerr=  -.000500"
BCERR:   "iteration 10: alpha=  2.734485 bcerr=  0.000000"
BCERR:   "iteration 11: alpha=  2.734485 bcerr=  -.003434"
BCERR:   "iteration 12: alpha=  2.734485 bcerr=  -.003434"
BCERR:   "iteration 13: alpha=  2.738262 bcerr=   .000000"

> PARAM;

> 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] );

November 1997  ICTCM X (Chicago)  17

> 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;
> PARAM, SOL := op(fBVPsolve({ODE,IC,BC}, y(x), D(y)(1)=0..15)):

BCERR: "iteration  1: alpha=  7.500000 bcerr= -2.262600"
BCERR: "iteration  2: alpha=  0.000000 bcerr= -5.962100"
BCERR: "iteration  3: alpha=   .100000 bcerr= -5.912900"
BCERR: "iteration  4: alpha=  -.010000 bcerr= -5.967200"
BCERR: "iteration  5: alpha= 11.728022 bcerr=  -.177000"
BCERR: "iteration  6: alpha= 12.085920 bcerr=   .000500"
BCERR: "iteration  7: alpha= 12.084912 bcerr=  -.002450"
BCERR: "iteration  8: alpha= 12.084914 bcerr=  -.002449"
BCERR: "iteration  9: alpha= 12.089880 bcerr=  -.000000"

> PARAM;

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

> dsolve( { ODE, IC, BC }, y(x) );
> exsol := BVPsolve( { ODE, IC, BC }, y(x));

exsol := 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;
> PARAM, SOL := op(fBVPsolve({ODE,IC,BC}, y(x), D(y)(1)=0..5 )):

BCERR: "iteration  1: alpha=  2.500000 bcerr=  -.313240"
BCERR: "iteration  2: alpha=  0.000000 bcerr=   .609290"
BCERR: "iteration  3: alpha=   .100000 bcerr=   .605030"
BCERR: "iteration  4: alpha=  -.010000 bcerr=   .609630"
BCERR: "iteration  5: alpha=  1.250000 bcerr=   .313420"
BCERR: "iteration  6: alpha=  1.375000 bcerr=   .262170"
BCERR: "iteration  7: alpha=  2.044988 bcerr=  -.059830"
BCERR: "iteration  8: alpha=  1.934214 bcerr=  -.002030"
BCERR: "iteration  9: alpha=  1.930259 bcerr=  -.000030"
BCERR: "iteration 10: alpha=  1.930200 bcerr=  -.005350"
BCERR: "iteration 11: alpha=  1.930200 bcerr=  -.005351"
BCERR: "iteration 12: alpha=  1.919758 bcerr=  -.000008"
BCERR: "iteration 13: alpha=  1.919743 bcerr=  -.000000"

> PARAM;

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

Application — Infinite Rotating Disk

• Physical Problem determine the steady fluid flow generated when the infinite plane z = 0, immersed in a Newtonian viscous fluid, rotates about the axis r = 0 with constant angular velocity ω

November 1997  ICTCM X (Chicago)  24

Application — Infinite Rotating Disk

• Notation

p  pressure
(v^{r}, v^{θ}, v^{z}) velocity
ρ  fluid density
ν  kinematic viscosity

• dimensionless variables

z^{∗} = \sqrt{\frac{ω}{ν}} z
v^{r} = rωF (z^{∗})
v^{θ} = rωG(z^{∗})
v^{z} = \sqrt{ων}H(z^{∗})
p = ρωνP (z^{∗})

November 1997  ICTCM X (Chicago)  25

Application — Infinite Rotating Disk

• boundary value problem

H′ = −2F
F′′ = F^{2} − G^{2} + F′H
G′′ = 2F G + HG′
P′ = −HH′ + H′′
F (0) = H(0) = P (0) = 0
G(0) = 1
lim_{z^{∗}→∞} F (z^{∗}) = lim_{z^{∗}→∞} G(z^{∗}) = 0

Maple Solution via Shooting Method

> restart;
> read shoot;
> with(plots):
> FNS:={ F(Z), G(Z), H(Z), Fp(Z), Gp(Z) }:
> ODE:={diff(H(Z),Z)  = -2*F(Z),
>  diff(F(Z),Z)  = Fp(Z),
>  diff(Fp(Z),Z) = -G(Z)^2+F(Z)^2+Fp(Z)*H(Z),
>  diff(G(Z),Z)  = Gp(Z),
>  diff(Gp(Z),Z) = 2*F(Z)*G(Z)+H(Z)*Gp(Z)}:
> IC:={ F(0)=0, G(0)=1, H(0)=0, Fp(0)=alpha, Gp(0)=beta }:
> BC:={ F(10)=0, G(10)=0 }:

Maple Solution via Shooting Method

> infolevel[shoot]:=1:
> S:=shoot( ODE union IC, FNS, BC, [alpha=0.51, beta=-0.62] ):

shoot:   alpha = .51, beta = -.62
shoot:   alpha = .5100033193, beta = -.6161487475
shoot:   alpha = .5102117355, beta = -.6159070527
shoot:   alpha = .5102135910, beta = -.6159096496

Maple Solution via Shooting Method

> P:=Z->-H(Z)^2/2-2*F(Z):
> C:=odeplot(S,[ [Z,F(Z)], [Z,G(Z)],
>  [Z,H(Z)], [Z,-P(Z)] ], 0..10 ):

>

> subs( [seq( op(0,lhs(f))=rhs(f(3/4)),
>  f=select(has,S,{F,G,H}) )] ):
> [op(%),P=subs(%, H^2+2*F)]:

>

> L:=textplot( subs(z=3/4, map(f->[z,rhs(f)+1/45,lhs(f)],%)),
>  font=[TIMES,ITALIC,18], align=ABOVE ):
> display( {C,L}, labels=[`z*`,``],
>  title=`Infinite Disk: alpha=0.51, beta=-0.62` );

Maple Solution via Shooting Method

> S2:=shoot( ODE union IC, FNS, BC, [alpha=0.5, beta=-0.61] ):

shoot:   alpha = .5, beta = -.61
shoot:   alpha = .5181506133, beta = -.5834922887
shoot:   alpha = .5178653168, beta = -.5900697906
shoot:   alpha = .5019889119, beta = -.5940268791
shoot:   alpha = .5040798081, beta = -.5927957583
shoot:   alpha = .5042719718, beta = -.5927086067
shoot:   alpha = .5042737786, beta = -.5927079110

Maple Solution via Shooting Method

> C:=odeplot(S2,[ [Z,F(Z)], [Z,G(Z)],
>  [Z,H(Z)], [Z,-P(Z)] ], 0..12 ):

>

> subs( [seq( op(0,lhs(f))=rhs(f(3/4)),
>  f=select(has,S2,{F,G,H}) )] ):
> [op(%),P=subs(%, H^2+2*F)]:

>

> L:=textplot( subs(z=3/4, map(f->[z,rhs(f)+1/45,lhs(f)],%)),
>  font=[TIMES,ITALIC,18], align=ABOVE ):
> display( {C,L}, labels=[`z*`,``],
>  title=`Infinite Disk: alpha=0.50, beta=-0.61` );

Ongoing Efforts

• 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 fBVPsolve (and shoot)

• 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?

November 1997  ICTCM X (Chicago)  34