ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 17 -- Shooting Method

Outline of Unit 17

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

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

>

The basic idea behind the shooting method is to convert a boundary value problem (BVP) into an initial value problem (IVP). One parameter is introduced for each missing initial condition. The boundary conditions that are not initial conditions serve as the constraints to determine appropriate values for the parameters. That is, given an initial guess for the parameters, an iterative solver is used to find values of the parameters that produce (approximate) solutions that satisfy the boundary conditions.

The shooting method is quite general. Here it will be demonstrated only for two second-order linear two-point boundary value problems of the form

.

The IVP that will be used as the basis for the shooting method is

where is the unknown parameter. Call the solution to this IVP . The goal is to find such that ( ) = . To accomplish this, define a Maple procedure that accepts the IVP with a specific value for , the unknown function, and the value of and returns the value of ( )

> Yb := ( ivp, fn, b ) ->

> eval( fn, dsolve( ivp, fn, type=numeric )(b) );

>

and then use, e.g. , the Secant Method to solve for the value of for which ( ) - = 0.

Consider the BVP

> ode1 := diff( y(x), x$2 ) + 2 * diff( y(x), x ) - 10 * y(x) = 7*exp(-x) + 4;

> bc1 := y(0)=2, y(1)=-5;

>

with boundary conditions specified at the points

> a := 0;

> b := 1;

>

The corresponding IVP has initial conditions

> ic1 := bc1[1], D(y)(a)=alpha;

>

where is to be chosen to make

> constraint := lhs(bc1[2])-rhs(bc1[2]):

> constraint = 0;

>

Here is an implementation of the Secant Method :

> a0 := 0:

> y0 := eval( constraint, y(b)=Yb( eval( {ode1,ic1}, alpha=a0 ), y(x), b ) ):

> a1 := 1:

> y1 := eval( constraint, y(b)=Yb( eval( {ode1,ic1}, alpha=a1 ), y(x), b ) ):

> while fnormal(y1-y0)<>0 do

> z := solve( y1 = (y1-y0)/(a1-a0)*(a1-x), x );

> a0, a1 := a1, z;

> y0, y1 := y1, eval( constraint, y(b)=Yb( eval( {ode1,ic1}, alpha=z ), y(x), b ) );

> end do:

> alpha_opt1 := z;

>

The shooting method solution is the solution to the IVP corresponding to this optimal choice of the parameter :

> shoot_sol1 := dsolve( eval( {ode1,ic1}, alpha=alpha_opt1 ), y(x) );

>

Not surprisingly, because of the simple nature of this example, Maple is able to determine the exact solution to this problem.

> infolevel[dsolve] := 3:

> exact_sol1 := combine(dsolve( { ode1, bc1 }, y(x) ));

> infolevel[dsolve] := 0:

Methods for second order ODEs:

Trying to isolate the derivative d^2y/dx^2...

Successful isolation of d^2y/dx^2

-> Trying classification methods

trying a quadrature

trying high order exact linear fully integrable

trying differential order: 2; linear nonhomogeneous with symmetry [0,1]

trying 2nd order linear exact nonhomogeneous

trying a double symmetry of the form [xi=0, eta=F(x)]

trying linear constant coefficient

linear constant coefficient successful

>

To compare the solutions, plot the exact and shooting solutions

> plot( [rhs(exact_sol1),rhs(shoot_sol1)], x=0..1, color=[BLACK,RED],

> thickness=[1,3], legend=['exact','shooting'] );

>

Note that the value of the shooting method solution at the right-hand boundary point is

> evalf( eval( rhs(shoot_sol1), x=1 ) );

>

This also shows good agreement between the shooting method and exact solutions.

>

Consider the BVP

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

> bc2 := y(1)=7, y(4)=-1;

>

The boundary conditions are given at

> a := 1;

> b := 4;

>

The corresponding IVP has initial conditions

> ic2 := bc2[1], D(y)(a)=alpha;

>

where is to be chosen to make

> constraint := lhs(bc2[2])-rhs(bc2[2]):

> constraint = 0;

>

Here is an implementation of the Secant Method :

> a0 := 0:

> y0 := eval( constraint, y(b)=Yb( eval( {ode2,ic2}, alpha=a0 ), y(x), b ) ):

> a1 := 1:

> y1 := eval( constraint, y(b)=Yb( eval( {ode2,ic2}, alpha=a1 ), y(x), b ) ):

> while fnormal(y1-y0)<>0 do

> z := solve( y1 = (y1-y0)/(a1-a0)*(a1-x), x );

> a0, a1 := a1, z;

> y0, y1 := y1, eval( constraint, y(b)=Yb( eval( {ode2,ic2}, alpha=z ), y(x), b ) );

> end do:

> alpha_opt2 := z;

>

The shooting method solution is the solution to the IVP corresponding to this optimal choice of the parameter.

> shoot_sol2 := dsolve( eval( {ode2,ic2}, alpha=alpha_opt2 ), y(x) );

>

Again, because of the simple nature of this example, Maple is able to determine the exact solution to this problem.

> exact_sol2 := combine(dsolve( { ode2, bc2 }, y(x) ));

>

To compare the solutions, plot the exact and shooting solutions

> plot( [rhs(exact_sol2),rhs(shoot_sol2)], x=1..4,

> color=[BLACK,RED],thickness=[1,3],

> legend=['exact','shooting'] );

>

Note also that the value of the shooting method solution at the right-hand boundary point is

> evalf( eval( rhs(shoot_sol2), x=1 ) );

>

This also shows good agreement between the shooting method and exact solutions.

>