Scilab Function
Last update : 00/00/0000
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
-
zThe solution of the ode evaluated on the mesh given by points
-
pointsan array which gives the points for which we want the solution
-
ncompnumber of differential equations (ncomp <= 20)
-
ma vector of size ncomp. m(j) gives the order of the j-th differential equation
-
aleftleft end of interval
-
arightright end of interval
-
zetazeta(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.
-
iparan 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 )
-
= -1for full diagnostic printout
-
= 0for selected printout
-
= 1for no printout
-
ipar(8) ( = iread )
-
= 0causes bvode to generate a uniform initial mesh.
-
= xxOther values are not implemented yet in Scilab
-
= 1if 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 useras with ipar(8)=1, and in addition no adaptive mesh selection is to be done.
-
ipar(9) ( = iguess )
-
= 0if no initial guess for the solution is provided.
-
= 1if an initial guess is provided by the user in subroutine guess.
-
= 2if an initial mesh and approximate solution coefficients are provided by the user in fspace. (the former and new mesh are the same).
-
= 3if 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.
-
= 4if 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)
-
= 0if the problem is regular
-
= 1if the first relax factor is =rstart, and the nonlinear iteration does not rely on past covergence (use for an extra sensitive nonlinear problem only).
-
= 2if 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.
-
ltolan 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
-
tolan 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.
-
fixpntan array of dimension ipar(11). it contains the points, other than aleft and aright, which are to be included in every mesh.
-
externalsThe 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.
-
fsubname 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 .
-
dfsubname 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.
-
gsubname 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.
-
dgsubname 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
-
guessname 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