Solution of a two point boundary problem: shooting (BVPSOL)

directly to the input form: predefined cases

directly to the input form: user defined case

• In the following, y, F and R denote functions with n components and y' denotes componentwise differentiation. Typically where will be a ''free'' real parameter x with the meaning ''space'' (or sometimes ''time''). We search for a solution y as function of x, given some relation between y, y' and x and some conditions on y at two points x=a, x=b, the so called boundary conditions.
  Specifically the problem dealt with here is:
  Find an (at least) once differentiable function y defined on the interval [a,b] such that there holds
  y'(x) = F(x,y(x)) for a <= x <= b and R(y(a),y(b)) = 0 .
  
• The essential difference opposed to initial value problems is the fact that the length of the interval of existence of the solution is part of the problem setting and that conditions on the solution are given at two different points. This type of problem is called a boundary value problem, BVP for short. Even for a F linear in y, i.e. F(x,y(x)) = A(x)y(x) + b(x), the solvability of the problem is not guaranteed, as can be seen by from the simple case
  y₁' = y₂, y₂' = -y₁
  with the general solution
  y₁(x)=A sin(x)+ B cos(x)
  if one requires
  y₁(0) = 0 and y₁(π) = 1
  as boundary conditions.
• In practice the differential equations appearing here are often of higher order, but since every such ODE can be reformulated as a first order one and the software used here requires this, we stick on the special form given above.
• There are quite little general results on the solvability and even less on the unique solvability of the problem and one must be aware of the possibility of several or no solutions at all. A simple example with two well defined isolated solutions is given in the prepared examples here.
• Since initial value problems can be solved much easier it is an old idea to reformulate the problem as follows:
  Find s such that the solution of the initial value problem
  y'(x) = F(x,y(x)) with y(a) = s exists for a <= x <= b and satisfies R(s,y(b)) = 0 .
• This is most simply explained by the case of the oblique shot, which gave the method its name. The physical problem formulation, neglecting air resistance, is given by
  
  1. (d/dt)² y(t) = -g
  2. x(t)= v₀ t cos(φ), v₀ = (((d/dt)x)(0)²+((d/dt)y)(0)²)^1/2
  3. x(0) = y(0) = 0 , t >= 0
  where v₀ is initial speed, φ the angle of the flight curve with the x-axis at the origin and (x(t),y(t)) the parameter representation of the curve (with t as time). Solving the differential equation for y and expressing y as a function of x we arrive at
  y(x) = x tan(φ) - g x² / (2 v₀² cos²(φ)) .
  With φ in the range [0,π/2] y will be positive in the interval
  0 <= x <= v₀² sin( 2 φ )/g .
  Hence we cannot presribe arbitrarily large b such that y(b) = 0 for the fixed given v₀. But with b in the given range we might now determine s = y'(0) such that the solution of the initial value problem
  y'' = -C(φ) , y(0) = 0 , y'(0) = s , C(φ) = g/ (2 v₀² cos²(φ))
  solves the boundary problem
  y'' = -C(φ) , y(0) = 0 , y(b) = 0 .
  (Observe that the physical situation is a bit tricky since C depends on
  φ = tan^-1(s) ).
• The choice of s as initial value makes the solution y of the initial value problem also a function of s which we denote by y = y(x;s) . Hence the reformulation of the boundary value problem in these terms reads
  Find s such that the solution of
  1. y(x;s)' = F(x,y(x;s)) for a <= x <= b
  2. y(a;s) = s
  satisfies
  R(s,y(b;s)) = 0 .
  
• This approach has a fundamental and unavoidable problem: Although for smooth F the solution of the initial problem exists on an interval [a,a+δ] it may well be that δ < b-a because the solution manifold of the ODE has moving singularities. This may be overcome by
• Multiple shooting: we decompose the long interval [a,b] into sufficiently small subintervals, solve on each subinterval the initial value problem independently and finally fit the parts together: this runs as follows (in order to discern the solutions on the different subintervals, the solution is now denoted by y_[i](x;s) ).
  1. [a,b] = [x₁,x₂]U....U[x_m-1,x_m]
  2. solve y_[i](x;s_i)' = F(x,y_[i](x;s_i)), y_[i](x_i;s_i) = s_i
     for x_i <= x <= x_i+1 , i=1,...,m-1 .
     (Solve the ODE on subinterval i with initial value s_i)
  and require now continuity at the grid points and the boundary conditions:
  1. R(s₁,s_m) = 0 (original boundary conditions)
  2. y_[i](x_i+1;s_i ) = s_i+1 i=1,...m-2 (solution coming from the left at x_i+1 is identical to initial value used for the right interval, the continuity condition)
  3. s_m = y_[m-1](x_m;s_m-1 ) (definition of s_m )
  and this constitutes system of m nonlinear equations each of dimension n for the m × n unknowns s₁,...,s_m .
  
• This nonlinear system is solved using the damped Newton's method as it is used in the chapter dealing with nonlinear equations. The problem now is shifted to the appropriate choice of the initial values for Newton's method, which is also not an easy task. In real life this may be quite tricky as you might try using the examples prepared here. If these initial values are not appropriate, then either the Newton solver or the initial value solver will signal trouble and the process breaks down.
• The code we are using here is BVPSOL by P. Deuflhard, G. Bader and L. Weimann, available from the ELIB of ZIP. It uses the Gragg-Bulirsch-Stoer method as integrator. This is not suitable for stiff ordinary differential equations, hence be careful with the choice of your own testcases.
• Newton's method requires the computation of the Jacobian of the nonlinear system and this involves the computation of the derivatives of the solutions y_[i](x;s) with respect to s, the matrices (d/d s)y_[i](x;s) which in turn requires the solution of the so called variational equations
  (d/d s) y_[i](x;s)' = (d/dy) F(x,y_[i](x;s))(d/ds)y_[i](x;s) , (d/ds)y_[i](x;s) = I for x = x_i .
  BVPSOL avoids this and uses an internal numerical differentiation (of first order accuracy, using the square root of the integrators tolerance requirement as difference stepsize) for this, by solving n additional initial value problems with slightly perturbed initial value on each subinterval. This in turn requires a highly precise numerical integration of the ODE's. Hence one must be somewhat restraining in requiring final precision of the BVP solver (the parameter tol).

There are two different input forms: Predefined
and User defined

Input: Predefined cases

• There is a choice between 4 different problems.
• These problems are
  1. 1. A scalar linear problem of second order transformed to first order, i.e. n=2 with inhomogeneous Dirichlet conditions. The solution is unique. This is simple and we use m=2, i.e. simple shooting.
  2. 2. A nonlinear two point boundary value problem of second order, transformed to first order, uniquely solvable but with a critically moving singularity. We use 4 subintervals, m=5.
  3. 3. A nonlinear two point boundary value problem of second order with homogeneous Dirichlet boundary conditions which has two solutions. Which one is obtained depends on the initial value choice. The default obtains the smaller of the two.
  4. 4. The famous reentry-problem of an Apollo-like space vehicle as described in ''Introduction to Numerical analysis'' by J. Stoer and R. Bulirsch, also obtainable from the ELIB of ZIP. This problem has dimension 7 and we use 5 subintervals.
• You have the possibility to use the internally defined initial values or to choose the initial values and the x-grid yourself. The testcases 1,2 and 3 have separable linear boundary conditions: your input must satisfy these conditions.
• You have the choice of some internal parameters influencing the integrator and the desired final precision (in the nonlinear system). These are tol mentioned above already, the attainable order maxord of the integrator (even and in the range {2,..,16}), the maximal allowed stepsize hmax to be taken (which of course must be considerably smaller than b-a and which might internally reduced even more) , the maximum allowed number of steps itmax in Newton's method, the amount of output of the BVP-solver you want to see ( info ) and the components of y which should be shown graphically (index1, index2). For the cases 1, 2 and 3 the current solution is shown (red) together with the known analytic solution (green). The default values of these parameters are maxord =16 , itmax = 30 , info = -1 (no intermediate output), and

  case	| tol	| hmax	| index1	| index2	|
  1	| 1.0e-6	| 0.05	| 1	| 1	|
  2	| 1.0e-6	| 0.1	| 1	| 1	|
  3	| 1.0e-5	| 0.05	| 1	| 1	|
  4	| 1.0e-6	| 0.01	| 2	|3	|

• The required precision for the integrator is tol/100 by default.

Output: predefined case

• the final solution and the effort involved:
  1. the number of Newton steps
  2. the number of computed time steps of the integrator.
  3. the number of rejected time steps (local error too large, step size reduced).
  4. the number of evaluations of the ODE function
  5. computing time
• In addition you get a plot of:
  1. step sizes used (in log form) (last iteration).
  2. local error estimation(in log form) (last iteration).
  3. the actual order used (last iteration).
  4. the history of the development of the solution for the two selected components resp. the true and the computed solution.
• Should the integrator or the BVP-solver issue warning messages you get these as a protocol.

Questions: Predefined case?!

• How sensitive is the solver against low precision requirements and why?
• Behaves the local error of the integrator as expected (tol/100)?
• How is the relation between total effort and the increase of m?
• In which cases exhibits the method high sensitivity against bad initial guesses and in which not? why?

to the input form: predefined case

Input: user defined case

• Here you have the possibility to define a BVP yourself. This requires
  1. Specification of n. n is the number of components of y resp. y', R, F.
  2. You write a piece of code for computing the components of F and a second piece of code for computing the components of R . Details are given on the input form and in the help text.
  3. The number m of grid points x_i
  4. The grid points. In the above a = x₁, b = x_m then.
  5. The initial guess of the solution at these gridpoints, i.e. a list of n × m values, each block of n values representing one value of y.
• Then you specify the parameters tol, hmax, maxord, itmax, info explained above already.
• For graphical output you can specify two components for plotting.

Output: user defined case

• This is the same as for the predefined cases.

Questions: user defined case?!

• This case is meant for the solution of a BVP problem of your own interest.
• The questions to be asked here are of course the same as for the predefined cases.

to the input form: user defined case