6 Non-linear tools (optimization and simulation)
6.1 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
( mstar = m(1) + ... + m(ncomp) <= 40 )
- 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.
1
- ipar(1)
( = nonlin )
2
- = 0
if the problem is linear
- = 1 if the problem is nonlinear
1
- 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 )
2
- = -1
for full diagnostic printout
- = 0
for selected printout
- = 1
for no printout
1
- ipar(8)
( = iread )
2
- = 0
causes bvode to generate a uniform initial mesh.
- = xx
Other values are not implemented yet in Scilab
3
- = 1
if the initial mesh is provided by the user. it
is defined in fspace as follows: the mesh
aleft=x(1)<x(2)< ... <x(n)<x(n+1)=aright
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.
2
1
- ipar(9)
( = iguess )
2
- = 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.)
1
- ipar(10)
2
- = 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.
1
- 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.
0
- 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)) \le 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.
1
- 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))))
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
dg(i,j) = dgi / dzj, j=1,...,mstar.
- 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 .
0
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 ) is the exact solution vector
1 2 ncomp
(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)) is a (generally) nonlinear function of
i
z(u)=z(u(x)).
g (zeta(j);z(u)) is a (generally) nonlinear function
j
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.
function [z,z1]=col1()
fixpnt=0
m=[4]
ncomp=1
aleft=1
aright=2
zeta=[1,1,2,2]
ipar=0*ones(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;
function [df]=dfsub(x,z)
df=[0,0,-6/x**2,-6/x]
function [f]=fsub(x,z)
f=(1 -6*x**2*z(4)-6*x*z(3))/x**3
function [g]=gsub(i,z)
g=[z(1),z(3),z(1),z(3)]
g=g(i)
function [dg]=dgsub(i,z)
dg=[1,0,0,0;0,0,1,0;1,0,0,0;0,0,1,0]
dg=dg(i,:)
function [z,mpar]=guess(x)
// unused here
z=0
mpar=0
function [u]=trusol(x)
u=0*ones(4,1)
U(1) = .25* (10.*LOG(2.)-3.) * (1.-X) +0.5* (1./X+ (3.+X)*LOG(X) - X)
U(2) = -.25* (10.*LOG(2.) - 3.) + .5 * (-1./X/X + LOG(X) + (3.+X)/X - 1.)
U(3) = .5 * (2./X**3 + 1./X -3./X/X)
U(4) = .5 * (-6./X**4 - 1./X/X + 6./X**3)
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
Fotran subroutine colnew.f
See Also :
fort X, link X, external X, ode X, dassl X
6.2 colnew boundary value problems for ODE
CALLING SEQUENCE :
This function has been renamed bvode.
6.3 dasrt DAE solver with zero crossing
CALLING SEQUENCE :
[r,nn,[,hd]]=dasrt(x0,t0,t [,atol,[rtol]],res [,jac],ng, surf, info [,hd])
PARAMETERS :
- x0
: is either y0 (ydot0 is estimated by dassl with zero as first estimate) or the matrix [y0 ydot0].
g(t,y0,ydot0) must be equal to zero. If you only know an estimate
of ydot0 set info(7)=1
1
- y0
: real column vector of initial conditions.
- ydot0
: real column vector of the time derivative of y at t0 (may be an estimate).
0
- t0
: real number is the initial instant.
- t
: real scalar or vector. Gives instants for which you want the solution. Note that you can get solution at each dassl's step point by setting info(2)=1.
- nn
: a vector with two entries [times num] times is the value
of the time at which the surface is crossed, num is the number
of the crossed surface
- atol,rtol
: real scalars or column vectors of same size as y. atol,rtol give respectively absolute and relative error tolerances of solution.
If vectors the tolerances are specified for each component of y.
- res
: external (function or list or string). Computes the value of g(t,y,ydot).
1
- function
: Its calling sequence must be [r,ires]=res(t,y,ydot) and res must return the residue r=g(t,y,ydot) and error flag
ires. ires = 0 if res succeeds to compute r, =-1 if residue is locally not defined for (t,y,ydot), =-2 if
parameters are out of admissible range.
- list
: it must be as follows:
list(res,x1,x2,...)
where the calling sequence of the function res is now
r=res(t,y,ydot,x1,x2,...)
res still returns r=g(t,y,ydot) as a function of
(t,y,ydot,x1,x2,...).
- string
: it must refer to the name of
a fortran subroutine (see source code of fresd.f).
0
- jac
: external (function or list or string). Computes the value
of dg/dy+cj*dg/dydot for a given value of parameter cj
1
- function
: Its calling sequence must be r=jac(t,y,ydot,cj) and the jac function must return
r=dg(t,y,ydot)/dy+cj*dg(t,y,ydot)/dydot where cj is a real scalar
- list
: it must be as follows
list(jac,x1,x2,...)
where the calling sequence of the function jac is now
r=jac(t,y,ydot,x1,x2,...)
jac still returns dg/dy+cj*dg/dydot as a function of
(t,y,ydot,cj,x1,x2,...).
- character string
: it must refer to the name of a fortran subroutine
(see source code of jacdd.f).
0
- surf
: external (function or list or string). Computes the value
of the column vector surf(t,y) with ng components.
Each component defines a surface.
1
- function
: Its calling sequence must be surf(t,y)
- list
: it must be as follows
list(surf,x1,x2,...)
where the calling sequence of the function surf is now
r=surf(t,y,x1,x2,...)
- character string
: it must refer to the name of a fortran subroutine
(see source code of fsurfd.f) in directory SCDIR/default
0
- info
: list which contains 7 elements:
1
- info(1)
: real scalar which gives the maximum time for which g is allowed
to be evaluated or an empty matrix [] if no limits imposed for time.
- info(2)
: flag which indicates if dassl returns its intermediate
computed values (flag=1) or only the user specified time point
values (flag=0).
- info(3)
: 2 components vector which give the definition [ml,mu] of band
matrix computed by jac;
r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)".
If jac returns a full matrix set info(3)=[].
- info(4)
: real scalar which gives the maximum step size. Set info(4)=[] if no
limitation.
- info(5)
: real scalar which gives the initial step size. Set info(4)=[] if
not specified.
- info(6)
: set info(6)=1 if the solution is known to be non negative,
else set info(6)=0.
- info(7)
: set info(7)=1 if ydot0 is just an estimation, info(7)=0 if g(t0,y0,ydot0)=0.
0
- hd
: real vector which allows to store the dassl context and to
resume integration
- r
: real matrix . Each column is the vector [t;x(t);xdot(t)] where t is time
index for which the solution had been computed
DESCRIPTION :
Solution of the implicit differential equation
g(t,y,ydot)=0
y(t0)=y0 and ydot(t0)=ydot0
Returns the surface crossing instants and the number of the
surface reached in nn.
Detailed examples can be found in SCIDIR/tests/dassldasrt.tst
EXAMPLE :
//dy/dt = ((2*log(y)+8)/t -5)*y, y(1) = 1, 1<=t<=6
//g1 = ((2*log(y)+8)/t - 5)*y
//g2 = log(y) - 2.2491
y0=1;t=2:6;t0=1;y0d=3;
info=list([],0,[],[],[],0,0);
atol=1.d-6;rtol=0;ng=2;
deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2*log(y)+8)/t-5)*y')
deff('[rts]=gr1(t,y)','rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')
[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
//(Should return nn=[2.4698972 2])
See Also :
ode X, dassl X, impl X, fort X, link X, external X
6.4 dassl differential algebraic equation
CALLING SEQUENCE :
[r [,hd]]=dassl(x0,t0,t [,atol,[rtol]],res [,jac],info [,hd])
PARAMETERS :
- x0
: is either y0 (ydot0 is estimated by dassl with zero as first estimate) or the matrix [y0 ydot0].
g(t,y0,ydot0) must be equal to zero. If you only know an estimate
of ydot0 set info(7)=1
1
- y0
: real column vector of initial conditions.
- ydot0
: real column vector of the time derivative of y at t0 (may be an estimate).
0
- t0
: real number is the initial instant.
- t
: real scalar or vector. Gives instants for which you want the solution. Note that you can get solution at each dassl's step point by setting info(2)=1.
- atol,rtol
: real scalars or column vectors of same size as y. atol,rtol give respectively absolute and relative error tolerances of solution.
If vectors the tolerances are specified for each component of y.
- res
: external (function or list or string). Computes the value of g(t,y,ydot).
1
- function
: Its calling sequence must be [r,ires]=res(t,y,ydot) and res must return the residue r=g(t,y,ydot) and error flag
ires. ires = 0 if res succeeds to compute r, =-1 if residue is locally not defined for (t,y,ydot), =-2 if
parameters are out of admissible range.
- list
: it must be as follows:
list(res,x1,x2,...)
where the calling sequence of the function res is now
r=res(t,y,ydot,x1,x2,...)
res still returns r=g(t,y,ydot) as a function of
(t,y,ydot,x1,x2,...).
- string
: it must refer to the name of
a fortran subroutine (see source code of fresd.f).
0
- jac
: external (function or list or string). Computes the value
of dg/dy+cj*dg/dydot for a given value of parameter cj
1
- function
: Its calling sequence must be r=jac(t,y,ydot,cj) and the jac function must return
r=dg(t,y,ydot)/dy+cj*dg(t,y,ydot)/dydot where cj is a real scalar
- list
: it must be as follows
list(jac,x1,x2,...)
where the calling sequence of the function jac is now
r=jac(t,y,ydot,x1,x2,...)
jac still returns dg/dy+cj*dg/dydot as a function of
(t,y,ydot,cj,x1,x2,...).
- character string
: it must refer to the name of a fortran subroutine
(see source code of jacdd.f).
0
- info
: list which contains 7 elements:
1
- info(1)
: real scalar which gives the maximum time for which g is allowed
to be evaluated or an empty matrix [] if no limits imposed for time.
- info(2)
: flag which indicates if dassl returns its intermediate
computed values (flag=1) or only the user specified time point
values (flag=0).
- info(3)
: 2 components vector which give the definition [ml,mu] of band
matrix computed by jac;
r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)".
If jac returns a full matrix set info(3)=[].
- info(4)
: real scalar which gives the maximum step size. Set info(4)=[] if no
limitation.
- info(5)
: real scalar which gives the initial step size. Set info(4)=[] if
not specified.
- info(6)
: set info(6)=1 if the solution is known to be non negative,
else set info(6)=0.
- info(7)
: set info(7)=1 if ydot0 is just an estimation, info(7)=0 if g(t0,y0,ydot0)=0.
0
- hd
: real vector which allows to store the dassl context and to
resume integration
- r
: real matrix . Each column is the vector [t;x(t);xdot(t)] where t is time
index for which the solution had been computed
DESCRIPTION :
Solution of the implicit differential equation
g(t,y,ydot)=0
y(t0)=y0 and ydot(t0)=ydot0
Detailed examples are given in SCIDIR/tests/dassldasrt.tst
EXAMPLES :
deff('[r,ires]=chemres(t,y,yd)',[
'r(1)=-0.04*y(1)+1d4*y(2)*y(3)-yd(1);';
'r(2)=0.04*y(1)-1d4*y(2)*y(3)-3d7*y(2)*y(2)-yd(2);'
'r(3)=y(1)+y(2)+y(3)-1;'
'ires=0']);
deff('[pd]=chemjac(x,y,yd,cj)',[
'pd=[-0.04-cj , 1d4*y(3) , 1d4*y(2);';
'0.04 ,-1d4*y(3)-2*3d7*y(2)-cj ,-1d4*y(2);';
'1 , 1 , 1 ]'])
y0=[1;0;0];
yd0=[-0.04;0.04;0];
t=[1.d-5:0.02:.4,0.41:.1:4,40,400,4000,40000,4d5,4d6,4d7,4d8,4d9,4d10];
info=list([],0,[],[],[],0,0);
y=dassl([y0,yd0],0,t,chemres,info);
info(2)=1;
y=dassl([y0,yd0],0,4d10,chemres,info);
y=dassl([y0,yd0],0,4d10,chemres,chemjac,info);
See Also :
ode X, dasrt X, impl X, fort X, link X, external X
6.5 fit_dat Parameter identification based on measured data
CALLING SEQUENCE :
[p,err]=fit_dat(G,p0,Z [,W] [,pmin,pmax] [,DG])
PARAMETERS :
- G
: Scilab function (e=G(p,z), e: nex1, p: npx1, z: nzx1)
- p0
: initial guess (size npx1)
- Z
: matrix [z_1,z_2,...z_n] where z_i (nzx1) is the ith measurement
- W
: weighting matrix of size nexne (optional; default 1)
- pmin
: lower bound on p (optional; size npx1)
- pmax
: upper bound on p (optional; size npx1)
- DG
: partial of G wrt p (optional; S=DG(p,z), S: nexnp)
DESCRIPTION :
fit_dat is used for fitting data to a model.
For a given function G(p,z), this function finds the best vector
of parameters p for approximating G(p,z_i)=0 for a set of measurement
vectors z_i. Vector p is found by minimizing
G(p,z_1)'WG(p,z_1)+G(p,z_2)'WG(p,z_2)+...+G(p,z_n)'WG(p,z_n)
EXAMPLE :
deff('y=FF(x)','y=a*(x-b)+c*x.*x')
X=[];Y=[];
a=34;b=12;c=14;for x=0:.1:3, Y=[Y,FF(x)+100*(rand()-.5)];X=[X,x];end
Z=[Y;X];
deff('e=G(p,z)','a=p(1),b=p(2),c=p(3),y=z(1),x=z(2),e=y-FF(x)')
[p,err]=fit_dat(G,[3;5;10],Z)
xset('window',0)
xbasc();
plot2d(X',Y',-1)
plot2d(X',FF(X)',5,'002')
a=p(1),b=p(2),c=p(3);plot2d(X',FF(X)',12,'002')
a=34;b=12;c=14;
deff('s=DG(p,z)','y=z(1),x=z(2),s=-[x-p(2),-p(1),x*x]')
[p,err]=fit_dat(G,[3;5;10],Z,DG)
xset('window',1)
xbasc();
plot2d(X',Y',-1)
plot2d(X',FF(X)',5,'002')
a=p(1),b=p(2),c=p(3);plot2d(X',FF(X)',12,'002')
See Also :
optim X
6.6 fsolve find a zero of a system of n nonlinear functions
CALLING SEQUENCE :
[x [,v [,info]]]=fsolve(x0,fct [,fjac] [,tol])
PARAMETERS :
- x0
: real vector (initial value of function argument).
- fct
: external (i.e function or list or string).
- fjac
: external (i.e function or list or string).
- tol
: real scalar. precision tolerance: termination occurs
when the algorithm estimates that the relative error
between x and the solution is at most tol.
(tol=1.d-10 is the default value).
- x :
real vector (final value of function argument, estimated zero).
- v :
real vector (value of function at x).
- info
: termination indicator
1
- 0
: improper input parameters.
- 1
: algorithm estimates that the relative error between x and the solution
is at most tol.
- 2
: number of calls to fcn reached
- 3
: tol is too small. No further improvement in the approximate solution
x is possible.
- 4
: iteration is not making good progress.
0
DESCRIPTION :
find a zero of a system of
n nonlinear functions in n variables by a modification
of the powell hybrid method. Jacobian may be provided.
0 = fct(x) w.r.t x.
fct is an "external".
This external returns v=fct(x) given x.
The simplest calling sequence for fct is:
[v]=fct(x).
If fct is a character string, it refers to a fortran routine
which must be linked to Scilab (see fsolvf.f for the calling
sequence). Dynamic link is possible (help link).
jac is an "external".
This external returns v=d(fct)/dx (x) given x.
The simplest calling sequence for jac is:
[v]=jac(x).
If jac is a character string, it refers to a fortran routine
which must be linked to Scilab (see fsolvj.f for the calling
sequence). Dynamic link is possible (help link).
EXAMPLES :
// A simple example with fsolve
a=[1,7;2,8];b=[10;11];
deff('[y]=fsol1(x)','y=a*x+b');
deff('[y]=fsolj1(x)','y=a');
[xres]=fsolve([100;100],fsol1);
a*xres+b
[xres]=fsolve([100;100],fsol1,fsolj1);
a*xres+b
// See routines/default/Ex-fsolve.f
[xres]=fsolve([100;100],'fsol1','fsolj1',1.e-7);
a*xres+b
See Also :
external X, quapro X, linpro X, optim X
6.7 impl differential algebraic equation
DESCRIPTION :
y=impl([type],y0,ydot0,t0,t [,atol, [rtol]],res,adda [,jac])
PARAMETERS :
- y0,ydot0
: real vectors or matrix (initial conditions).
- t0
: real scalar (initial time).
- t
: real vector (times at which the solution is computed).
- res,adda
: externals (function or character string or list).
- type
: string 'adams' or 'stiff'
- atol,rtol
: real scalar or real vector of the same size as as y.
- jac
: external (function or character string or list).
DESCRIPTION :
Solution of the linear implicit differential equation
A(t,y) dy/dt=g(t,y), y(t0)=y0
t0 is the initial instant, y0 is the vector of initial conditions
Vector ydot0 of the time derivative of y at t0 must
also be given.
r
The input res is an external i.e. a function with
specified syntax, or the name a Fortran subroutine or a C function
(character string) with specified calling sequence or a list.
If res is a function, its syntax must be as follows:
r = res(t,y,ydot)
where t is a real scalar (time) and y and ydot are
real vector (state and derivative of the state).
This function must return r=g(t,y)-A(t,y)*ydot.
If res is a character string, it refers to the name of a Fortran
subroutine or a C function. See
SCIDIR/routines/default/Ex-impl.f for an example to do that.
res can also be a list: see the help of ode.
The input adda is also an external.
If adda is a function, its syntax must be as follows:
r = adda(t,y,p)
and it must return r=A(t,y)+p where p is a matrix to be
added to A(t,y).
If adda is a character string, it refers to the name of a Fortran
subroutine or a C function. See
SCIDIR/routines/default/Ex-impl.f for an example to do that.
adda can also be a list: see the help of ode.
The input jac is also an external.
If adda is a function, its syntax must be as follows:
j = jac(t,y,ydot)
and it must return the Jacobian of r=g(t,y)-A(t,y)*ydot with
respect to y.
If jac is a character string, it refers to the name of a Fortran
subroutine or a C function. See
SCIDIR/routines/default/Ex-impl.f for an example to do that.
jac can also be a list: see the help of ode.
EXAMPLES :
y=impl([1;0;0],[-0.04;0.04;0],0,0.4,'resid','aplusp');
// Using hot restart
//[x1,w,iw]=impl([1;0;0],[-0.04;0.04;0],0,0.2,'resid','aplusp');
// hot start from previous call
//[x1]=impl([1;0;0],[-0.04;0.04;0],0.2,0.4,'resid','aplusp',w,iw);
//maxi(abs(x1-x))
See Also :
dassl X, ode X, external X
6.8 int2d definite 2D integral by quadrature and cubature method
CALLING SEQUENCE :
[I,err]=int2d(X,Y,f [,params])
PARAMETERS :
- X
: a 3 by N array containing the abscissae of the vertices
of the N triangles.
- Y
: a 3 by N array containing the ordinates of the vertices
of the N triangles.
- f
: external (function or list or string) defining the integrand f(u,v);
- params
: real vector [tol, iclose, maxtri, mevals, iflag]. default
value is [1.d-10, 1, 50, 4000, 1].
1
- tol
:the desired bound on the error. If iflag=0,
tol is interpreted as a bound on the relative error;
if iflag=1, the bound is on the absolute error.
- iclose
:an integer parameter that determines the selection
of LQM0 or LQM1 methods. If iclose=1 then LQM1 is used.
Any other value of iclose causes LQM0 to be used.
LQM0 uses function values only at interior points of
the triangle. LQM1 is usually more accurate than LQM0
but involves evaluating the integrand at more points
including some on the boundary of the triangle. It
will usually be better to use LQM1 unless the integrand
has singularities on the boundary of the triangle.
- maxtri
:the maximum number of triangles in the final triangulation of the
region
- mevals
: the maximum number of function evaluations
to be allowed. This number will be effective in limiting
the computation only if it is less than 94*maxtri when LQM1
is specified or 56*maxtri when LQM0 is specified.
- iflag
:
0
- I
: the integral value
- err
: the estimated error
DESCRIPTION :
int2d computes the two-dimensional integral of a
function f over a region consisting of n triangles.
A total error estimate is obtained and compared with a
tolerance - tol - that is provided as input to the subroutine.
The error tolerance is treated as either relative or absolute
depending on the input value of iflag. A 'Local Quadrature
Module' is applied to each input triangle and estimates of the
total integral and the total error are computed. The local
quadrature module is either subroutine LQM0 or subroutine
LQM1 and the choice between them is determined by the
value of the input variable iclose.
If the total error estimate exceeds the tolerance, the triangle
with the largest absolute error is divided into two triangles
by a median to its longest side. The local quadrature module
is then applied to each of the subtriangles to obtain new
estimates of the integral and the error. This process is
repeated until either (1) the error tolerance is satisfied,
(2) the number of triangles generated exceeds the input
parameter maxtri, (3) the number of integrand evaluations
exceeds the input parameter mevals, or (4) the function
senses that roundoff error is beginning to contaminate
the result.
EXAMPLES :
X=[0,0;1,1;1,0];
Y=[0,0;0,1;1,1];
deff('z=f(x,y)','z=cos(x+y)')
[I,e]=int2d(X,Y,f)
// computes the integrand over the square [0 1]x[0 1]
See Also :
intc X, intl X, int3d X, intg X, mesh2d X
REFERENCES :
Fortran routine twodq,Authors: Kahaner,D.K.,N.B.S., Rechard,O.W.,N.B.S.,
Barnhill,Robert,Univ. of UTAH
6.9 int3d definite 3D integral by quadrature and cubature method
CALLING SEQUENCE :
[result,err]=int3d(X,Y,Z,f [,nf[,params]])
PARAMETERS :
- X
: a 4 by NUMTET array containing the abscissae of the vertices
of the NUMTET tetrahedrons.
- Y
: a 4 by NUMTET array containing the ordinates of the vertices
of the NUMTET tetrahedrons.
- Z
: a 4 by NUMTET array containing the third coordinates of the vertices
of the NUMTET tetrahedrons.
- f
: external (function or list or string) defining the integrand
f(xyz,nf), where xyz is the vector of a point coordinates
and nf the number functions
- nf
: the number of function to integate (default is 1)
- params
: real vector [minpts, maxpts, epsabs, epsrel]. default
value is [0, 1000, 0.0, 1.d-5].
1
- epsabs
: Desired bound on the absolute error.
- epsrel
: Desired bound on the relative error.
- minpts
: Minimum number of function evaluations.
- maxpts
: Maximum number of function evaluations. The number of
function evaluations over each subregion is 43
0
- result
: the integral value,or vector of the integral values.
- err
: Estimates of absolute errors.
DESCRIPTION :
The function calculates an approximation to a given
vector of definite integrals
I I I (F ,F ,...,F ) dx(3)dx(2)dx(1),
1 2 numfun
where the region of integration is a collection of
NUMTET tetrahedrons and where
F = F (X(1),X(2),X(3)), J = 1,2,...,NUMFUN.
J J
A globally adaptive strategy is applied in order to
compute approximations result(k) hopefully satisfying, for each component of I, the
following claim for accuracy:
ABS(I(K)-RESULT(K))<=MAX(EPSABS,EPSREL*ABS(I(K)))
int3d repeatedly subdivides the tetrahedrons with greatest
estimated errors and estimates the integrals and the errors over the
new subtetrahedrons until the error request is met or MAXPTS function evaluations have been used.
A 43 point integration rule with all evaluation points inside the
tetrahedron is applied. The rule has polynomial degree 8.
If the values of the input parameters EPSABS or EPSREL are
selected great enough, an integration rule is applied over each
tetrahedron and the results are added up to give the approximations
RESULT(K). No further subdivision of the tetrahedrons will then
be applied.
When int3d computes estimates to a vector of
integrals, all components of the vector are given
the same treatment. That is, I(Fj) and I(Fk) for
j not equal to k, are estimated with the same
subdivision of the region of integration.
For integrals with enough similarity, we may save
time by applying int3d to all integrands in one call.
For integrals that varies continuously as functions of
some parameter, the estimates produced by int3d will
also vary continuously when the same subdivision is
applied to all components. This will generally not be
the case when the different components are given
separate treatment.
On the other hand this feature should be used with
caution when the different components of the integrals
require clearly different subdivisions.
EXAMPLES :
X=[0;1;0;0];
Y=[0;0;1;0];
Z=[0;0;0;1];
deff('v=f(xyz,numfun)','v=exp(xyz''*xyz)')
[RESULT,ERROR]=int3d(X,Y,Z,'int3dex')
// computes the integrand exp(x*x+y*y+z*z) over the
//tetrahedron (0.,0.,0.),(1.,0.,0.),(0.,1.,0.),(0.,0.,1.)
See Also :
intc X, intl X, int3d X
REFERENCES :
Fortran routine dqtet.f
Authors:
Jarle Berntsen, The Computing Centre,
University of Bergen, Thormohlens gt. 55,
N-5008 Bergen, Norway
Phone.. 47-5-544055
Email.. jarle@eik.ii.uib.no
Ronald Cools, Dept. of Computer Science,
Katholieke Universiteit Leuven, Celestijnenlaan 200A,
B-3030 Heverlee, Belgium
Phone.. 32-16-201015 (3562)
Email.. ronald@cs.kuleuven.ac.be
Terje O. Espelid, Department of Informatics,
University of Bergen, Thormohlens gt. 55,
N-5008 Bergen, Norway
Phone.. 47-5-544180
Email.. terje@eik.ii.uib.no
6.10 intc Cauchy integral
CALING SEQUENCE :
[y]=intc(a,b,f)
PARAMETERS :
- a,b
: two complex numbers
- f
: "external" function
DESCRIPTION :
If f is a complex-valued function, intc(a,b,f) computes the integral from a to b of f(z)dz along the
straight line a b of the complex plane.
See Also :
intg X, intl X
Author :
F. D.
6.11 intg definite integral
CALLING SEQUENCE :
[v,err]=intg(a,b,f [,ea [,er])
PARAMETERS :
- a,b
: real numbers
- f
: external (function or list or string)
- ea, er
: real numbers
- ea
: absolute error required on the result. Default value: 0
- er
: relative error required on the result. Default value: 1.d-8
- err
: estimated absolute error on the result.
DESCRIPTION :
intg(a,b,f) evaluates the definite integral from a to
b of f(t)dt.
The evaluation hopefully satisfies following claim for accuracy:
abs(I-v)<= max(ea,er*abs(I)) where I stands for the
exact value of the integral.
f is an external :
If f is function its definition must be as follows
y = f(t)
If f is a list the list must be as follows:
list(f,x1,x2,...) where f is a function with calling
sequence f(t,x1,x2,...).
If f is a string it refers to a the name
of a Fortran subroutine (see source code of fintg.f)
EXAMPLE :
deff('[y]=f(x)','y=x*sin(30*x)/sqrt(1-((x/(2*%pi))^2))')
exact=-2.5432596188;
abs(exact-intg(0,2*%pi,f))
// See file routines/default/Ex-intg.f
abs(exact-intg(0,2*%pi,'intgex'))
See Also :
intc X, intl X
6.12 intl Cauchy integral
CALLING SEQUENCE :
[y]=intl(a,b,z0,r,f)
PARAMETERS :
- z0
: complex number
- a,b
: two real numbers
- r
: positive real number
- f
: "external" function
DESCRIPTION :
If f is a complex-valued function, intl(a,b,z0,r,f) computes
the integral of f(z)dz along the curve
in the complex plane defined by z0 + r.exp(%i*t) for
a<=t<=b .(part of the circle with center z0 and radius r with phase between a and b)
See Also :
intc X
Author :
F. D.
6.13 karmarkar karmarkar (test program)
CALLING SEQUENCE :
[x1]=karmarkar(a,b,c,x0)
PARAMETERS :
- a
: matrix (n,p)
- b
: n - vector
- c
: p - vector
- x0
: initial vector
- eps
: threshold (default value : 1.d-5)
- gamma
: descent step 0<gamma<1 , default value : 1/4
- x1
: solution
- crit
: value of c'*x1
DESCRIPTION :
Computes x which minimizes
c'*x
under constraints:
a*x = b
x>=0
EXAMPLE :
// n=10;p=20;
// a=rand(n,p);c=rand(p,1);x0=abs(rand(p,1));b=a*x0;x1=karmarkar(a,b,c,x0);
6.14 linpro linear programming solver
CALLING SEQUENCE :
[x,lagr,f]=linpro(p,C,b [,x0])
[x,lagr,f]=linpro(p,C,b,ci,cs [,x0])
[x,lagr,f]=linpro(p,C,b,ci,cs,mi [,x0])
[x,lagr,f]=linpro(p,C,b,ci,cs,mi,x0 [,imp])
PARAMETERS :
- p
: real (column) vector (dimension n)
- C
: real matrix (dimension (mi + md) x n)
(If no constraints are given, you can set C = [])
- b
: RHS vector (dimension 1 x (mi + md))
- ci
: (column) vector of lower-bounds (dimension n).
If there are no lower bound constraints, put ci = [].
If some components of x are bounded from below,
set the other (unconstrained) values of ci to a very
large negative number (e.g. ci(j) = -(% eps)^(-1).
- cs
: (column) vector of upper-bounds. (Same remarks as above).
- mi
: number of equality constraints (i.e. C(1:mi,:)*x = b(1:mi))
- x0
: either an initial guess for x or one of the character strings 'v' or 'g'.
If x0='v' the calculated initial feasible point is a vertex.
If x0='g' the calculated initial feasible point is arbitrary.
- imp
: verbose option (optional parameter) (Try imp=7,8,...)
- x
: optimal solution found.
- f
: optimal value of the cost function (i.e. f=p'*x).
- lagr
: vector of Lagrange multipliers.
If lower and upper-bounds ci,cs are provided, lagr has
n + mi + md components and lagr(1:n) is the Lagrange
vector associated with the bound constraints and
lagr (n+1 : n + mi + md) is the Lagrange vector associated
with the linear constraints.
(If an upper-bound (resp. lower-bound) constraint i is active
lagr(i) is > 0 (resp. <0).
If no bounds are provided, lagr has only mi + md components.
DESCRIPTION :
[x,lagr,f]=linpro(p,C,b [,x0])
Minimize p'*x
under the constraints
C*x <= b
[x,lagr,f]=linpro(p,C,b,ci,cs [,x0])
Minimize p'*x
under the constraints
C*x <= b ci <= x <= cs
[x,lagr,f]=linpro(p,C,b,ci,cs,mi [,x0])
Minimize p'*x
under the constraints
C(j,:) x = b(j), j=1,...,mi
C(j,:) x <= b(j), j=mi+1,...,mi+md
ci <= x <= cs
If no initial point is given the
program computes a feasible initial point
which is a vertex of the region of feasible points if
x0='v'.
If x0='g', the program computes a feasible initial
point which is not necessarily a vertex. This mode is
advisable when the quadratic form is positive
definite and there are a few constraints in
the problem or when there are large bounds
on the variables that are security bounds and
very likely not active at the optimal solution.
EXAMPLE :
//Find x in R^6 such that:
//C1*x = b1 (3 equality constraints i.e mi=3)
C1= [1,-1,1,0,3,1;
-1,0,-3,-4,5,6;
2,5,3,0,1,0];
b1=[1;2;3];
//C2*x <= b2 (2 inequality constraints)
C2=[0,1,0,1,2,-1;
-1,0,2,1,1,0];
b2=[-1;2.5];
//with x between ci and cs:
ci=[-1000;-10000;0;-1000;-1000;-1000];cs=[10000;100;1.5;100;100;1000];
//and minimize p'*x with
p=[1;2;3;4;5;6]
//No initial point is given: x0='v';
C=[C1;C2]; b=[b1;b2] ; mi=3; x0='v';
[x,lagr,f]=linpro(p,C,b,ci,cs,mi,x0)
// Lower bound constraints 3 and 4 are active and upper bound
// constraint 5 is active --> lagr(3:4) < 0 and lagr(5) > 0.
// Linear (equality) constraints 1 to 3 are active --> lagr(7:9) <> 0
See Also :
quapro X
Author :
E. Casas, C. Pola Mendez
6.15 lmisolver linear matrix inequation solver
CALLING SEQUENCE :
[XLISTF[,OPT]] = lmisolver(XLIST0,evalfunc [,options])
PARAMETERS :
- XLIST0
: a list of containing initial guess (e.g. XLIST0=list(X1,X2,..,Xn))
- evalfunc
: a Scilab function ("external" function with specific syntax)
- XLISTF
: a list of matrices (e.g. XLIST0=list(X1,X2,..,Xn))
- options
: optional parameter. If given, options is
a real row vector with 5 components
[Mbound,abstol,nu,maxiters,reltol]
The syntax the function evalfunc must be as follows:
[LME,LMI,OBJ]=evalfunct(X) where X is a list of matrices,
LME, LMI are lists and OBJ a real scalar.
DESCRIPTION :
lmisolver solves the following problem:
minimize f(X1,X2,...,Xn) a linear function of Xi's
under the linear constraints:
Gi(X1,X2,...,Xn)=0 for i=1,...,p and LMI (linear matrix
inequalities) constraints:
Hj(X1,X2,...,Xn) > 0 for j=1,...,q
The functions f, G, H are coded in the Scilab function evalfunc and the set of matrices Xi's in the list X (i.e.
X=list(X1,...,Xn)).
The function evalfun must return in the list LME the matrices
G1(X),...,Gp(X) (i.e. LME(i)=Gi(X1,...,Xn), i=1,...,p).
evalfun must return in the list LMI the matrices
H1(X0),...,Hq(X) (i.e. LMI(j)=Hj(X1,...,Xn), j=1,...,q).
evalfun must return in OBJ the value of f(X) (i.e. OBJ=f(X1,...,Xn)).
lmisolver returns in XLISTF, a list of real matrices,
i. e. XLIST=list(X1,X2,..,Xn) where the Xi's solve the LMI
problem:
Defining Y,Z and cost by:
[Y,Z,cost]=evalfunc(XLIST), Y is a list of zero matrices,
Y=list(Y1,...,Yp), Y1=0, Y2=0, ..., Yp=0.
Z is a list of square symmetric matrices,
Z=list(Z1,...,Zq) , which are semi positive definite
Z1>0, Z2>0, ..., Zq>0 (i.e. spec(Z(j)) > 0),
cost is minimized.
lmisolver can also solve LMI problems in which the Xi's are not matrices but lists of matrices. More details are given in the
documentation of LMITOOL.
EXAMPLE :
//Find diagonal matrix X (i.e. X=diag(diag(X), p=1) such that
//A1'*X+X*A1+Q1 < 0, A2'*X+X*A2+Q2 < 0 (q=2) and trace(X) is maximized
n=2;A1=rand(n,n);A2=rand(n,n);
Xs=diag(1:n);Q1=-(A1'*Xs+Xs*A1+0.1*eye());
Q2=-(A2'*Xs+Xs*A2+0.2*eye());
deff('[LME,LMI,OBJ]=evalf(Xlist)','X=Xlist(1),LME=X-diag(diag(X));...
LMI=list(-(A1''*X+X*A1+Q1),-(A2''*X+X*A2+Q2)),OBJ= -sum(diag(X)) ');
X=lmisolver(list(zeros(A1)),evalf);X=X(1)
[Y,Z,c]=evalf(X)
See Also :
lmitool X
6.16 lmitool tool for solving linear matrix inequations
CALLING SEQUENCE :
lmitool()
lmitool(filename)
txt=lmitool(probname,varlist,datalist)
PARAMETERS :
- filename
: a string referring to a .sci function
- probname
: a string containing the name of the problem
- varlist
: a string containing the names of the unknown matrices
(separated by commas if there are more than one)
- datalist
: a string containing the names of data matrices (separated
by commas if there are more than one)
- txt
: a string providing information on what the user should do
next
DESCRIPTION :
lmitool() or lmitool(filename) is used to define
interactively a LMI problem.
In the non interactive mode,
txt=lmitool(probname,varlist,datalist) generates a file in the current directory. The name of this file
is obtained by adding .sci to the end of probname.
This file is the skeleton of a solver function and the corresponding
evaluation function needed by lmisolver.
See Also :
lmisolver X
6.17 ode ordinary differential equation solver
CALLING SEQUENCE :
[y]=ode(y0,t0,t,f)
[y,w,iw]=ode([type],y0,t0,t [,rtol [,atol]],f [,jac] [,w,iw])
[y,rd,w,iw]=ode('root',y0,t0,t [,rtol [,atol]],f [,jac],ng,g [,w,iw])
[y]=ode('discrete',y0,k0,kvect,f)
PARAMETERS :
- y0
: real vector or matrix (initial conditions).
- t0,
: real scalar (initial time).
- t
: real vector (times at which the solution is computed).
- f
: external (function or character string or list).
- type
: one of the following character string:
'adams' 'stiff' 'rk' 'rkf' 'fix' 'discrete' 'roots'
- rtol,atol
: real constants or real vector of the same size as y.
- jac
: external (function or character string or list).
- w,iw
: real vectors.
- ng
: integer.
- g
: external (function or character string or list).
- k0
: integer (initial time).
kvect
: integer vector.
DESCRIPTION :
ode is the standard function for solving explicit ODE systems
defined by:
dy/dt=f(t,y) , y(t0)=y0.
It is an interface to various solvers, in particular to ODEPACK.
The type of problem solved and the method used depends on the value of
the first optional argument type which can be one of the
following strings:
- <not given>:
lsoda solver of package ODEPACK is called by default. It automatically
selects between nonstiff predictor-corrector Adams method and
stiff Backward Differentiation Formula (BDF) method. It uses
nonstiff method initially and dynamically
monitors data in order to decide which method to use.
- "adams":
This is for nonstiff problems. lsode solver of package ODEPACK is called and it uses the Adams method.
- "stiff":
This is for stiff problems. lsode solver of package ODEPACK is called and it uses the BDF method.
- "rk":
Adaptive Runge-Kutta of order 4 (RK4) method.
- "rkf":
The Shampine and Watts program based on Fehlberg's
Runge-Kutta pair of order 4 and 5 (RKF45) method is
used. This is for non-stiff and mildly stiff problems when
derivative evaluations are inexpensive. This method should
generally not be used when the user is demanding high accuracy.
- "fix":
Same solver as 'rkf', but the
user interface is very simple, i.e. only rtol and atol parameters can be passed to the solver. This is the
simplest method to try.
- "root":
ODE solver with rootfinding capabilities.
The lsodar solver of package ODEPACK is used. It is a variant of the
lab{lsoda}solver where it finds the roots of a given
vector function.
- "discrete":
Discrete time simulation.
In this help we only describe the use of ode for standard
explicit ODE systems.
Other helps are available for root finding (ode_root help) and
discrete time simulation (ode_discrete help).
The simplest call of ode is:
y=ode(y0,t0,t,f) where y0 is the vector of initial conditions, t0 is the
initial time, t is the vector of times at which the solution
y is computed and y is the solution vector
y=[y(t(1)),y(t(2)),...].
The input f to ode is an external i.e. a function with
specified syntax, or the name a Fortran subroutine or a C function
(character string) with specified calling sequence or a list.
If f is a function, its syntax must be as follows:
ydot = f(t,y)
where t is a real scalar (time) and y a real vector (state).
This function is the RHS of the differential equation dy/dt=f(t,y).
If f is a character string, it refers to the name of a Fortran
subroutine or a C function, i.e. if ode(y0,t0,t,'fex') is the
command, then the subroutine fex is called.
This routine must have the following calling
sequence:f(n,t,y,ydot). It can be dynamically linked to Scilab
by the link function. Examples of such programs can be seen in
the files SCIDIR/routines/default/README and
SCIDIR/routines/default/Ex-ode.f.
The f argument can also be a list:
if ode(y0,t0,t,lst) is the command, then lst must be a list with the following
structure:
lst=list(f,u1,u2,...un)
where f is a function with syntax:
ydot = f(t,y,u1,u2,...,un)
this allows to use parameters as the arguments of f.
The function f can return a p X q matrix instead of a vector.
With this matrix notation, we solve the n=p+q ODE's
system Y/dt=F(t,Y) where Y is a p X q matrix.
Then initial conditions, Y0, must also be
a p X q matrix and the result of ode is the
p X q(T+1) matrix [Y(t_0),Y(t_1),...,Y(t_T)].
Optional parameters can be given for the error of the solution:
rtol and atol are threshold for relative and absolute estimated errors.
The estimated error on y(i) is rtol(i)*abs(y(i))+*atol(i) and integration is carried out as far as this error is small
for all components of the state.
If rtol and/or atol is a constant rtol(i) and/or
atol(i) are
set to this constant value. Default values for rtol and atol are respectively rtol=1.d-5 and atol=1.d-7 for most
solvers and rtol=1.d-3 and atol=1.d-4 for 'rfk' and
'fix'.
For stiff problem, it is better to give the Jacobian of the RHS
function as the optional argument jac.
It is an external i.e. a function with
specified syntax, or the name a Fortran subroutine or a C function
(character string) with specified calling sequence or a list.
If jac is function the syntax should be as follows:
J=jac(t,y)
where t is a real scalar (time) and y a real vector (state).
The result matrix J must evaluate to df/dx i.e.
J(k,i) = d fk /d xi with fk = kth component of f.
If jac is a character string it refers to the name of a Fortran
subroutine or a C function, with the following calling sequence:
jac(n,t,y,ml,mu,J,nrpd)). In most cases you have not to refer ml,
mu and nrpd (see source code in
SCIDIR/routines/default/Ex-ode.f for an example) .
If jac is a list the same conventions as for f apply.
Optional arguments w and iw are
vectors for storing information returned by the
integration routine. When these vectors are provided in RHS
of ode the integration re-starts with the same parameters as
in its previous stop.
More options can be given to ODEPACK solvers by using
%ODEOPTIONS variable. See odeoptions help.
EXAMPLE :
// Simple one dimension ODE
// dy/dt=y^2-y sin(t)+cos(t), y(0)=0
deff('[ydot]=f(t,y)','ydot=y^2-y*sin(t)+cos(t)')
y0=0;t0=0;t=0:0.1:%pi;
y=ode(y0,t0,t,f)
plot(t,y)
// Simulation of dx/dt = A x(t) + B u(t) with u(t)=sin(omega*t),
// x0=[1;0] ;
// solution x(t) desired at t=0.1, 0.2, 0.5 ,1.
// A and u function are passed to RHS function in a list.
// B and omega are passed as global variables
deff('[xdot]=linear(t,x,A,u)','xdot=A*x+B*u(t)')
deff('[ut]=u(t)','ut=sin(omega*t)')
A=[1 1;0 2];B=[1;1];omega=5;
ode([1;0],0,[0.1,0.2,0.5,1],list(linear,A,u))
//
// Matrix notation
// Integration of the Riccati differential equation
// Xdot=A'*X + X*A - X'*B*X + C , X(0)=Identity
// Solution at t=[1,2]
deff('[Xdot]=ric(t,X)','Xdot=A''*X+X*A-X''*B*X+C')
A=[1,1;0,2]; B=[1,0;0,1]; C=[1,0;0,1];
X=ode(eye(A),0,t,ric)
//
// Computation of exp(A)
A=[1,1;0,2];
deff('[xdot]=f(t,x)','xdot=A*x');
ode(eye(A),0,1,f)
ode('adams',eye(A),0,1,f)
// with stiff matrix, Jacobian given
A=[10,0;0,-1];
deff('[J]=Jacobian(t,y)','J=A')
ode('stiff',[0;1],0,1,f,Jacobian)
See Also :
odedc X, dassl X, csim X, expm X, external X, impl X, ltitr X, rtitr X
6.18 ode_discrete ordinary differential equation solver, discrete time simulation
CALLING SEQUENCE :
[y]=ode('discrete',y0,k0,kvect,f)
PARAMETERS :
- y0
: real vector or matrix (initial conditions).
- t0,
: real scalar (initial time).
- f
: external i.e. function or character string or list.
k0
: integer (initial time).
kvect
: integer vector.
DESCRIPTION :
With this syntax (first argument equal to 'discrete')
ode computes recursively y(k+1)=f(k,y(k))
from an initial state
y(k0) and returns y(k) for k in kvect.
kvect(1) must be greater or equal to k0.
Other arguments and other options are the same as for ode, see the
ode help.
EXAMPLE :
y1=[1;2;3];deff('yp=a_function(k,y)','yp=A*y+B*u(k)')
A=diag([0.2,0.5,0.9]);B=[1;1;1];u=1:10;n=5;
y=ode('discrete',y1,1,1:n,a_function);
y(:,2)-(A*y1+B*u(1))
// Now y evaluates at [y3,y5,y7,y9]
y=ode('discrete',y1,1,3:2:9,a_function);
See Also :
ode X
6.19 ode_root ordinary differential equation solver with root finding
CALLING SEQUENCE :
[y,rd,w,iw]=ode('root',y0,t0,t [,rtol [,atol]],f [,jac],ng,g [,w,iw])
PARAMETERS :
- y0
: real vector or matrix (initial conditions).
- t0,
: real scalar (initial time).
- t
: real vector (times at which the solution is computed).
- f
: external i.e. function or character string or list.
- rtol,atol
: real constants or real vector of the same size as y.
as y.
- jac
: external i.e. function or character string or list.
- w,iw
: real vectors.
- ng
: integer.
- g
: external i.e. function or character string or list.
DESCRIPTION :
With this syntax (first argument equal to 'root')
ode computes the solution
of the differential equation dy/dt=f(t,y) until the state
y(t) crosses the surface g(t,y)=0.
g should give the equation of the surface.
It is an external i.e. a function with
specified syntax, or the name a Fortran subroutine or a C function
(character string) with specified calling sequence or a list.
If g is function the syntax should be as follows:
z=g(t,y)
where t is a real scalar (time) and y a real vector (state).
It returns a vector of size ng which corresponds to
the ng constraints.
If g is a character string it refers to the name of a Fortran
subroutine or a C function, with the following calling sequence:
g(n,t,y,ng,gout).
where ng is the number of constraints and
gout is the value of g (output of the program).
If is a list the same conventions as for f apply (see ode
help).
Ouput rd is a 1xk vector. The first entry contains the stopping
time. Other entries indicate which components of g have changed
sign. k larger than two indicates that more than one surface
(k-1) has been simultaneously traversed.
Other arguments and other options are the same as for ode, see the
ode help.
EXAMPLE :
// Integration of the differential equation
// dy/dt=y , y(0)=1, and finds the minimum time t such that y(t)=2
deff('[ydot]=f(t,y)','ydot=y')
deff('[z]=g(t,y)','z=y-2')
y0=1;ng=1;
[y,rd]=ode('roots',y0,0,2,f,ng,g)
See Also :
dasrt X, ode X
6.20 odedc discrete/continuous ode solver
CALLING SEQUENCE :
yt=odedc(y0,nd,stdel,t0,t,f)
PARAMETERS :
- y0
: real column vector (initial conditions),y0=[y0c;y0d] where
y0d has nd components.
- nd
: integer, dimension of y0d
- stdel
: real vector with one or two entries, stdel=[h, delta] (with
delta=0 as default value).
- t0
: real scalar (initial time).
- t
: real (row) vector, instants where yt is calculated
- f
: Scilab "external" i.e. function or character string or list with
calling sequence: yp=f(t,yc,yd,flag)
DESCRIPTION :
y=odedc([y0c;y0d],nd,[h,delta],t0,t,f) computes the solution of a mixed discrete/continuous system.
The discrete system state yd_k is embedded into a piecewise constant yd(t) time function as follows: yd(t)= yd_k for t in [t_k=delay+k*h,t_(k+1)=delay+(k+1)*h[ (with delay=h*delta).
The simulated equations are now:
dyc/dt=f(t,yc(t),yd(t),0) , for t in [t_k,t_(k+1)[
yc(t0)=y0c
and at instants t_k the discrete variable yd is updated by:
yd(t_k+)=f(yc(t_k-),yd(t_k-),1)
Note that, using the definition of yd(t) the last equation gives
yd_k = f (t_k,yc(t_k-),yd(t_(k-1)),1) (yc is time-continuous: yc(t_k-)=yc(tk))
The calling parameters of f are fixed: ycd=f(t,yc,yd,flag);
this function must return either the derivative of the vector yc if
flag=0 or the update of yd if flag=1.
ycd=dot(yc) must be a vector with same dimension as yc if flag=0 and ycd=update(yd) must be a vector with same
dimension as yd if flag=1
t is a vector of instants where the solution y is computed.
y is the vector y=[y(t(1)),y(t(2)),...].
This function can be called with the same optional parameters as the
ode function (provided nd and stdel are given
in the calling sequence as second and third parameters).
It particular integration flags, tolerances can be set. Optional
parameters can be set by the odeoptions function.
An example for calling an external routine is given in directory
default/fydot2.f External routines can be dynamically linked (see link).
EXAMPLE :
//Linear system with switching input
deff('xdu=phis(t,x,u,flag)','if flag==0 then xdu=A*x+B*u; else xdu=1-u;end');
x0=[1;1];A=[-1,2;-2,-1];B=[1;2];u=0;nu=1;stdel=[1,0];u0=0;t=0:0.05:10;
xu=odedc([x0;u0],nu,stdel,0,t,phis);x=xu(1:2,:);u=xu(3,:);
nx=2;
plot2d1('onn',t',x',[1:nx],'161');
plot2d2('onn',t',u',[nx+1:nx+nu],'000');
//Fortran external( see fydot2.f):
norm(xu-odedc([x0;u0],nu,stdel,0,t,'phis'),1)
//Sampled feedback
//
// | xcdot=fc(t,xc,u)
// (system) |
// | y=hc(t,xc)
//
//
// | xd+=fd(xd,y)
// (feedback) |
// | u=hd(t,xd)
//
deff('xcd=f(t,xc,xd,iflag)',...
['if iflag==0 then '
' xcd=fc(t,xc,e(t)-hd(t,xd));'
'else '
' xcd=fd(xd,hc(t,xc));'
'end']);
A=[-10,2,3;4,-10,6;7,8,-10];B=[1;1;1];C=[1,1,1];
Ad=[1/2,1;0,1/20];Bd=[1;1];Cd=[1,1];
deff('st=e(t)','st=sin(3*t)')
deff('xdot=fc(t,x,u)','xdot=A*x+B*u')
deff('y=hc(t,x)','y=C*x')
deff('xp=fd(x,y)','xp=Ad*x + Bd*y')
deff('u=hd(t,x)','u=Cd*x')
h=0.1;t0=0;t=0:0.1:2;
x0c=[0;0;0];x0d=[0;0];nd=2;
xcd=odedc([x0c;x0d],nd,h,t0,t,f);
norm(xcd-odedc([x0c;x0d],nd,h,t0,t,'fcd1')) // Fast calculation (see fydot2.f)
plot2d([t',t',t'],xcd(1:3,:)');
xset("window",2);plot2d2("gnn",[t',t'],xcd(4:5,:)');
xset("window",0);
See Also :
ode X, odeoptions X, csim X, external X
6.21 odeoptions sets options for ode solvers
CALLING SEQUENCE :
odeoptions()
DESCRIPTION :
This functions interactively displays a command which
should be executed to set various options of ode solvers.
The global variable %ODEOPTIONS sets the options.
CAUTION: the ode function checks if this variable
exists and in this case uses it. For using default
values you should clear this variable. Note that odeoptions does not create this variable. To create it you must execute
the command line displayed by odeoptions.
The variable %ODEOPTIONS is a vector with the following elements:
[itask,tcrit,h0,hmax,hmin,jactyp,mxstep,maxordn,maxords,ixpr, ml,mu]
The default value is: [1,0,0,%inf,0,2,500,12,5,0,-1,-1].
The meaning of the elements is described below.
itask 1 : normal computation at specified times
2 : computation at mesh points (given in first row of output of ode)
3 : one step at one internal mesh point and return
4 : normal computation without overshooting tcrit 5 : one step, without passing tcrit, and return
tcrit assumes itask equals 4 or 5, described above
h0 first step tried
hmax max step size
hmin min step size
jactype 0 : functional iterations, no jacobian used ('adams' or
'stiff' only)
1 : user-supplied full jacobian
2 : internally generated full jacobian
3 : internally generated diagonal jacobian ('adams' or
'stiff' only))
4 : user-supplied banded jacobian (see ml and mu below)
5 : internally generated banded jacobian (see ml and mu below)
maxordn maximum non-stiff order allowed, at most 12
maxords maximum stiff order allowed, at most 5
ixpr print level, 0 or 1
ml, mu If jactype equals 4 or 5, ml and mu are the lower
and upper half-banwidths of the banded jacobian: the band is the
i,j's with i-ml <= j <= ny-1
If jactype equals 4 the jacobian function must return
a matrix J which is ml+mu+1 x ny (where ny=dim of y in ydot=f(t,y))
such that column 1 of J is made of mu zeros followed by
df1/dy1, df2/dy1, df3/dy1,... (1+ml possibly non-zero entries)
column 2 is made of mu-1 zeros followed by df1/dx2, df2/dx2,etc'
See Also :
ode X
6.22 optim non-linear optimization routine
CALLING SEQUENCE :
[f,xopt]=optim(costf,x0)
[f,[xopt,[gradopt,[work]]]]=optim(costf,[contr],x0,['algo'],[df0,[mem]],
[work],[stop],['in'],[imp=iflag])
PARAMETERS :
- costf
: external, i.e Scilab function or string (costf is the cost function: see below its
calling sequence (Scilab or Fortran)).
- x0
: real vector (initial value of variable to be minimized).
- f
: value of optimal cost (f=costf(xopt))
- xopt
: best value of x found.
- contr
: 'b',binf,bsup with binf and bsup real vectors with same
dimension as x0. binf and bsup are lower and upper bounds on x.
- "algo"
: 'qn' or 'gc' or 'nd' . This string stands for quasi-Newton (default),
conjugate gradient or non-differentiable respectively.
Note that 'nd' does not accept bounds on x ).
- df0
: real scalar. Guessed decreasing of f at first iteration.
(df0=1 is the default value).
- mem :
integer, number of variables used to approximate the
Hessian, (algo='gc' or 'nd'). Default value is around 6.
- stop
: sequence of optional parameters controlling the
convergence of the algorithm.
top=top-1stop= 'ar',nap, [iter [,epsg [,epsf [,epsx]]]]
1
- "ar"
: reserved keyword for stopping rule selection defined as follows:
- nap
: maximum number of calls to costf allowed.
- iter
: maximum number of iterations allowed.
- epsg
: threshold on gradient norm.
- epsf
: threshold controlling decreasing of f
- epsx
: threshold controlling variation of x.
This vector (possibly matrix) of same size as x0 can be used
to scale x.
0
- "in"
: reserved keyword for initialization of parameters
used when costf in given as a Fortran routine (see below).
- "imp=iflag"
: named argument used to set the trace mode. iflag=0 nothing
(execpt errors) is reported, iflag=1 initial and final reports,
iflag=2 adds a report per iteration, iflag>2 add reports on
linear search. Warning, most of these reports are written on the
Scilab standard output.
- gradopt
: gradient of costf at xopt
- work
: working array for hot restart for quasi-Newton method.
This array is automatically initialized by optim when
optim is invoked. It can be used as input parameter to
speed-up the calculations.
DESCRIPTION :
Non-linear optimization routine for programs without constraints or with bound constraints:
min costf(x) w.r.t x.
costf is an "external" i.e function, or list or Fortran routine (see "external").
This external must return f (costf(x)) and g (gradient of costf)
given x.
If costf is a function, the calling sequence for costf must be:
[f,g,ind]=costf(x,ind).
Here, costf is a function which returns f, value (real number) of
cost function at x, and g, gradient vector of cost function at x.
The variable ind is used by optim and is described below.
If ind=2 (resp. 3, 4), costf must provide f (resp. g, f and g).
If ind=1 nothing is computed (used for display purposes only).
On output, ind<0 means that f cannot be evaluated at x and
ind=0 interrupts the optimization.
If costf is a character string, it refers to the name of a Fortran routine
which must be linked to Scilab (see examples in the routines foptim.f and e.g. genros.f in the directory SCIDIR/default)
Dynamic link of Fortran routine is also possible (help link).
Here, the generic calling sequence for the Fortran subroutine is:
function costf(ind,n,x,f,g,ti,tr,td)
ind has the same meaning as above if set to 0,1,2 but the
values ind=10 and ind=11 are now valid.
These values are used for initializations (see below).
n is the dimension of x, x is an n vector,
ti,tr,td are working arrays.
The Fortran function costf must return f and the vector g,
given x, ind, n, ti, tr, td.
If costf is given as a Fortran routine, it is possible to initialize
parameters or to send Scilab variables to this routine.
This facility is managed by the parameter 'in.
If the string 'in' is present, initialization is done by Fortran:
optim makes two calls to the Fortran function costf, once with
ind=10 and once with ind=11. In this case, for ind=10,
costf must set the dimensions nti, ntr, ntd of ti, tr, td in the common/nird/nti, ntr, ntd and, for ind=11, costf must initialize the vectors ti , tr, td (integer vector, real vector,
double precision vector respectively).
In the calling sequence of optim, the string 'in' can be
replaced by 'ti', valti ,'td' , valtd. Then, the Fortran function
costf(ind, x, f, g, ti, tr, td) is evaluated with ti=valti and td=valtd whatever the value of ind.
Thus, the Scilab variables valti and valtd (integer vector and
real vector) are sent to the Fortran function costf.
It is also possible to save the content of of the working arrays
ti and td. This is possible by adding the strings 'si' and/or 'sd'
at the ned of the calling sequence of optim.
Then, the output variables must be: [f,[x,[g],[to]]],[ti],[td]].
EXAMPLES :
xref=[1;2;3];x0=[1;-1;1]
deff('[f,g,ind]=cost(x,ind)','f=0.5*norm(x-xref)^2,g=x-xref');
[f,xopt]=optim(cost,x0) //Simplest call
[f,xopt,gopt]=optim(cost,x0,'gc') // By conjugate gradient
[f,xopt,gopt]=optim(cost,x0,'nd') //Seen as non differentiable
[f,xopt,gopt]=optim(cost,'b',[-1;0;2],[0.5;1;4],x0) // Bounds on x
[f,xopt,gopt]=optim(cost,'b',[-1;0;2],[0.5;1;4],x0,'gc') // Bounds on x
[f,xopt,gopt]=optim(cost,'b',[-1;0;2],[0.5;1;4],x0,'gc','ar',3)
// Here, 3 calls to cost are allowed.
// Now calling the Fortran subroutine "genros" in SCIDIR/default/Ex-optim.f
// See also the link function for dynamically linking an objective function
[f,xopt,gopt]=optim('genros',[1;2;3]) //Rosenbrock's function
See Also :
external X, quapro X, linpro X
6.23 quapro linear quadratic programming solver
CALLING SEQUENCE :
[x,lagr,f]=quapro(Q,p,C,b [,x0])
[x,lagr,f]=quapro(Q,p,C,b,ci,cs [,x0])
[x,lagr,f]=quapro(Q,p,C,b,ci,cs,mi [,x0])
[x,lagr,f]=quapro(Q,p,C,b,ci,cs,mi,x0 [,imp])
PARAMETERS :
- Q
: real symmetric matrix (dimension n x n).
- p
: real (column) vector (dimension n)
- C
: real matrix (dimension (mi + md) x n)
(If no constraints are given, you can set C = [])
- b
: RHS vector (dimension 1 x (mi + md))
- ci
: (column) vector of lower-bounds (dimension 1 x n).
If there are no lower bound constraints, put ci = [].
If some components of x are bounded from below,
set the other (unconstrained) values of ci to a very
large negative number (e.g. ci(j) = -(% eps)^(-1).
- cs
: (column) vector of upper-bounds. (Same remarks as above).
- mi
: number of equality constraints (i.e. C(1:mi,:)*x = b(1:mi))
- x0
: either an initial guess for x or one of the character strings 'v' or 'g'.
If x0='v' the calculated initial feasible point is a vertex.
If x0='g' the calculated initial feasible point is arbitrary.
- imp
: verbose option (optional parameter) (Try imp=7,8,...).
- x
: optimal solution found.
- f
: optimal value of the cost function (i.e. f=p'*x).
- lagr
: vector of Lagrange multipliers.
If lower and upper-bounds ci,cs are provided, lagr has
n + mi + md components and lagr(1:n) is the Lagrange
vector associated with the bound constraints and
lagr (n+1 : n + mi + md) is the Lagrange vector associated
with the linear constraints.
(If an upper-bound (resp. lower-bound) constraint i is active
lagr(i) is > 0 (resp. <0).
If no bounds are provided, lagr has only mi + md components.
DESCRIPTION :
[x,lagr,f]=quapro(Q,p,C,b [,x0])
Minimize 0.5*x'*Q*x + p'*x
under the constraint
C*x <= b
[x,lagr,f]=quapro(Q,p,C,b,ci,cs [,x0])
Minimize 0.5*x'*Q*x + p'*x
under the constraints
C*x <= b ci <= x <= cs
[x,lagr,f]=quapro(Q,p,C,b,ci,cs,mi [,x0])
Minimize 0.5*x'*Q*x + p'*x
under the constraints
C(j,:) x = b(j), j=1,...,mi
C(j,:) x <= b(j), j=mi+1,...,mi+md
ci <= x <= cs
If no initial point is given the
program computes a feasible initial point
which is a vertex of the region of feasible points if
x0='v'.
If x0='g', the program computes a feasible initial
point which is not necessarily a vertex. This mode is
advisable when the quadratic form is positive
definite and there are few constraints in
the problem or when there are large bounds
on the variables that are just security bounds and
very likely not active at the optimal solution.
Note that Q is not necessarily non-negative, i.e.
Q may have negative eigenvalues.
EXAMPLE :
//Find x in R^6 such that:
//C1*x = b1 (3 equality constraints i.e mi=3)
C1= [1,-1,1,0,3,1;
-1,0,-3,-4,5,6;
2,5,3,0,1,0];
b1=[1;2;3];
//C2*x <= b2 (2 inequality constraints)
C2=[0,1,0,1,2,-1;
-1,0,2,1,1,0];
b2=[-1;2.5];
//with x between ci and cs:
ci=[-1000;-10000;0;-1000;-1000;-1000];cs=[10000;100;1.5;100;100;1000];
//and minimize 0.5*x'*Q*x + p'*x with
p=[1;2;3;4;5;6]; Q=eye(6,6);
//No initial point is given;
C=[C1;C2] ; //
b=[b1;b2] ; //
mi=3;
[x,lagr,f]=quapro(Q,p,C,b,ci,cs,mi)
//Only linear constraints (1 to 4) are active (lagr(1:6)=0):
[x,lagr,f]=quapro(Q,p,C,b,[],[],mi) //Same result as above
See Also :
linpro X
Author :
E. Casas, C. Pola Mendez
6.24 semidef semidefinite programming
CALLING SEQUENCE :
[x,Z,ul,info]=semidef(x0,Z0,F,blck_szs,c,options)
PARAMETERS :
- x0
: m x 1 real column vector (must be strictly primal feasible, see below)
- Z0
: L x 1 real vector (compressed form of a strictly feasible dual
matrix, see below)
- F
: L x (m+1) real matrix
- blck_szs
: p x 2 integer matrix (sizes of the blocks) defining the dimensions
of the (square) diagonal blocks size(Fi(j)=blck_szs(j) j=1,...,m+1.
- c
: m x 1 real vector
- options
: row vector with five entries [nu,abstol,reltol,0,maxiters]
- ul
: row vector with two entries
DESCRIPTION :
[x,Z,ul,info]=semidef(x0,Z0,F,blck_szs,c,options) solves semidefinite program:
minimize c'*x
subject to F_0 + x_1*F_1 + ... + x_m*F_m >= 0
and its dual
maximize -Tr F_0 Z
subject to Tr F_i Z = c_i, i=1,...,m
Z >= 0
exploiting block structure in the matrices F_i.
It interfaces L. Vandenberghe and S. Boyd sp.c program.
The Fi's matrices are stored columnwise in F in
compressed format: if F_i^j, i=0,..,m, j=1,...,L denote the jth
(symmetric) diagonal block of F_i, then
[ pack(F_0^1) pack(F_1^1) ... pack(F_m^1) ]
[ pack(F_0^2) pack(F_1^2) ... pack(F_m^2) ]
F= [ ... ... ... ]
[ pack(F_0^L) pack(F_1^L) ... pack(F_m^L) ]
where pack(M), for symmetric M, is the vector
[M(1,1);M(1,2);...;M(1,n);M(2,2);M(2,3);...;M(2,n);...;M(n,n)] (obtained by scanning columnwise the lower triangular part of M).
blck_szs gives the size of block j, ie,
size(F_i^j)=blck_szs(j).
Z is a block diagonal matrix with L blocks Z^0, ..., Z^{L-1}.
Z^j has size blck_szs[j] times blck_szs[j].
Every block is stored using packed storage of the lower triangular part.
The 2 vector ul contains the primal objective value c'*x and the dual objective value -Tr F_0*Z.
The entries of options are respectively:
nu = a real parameter which ntrols the rate of convergence.
abstol = absolute tolerance.
reltol = relative tolerance (has a special meaning when negative).
tv target value, only referenced if reltol < 0.
iters = on entry: maximum number of iterations >= 0, on exit:
the number of iterations taken.
info returns 1 if maxiters exceeded, 2 if absolute accuracy
is reached, 3 if relative accuracy is reached, 4 if target value is
reached, 5 if target value is not achievable; negative values
indicate errors.
Convergence criterion:
(1) maxiters is exceeded
(2) duality gap is less than abstol
(3) primal and dual objective are both positive and
duality gap is less than (reltol * dual objective)
or primal and dual objective are both negative and
duality gap is less than (reltol * minus the primal objective)
(4) reltol is negative and
primal objective is less than tv or dual objective is greater
than tv
EXAMPLE :
F0=[2,1,0,0;
1,2,0,0;
0,0,3,1
0,0,1,3];
F1=[1,2,0,0;
2,1,0,0;
0,0,1,3;
0,0,3,1]
F2=[2,2,0,0;
2,2,0,0;
0,0,3,4;
0,0,4,4];
blck_szs=[2,2];
F01=F0(1:2,1:2);F02=F0(3:4,3:4);
F11=F1(1:2,1:2);F12=F1(3:4,3:4);
F21=F2(1:2,1:2);F22=F2(3:4,3:4);
x0=[0;0]
Z0=2*F0;
Z01=Z0(1:2,1:2);Z02=Z0(3:4,3:4);
FF=[[F01(:);F02(:)],[F11(:);F12(:)],[F21(:);F22(:)]]
ZZ0=[[Z01(:);Z02(:)]];
c=[trace(F1*Z0);trace(F2*Z0)];
options=[10,1.d-10,1.d-10,0,50];
[x,Z,ul,info]=semidef(x0,pack(ZZ0),pack(FF),blck_szs,c,options)
w=vec2list(unpack(Z,blck_szs),[blck_szs;blck_szs]);Z=sysdiag(w(1),w(2))
c'*x+trace(F0*Z)
spec(F0+F1*x(1)+F2*x(2))
trace(F1*Z)-c(1)
trace(F2*Z)-c(2)