7.22 boundary value problem, to solve (13.10.95)

7.22.1 Ulrich Klein

We have a problem in solving a differential equation. We have a differential 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 
 \  MAPLE  /  University of Waterloo. All rights reserved. Maple and Maple V 
 <____ ____>  are registered trademarks of Waterloo Maple Software. 
      |       Type ? for help. 
> read shoot; 
> ?shoot 
FUNCTION: shoot - shooting method for two-point boundary value 
CALLING SEQUENCE: shoot(deqns, vars, bc, init, options) 
           deqns   - ordinary differential equations in vars 
           vars    - variables to be solved for 
           bc      - boundary conditions 
           init    - initial values for shooting method 
           options - optional arguments 
- 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 
- 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') 
   - 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 
- 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. 
- Additional information about each iteration can be obtained 
  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:. 
- For more information on the shooting method, see Stoer and 
  Bulirsch (p. 469) or other standard numerical analysis texts. 
> restart; 
> with(plots): 
> read shoot; 
# 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