#### 7.22 boundary value problem, to solve (13.10.95)

##### 7.22.1 Ulrich Klein

We have a problem in solving a diﬀerential equation. We have a diﬀerential equation, say de and an initial value condition: $$y(0)=0$$ and a boundary value condition: D(y)(20)=0, eg.

If we try to solve this equation using the dsolve command, Maple complains about the second condition. Besides, we use the numeric option of the dsolve command.

Do you know a way to solve such an equation with Maple?

##### 7.22.2 Douglas B. Meade (16.10.95)

Maple’s numeric dsolve is, as the on-line help describes, of use only for initial value problems. You are asking for the numerical solution of a two-point boundary value problem.

While you do not provide enough information to know exactly what type of ODE you are attempting to solve, I might have something that will be of interest.

I have recently developed a Maple implementation of the simple shooting method.

This procedure works very much like dsolve[numeric]. It must be used with caution, as the numerical method is not always stable. A paper introducing this procedure (shoot) and demonstrating its usage on three problems will appear in the upcoming Special Issue of MapleTech. (These examples are nonlinear and some involve systems.)

The on-line help for shoot, which includes two (simple) examples, is attached to the end of this message.

Anyone interested in this paper and/or procedure should send mail to: Douglas B. Meade

|\^/|     Maple V Release 3 (University of South Carolina)
._|\|   |/|_. Copyright (c) 1981-1994 by Waterloo Maple Software and the
<____ ____>  are registered trademarks of Waterloo Maple Software.
|       Type ? for help.
> ?shoot
FUNCTION: shoot - shooting method for two-point boundary value
problems

CALLING SEQUENCE: shoot(deqns, vars, bc, init, options)

PARAMETERS:
deqns   - ordinary differential equations in vars
vars    - variables to be solved for
bc      - boundary conditions
init    - initial values for shooting method
options - optional arguments

SYNOPSIS:
- shoot uses the shooting method to compute approximate solutions
to two-point boundary value problems for a first-order system
of ordinary differential equations using the shooting method.
The default output of this procedure is a list of procedures,
as generated by numeric dsolve (see ?dsolve,numeric)

- Each of the first four arguments must be specified as a single
equation or a list or set of equations (very much like the
arguments to dsolve).

- The first two arguments are exactly the same as for dsolve.
That is, the first argument contains the differential equations
and initial conditions. The second argument contains the
dependent variable(s) in the problem.

- The initial conditions will include both the boundary values
at one endpoint and the auxiliary initial conditions needed
for the shooting method (but not the initial conditions for
the sensitivity equations). The initial values for the auxiliary
conditions will be given in terms of the shooting method
parameters. There should be an equal number of differential
equations and initial conditions.

- The the third argument contains the boundary condition(s) -- as
a single equation or as a list or set of equations -- at the
second boundary point. These are the control equations for the
iterations.

- The final required argument is a single equation or list or
set of equations containing the parameter values for the
first iteration of the shooting method.

- Any optional arguments must be in the form of keyword=value.
These can be any valid optional argument to dsolve, or one of
the following (with default settings listed in parentheses):
'bcerr'    = bctol           (Float(1,2-Digits))
'iterr'    = ittol           (Float(1,2-Digits))
'maxiter'  = maxit           (10)
'itmethod' = itmeth          ('newton')
where:
- bctol is the maximum acceptable L2-error between the computed
and exact boundary conditions;
- ittol is the maximum acceptable L2-error between successive
parameter estimates;
- maxiter is the maximum number of iterations
- itmeth is the root-finding method used to compute new
parameter estimates (only Newton-Raphson is available)

- The only known option to dsolve that is not supported is
output=listprocedure. However, THERE IS NO GUARANTEE THAT ALL
COMBINATIONS OF OPTIONS WILL WORK AS EXPECTED.

- Each iteration of the shooting method begins by solving (by
numeric dsolve) the initial value problem using the most recent
set of shooting method parameters. Each boundary condition is
then evaluated. If these conditions are not satisfied to the
requested tolerance, one iteration of the Newton-Raphson method
is used to generate a new set of shooting method parameters.

via the infolevel command. To see each set of parameters,
use infolevel[shoot]:=1:. To include the error at the
boundary, use infolevel[shoot]:=2:. Use infolevel[shoot]:=3:
to see the full system of differential equations involved in
the iterations. To reset to normal usage, use infolevle[shoot]:=0:.

Bulirsch (p. 469) or other standard numerical analysis texts.

EXAMPLES:
> restart;
> with(plots):
--------------------------------------------------------------------------------
# Example 1: x'' + x = cos(t), x(1)=3, x'(3)=1
> ODE := { diff( x(t), t ) = y(t), diff( y(t), t ) = -x(t)+cos(t) };

d                 d
ODE := {---- x(t) = y(t), ---- y(t) = - x(t) + cos(t)}
dt                dt
> FNS := { x(t), y(t) };

FNS := {x(t), y(t)}
> IC := { x(1)=3, y(1)=alpha };

IC := {x(1) = 3, y(1) = alpha}
> BC := y(3)=1;

BC := y(3) = 1
> S := shoot( ODE union IC, FNS, BC, alpha=0 );

S := [t = proc(t) ... end, x(t) = proc(t) ... end, y(t) = proc(t) ... end]

> odeplot( S, [t,x(t)], 1..3 );
--------------------------------------------------------------------------------
# Example 2: x''+x = cos(t), x(1)=3, x(3)+x'(3)=1
> infolevel[shoot]:=1:
> S := shoot( ODE union IC, FNS, x(3)+y(3)=1, alpha=0, bcerr=Float(5,-7) ):
shoot:   New parameter values:    alpha = 0
shoot:   New parameter values:    alpha = 12.08987391
shoot:   New parameter values:    alpha = 12.08987493

> infolevel[shoot]:=0:
> odeplot( S, [ [t,x(t)], [t, x(t)+y(t)] ], 1..3, 0..20, labels=[x,t] );
# The exact solution is not hard to find for this problem:
> X := rhs( dsolve( { diff(x(t),t\$2)+x(t)=cos(t), x(1)=3 }, x(t), method=laplace ) ):
> Xp := diff(X,t):
> c := solve( subs(t=3,X+Xp)=1, {D(x)(1)} ):
# Note that the final value of the shooting parameter is correct
# to 5 digits to the right of the decimal point
> evalf(c);

{D(x)(1) = 12.08987954}
> X := subs(c,X): evalf(X);

.4207354924 (t - 1.) cos(t - 1.) + .2701511530 (t - 1.) sin(t - 1.)

+ 11.66914405 sin(t - 1.) + 3. cos(t - 1.)
# The results are very comparable when viewed graphically
> p1 := odeplot( S, [t,x(t)], 1..3, style=point ):
> p2 := plot( X, t=1..3, color=red ):
> display( {p1,p2} );

SEE ALSO: dsolve[numeric], userinfo PS: It worked very well for our problem. U. Klein