Effect of the Guess Function & Continuation Method on the Run Time of MATLAB BVP Solvers

A well-known statement says that the PID controller is the â€œbread and butterâ€ ​ of the control engineer. This is indeed true, from a scientific standpoint. However, nowadays, in the era of computer science, when the paper and pencil have been replaced by the keyboard and the display of computers, one may equally say that MATLAB is the â€œbreadâ€ ​ in the above statement. MATLAB has became a de facto tool for the modern system engineer. This book is written for both engineering students, as well as for practicing engineers. The wide range of applications in which MATLAB is the working framework, shows that it is a powerful, comprehensive and easy-to-use environment for performing technical computations. The book includes various excellent applications in which MATLAB is employed: from pure algebraic computations to data acquisition in real-life experiments, from control strategies to image processing algorithms, from graphical user interface design for educational purposes to Simulink embedded systems

but also whether any solution is achieved or not depends strongly on the initial guess.Therefore, depending of the guess function, BVPs may have no solution or a single solution, or multiple solutions.Moreover, the quality of the initial guess can be critical to the solver performance, which reduces or augments the run time.However, coming up with a sufficiently good guess can be the most challenging part of solving a BVP.Certainly, the user should apply the knowledge of the problem's physical origin (MATLAB Documentation).In MATLAB, when solving BVPs the user must provide a guess to assist the solver in computing the desired solution (Kierzenka & Shampine, 2001).MATLAB BVP solvers call for users to provide guesses for the mesh and solution.Although MATLAB BVP solvers take an unusual approach to the control of error in case of having poor guesses for the mesh and solution, especially for the nonlinear BVP, a good guess is necessary to obtain convergence (Shampine et al., 2003).Whatever intuitive guess values/functions are imposed, eventually the BVP solver fails for some parameters or for some lengths.If any guess values works for the range of length, the rest of the length may be extended using continuation.The method of continuation exploits the fact that the solution obtained for one input will serve as the initial guess for the next value tried.In case of any difficulty in finding a guess for the interval of interest, generally it will be easier to solve the problem on a shorter interval.Then the solution of the sequence of BVPs on the shorter interval will be used as a guess for the next section.With modest increases in the interval, this will continue until the interval of interest is spanned (Shampine et al., 2003).The cost of the continuation method is the increased run time.How the guess value good is, the less computation time it takes with the continuation method.This is due the fact that, the remaining length depends of the convergence length (based on the guess value) which its higher value reduces the computation time.

Initial setup
The first step in solving a problem is defining it in a way the software can understand.The bvp4c framework uses a number of subfunctions which make it as simple as possible for the user to enter the ODE function, initial data and parameters for a given problem.By way of the following example we see exactly how a problem is supplied and solved by bvp4c.For the evaluation of the guess value /function, the steady-state Brillouin equation is exploited.The coupled ODEs for the evolution of the intensities of pump I p and Stokes I s can be written as (Agrawal, 2001), where 0 ≤ z ≤ L is the propagation distance along the optical fiber of the total length L, α is the fiber loss coefficient, g B is the Brillouin gain coefficient, respectively.Here, it is assumed that, Stokes wave is launched from the rear end of the fiber.Then the known values of the input pump power I p (0) and the Stokes wave power I s (L) are referred as the boundary values.
The first task is to define the ODEs in MATLAB as a function to return these equations.
Similarly the user then rewrites the boundary conditions to correspond to this form of the problem.We may code the ODEs for scalar evaluation and boundary conditions, respectively as, The next step is to create an initial guess for the form of the solution using a specific MATLAB subroutine called bvpinit .The user passes a vector x and an initial guess on this mesh in the form bvpinit (x, Yinit), which is then converted into a structure useable by bvp4c.Aside from a sensible guess being necessary for a convergent solution the mesh vector passed to bvpinit will also define the boundary points of the problem, i.e.
The initial guess for the solution may take one of two forms.One option is a vector where Yinit(i) is a constant guess for the i-th component y(i,:) of the solution at all the mesh points in x.The other is as a function of a scalar x, for example bvpinit(x,@yfun) where for any x in [a, b], yfun(x) returns a guess for the solution y(x).It must be pointed out that even when a function is supplied that approximates the solution everywhere in the interval; the solver uses its values only on the initial mesh.The guess can be coded as a function of a scalar x as, The next subroutine to look at is bvpset, that specifies which options bvp4c should be use in solving it.The function is called options = bvpset('name1',value1,...) and since MATLAB documentation gives an in depth account of each of the options only a brief outline of those notable is given here (Hale, 2006).options = []; % default %options = bvpset ('Stats','on','RelTol','abstol',; %options = bvpset(options,'Vectorized','on'); %options = bvpset(options,'FJacobian',@odeJac); %options = bvpset(options,'BCJacobian',@bcJac); RelTol -Relative tolerance for the residual [ positive scalar 1e-3 ] The computed solution S(x) is the exact solution of S ' (x) = F(x, S(x)) + res(x).On each subinterval of the mesh, component i of the residual must satisfy norm AbsTol -Absolute tolerance for the residual [positive scalar or vector 1e-6] Elements of a vector of tolerances apply to corresponding components of the residual vector.AbsTol defaults to 1e-6.FJacobian \ BCJacobian -Analytical partial derivatives of ODEFUN \ BCFUN Computation of the Jacobian matrix at each mesh point can be a very expensive process.By passing an analytic derivative of the ODE and BC functions the user can greatly reduce computational time.For example when solving y' = f(x, y), setting FJacobian to @FJAC where ∂f/∂y = FJAC(x, y) evaluates the Jacobian of f with respect to y. Stats -Display computational cost statistics [ onoff ] Vectorized -Vectorized ODE function [ on -off ] As will be discussed in section 6, bvp4c is able to accept a vectorised function which can markedly increase the efficiency of calculating local Jacobians over using finite differences with the odenumjac subroutine.Hence in the following programs, we will define options = bvpset ('Stats','on','RelTol','abstol',; solinit = bvpinit(linspace(0,L,2), @guess); And call the bvp4c routine with: sol = bvp4c(@ode,@bc,solinit,options); The above essentially ends the user input in solving the BVP system and the rest is left to bvp4c.Within the framework there are several notable steps which should be expounded.

Derivation of the guess
In this chapter, four guess functions are derived for the assistance of the MATLAB BVP solvers with the help of MATLAB symbolic toolbox.

4 th Guess
Highly intuitive guess function may be derived the using the solution of lossless system, i.e., with eliminating the α coefficient in the Eq (3) and Eq. ( 4), With neglecting the attenuation coefficient, the solution of the Eq.( 5) and Eq.( 6) is found as (Kobyakov et al., 2006), where, where, Exploiting the solution of Eq.( 7) and Eq. ( 8), general expression of P p (z) and P S (z) can be derived as, If a→0, then () () We are interested in the behavior as z → 0 and so, the higher the power of x, the less effect it has in these expansions.Our goal is to satisfy the equations as well as possible, so we want to choose coefficients that make as many successive terms zero as possible, starting with the lowest power.To eliminate the constant terms, we see from the expansions that we must take  As can be seen from Table 1 and Fig. 4 the best estimation is the 4 th guess.Because its' convergence length (30000) is more than the others (15000, 9000, 8000, respectively).The performance of the 2 nd guess is approximately same as the first one.Because it hardly converge the solution using 453 points at 9000 meter.However, its performance is same as the first one with 40 points at 8000 meter.

Continuation
The method of continuation exploits the fact that generally the solution of one BVP is a good guess for the solution of another with slightly different parameters.If you have difficulty in finding a guess for the solution that is good enough to achieve convergence for the interval of interest, it is frequently the case that the problem is easier to solve on a shorter interval.
The idea is then to solve a sequence of BVPs with the solution on one interval being used as a guess for the problem posed on a longer interval.Of course, this technique cannot be used to extend the interval ad infinitum; no matter how good the guess is, eventually the solver will not be able to distinguish the different kinds of solutions (Shampine et al., 2003).

Effect of the step size on the run time
Using the snippets 5 and the 1 st Guess of Table 1, the performance of the continuation over step size is illustrated in Table 2 Warning: Unable to meet the tolerance without using more than 5000 mesh points.
The last mesh of 4106 points and the solution are available in the output argument.
The maximum residual is 0.00018931, while requested accuracy is 1e-005.
As can be seen on Table 2, for some step size over the modest increment, computation blows up (i.e.600 m).When the step size is between 50 and 400, the number of used mesh is slightly different from each other.However, when it is 500 m, abruptly increase in the number of used mesh is a sign of lack of confidence.In Fig. 5, it can be seen that the distance between some mesh points, especially near the boundaries are denser than the others.This is because the solver tries to control the residual of the interpolating polynomial: r(x) = S'(x) -f(x,S(x)).The behavior of this residual depends on the behavior of some high derivatives of the solution (that the solver does not have access to).In the solver, the residual is estimated at each mesh subinterval, and additional mesh points are introduced if the estimate is bigger than the tolerance.The mesh selection algorithm is 'localized', which means that if the residual is just above the tolerance, the interval will be divided into two (and likely on each of those subintervals, the residual will be much smaller than the tolerance).Also, the algorithm for removing mesh points is quite conservative, so there could be regions where the residual will be quite a bit smaller that the tolerance (i.e., the mesh could be quite a bit denser than necessary) ( Kierzenka & Shampine, 2001).

Effect of the constructing an initial guess with bvpinit function on the run time
As can be seen from Table 3, when constructing an initial guess using bvpinit(linspace(0,L,N), starting with a mesh of 5-10 points could often result in a more efficient run.It must be pointed out that with adaptive collocation solvers, using that many points (N=50,100) with a poor guess could often be counterproductive.In the case of N=100, the solver still achieved the sufficient accuracy as it is between 2 and 10.

Speeding up the run time of BVP solvers
The first technique which is used to reduce run time is vectorizing the evaluation of the differential equations.Vectorization is a valuable tool for speeding up MATLAB programs and this greatly reduces the run time (Shampine et al., 2003).By vectorization, the function f(x,y) is coded so that when given a vector x=[x 1 ,x 2 ,...] and a corresponding array of column vectors y=[y 1 ,y 2 ,...], it returns an array of column vectors [f(x 1 ,y 1 ),f(x 2 ,y 2 ),...]).By default, bvp4c and bvp4c approximate a Jacobian using finite differences.The evaluation of the ODEs is vectorized by changing the vectors to arrays and changing the multiplication to an array multiplication.It can be coded by changing scalar quantities like y(1) into arrays like y(1,:) and changing from scalar operations to array operations by replacing * and ˆ with .*and .ˆ,respectively.When vectorizing the ODEs, the solver must be informed about the presence of vectorization by means of the option 'Vectorized','on'.options = bvpset('Stats','on','RelTol',1e-3,'Vectorized','on'); The second technique is that of supplying analytical partial derivatives or to supply a function for evaluating the Jacobian matrix.This is because, in general, BVPs are solved much faster with analytical partial derivatives.However, this is not an easy task since it is too much trouble and inconvenient, although MATLAB Symbolic Toolbox can be exploited when obtaining analytical Jacobians.The third technique is to supply analytical partial derivatives for the boundary conditions.However, it has less effect on the computation time compared with supplying analytical Jacobians and vectorization.The solver permits the user to supply as much information as possible.It must be emphasized that supplying more information for the solvers results in a shorter computation run time (Gokhan & Yilmaz, 2011b).
The set of equations ( 3) and ( 4) is vectorized by changing the vectors to arrays and changing the multiplication to an array multiplication as seen below, The code bvp4c permits you to supply analytical partial derivatives for either the ODEs or the boundary conditions or both.It is far more important to provide partial derivatives for the ODEs than the boundary conditions.The solver is informed that a function is written for evaluating ∂f/∂y by providing its handle as the value of the FJacobian option.Similarly, the solver can be informed of a function for evaluating analytical partial derivatives of the boundary conditions with the option BCJacobian (Shampine et al., 2003).FJacobian and BCJacobian can be introduced as with the below codes, %options = bvpset(options,'FJacobian',@odeJac); %options = bvpset(options,'BCJacobian',@bcJac); The MATLAB Symbolic Toolbox has a function jacobian that can be very helpful when working out partial derivatives for complicated functions.Its use is illustrated with a script for the partial derivatives of the ODEs of this example.The performance of the insertion of analytical partial derivatives and vectorization, and both are compared in Table 5.As can be seen from If the comparison among two solvers has made, it could be expressed that bvp5c "looks" exactly like bvp4c.However, bvp5c controls scaled residual and true error but bvp4c controls residual in a different norm.And, bvp5c is more efficient at stringent tolerances.
Also, bvp5c solves singular BVPs, but not multipoint BVPs.Moreover, bvp5c handles unknown parameters in a different way.And also, bvp5c was added to MATLAB at R2007b (Shampine, 2008)

Conclusion
Within the chapter, in order to analyze the effect of guess functions on the computation time, four guess functions are derived.For better understanding, while exploiting physical origin, guess functions are derived with the help of MATLAB Symbolic toolbox.Continuation method with functional snippets is presented to cope with poor guesses.Effect of the step size and bvpinit function on the computation time is analyzed.Speeding up the run time with vectorization and analytical partial derivatives are discussed and the comparison between bvp4c and bvp5c has been made.As a conclusion, it is illustrated that, intuitive guess values/functions improves the convergence length, leads the computation with fewer mesh points and consequently lessens the computation time.On the other hand, regarding with the continuation, adjusting the step size is important for the reduction of run time.It is illustrated that, over the modest step size, the solver fails and below the optimum step size, the computation time is increased.Moreover, it is showed that when constructing an initial guess using bvpinit(linspace(0,L,N), starting with a mesh of 5-10 points could often result in a more efficient run.Another outcome of the chapter is the illustration of the efficiency of the vectorization and analytical partial derivatives.It is showed specifically with an example and with bvp4c that, with the application of vectorization and analytical partial derivatives, the computation time is reduced approximately 15 %.The performance of the bvp4c and bvp5c is also compared.In terms of scalar evaluation, the performance of bvp5c solver is better than bvpc4 and it is evident as the computation length is increased.Compared with the scalar evaluation, for the bvp5c, only with vectorization and only with analytical partial derivatives this improvement is 8% and 13 % respectively.If both is used this improvement is about 24 %.

Table 1 .
Guess values/functions versus convergence length/mesh and run time.

Table 2 .
. The performance of continuation method versus step size If the step size is increased over the 600 meter, the computation fails with the below warning message;

Table 3 .
Performance of equally spaced N points for the mesh of a guess

Table 4 .
Comparison of the computation time with scalar evaluation and with vectorization

Table 5 .
Comparison of the computation time of bvp4c with vectorization, with analytical partial derivatives and with bothIn Table6, the performance of the bvp5c is illustrated.In terms of scalar evaluation, the performance of bvp5c solver is better than bvpc4 and it is evident as the length is increased.This improvement is about 47 % at 40 km.As in the case of bvp4c, the performance can be increased with vectorization and analytical partial derivatives or with both.Compared with the scalar evaluation, only with vectorization and only with analytical partial derivatives this improvement is 8% and 13 %, respectively.If both is used this improvement is about 24 %.

Table 6 .
Comparison of the computation time of bvp5c with vectorization, with analytical partial derivatives and with both