Scilab Function
Last update : 28/6/2005

bvode - boundary value problems for ODE

### Calling Sequence

  `[z]=bvode(points,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...`
  `fsub1,dfsub1,gsub1,dgsub1,guess1)`

### Parameters

- `z`The solution of the ode evaluated on the mesh given by points
- `points`an array which gives the points for which we want the solution
- `ncomp`number of differential equations (ncomp <= 20)
- `m`a vector of size `ncomp`. `m(j)` gives the order of the j-th differential equation
- `aleft`left end of interval
- `aright`right end of interval
- `zeta` `zeta(j)` gives j-th side condition point (boundary point). must have

  `zeta(j) <= zeta(j+1)`

  all side condition points must be mesh points in all meshes used, see description of `ipar(11)` and `fixpnt` below.

- `ipar`an integer array dimensioned at least 11. a list of the parameters in `ipar` and their meaning follows some parameters are renamed in bvode; these new names are given in parentheses.
  - `ipar(1)` 0 if the problem is linear, 1 if the problem is nonlinear
  - `ipar(2)` = number of collocation points per subinterval (= k ) where

    `max m(i) <= k <= 7 .`

    if `ipar(2)=0` then bvode sets

    `k = max ( max m(i)+1, 5-max m(i) )`

  - `ipar(3)` = number of subintervals in the initial mesh ( = n ). if `ipar(3) = 0` then bvode arbitrarily sets `n = 5`.
  - `ipar(4)`= number of solution and derivative tolerances. ( = ntol ) we require

    `0 < ntol <= mstar.`

  - `ipar(5)`= dimension of `fspace` ( = ndimf ) a real work array. its size provides a constraint on nmax. choose ipar(5) according to the formula:

    `ipar(5)>=nmax*nsizef`

    where

    `nsizef=4+3*mstar+(5+kd)*kdm+(2*mstar-nrec)*2*mstar`.

  - `ipar(6)`= dimension of ispace ( = ndimi )an integer work array. its size provides a constraint on nmax, the maximum number of subintervals. choose `ipar(6)` according to the formula:

    `ipar(6)>=nmax*nsizei`

    where

    `nsizei=3 + kdm` with `kdm=kd+mstar` ; `kd=k*ncomp` ; `nrec`=number of right end boundary conditions.

  - `ipar(7)` output control ( = iprint )
    - `= -1`for full diagnostic printout
    - `= 0`for selected printout
    - `= 1`for no printout
  - `ipar(8)` ( = iread )
    - `= 0`causes bvode to generate a uniform initial mesh.
    - `= xx`Other values are not implemented yet in Scilab
      - `= 1`if the initial mesh is provided by the user. it is defined in fspace as follows: the mesh

        will occupy `fspace(1), ..., fspace(n+1)`. the user needs to supply only the interior mesh points `fspace(j) = x(j), j = 2, ..., n.`

      - `= 2 if the initial mesh is supplied by the user`as with `ipar(8)=1`, and in addition no adaptive mesh selection is to be done.
  - `ipar(9)` ( = iguess )
    - `= 0`if no initial guess for the solution is provided.
    - `= 1`if an initial guess is provided by the user in subroutine `guess`.
    - `= 2`if an initial mesh and approximate solution coefficients are provided by the user in fspace. (the former and new mesh are the same).
    - `= 3`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.
    - `= 4`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.)
  - `ipar(10)`
    - `= 0`if the problem is regular
    - `= 1`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).
    - `= 2`if we are to return immediately upon (a) two successive nonconvergences, or (b) after obtaining error estimate for the first time.
  - `ipar(11)` = number of fixed points in the mesh other than `aleft` and `aright`. ( = nfxpnt , the dimension of `fixpnt`) the code requires that all side condition points other than `aleft` and `aright` (see description of zeta ) be included as fixed points in `fixpnt`.
- `ltol`an array of dimension `ipar(4)`. `ltol(j) = l` specifies that the j-th tolerance in tol controls the error in the l-th component of `z(u)`. also require that:

  `1 <= ltol(1) < ltol(2) < ... < ltol(ntol) <= mstar`

- `tol`an array of dimension `ipar(4)`. `tol(j)` is the error tolerance on the `ltol(j)` -th component of `z(u)`. thus, the code attempts to satisfy for `j=1:ntol` on each subinterval

              abs(z(v)-z(u))       <=     tol(j)*abs(z(u))     +tol(j)
                           ltol(j)                       ltol(j)

  if `v(x)` is the approximate solution vector.

- `fixpnt`an array of dimension `ipar(11)`. it contains the points, other than `aleft` and `aright`, which are to be included in every mesh.
- `externals`The function `fsub,dfsub,gsub,dgsub,guess` 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 `fcol.f`.
  - `fsub`name of subroutine for evaluating

                                                        t
                        f(x,z(u(x))) =    (f ,...,f     )
                                            1      ncomp

    at a point x in `(aleft,aright)`. it should have the heading `[f]=fsub(x,z)` where `f` is the vector containing the value of `fi(x,z(u))` in the i-th component and

                                                            t
                                 z(u(x))=(z(1),...,z(mstar))

    is defined as above under purpose .

  - `dfsub`name of subroutine for evaluating the Jacobian of `f(x,z(u))` at a point x. it should have the heading `[df]=dfsub (x , z )` where `z(u(x))` is defined as for `fsub` and the (`ncomp`) by (`mstar`) array df should be filled by the partial derivatives of f, viz, for a particular call one calculates

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

  - `gsub`name of subroutine for evaluating the i-th component of `g(x,z(u(x))) = g (zeta(i),z(u(zeta(i))))` at a point `x = zeta(i)` where `1<=i<=mstar.`

    it should have the heading`[g]=gsub (i , z)` where `z(u)` is as for `fsub`, and `i` and `g=gi` are as above. Note that in contrast to `f` in `fsub` , here only one value per call is returned in `g`.

  - `dgsub`name of subroutine for evaluating the i-th row of the Jacobian of `g(x,u(x))`. it should have the heading `[dg]=dgsub (i , z )` where `z(u)` is as for fsub, i as for gsub and the mstar-vector `dg` should be filled with the partial derivatives of g, viz, for a particular call one calculates
  - `guess`name of subroutine to evaluate the initial approximation for `z(u(x))` and for `dmval(u(x))`= vector of the mj-th derivatives of `u(x)`. it should have the heading `[z,dmval]= guess (x )` note that this subroutine is used only if `ipar(9) = 1`, and then all `mstar` components of z and ncomp components of dmval should be specified for any x,

    `aleft <= x <= aright .`

### Description

this package solves a multi-point boundary value problem for a mixed order system of ode-s given by

           (m(i))
           u       =  f  ( x; z(u(x)) )      i = 1, ... ,ncomp
            i          i                     aleft < x  < aright,

           g  ( zeta(j); z(u(zeta(j))) ) = 0   j = 1, ... ,mstar
            j
          mstar = m(1)+m(2)+...+m(ncomp),

where

                                            t
                 u = (u , u , ... ,u     )
                       1   2        ncomp

is the exact solution vector

                  (mi)
                 u     is the mi=m(i) th  derivative of u
                  i                                      i

                                    (1)        (m1-1)       (mncomp-1)
                 z(u(x)) = ( u (x),u  (x),...,u    (x),...,u      (x) )
                              1     1          1            ncomp

                  f (x,z(u))
                   i

is a (generally) nonlinear function of `z(u)=z(u(x))`.

                  g (zeta(j);z(u))
                   j

is a (generally) nonlinear function used to represent a boundary condition.

the boundary points satisfy

`aleft <= zeta(1) <= .. <= zeta(mstar) <= aright`.

the orders `mi` of the differential equations satisfy

`1<=m(i)<=4`.

### Examples

    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

### See Also

`fort`, `link`, `external`, `ode`, `dassl`,

### 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