Address issues raised in #54 (mainly documentation).

RK-coeff-opt/oc_albrecht.m

@@ -1,9 +1,9 @@
 function coneq=oc_albrecht(A,b,c,p)
 % function coneq=oc_albrecht(A,b,c,p)
 %
-% Order conditions for SSP RK Methods.
+% Order conditions for SSP RK methods.
 %
-% This version is based on Albrecht's approach
+% This version is based on Albrecht's approach.
 
 % ..warning:: 
 %

RK-coeff-opt/oc_butcher.m

@@ -1,7 +1,7 @@
 function coneq=oc_butcher(A,b,c,p)
 % function coneq=oc_butcher(A,b,c,p)
 %
-% Order conditions for RKMs.
+% Order conditions for Runge-Kutta methods.
 % This version is based on Butcher's approach.
 %
 % Assumes `p>1`.

RK-coeff-opt/rk_obj.m

@@ -1,17 +1,18 @@
 function [r,g]=rk_obj(x,class,s,p,objective)
 % function [r,g]=rk_obj(x,class,s,p,objective)
+%
 % Objective function for RK optimization.
 %
-% The meaning of the input arguments is as follow:
+% The meaning of the input arguments is as follows:
 %     * :math:`x`: vector of the unknowns.
 %     * class: class of method to search ('erk' = explicit RK; 'irk' = implicit RK; 'dirk' = diagonally implicit RK; 'sdirk' = singly diagonally implicit RK; '2S', '3S', '2S*', '3S*' = low-storage formulations).
 %     * :math:`s`:number of stages.
 %     * :math:`p`: order of the RK scheme.
 %     * objective: objective function ('ssp' = maximize SSP coefficient; 'acc' = minimize leading truncation error coefficient).
 %
-% The meaning of the output arguments is as follow:
+% The meaning of the output arguments is as follows:
 %     * r: it is a scalar containing the radius of absolute monotonicity if objective = 'ssp' or the value of the leading truncation error coefficient if objective = 'acc'.
-%     * g: it is a vector and contains the gradient of the objective function respect to the unknowns.  It is an array with all zero elements except for the last component which is equal to one if objective = 'ssp' or it is an empty array if objective = 'acc'. 
+%     * g: a vector containing the gradient of the objective function respect to the unknowns.  It is an array with all zero elements except for the last component which is equal to one if objective = 'ssp' or it is an empty array if objective = 'acc'. 
 
 
 if strcmp(objective,'ssp')

RK-coeff-opt/rk_opt.m

@@ -22,7 +22,7 @@
 %       Accuracy optimization is not currently supported for multistep RK methods
 %     * poly_coeff_ind: index of the polynomial coefficients to constrain (`\beta_j`) for `j > p`  (j denotes the index of the stage). The default value is an empty array.  Note that one should not include any indices `i \le p`, since those are determined by the order conditions.
 %     * poly_coeff_val: constrained values of the polynomial coefficients (`\beta_j`) for `j > p` (tall-tree elementary weights). The default value is an empty array.
-%     * startvec: vector of the initial guess ('random' = random approach; 'smart' = smart approach; alternatively, the user can provide the startvec array. By default startvec is initialize with random numbers.
+%     * startvec: vector of the initial guess ('random' = random approach; 'smart' = smart approach; alternatively, the user can provide the startvec array. By default startvec is initialized with random numbers.
 %     * solveorderconditions: if set to 1, solve the order conditions first before trying to optimize. The default value is 0.
 %     * np: number of processor to use. If np `> 1` the MATLAB global optimization toolbox *Multistart* is used. The default value is 1 (just one core).
 %     * num_starting_points: Number of starting points for the global optimization per processor. The default value is 10.
@@ -36,7 +36,7 @@
 %        **numerical experiments have shown that when the objective function is the minimization of the leading truncation error coefficient, the interior-point algorithm performs much better than the sqp one.**
 %
 %     * display: level of display of fmincon solver ('off', 'iter', 'notify' or 'final'). The default value is 'notify'.
-%     * problem_class: class of problems for which the RK is designed ('linear' or 'nonlinear' problems). This option changes the type of order conditions check, i.e. linear or nonlinear order conditions controll. The default value is 'nonlinear'.
+%     * problem_class: class of problems for which the RK is designed ('linear' or 'nonlinear' problems). This option changes the type of order conditions check, i.e. linear or nonlinear order conditions control. The default value is 'nonlinear'.
 %
 %
 %     .. note::

RK-coeff-opt/write_field.m

@@ -1,6 +1,8 @@
-function writeField(writeFid,name,value)
-% function writeField(writeFid,name,value)
+function write_field(writeFid,name,value)
+% function write_field(writeFid,name,value)
+%
+%
+% Utility function to write a single parameter and value.
 
 fprintf(writeFid,'\n%s\n',name);
 fprintf(writeFid, [repmat('%5.16E\t', 1, size(value,2)),'\n'], value'); 
-

RKtools/README.rst

@@ -1 +1,8 @@
 Some general utilities for analyzing Runge-Kutta methods.
+
+Some of the routines expect as input a structured array `rk`.
+This structure must have the fields `A, b, c`, containing its
+Butcher coefficients.  Optionally, it may represent an additive
+Runge-Kutta method or an embedded pair in which case it should also have
+`Ahat`, `bhat`, `chat` containing the coefficients of the secondary
+method.

RKtools/am_radius.m

@@ -5,23 +5,21 @@
 % of a Runge-Kutta method, given the Butcher array.
 %
 % For an `m`-stage method, `A` should be an `m \times m` matrix
-% and `b` should be a column vector of length m.
+% and `b` and `c` should be column vectors of length m.
 %
-% Accuracy can be changed by modifying the value of eps (default `10^-{10}`)
+% Accuracy can be changed by modifying the value of eps (default `10^{-10}`)
 % Methods with very large radii of a.m. (>50) will require
 % the default value of rmax to be increased.
 %
 % The radius of absolute monotonicity is the largest value of `r`
 % such that
 %
-% .. raw:: latex
+%    \\begin{align*}
+%    K(I+rA)^{-1} \\ge & 0    \\\\\\
+%    rK(I+rA)^{-1}e_m \\le & e_{m+1}
+%    \\end{align*}
 %
-%    \begin{eqnarray}
-%    K(I+rA)^{-1} \ge & 0    \\
-%    rK(I+rA)^{-1}e_m \le & e_{m+1}
-%    \end{eqnarray}
-%
-%    where $$ K = \left(\begin{array}{c} A \\ b^T \end{array}\right) $$
+%    where $$ K = \\left(\\begin{array}{c} A \\\\\\ b^T \\end{array}\\right) $$
 
 
 if nargin<5 rmax=50; end

am_radius-opt/README.rst

@@ -3,6 +3,6 @@ This includes codes for optimizing stability functions of
 multistep, multistage methods and even methods with downwinding.
 
 Generally, the optimization problem is phrased as a sequence of linear 
-programming feasibility problems.  For details, see [ketcheson2009]_.
+programming feasibility problems.  For details, see :cite:`2009_monotonicity`.
 
 The optimization of rational functions is experimental.

am_radius-opt/Rkp.m

@@ -3,12 +3,13 @@
 %
 % Find the optimal SSP k-step explicit LMM with order of accuracy p.
 %
-% Inputs: 
-%       * `k` = # of steps
-%       * `p` = order of accuracy
+% Inputs:
+%       * k = # of steps
+%       * p = order of accuracy
 %
-% Outputs: 
-%       * `\alpha, \beta` = the coefficients of the method
+% Outputs:
+%       * R = the SSP coefficient
+%       * alpha, beta = the coefficients of the method
 %
 % Requires MATLAB's optimization toolbox for the LP solver.
 

am_radius-opt/Rkp_dw.m

@@ -4,13 +4,16 @@
 % Finds the optimal SSP k-step explicit LMM with order of accuracy p
 % allowing downwind operators
 %
-% Inputs: k = # of steps
-%         p = order of accuracy
+% Inputs:
+%       * k = # of steps
+%       * p = order of accuracy
 %
-% Outputs: alpha, beta, tbeta = the coefficients of the method
+% Outputs:
+%       * R = the SSP coefficient
+%       * alpha, beta, tbeta = the coefficients of the method
 %
 % The method is given by
-% `u_n = \sum_{j=0}^{k-1} (\alpha[j] + \beta[j] F(u_{n-k+j} + tbeta[j] tF(u_{n-k+j}))`
+% `u_n = \sum_{j=0}^{k-1} (\alpha[j] + \beta[j] F(u_{n-k+j}) + tbeta[j] tF(u_{n-k+j}))`
 % where tF(u) is the negated downwind operator.
 %
 % Depends on MATLAB's optimization toolbox for the LP solver

am_radius-opt/Rkp_imp.m

@@ -3,11 +3,13 @@
 %
 % Find the optimal SSP k-step implicit LMM with order of accuracy p
 %
-% Inputs: 
+% Inputs:
 %       * k = # of steps
 %       * p = order of accuracy
 %
-% Outputs: alpha, beta = the coefficients of the method
+% Outputs:
+%       * R = the SSP coefficient
+%       * alpha, beta = the coefficients of the method
 %
 % Depends on MATLAB's optimization toolbox for the LP solver
 

am_radius-opt/Rkp_imp_dw.m

@@ -4,10 +4,13 @@
 % Finds the optimal k-step implicit LMM with order of accuracy p
 % allowing downwinding
 %
-% Inputs: k = # of steps
-%        p = order of accuracy
+% Inputs:
+%       * k = # of steps
+%       * p = order of accuracy
 %
-% Outputs: alpha, beta, tbeta = the coefficients of the method
+% Outputs:
+%       * R = the SSP coefficient
+%       * alpha, beta, tbeta = the coefficients of the method
 %    
 % Depends on MATLAB's optimization toolbox for the LP solver
 

am_radius-opt/Rskp.m

@@ -4,21 +4,21 @@
 % Finds the optimal contractive k-step, s-stage GLM with order of accuracy p
 % for linear problems
 %
-% Inputs: s = # of stages
-%         k = # of steps
-%         p = order of accuracy
+% Inputs:
+%       * s = # of stages
+%       * k = # of steps
+%       * p = order of accuracy
 %
-% Outputs: 
-%        R = threshold factor
-%        gamma = coefficients of the polynomials
+% Outputs:
+%       * R = threshold factor
+%       * gamma = coefficients of the polynomials
 %         
 %        for k=1, the resulting polynomial is
 %        `\sum_{j=0}^m (1+z/R)^j`
 %
-%        in general, the resulting stability function is
-%        (Fill in)
+%        For details on the general case, see :cite:`2009_monotonicity`.
 %
-%Depends on MATLAB's optimization toolbox for the LP solver
+% This routine requires MATLAB's optimization toolbox for the LP solver.
 
 %=========================================================
 %Initialize

am_radius-opt/multi_R_opt.m

@@ -3,28 +3,21 @@
 %
 % This function is a script to run the routines Rskp, Rkp_dw, Rkp_imp, or
 % Rkp_imp_dw several times with different inputs, in order to construct tables
-% of optimal values like those that appear in [ketcheson2009]_.
-% different values of the input parameters, i.e.: 
-%
-% k = [k1, k2, ..., kK]^T, K = length(k),  ith-element = # of steps
-% p = [p1, p2, ..., pP]^T, P = length(p),  ith-element = order of accuracy
+% of optimal values like those that appear in :cite:`2009_monotonicity`.
 % 
-% and 
-%
-% s = [s1, s2, ..., sS]^T, S = length(s),  ith-element = # of stages
-%
-% when optimal contractive k-step, s-stage GLM are investigated.
-%
-% The family of method to be considered is specified in the string 'class'.
+% The inputs k, p, and (optionally) s should be vectors containing
+% the numbers of steps, orders of accuracy, and numbers of stages
+% to be considered, respectively.  The output includes results for
+% all combinations of values from the input vectors.
+% 
+% The family of methods to be considered is specified in the string 'class'.
+% Valid values are:
 %
-% Note that in general `S\ne K\ne P`. Fixed the order of accuracy of the time 
-%    integration scheme, one is usually interested in understanding the
-%    behavior of the threshold factor R as a function of the number of
-%    stages. Therefore, for a fixed element of the array "p", this function
-%    loops over the elements of the array "s". Thus, min(s) => max(p). The
-%    equality holds for any order of accuracy because the number of 
-%    linear order conditions that will be imposed to construct the 
-%    GLM coefficients is p. 
+%  * 'skp': find optimal general linear methods (multistep, multistage).
+%           In this case the vector s must be included in the inputs.
+%  * 'kp_imp': find optimal implicit linear multistep methods.
+%  * 'kp_dw':  find optimal explicit downwind linear multistep methods.
+%  * 'kp_imp_dw':  find optimal implicit downwind linear multistep methods.
 
 
 % Parse optional input arguments

am_radius-opt/radimpfast.m

@@ -1,12 +1,14 @@
 function rad=radimpfast(p,q)
 % function rad=radimpfast(p,q)
 %
-% Compute the radius of absolute monotonicity of a rational function.
+% Compute the radius of absolute monotonicity of the rational function
+% whose numerator has coefficients p and denominator has coefficients q.
+% The coefficients are ordered in ascending powers.
 %
 % This function is outdated and needs to be fixed.
 %
-% Uses van de Griend's algorithm, assuming multiplicity one for all roots.
-% Uses high precision arithmetic.
+% Uses van de Griend's algorithm :cite:`vandegriend1986`, assuming multiplicity
+% one for all roots.  Uses high precision arithmetic.
 
 syms x P Q phi