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