%
% BVP2SH.M This function file solves a 2nd order BVP with linear BCs
% via the Shooting Method
%
% This method solves the BVP using an IVP solver with an initial guess
% for the unknown initial conditions. An iterative scheme automatically
% modifies this guess until the boundary conditions at the second boundary
% point are satisfied. This iterative technique is called the Shooting Method.
%
% A general 2nd order ODE can be written in the form
% y'' + a*y' + b*y = f(t)
% where, in general, a and b can be functions of y and t.
%
% We can also write this as a set of two 1st order equations in so-called
% state space form as
% dz1/dt = z2 and dz2/dt = f - a*z2 - b*z1
% where z1 = y and z2 = y'
%
% Now, if we had a set of initial conditions, zo = [y(to); y'(to)], then
% this state space system could be easily solved as an IVP.
%
% However, since this is a BVP, we do not have a complete set of initial
% conditions. Thus, the basic idea of the shooting method is to guess at the
% unknown variable at the initial state so that the solution of the resulting
% IVP gives the desired final state for the BVP. This is done iteratively by
% converging on the free parameter, alf, that is included as part of the
% initial condition vector such that the desired final state is achieved.
%
% For the case of a 2nd order 2-point BVP with linear BCs, a general way to
% write the boundary conditions is as follows:
% at p1: a1*y(p1) + b1*y'(p1) - w1 = 0
% at p2: a2*y(p2) + b2*y'(p2) - w2 = 0
% The coefficients in these expressions are input to this function in the
% zbc matrix as: zbc = [a1 b1 w1; a2 b2 w2]. This format essentially
% defines the BCs for the 2nd order BVP with linear BCs (note that the equation
% can be nonlinear -- but not the BCs).
%
% At the initial boundary point, p1, a check is made to determine the
% kind of boundary condition and thus the form of the initial condition
% vector for the iterative IVP solution:
% Dirichlet BC: if b1 = 0, y(p1) = w1/a1
% then zo = [w1/a1; alf] (guess initial slope)
% Neumann BC: if a1 = 0, y'(p1) = w1/b1
% then zo = [alf; w1/b1] (guess initial value)
% Robin (mixed) BC: if both a1 and b1 are nonzero,
% then zo = [alf; (w1 - a1*alf)/b1] (guess initial value)
% where alf is the unknown parameter which will be computed as part
% of the iterative solution. alf is determined such that the second
% BC is satisfied to within some small tolerance. Note that the Neumann
% BC is just a subset of the Robin BC with a1 = 0, and it is implemented as
% such in this routine.
%
% The inputs to this function are similar to those used with Matlab's ode23
% IVP solver (in its simplest implementation) except that the zbc matrix
% replaces the initial state vector. With zbc known, the zo vector is computed
% internally as described above. Also, there is an option to input a guess
% for alf, if desired
% func - name of function file that evaluates the RHS of the state equations
% tspan - vector that defines the range for the independent variable
% zbc - matrix containing the BC coefficients (as described above)
% options - standard options structure for use in ode23 (optional)
% alfo - inital guess for alf (optional)
%
% The outputs are:
% t - column vector containing the independent variable
% z - solution matrix with two columns with column 1 containing z1(t)
% and column 2 containing z2(t)
%
% Note #1: Extra parameters that may be needed in the ode function file must be
% passed to the function file via a global statement in the main program and
% the function file. The current, relatively simple implementation of the
% Shooting Method does not use the variable input arguments (varargin) that are
% often used in the ode23 routine. This is not overly restrictive, however, since
% it is just as easy to pass needed variables into the function via the global
% statement.
%
% Note #2: This routine assumes that you are integrating from boundary point
% p1 to boundary point p2. This normally implies that p2 > p1, but this is not
% absolutely essential. In some cases, it may be appropriate to integrate
% 'backwards' -- just be careful with the ordering of the BCs if you try this...
%
% Note #3: Nonlinear problems are sometimes sensitive to the initial guess for
% the free variable, alf, in the initial condition vector. Here we have simply
% set alfo = 1 as default. The user can override this as needed in particularly
% difficult problems. A poor initial guess usually just means that a few more
% iterations may be needed. If the system does not converge within 10 - 20
% iterations, you may want to change alfo via the last entry in the input argument
% list to see if this helps. This is not usually necessary, however...
%
% File prepared by J. R. White, UMass-Lowell (March 2003)
%
function [t,z] = bvp2sh(func,tspan,zbc,options,alfo)
%
% check to see if the options structure is available
if nargin < 4, options = []; end
%
% check to see if there is an initial guess for alfo
if nargin < 5, alfo = 1; end
%
% extract BC coefficients (for convenience)
a1 = zbc(1,1); b1 = zbc(1,2); w1 = zbc(1,3);
a2 = zbc(2,1); b2 = zbc(2,2); w2 = zbc(2,3);
%
% set some other iterative parameters
err = 1e10; tol = 1e-6; icnt = 1; icntmax = 25;
%
% start iteration
while abs(err) > tol & icnt <= icntmax
%
% determine initial condition vector based on structure of zbc matrix
if b1 == 0
zo = [w1/a1; alfo];
else
zo = [alfo; (w1 - a1*alfo)/b1];
end
%
% solve IVP with current guess for initial conditions
[t,z] = ode23(func,tspan,zo,options);
%
% evaluate error in second BC
e1 = a2*z(end,1) + b2*z(end,2) - w2; err = e1;
fprintf(1,' \n')
fprintf(1,'For iteration # %2i: \n',icnt)
fprintf(1,' Initial conditions are: %13.5e %13.5e \n',zo)
fprintf(1,' Error in 2nd BC is: %13.5e \n',err)
fprintf(1,' \n')
%
% estimate new values for alf (if necessary)
if abs(err) > tol
alfp = 1.01*alfo;
if b1 == 0
zo = [w1/a1; alfp];
else
zo = [alfp; (w1 - a1*alfp)/b1];
end
[t,z] = ode23(func,tspan,zo,options);
e2 = a2*z(end,1) + b2*z(end,2) - w2;
deda = (e2-e1)/(0.01*alfo); alfn = alfo-e1/deda;
icnt = icnt+1; alfo = alfn;
end
end % end of iteration loop
if icnt >= icntmax
fprintf(1,' WARNING -- Hit maximum iteration limit!!! \n');
end
%
% end of function