<?xml version="1.0" encoding="ISO-8859-1" standalone="no"?>
<!DOCTYPE MAN SYSTEM "../../manrev.dtd">
<MAN>
  <LANGUAGE>eng</LANGUAGE>
  <TITLE>bvode</TITLE>
  <TYPE>Scilab Function</TYPE>
  <DATE>June 1993</DATE>
  <SHORT_DESCRIPTION name="bvode"> boundary value problems for ODE</SHORT_DESCRIPTION>
  <CALLING_SEQUENCE>
    <CALLING_SEQUENCE_ITEM>   [z]=bvode(points,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...  </CALLING_SEQUENCE_ITEM>
    <CALLING_SEQUENCE_ITEM> fsub1,dfsub1,gsub1,dgsub1,guess1)  </CALLING_SEQUENCE_ITEM>
  </CALLING_SEQUENCE>
  <PARAM>
    <PARAM_INDENT>
      <PARAM_ITEM>
        <PARAM_NAME>z</PARAM_NAME>
        <PARAM_DESCRIPTION>
          <SP>The solution of the ode evaluated on the mesh given by points</SP>
        </PARAM_DESCRIPTION>
      </PARAM_ITEM>
      <PARAM_ITEM>
        <PARAM_NAME>points</PARAM_NAME>
        <PARAM_DESCRIPTION>
          <SP>an array which gives the points for which we want the solution</SP>
        </PARAM_DESCRIPTION>
      </PARAM_ITEM>
      <PARAM_ITEM>
        <PARAM_NAME>ncomp</PARAM_NAME>
        <PARAM_DESCRIPTION>
          <SP>number of differential equations   (ncomp &lt;= 20)</SP>
        </PARAM_DESCRIPTION>
      </PARAM_ITEM>
      <PARAM_ITEM>
        <PARAM_NAME>m</PARAM_NAME>
        <PARAM_DESCRIPTION>
          <SP>a vector of size <VERB>ncomp</VERB>. <VERB>m(j)</VERB> gives the  order of the j-th differential equation</SP>
        </PARAM_DESCRIPTION>
      </PARAM_ITEM>
      <PARAM_ITEM>
        <PARAM_NAME>aleft</PARAM_NAME>
        <PARAM_DESCRIPTION>
          <SP>left end of interval</SP>
        </PARAM_DESCRIPTION>
      </PARAM_ITEM>
      <PARAM_ITEM>
        <PARAM_NAME>aright</PARAM_NAME>
        <PARAM_DESCRIPTION>
          <SP>right end of interval</SP>
        </PARAM_DESCRIPTION>
      </PARAM_ITEM>
      <PARAM_ITEM>
        <PARAM_NAME>zeta</PARAM_NAME>
        <PARAM_DESCRIPTION>
          <SP><VERB>zeta(j)</VERB> gives j-th side condition point (boundary point). must have</SP>
          <P>
            <VERB>zeta(j) &lt;= zeta(j+1)</VERB>
          </P>
          <P>
    all side condition points must be mesh points in all meshes used, see description of <VERB>ipar(11)</VERB> and <VERB>fixpnt</VERB> below.
  </P>
        </PARAM_DESCRIPTION>
      </PARAM_ITEM>
      <PARAM_ITEM>
        <PARAM_NAME>ipar</PARAM_NAME>
        <PARAM_DESCRIPTION>
          <SP>an integer array dimensioned at least 11. a list of the parameters in <VERB>ipar</VERB> and their meaning follows some parameters are renamed in bvode; these new names are given in parentheses.</SP>
          <PARAM_INDENT>
            <PARAM_ITEM>
              <PARAM_NAME>ipar(1)  </PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>0  if the problem is linear, 1 if the problem is nonlinear</SP>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
            <PARAM_ITEM>
              <PARAM_NAME>ipar(2)  </PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>= number of collocation points per subinterval  (= k ) where</SP>
                <P>
                  <VERB>max m(i) &lt;=  k &lt;= 7 .</VERB>
                </P>
                <P>
    if <VERB>ipar(2)=0</VERB> then bvode sets  
  </P>
                <P>
                  <VERB>k = max ( max m(i)+1, 5-max m(i) )</VERB>
                </P>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>

            <PARAM_ITEM>
              <PARAM_NAME>ipar(3)  </PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>= number of subintervals in the initial mesh  ( = n ). if <VERB>ipar(3) = 0</VERB> then bvode arbitrarily sets <VERB>n = 5</VERB>.</SP>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
            <PARAM_ITEM>
              <PARAM_NAME>ipar(4)</PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>= number of solution and derivative tolerances.  ( = ntol ) we require</SP>
                <P>
                  <VERB>0 &lt; ntol &lt;= mstar.</VERB>
                </P>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>

            <PARAM_ITEM>
              <PARAM_NAME>ipar(5)</PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>= dimension of <VERB>fspace</VERB> ( = ndimf ) a
                real work array. its size provides a constraint on
                nmax. choose ipar(5) according to the formula:</SP>
                <P> <VERB>ipar(5)>=nmax*nsizef</VERB> </P>
                <P>where</P>
                <P> <VERB>nsizef=4+3*mstar+(5+kd)*kdm+(2*mstar-nrec)*2*mstar</VERB>.</P>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>

            <PARAM_ITEM>
              <PARAM_NAME>ipar(6)</PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>= dimension of ispace ( = ndimi )an integer work
                 array. its size provides a constraint on nmax, the
                 maximum number of subintervals. choose
                 <VERB>ipar(6)</VERB> according to the formula:</SP>
                 <P> <VERB>ipar(6)>=nmax*nsizei</VERB></P>
                  <P>where</P>
                 <P><VERB>nsizei=3 + kdm</VERB> with
                 <VERB>kdm=kd+mstar</VERB> ; <VERB>kd=k*ncomp </VERB>;
                 <VERB>nrec</VERB>=number of right end boundary
                 conditions.</P>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
            <PARAM_ITEM>
              <PARAM_NAME>ipar(7) </PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>output control ( = iprint )</SP>
                <PARAM_INDENT>
                  <PARAM_ITEM>
                    <PARAM_NAME>= -1</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>for full diagnostic printout</SP>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                  <PARAM_ITEM>
                    <PARAM_NAME>= 0</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>for selected printout</SP>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                  <PARAM_ITEM>
                    <PARAM_NAME>= 1</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>for no printout</SP>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                </PARAM_INDENT>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
            <PARAM_ITEM>
              <PARAM_NAME>ipar(8)  </PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>( = iread )</SP>
                <PARAM_INDENT>
                  <PARAM_ITEM>
                    <PARAM_NAME>= 0</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>causes bvode to generate a uniform initial mesh.</SP>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                  <PARAM_ITEM>
                    <PARAM_NAME>= xx</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>Other values are not implemented yet in Scilab</SP>
                      <PARAM_INDENT>
                        <PARAM_ITEM>
                          <PARAM_NAME>= 1</PARAM_NAME>
                          <PARAM_DESCRIPTION>
                            <SP>if the initial mesh is provided by the user.  it is defined in fspace as follows:  the mesh</SP>
                            <P>
    will occupy  <VERB>fspace(1), ..., fspace(n+1)</VERB>. the user needs to supply only the interior mesh points  <VERB>fspace(j) = x(j), j = 2, ..., n.</VERB>
  </P>
                          </PARAM_DESCRIPTION>
                        </PARAM_ITEM>
                        <PARAM_ITEM>
                          <PARAM_NAME>= 2 if the initial mesh is supplied by the user</PARAM_NAME>
                          <PARAM_DESCRIPTION>
                            <SP>as with <VERB>ipar(8)=1</VERB>, and in addition no adaptive mesh selection is to be done.</SP>
                          </PARAM_DESCRIPTION>
                        </PARAM_ITEM>
                      </PARAM_INDENT>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                </PARAM_INDENT>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
            <PARAM_ITEM>
              <PARAM_NAME>ipar(9)  </PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>( = iguess )</SP>
                <PARAM_INDENT>
                  <PARAM_ITEM>
                    <PARAM_NAME>= 0</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>if no initial guess for the solution is provided.</SP>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                  <PARAM_ITEM>
                    <PARAM_NAME>= 1</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>if an initial guess is provided by the user in subroutine  <VERB>guess</VERB>.</SP>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                  <PARAM_ITEM>
                    <PARAM_NAME>= 2</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>if an initial mesh and approximate solution coefficients are provided by the user in fspace. (the former and new mesh are the same).</SP>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                  <PARAM_ITEM>
                    <PARAM_NAME>= 3</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>if a former mesh and approximate solution coefficients are provided by the user in fspace, and the new mesh is to be taken twice as coarse; i.e.,every second point from the former mesh.</SP>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                  <PARAM_ITEM>
                    <PARAM_NAME>= 4</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>if in addition to a former initial mesh and approximate solution coefficients, a new mesh is provided in fspace as well. (see description of output for further details on iguess = 2, 3, and 4.)</SP>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                </PARAM_INDENT>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
            <PARAM_ITEM>
              <PARAM_NAME>ipar(10)  </PARAM_NAME>
              <PARAM_DESCRIPTION>
                <PARAM_INDENT>
                  <PARAM_ITEM>
                    <PARAM_NAME>= 0</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>if the problem is regular</SP>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                  <PARAM_ITEM>
                    <PARAM_NAME>= 1</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>if the first relax factor is =rstart, and the nonlinear iteration does not rely on past covergence (use for an extra sensitive nonlinear problem only).</SP>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                  <PARAM_ITEM>
                    <PARAM_NAME>= 2</PARAM_NAME>
                    <PARAM_DESCRIPTION>
                      <SP>if we are to return immediately upon  (a) two successive nonconvergences, or  (b) after obtaining error estimate for the first time.</SP>
                    </PARAM_DESCRIPTION>
                  </PARAM_ITEM>
                </PARAM_INDENT>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
            <PARAM_ITEM>
              <PARAM_NAME>ipar(11)  </PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>= number of fixed points in the mesh other than <VERB>aleft</VERB>
  and <VERB>aright</VERB>. ( = nfxpnt , the dimension of
  <VERB>fixpnt</VERB>) the code requires that all side condition
  points other than <VERB>aleft</VERB> and <VERB>aright</VERB> (see
  description of zeta ) be included as fixed points in
  <VERB>fixpnt</VERB>.</SP>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
          </PARAM_INDENT>
        </PARAM_DESCRIPTION>
      </PARAM_ITEM>
      <PARAM_ITEM>
        <PARAM_NAME>ltol</PARAM_NAME>
        <PARAM_DESCRIPTION>
          <SP>an array of dimension <VERB>ipar(4)</VERB>.
          <VERB>ltol(j) = l</VERB> specifies that the j-th tolerance
          in tol controls the error in the l-th component of
          <VERB>z(u)</VERB>.  also require that:</SP>
          <P>
          <VERB>1 &lt;= ltol(1) &lt; ltol(2) &lt; ... &lt; ltol(ntol)
          &lt;= mstar</VERB> 
          </P>
          
        </PARAM_DESCRIPTION>
      </PARAM_ITEM>
      <PARAM_ITEM>
        <PARAM_NAME>tol</PARAM_NAME>
        <PARAM_DESCRIPTION>
          <SP>an array of dimension
          <VERB>ipar(4)</VERB>. <VERB>tol(j)</VERB> is the error
          tolerance on the <VERB>ltol(j)</VERB> -th component of
          <VERB>z(u)</VERB>. thus, the code attempts to satisfy for
          <VERB>j=1:ntol</VERB>  on each subinterval</SP>
<VERBATIM>
        abs(z(v)-z(u))       &lt;=     tol(j)*abs(z(u))     +tol(j)
                     ltol(j)                       ltol(j)
</VERBATIM>
          
          <P>
    if <VERB>v(x)</VERB> is the approximate solution vector.
  </P>
        </PARAM_DESCRIPTION>
      </PARAM_ITEM>
      <PARAM_ITEM>
        <PARAM_NAME>fixpnt</PARAM_NAME>
        <PARAM_DESCRIPTION>
          <SP>an array of dimension <VERB>ipar(11)</VERB>. it contains the points, other than <VERB>aleft</VERB> and <VERB>aright</VERB>, which are to be included in every mesh.</SP>
        </PARAM_DESCRIPTION>
      </PARAM_ITEM>
      <PARAM_ITEM>
        <PARAM_NAME>externals</PARAM_NAME>
        <PARAM_DESCRIPTION>
          <SP>The function <VERB>fsub,dfsub,gsub,dgsub,guess</VERB> are Scilab
  externals i.e. functions (see syntax below) or the name of a Fortran
  subroutine (character string) with specified calling sequence or a
  list. An external as a  character string refers to the name of a
  Fortran subroutine. The Fortran coded function interface to bvode
  are specified in the file <VERB>fcol.f</VERB>.</SP>
          <PARAM_INDENT>
            <PARAM_ITEM>
              <PARAM_NAME>fsub</PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>name of subroutine for evaluating</SP>
<VERBATIM>
                                                t
                f(x,z(u(x))) =    (f ,...,f     )  
                                    1      ncomp
</VERBATIM>
                <P>
    at a point x in <VERB>(aleft,aright)</VERB>. it should have the heading  <VERB>[f]=fsub(x,z)</VERB>  where <VERB>f</VERB> is the vector containing the value of <VERB>fi(x,z(u))</VERB> in the i-th component and 
  </P>
 <VERBATIM>
                                                    t
                         z(u(x))=(z(1),...,z(mstar))
</VERBATIM>
                <P>
    is defined as above under  purpose .
  </P>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
            <PARAM_ITEM>
              <PARAM_NAME>dfsub</PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>name of subroutine for evaluating the Jacobian of
                <VERB>f(x,z(u))</VERB> at a point x.  it should have
                the heading <VERB>[df]=dfsub (x , z )</VERB> where
                <VERB>z(u(x))</VERB> is defined as for
                <VERB>fsub</VERB> and the (<VERB>ncomp</VERB>) by
                (<VERB>mstar</VERB>) array df should be filled by the
                partial derivatives of  f, viz, for a particular call
                one calculates</SP>

 <VERBATIM>
                    df(i,j) = dfi / dzj, i=1,...,ncomp
                                         j=1,...,mstar.
</VERBATIM>

              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
            <PARAM_ITEM>
              <PARAM_NAME>gsub</PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>name of subroutine for evaluating the i-th
                component of
                 <VERB>g(x,z(u(x))) = g (zeta(i),z(u(zeta(i)))) </VERB>
                 at a point <VERB>x = zeta(i)</VERB> where
                  <VERB>1&lt;=i&lt;=mstar.</VERB>
                </SP>
                <P>
                it should have the heading<VERB>[g]=gsub (i ,
                 z)</VERB> where <VERB>z(u)</VERB> is as for
                 <VERB>fsub</VERB>, and <VERB>i</VERB> and
                 <VERB>g=gi</VERB> are as above. Note that in contrast
                 to <VERB>f</VERB> in <VERB>fsub</VERB> , here only one value per
                 call is returned in <VERB>g</VERB>.
  </P>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
            <PARAM_ITEM>
              <PARAM_NAME>dgsub</PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>name of subroutine for evaluating the i-th row of the Jacobian of
  <VERB>g(x,u(x))</VERB>.  it should have the heading <VERB>[dg]=dgsub
  (i , z )</VERB> where <VERB>z(u)</VERB> is as for fsub, i as for
  gsub and the mstar-vector  <VERB>dg</VERB> should be filled with the
  partial derivatives of g, viz, for a particular call one calculates</SP>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
            <PARAM_ITEM>
              <PARAM_NAME>guess</PARAM_NAME>
              <PARAM_DESCRIPTION>
                <SP>name of subroutine to evaluate the initial approximation for
  <VERB>z(u(x))</VERB> and for <VERB>dmval(u(x))</VERB>= vector of the
  mj-th derivatives of <VERB>u(x)</VERB>. it should have the heading
  <VERB>[z,dmval]= guess (x )</VERB> note that this subroutine is used
  only if  <VERB>ipar(9) = 1</VERB>, and then all  <VERB>mstar</VERB>
  components of z and  ncomp  components of  dmval  should be
  specified for any x,</SP>
                <P>
                  <VERB>aleft &lt;= x &lt;= aright .</VERB>
                </P>
              </PARAM_DESCRIPTION>
            </PARAM_ITEM>
          </PARAM_INDENT>
        </PARAM_DESCRIPTION>
      </PARAM_ITEM>
    </PARAM_INDENT>
  </PARAM>
  <DESCRIPTION>
    <P>
    this package solves a multi-point boundary value
    problem for a mixed order system of ode-s given by</P>
<VERBATIM>
       (m(i))
       u       =  f  ( x; z(u(x)) )      i = 1, ... ,ncomp
        i          i                     aleft &lt; x  &lt; aright,
                                        
</VERBATIM>

<VERBATIM>
       g  ( zeta(j); z(u(zeta(j))) ) = 0   j = 1, ... ,mstar
        j
      mstar = m(1)+m(2)+...+m(ncomp),
</VERBATIM>
<P> where</P>
<VERBATIM>
                                        t
             u = (u , u , ... ,u     )   
                   1   2        ncomp    
</VERBATIM>
<P>is the exact solution vector</P>
<VERBATIM>
              (mi)
             u     is the mi=m(i) th  derivative of u
              i                                      i
</VERBATIM> 
<VERBATIM>                                     

                                (1)        (m1-1)       (mncomp-1)
             z(u(x)) = ( u (x),u  (x),...,u    (x),...,u      (x) )
                          1     1          1            ncomp
</VERBATIM> 
<VERBATIM>  
              f (x,z(u))   
               i
</VERBATIM>
    <P> is a (generally) nonlinear function of <VERB>z(u)=z(u(x))</VERB>.</P>
<VERBATIM> 
              g (zeta(j);z(u))  
               j
</VERBATIM>
      <P>  is a (generally) nonlinear function   used to represent a boundary condition.</P>

    <P> the boundary points satisfy</P>

     <P><VERB>aleft &lt;= zeta(1)  &lt;= ..  &lt;= zeta(mstar)  &lt;= aright</VERB>. </P>
    <P> the orders <VERB>mi</VERB> of the differential equations satisfy</P>
    <P><VERB>1&lt;=m(i)&lt;=4</VERB>.</P>
  </DESCRIPTION>
  <EXAMPLE>
<![CDATA[
deff('df=dfsub(x,z)','df=[0,0,-6/x**2,-6/x]')
deff('f=fsub(x,z)','f=(1 -6*x**2*z(4)-6*x*z(3))/x**3')
deff('g=gsub(i,z)','g=[z(1),z(3),z(1),z(3)];g=g(i)')
deff('dg=dgsub(i,z)',['dg=[1,0,0,0;0,0,1,0;1,0,0,0;0,0,1,0]';
                      'dg=dg(i,:)'])
deff('[z,mpar]=guess(x)','z=0;mpar=0')// unused here

 //define trusol for testing purposes
deff('u=trusol(x)',[ 
   'u=0*ones(4,1)';
   'u(1) =  0.25*(10*log(2)-3)*(1-x) + 0.5 *( 1/x   + (3+x)*log(x) - x)'
   'u(2) = -0.25*(10*log(2)-3)       + 0.5 *(-1/x^2 + (3+x)/x      + log(x) - 1)'
   'u(3) = 0.5*( 2/x^3 + 1/x   - 3/x^2)'
   'u(4) = 0.5*(-6/x^4 - 1/x/x + 6/x^3)'])

fixpnt=0;m=4;
ncomp=1;aleft=1;aright=2;
zeta=[1,1,2,2];
ipar=zeros(1,11);
ipar(3)=1;ipar(4)=2;ipar(5)=2000;ipar(6)=200;ipar(7)=1;
ltol=[1,3];tol=[1.e-11,1.e-11];
res=aleft:0.1:aright;
z=bvode(res,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...
 fsub,dfsub,gsub,dgsub,guess)
z1=[];for x=res,z1=[z1,trusol(x)]; end;  
z-z1
 ]]>
  </EXAMPLE>
  <SEE_ALSO>
    <SEE_ALSO_ITEM>
      <LINK>fort</LINK>
    </SEE_ALSO_ITEM>
    <SEE_ALSO_ITEM>
      <LINK>link</LINK>
    </SEE_ALSO_ITEM>
    <SEE_ALSO_ITEM>
      <LINK>external</LINK>
    </SEE_ALSO_ITEM>
    <SEE_ALSO_ITEM>
      <LINK>ode</LINK>
    </SEE_ALSO_ITEM>
    <SEE_ALSO_ITEM>
      <LINK>dassl</LINK>
    </SEE_ALSO_ITEM>
  </SEE_ALSO>
  <AUTHOR>u. ascher, department of computer science, university of british; columbia, vancouver, b. c., canada   v6t 1w5; g. bader, institut f. angewandte mathematik university of heidelberg; im neuenheimer feld 294d-6900 heidelberg 1 ; ; Fortran subroutine colnew.f</AUTHOR>
</MAN>
