Similarly to scicos.c's 'jacpsol' and 'psol', doing preconditioner factorization in the evaluation function instead of the solving one.
Change-Id: I332594a2892348ad4f7f102852edb0697429cb9a
<?xml version="1.0" encoding="UTF-8"?>
<!--
* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2013 - Scilab Enterprises
+ * Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
*
* This file must be used under the terms of the CeCILL.
* This source file is licensed as described in the file COPYING, which
<listitem>
<para>
is either <literal>y0</literal> (<literal>ydot0</literal> is
- estimated by <literal>dassl</literal> with zero as first estimate)
+ estimated by <literal>daskr</literal> with zero as first estimate)
or the matrix <literal>[y0 ydot0]</literal>.
- <literal>g(t,y0,ydot0)</literal> must be equal to zero. If you only
- know an estimate of <literal>ydot0</literal> set
+ <literal>g(t, y0, ydot0)</literal> must be equal to zero.
+ If you only know an estimate of <literal>ydot0</literal>, set
<literal>info(7)=1</literal>.
</para>
<variablelist>
<term>ydot0</term>
<listitem>
<para>a real column vector of the time derivative of
- <literal>y</literal> at <literal>t0</literal> (may be an
- estimate).
+ <literal>y</literal> at <literal>t0</literal> (may be an estimate).
</para>
</listitem>
</varlistentry>
<term>t</term>
<listitem>
<para>a real scalar or vector. Gives instants for which you want the
- solution. Note that you can get solution at each dassl's step point
+ solution. Note that you can get solution at each daskr's step point
by setting <literal>info(2)=1</literal>.
</para>
</listitem>
</varlistentry>
<varlistentry>
- <term>nn</term>
- <listitem>
- <para>
- a vector with two entries <literal>[times num]</literal>
- <literal>times</literal> is the value of the time at which the
- surface is crossed, <literal>num</literal> is the number of the
- crossed surface.
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
<term>atol, rtol</term>
<listitem>
<para>real scalars or column vectors of same size as
- <literal>y</literal>. <literal>atol, rtol</literal> give respectively
- absolute and relative error tolerances of solution. If vectors the
- tolerances are specified for each component of
+ <literal>y</literal>. <literal>atol</literal> and <literal>rtol</literal> give respectively
+ absolute and relative error tolerances of solution.
+ If vectors, then the tolerances are specified for each component of
<literal>y</literal>.
</para>
</listitem>
<term>res</term>
<listitem>
<para>
- <link linkend="external">external</link> (function or list or string). Computes the value of
- <literal>g(t,y,ydot)</literal>. It may be :
+ <link linkend="external">external</link> (function, list or string). Computes the value of
+ <literal>g(t, y, ydot)</literal>. It may be :
</para>
<itemizedlist>
<listitem>
<para>A Scilab function.</para>
<para>Its calling sequence must be
- <literal>[r,ires]=res(t,y,ydot)</literal> and
- <literal>res</literal> must return the residual
- <literal>r=g(t,y,ydot)</literal> and error flag
+ <literal>[r, ires] = res(t, y, ydot)</literal> and must return the residual
+ <literal>r = g(t, y, ydot)</literal> and error flag
<literal>ires</literal>. <literal>ires = 0</literal> if
<literal>res</literal> succeeds to compute <literal>r</literal>,
- <literal>=-1</literal> if residual is locally not defined for
- <literal>(t,y,ydot)</literal>, <literal>=-2</literal> if
+ <literal>ires = -1</literal> if residual is locally not defined for
+ <literal>(t, y, ydot)</literal>, <literal>ires = -2</literal> if
parameters are out of admissible range.
</para>
</listitem>
<listitem>
<para>A list.</para>
- <para>This form allows to pass parameters other than t,y,ydot to
- the function. It must be as follows:
+ <para>This form allows to pass parameters other than t, y, ydot to the function.
+ It must be as follows:
</para>
<programlisting role="no-scilab-exec"><![CDATA[
-list(res,x1,x2,...)
+list(res, x1, x2, ...)
]]></programlisting>
<para>where the calling sequence of the function
<literal>res</literal> is now
</para>
<programlisting role="no-scilab-exec"><![CDATA[
-r=res(t,y,ydot,x1,x2,...)
+r = res(t, y, ydot, x1, x2, ...)
]]></programlisting>
<para>
<literal>res</literal> still returns
- <literal>r=g(t,y,ydot)</literal> as a function of
- <literal>(t,y,ydot,x1,x2,...)</literal>.
+ <literal>r = g(t, y, ydot)</literal> as a function of
+ <literal>(t, y, ydot, x1, x2, ...)</literal>.
</para>
<para>Warning: this form must not be used if there is no extra
- argument to pass to <literal>the function.</literal>
+ argument to pass to the function.
</para>
</listitem>
<listitem>
<para>A string.</para>
- <para>It must refer to the name of a C or Fortran subroutine
+ <para>It must refer to the name of a C function or a Fortran subroutine
linked with Scilab.
</para>
- <para>In C The calling sequence must be:</para>
+ <para>In C, the calling sequence must be:</para>
<programlisting role="no-scilab-exec"><![CDATA[
void res(double *t, double y[], double yd[], double r[],
int *ires, double rpar[], int ipar[])
]]></programlisting>
- <para>In Fortran it must be:</para>
+ <para>In Fortran, it must be:</para>
<programlisting role="no-scilab-exec"><![CDATA[
-subroutine res(t,y,yd,r,ires,rpar,ipar)
-double precision t, y(*),yd(*),r(*),rpar(*)
-integer ires,ipar(*)
+subroutine res(t, y, yd, r, ires, rpar, ipar)
+double precision t, y(*), yd(*),r(*),rpar(*)
+integer ires, ipar(*)
]]></programlisting>
<para>
- The <literal>rpar</literal> and <literal>ipar</literal> arrays must be present but cannot be
- used.
+ The <literal>rpar</literal> and <literal>ipar</literal>
+ arrays must be present but cannot be used.
</para>
</listitem>
</itemizedlist>
<term>jac</term>
<listitem>
<para>
- <link linkend="external">external</link> (function or list or string). Computes the value of
- <literal>dg/dy + cj*dg/dydot</literal> for a given value of parameter
+ <link linkend="external">external</link> (function, list or string).
+ Computes the value of <literal>dg/dy + cj*dg/dydot</literal> for a given value of parameter
<literal>cj</literal>.
</para>
<itemizedlist>
<listitem>
<para>A Scilab function.</para>
<para>Its calling sequence must be
- <literal>r=jac(t,y,ydot,cj)</literal> and the
- <literal>jac</literal> function must return
- <literal>r=dg(t,y,ydot)/dy+cj*dg(t,y,ydot)/dydot</literal> where
+ <literal>r = jac(t, y, ydot, cj)</literal> and must return
+ <literal>r = dg(t, y, ydot)/dy + cj*dg(t, y, ydot)/dydot</literal> where
<literal>cj</literal> is a real scalar.
</para>
</listitem>
<para>A list.</para>
<para>It must be as follows</para>
<programlisting role="no-scilab-exec"><![CDATA[
-list(jac,x1,x2,...)
+list(jac, x1, x2, ...)
]]></programlisting>
<para>where the calling sequence of the function
<literal>jac</literal> is now
</para>
<programlisting role="no-scilab-exec"><![CDATA[
-r=jac(t,y,ydot,cj,x1,x2,...)
+r=jac(t, y, ydot, cj, x1, x2,...)
]]></programlisting>
<para>
<literal>jac</literal> still returns
<literal>dg/dy + cj*dg/dydot</literal> as a function of
- <literal>(t,y,ydot,cj,x1,x2,...)</literal>.
+ <literal>(t, y, ydot, cj, x1, x2, ...)</literal>.
</para>
</listitem>
<listitem>
<para>A character string.</para>
- <para>It must refer to the name of a Fortran subroutine linked
- with scilab
+ <para>It must refer to the name of a C function or a Fortran subroutine linked with Scilab
</para>
- <para>In C The calling sequence must be:</para>
+ <para>In C, the calling sequence must be:</para>
<programlisting role="no-scilab-exec"><![CDATA[
void jac(double *t, double y[], double yd[], double pd[],
double *cj, double rpar[], int ipar[])
]]></programlisting>
- <para>In Fortran it must be:</para>
+ <para>In Fortran, it must be:</para>
<programlisting role="no-scilab-exec"><![CDATA[
-subroutine jac(t,y,yd,pd,cj,rpar,ipar)
-double precision t, y(*),yd(*),pd(*),cj,rpar(*)
+subroutine jac(t, y, yd, pd, cj, rpar, ipar)
+double precision t, y(*), yd(*), pd(*), cj, rpar(*)
integer ipar(*)
]]></programlisting>
</listitem>
<term>surf</term>
<listitem>
<para>
- <link linkend="external">external</link> (function or list or string). Computes the value of
- the column vector <literal>surf(t,y)</literal> with
+ <link linkend="external">external</link> (function, list or string).
+ Computes the value of the column vector <literal>surf(t, y)</literal> with
<literal>ng</literal> components. Each component defines a surface.
It may be defined by:
</para>
<listitem>
<para>A Scilab function.</para>
<para>Its calling sequence must be
- <literal>surf(t,y)</literal>
+ <literal>surf(t, y)</literal>
</para>
</listitem>
<listitem>
<para>A list.</para>
<para>It must be as follows</para>
<programlisting role="no-scilab-exec"><![CDATA[
-list(surf,x1,x2,...)
+list(surf, x1, x2, ...)
]]></programlisting>
<para>where the calling sequence of the function
<literal>surf</literal> is now
</para>
<programlisting role="no-scilab-exec"><![CDATA[
-r=surf(t,y,x1,x2,...)
+r = surf(t, y, x1, x2, ...)
]]></programlisting>
</listitem>
<listitem>
<para>A character string.</para>
- <para>It must refer to the name of a Fortran subroutine linked
- with scilab.
+ <para>It must refer to the name of a C function or a Fortran subroutine linked with Scilab.
</para>
- <para>In C the calling sequence must be:</para>
+ <para>In C, the calling sequence must be:</para>
<programlisting role="no-scilab-exec"><![CDATA[
void surf(int *ny, double *t, double y[], int *ng, double gout[])
]]></programlisting>
- <para>In Fortran it must be:</para>
+ <para>In Fortran, it must be:</para>
<programlisting role="no-scilab-exec"><![CDATA[
-subroutine surf(ny,t,y,ng,gout)
-double precision t, y(*),gout(*)
-integer ny,ng
+subroutine surf(ny, t, y, ng, gout)
+double precision t, y(*), gout(*)
+integer ny, ng
]]></programlisting>
</listitem>
</itemizedlist>
<term>info</term>
<listitem>
<para>
- list which contains <literal>14</literal> elements. Default
- value is <literal>list([],0,[],[],[],0,[],0,[],0,[],[],[],1)</literal>.
+ list which contains <literal>14</literal> elements. Default value is
+ <literal>list([], 0, [], [], [], 0, [], 0, [], 0, [], [], [], 1)</literal>.
</para>
<variablelist>
<varlistentry>
<term>info(1)</term>
<listitem>
<para>real scalar which gives the maximum time for which
- <literal>g</literal> is allowed to be evaluated or an empty
- matrix <literal>[]</literal> if no limits imposed for
- time.
+ <literal>g</literal> is allowed to be evaluated or an empty matrix
+ <literal>[]</literal> if no limits imposed for time.
</para>
</listitem>
</varlistentry>
<term>info(2)</term>
<listitem>
<para>
- flag which indicates if <literal>dassl</literal> returns
- its intermediate computed values (<literal>flag=1</literal>)
+ flag which indicates if <literal>daskr</literal> returns
+ its intermediate computed values (<literal>= 1</literal>)
or only the user specified time point values
- (<literal>flag=0</literal>).
+ (<literal>= 0</literal>).
</para>
</listitem>
</varlistentry>
<para>
<literal>2</literal> components vector which give the
definition <literal>[ml,mu]</literal> of band matrix computed
- by <literal>jac</literal>; <literal>r(i - j + ml + mu + 1,j) =
- "dg(i)/dy(j)+cj*dg(i)/dydot(j)"
- </literal>
- . If <literal>jac</literal> returns a full matrix set
+ by <literal>jac</literal>; <literal>r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)"</literal>.
+ If <literal>jac</literal> returns a full matrix set
<literal>info(3)=[]</literal>.
Treat as dummmy if
<literal>info(8)=1</literal>.
<listitem>
<para>
if ydot0 is set so that
- <literal>g(t0,y0,ydot0)=0</literal> then set
+ <literal>g(t0, y0, ydot0) = 0</literal> then set
<literal>info(7)=[]</literal>. Otherwise, set
<literal>info(7)=[+-1, ..., +-1]</literal>, with
- <literal>info(7)(i)=1</literal> if y(i) is a differential variable and
- <literal>info(7)(i)=-1</literal> if y(i) is an algebraic variable
- (if its derivatives do not appear explicitly in the function g(t,y,ydot)).
+ <literal>info(7)(i) = 1</literal> if y(i) is a differential variable and
+ <literal>info(7)(i) = -1</literal> if y(i) is an algebraic variable
+ (if its derivatives do not appear explicitly in the function g(t, y, ydot)).
</para>
</listitem>
</varlistentry>
<term>info(8)</term>
<listitem>
<para>
- direct / Krylov. Set <literal>info(8)=1</literal> and provide a routine in <literal>psol</literal>
- if you want the solver to use Krylov iterations, else (dassl's direct method) set
+ direct / Krylov method. Set <literal>info(8)=1</literal> and provide a routine in <literal>psol</literal>
+ if you want the solver to use Krylov iterations, else (daskr's direct method) set
<literal>info(8)=0</literal>.
</para>
</listitem>
<para>
Krylov parameters. Treat as dummy argument if you have set
<literal>info(8)=0</literal>. Otherwise, set
+ <literal>info(9)=[]</literal> or
<literal>info(9)=[maxl kmp nrmax epli]</literal>, where:
</para>
<para>
- - maxl = maximum number of iterations in the SPIGMR algorithm (default
- <literal>min(5,neq)</literal>),
+ - maxl = maximum number of iterations in the GMRes algorithm (default
+ <literal>min(5, neq)</literal>),
</para>
<para>
- - kmp = number of vectors on which orthogonalization is done in the SPIGMR algorithm
+ - kmp = number of vectors on which orthogonalization is done in the GMRes algorithm
(default maxl),
</para>
<para>
- - nrmax = maximum number of restarts of the SPIGMR algorithm per nonlinear iteration
+ - nrmax = maximum number of restarts of the GMRes algorithm per nonlinear iteration
(default <literal>5</literal>),
</para>
<para>
- - epli = convergence test constant in SPIGMR algorithm (default <literal>0.05</literal>).
+ - epli = convergence test constant in GMRes algorithm (default <literal>0.05</literal>).
</para>
</listitem>
</varlistentry>
if you wish to control errors locally on all the variables then set
<literal>info(12)=[]</literal>. Otherwise, set
<literal>info(12)=[+-1, ..., +-1]</literal>, with
- <literal>info(12)(i)=1</literal> if y(i) is a differential variable and
- <literal>info(12)(i)=-1</literal> if y(i) is an algebraic variable
- (if its derivatives do not appear explicitly in the function g(t,y,ydot)).
+ <literal>info(12)(i) = 1</literal> if y(i) is a differential variable and
+ <literal>info(12)(i) = -1</literal> if y(i) is an algebraic variable
+ (if its derivatives do not appear explicitly in the function g(t, y, ydot)).
</para>
</listitem>
</varlistentry>
<para>
heuristic parameters. Treat as dummy argument if
<literal>info(7)=[]</literal>. Otherwise, set
+ <literal>info(13)=[]</literal> or
<literal>info(13)=[mxnit mxnj mxnh lsoff stptol epinit]</literal>, where:
</para>
<para>
<literal>6</literal> if <literal>info(8)=0</literal>, <literal>2</literal> otherwise),
</para>
<para>
- - mxnh = maximum number of values of the artificial stepsize parameter H to be tried if info(7) ≠ [] (default
+ - mxnh = maximum number of values of the artificial stepsize parameter <literal>h</literal> to be tried if info(7) ≠ [] (default
<literal>5</literal>),
</para>
<para>
<term>psol</term>
<listitem>
<para>
- <link linkend="external">external</link> (function or list or string). Solves a linear system
- <literal>P*x = b</literal>, with P being the preconditioner that routine <literal>pjac</literal>
- computed beforehand and stocked in wp and iwp.
+ <link linkend="external">external</link> (function, list or string). Solves a linear system
+ <literal>P*x = b</literal>, with P being the factored preconditioner that routine <literal>pjac</literal>
+ computed beforehand and stored in <literal>wp</literal> and <literal>iwp</literal>.
</para>
<itemizedlist>
<listitem>
<para>A Scilab function.</para>
<para>Its calling sequence must be
- <literal>[r, ier] = psol(wp, iwp, b)</literal> and the
- <literal>psol</literal> function must return the solution of that system in
+ <literal>[r, ier] = psol(wp, iwp, b)</literal> and must return the solution of the system in
<literal>r</literal> and an error flag <literal>ier</literal>.
</para>
</listitem>
<listitem>
<para>A list.</para>
- <para>It must be as follows</para>
+ <para>It must be as follows:</para>
<programlisting role="no-scilab-exec"><![CDATA[
-list(psol,x1,x2,...)
+list(psol, x1, x2, ...)
]]></programlisting>
- <para>where the calling sequence of the function
- <literal>psol</literal> is now
+ <para>
+ where the calling sequence of <literal>psol</literal> is now
</para>
<programlisting role="no-scilab-exec"><![CDATA[
-psol(wp,iwp,b,x1,x2,...)
+psol(wp, iwp, b, x1, x2, ...)
]]></programlisting>
<para>
<literal>psol</literal> still returns the solution in <literal>r</literal>.
</listitem>
<listitem>
<para>A character string.</para>
- <para>It must refer to the name of a Fortran subroutine linked
- with scilab
+ <para>It must refer to the name of a C function or a Fortran subroutine linked with Scilab
</para>
- <para>In C The calling sequence must be:</para>
+ <para>In C, the calling sequence must be:</para>
<programlisting role="no-scilab-exec"><![CDATA[
void psol (int*neq, double*t, double*y, double*ydot, double*savr,
double*wk, double*cj, double*wght, double*wp, int*iwp, double*b, double*eplin, int*ier, double*rpar, int*ipar)
]]></programlisting>
- where the arrays <literal>wp</literal> and <literal>iwp</literal> contain matrix elements of preconditioner
- <literal>P</literal> in compressed row format.
- <para>In Fortran it must be:</para>
+ where the arrays <literal>wp</literal> and <literal>iwp</literal> contain matrix elements of LU-factored preconditioner
+ <literal>P</literal>, <literal>wp</literal> being the values and
+ <literal>iwp</literal> the pivots used in the factorization.
+ <para>In Fortran, it must be:</para>
<programlisting role="no-scilab-exec"><![CDATA[
subroutine psol (neq, t, y, ydot, savr, wk, cj, wght,
wp, iwp, b, eplin, ier, rpar, ipar)
-double precision t,y(*),ydot(*),savr(*),wk(*),cj,wght(*),wp(*),
- b(*),eplin,rpar(*)
-integer neq,iwp(*),ier,ipar(*)
+double precision t,y(*), ydot(*), savr(*), wk(*), cj, wght(*), wp(*),
+ b(*), eplin, rpar(*)
+integer neq, iwp(*), ier, ipar(*)
]]></programlisting>
</listitem>
</itemizedlist>
<term>pjac</term>
<listitem>
<para>
- <link linkend="external">external</link> (function or list or string). Computes the value of
+ <link linkend="external">external</link> (function, list or string). Computes the value of
<literal>dg/dy + cj*dg/dydot</literal> for a given value of parameter
- <literal>cj</literal> and LU-factorizes it in double and int arrays.
+ <literal>cj</literal> and LU-factorizes it in two arrays, real and integer.
</para>
<itemizedlist>
<listitem>
<para>A Scilab function.</para>
<para>Its calling sequence must be
<literal>[wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)</literal> and in return,
- the arrays wp and iwp must contain all preconditioner information in compressed sparse row format.
+ the arrays <literal>wp</literal> and <literal>iwp</literal> must contain all factored preconditioner information.
</para>
</listitem>
<listitem>
<para>A list.</para>
<para>It must be as follows</para>
<programlisting role="no-scilab-exec"><![CDATA[
-list(pjac,x1,x2,...)
+list(pjac, x1, x2, ...)
]]></programlisting>
- <para>where the calling sequence of the function
- <literal>pjac</literal> is now
+ <para>
+ where the calling sequence of <literal>pjac</literal> is
</para>
<programlisting role="no-scilab-exec"><![CDATA[
-pjac(neq,t,y,ydot,h,cj,rewt,savr,x1,x2,...)
+pjac(neq, t, y, ydot, h, cj, rewt, savr, x1, x2,...)
]]></programlisting>
<para>
<literal>pjac</literal> still returns factorized
<literal>dg/dy + cj*dg/dydot</literal> as a function of
- <literal>(neq,t,y,ydot,h,cj,rewt,savr,x1,x2,...)</literal>.
+ <literal>(neq, t, y, ydot, h, cj, rewt, savr, x1, x2, ...)</literal>.
</para>
</listitem>
<listitem>
<para>A character string.</para>
- <para>It must refer to the name of a Fortran subroutine linked
- with scilab
+ <para>It must refer to the name of a C function or a Fortran subroutine linked with Scilab
</para>
- <para>In C The calling sequence must be:</para>
+ <para>In C, the calling sequence must be:</para>
<programlisting role="no-scilab-exec"><![CDATA[
void pjac (double*res, int*ires, int*neq, double*t, double*y, double*ydot, double*rewt, double*savr,
double*wk, double*h, double*cj, double*wp, int*iwp, int*ier, double*rpar, int*ipar)
]]></programlisting>
- <para>In Fortran it must be:</para>
+ <para>In Fortran, it must be:</para>
<programlisting role="no-scilab-exec"><![CDATA[
subroutine pjac (res, ires, neq, t, y, ydot, rewt, savr,
wk, h, cj, wp, iwp, ier, rpar, ipar)
<term>hd</term>
<listitem>
<para>
- real vector which allows to store the <literal>dassl</literal>
+ real vector which allows to store the <literal>daskr</literal>
context and to resume integration.
</para>
</listitem>
<term>r</term>
<listitem>
<para>
- real matrix . Each column is the vector <literal>[t;x(t);xdot(t)]</literal> where
- <literal>t</literal> is time index for which the solution had been computed.
+ real matrix. Each column is the vector <literal>[t; x(t); xdot(t)]</literal> where
+ <literal>t</literal> is the time index for which the solution has been computed,
+ <literal>x(t)</literal> is the value of the computed solution,
+ <literal>xdot(t)</literal> is the derivative of the computed solution.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>nn</term>
+ <listitem>
+ <para>
+ a vector with two entries <literal>[times num]</literal>,
+ <literal>times</literal> is the value of the time at which the surface is crossed,
+ <literal>num</literal> is the number of the crossed surface.
</para>
</listitem>
</varlistentry>
</refsection>
<refsection>
<title>Description</title>
- <para>Solution of the implicit differential equation.</para>
+ <para>Solution of the implicit differential equation:</para>
<programlisting role="no-scilab-exec"><![CDATA[
-g(t,y,ydot) = 0
-y(t0) = y0 and ydot(t0) = ydot0
+g(t, y, ydot) = 0
+y(t0) = y0 and ydot(t0) = ydot0
]]></programlisting>
<para>Returns the surface crossing instants and the number of the surface
reached in <literal>nn</literal>.
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
-// dy/dt = ((2*log(y)+8)/t -5)*y, y(1) = 1, 1 <= t <=6
-// g1 = ((2*log(y)+8)/t - 5)*y
+// 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;
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]')
+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] = daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1);
+[yy, nn] = daskr([y0, y0d], t0, t, atol, rtol, res1, ng, gr1);
nn
// Should return nn = [2.4698972 2]
--- /dev/null
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
+ *
+ * This file must be used under the terms of the CeCILL.
+ * This source file is licensed as described in the file COPYING, which
+ * you should have received as part of this distribution.
+ * The terms are also available at
+ * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ *
+ -->
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="daskr" xml:lang="fr">
+ <refnamediv>
+ <refname>daskr</refname>
+ <refpurpose>solveur de DAE avec traversées de zéros</refpurpose>
+ </refnamediv>
+ <refsynopsisdiv>
+ <title>Séquence d'appel</title>
+ <synopsis>[r, nn [, hd]] = daskr(x0, t0, t [, atol [, rtol]], res [, jac], ng, surf [, info [, psol] [, pjac]] [, hd])</synopsis>
+ </refsynopsisdiv>
+ <refsection>
+ <title>Paramètres</title>
+ <variablelist>
+ <varlistentry>
+ <term>x0</term>
+ <listitem>
+ <para>
+ représente soit <literal>y0</literal> (<literal>ydot0</literal> sera
+ estimé par <literal>daskr</literal> avec zéro comme première estimation),
+ soir la matrice <literal>[y0 ydot0]</literal>.
+ <literal>g(t, y0, ydot0)</literal> doit être égal à zéro. Si vous ne connaissez
+ qu'une estimation de <literal>ydot0</literal>, assignez
+ <literal>info(7)=1</literal>.
+ </para>
+ <variablelist>
+ <varlistentry>
+ <term>y0</term>
+ <listitem>
+ <para>vecteur colonne réel des conditions initiales.</para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>ydot0</term>
+ <listitem>
+ <para>vecteur colonne réel de la dérivée en temps de
+ <literal>y</literal> à <literal>t0</literal> (peut être une estimation.
+ </para>
+ </listitem>
+ </varlistentry>
+ </variablelist>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>t0</term>
+ <listitem>
+ <para>réel, temps initial.</para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>t</term>
+ <listitem>
+ <para>réel, scalaire ou vecteur. Temps auxquels la solution est désirée.
+ Notez que vous pouvez obtenir la solution à chaque étape de daskr en fixant
+ <literal>info(2)=1</literal>.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>atol, rtol</term>
+ <listitem>
+ <para>réels scalaires ou vecteurs colonnes de même taille que
+ <literal>y</literal>. <literal>atol</literal> et <literal>rtol</literal> représentent
+ les tolérances d'erreur absolue et relative de la solution.
+ Si ce sont des vecteurs, alors les tolérances sont spécifiées pour chaque composante de
+ <literal>y</literal>.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>res</term>
+ <listitem>
+ <para>
+ <link linkend="external">external</link> (fonction, liste ou chaîne de caractères).
+ Calcule la valeur de
+ <literal>g(t, y, ydot)</literal>. Elle peut être :
+ </para>
+ <itemizedlist>
+ <listitem>
+ <para>Une fonction Scilab.</para>
+ <para>Sa séquence d'appel doit être
+ <literal>[r, ires] = res(t, y, ydot)</literal> et doit retourner le résidu
+ <literal>r = g(t, y, ydot)</literal> et un drapeau d'erreur
+ <literal>ires</literal>. <literal>ires = 0</literal> si
+ <literal>res</literal> a correctement calculé <literal>r</literal>,
+ <literal>ires = -1</literal> si le résidu est localement non défini pour
+ <literal>(t, y, ydot)</literal>, <literal>ires = -2</literal> si
+ des paramètres sont hors du champ admissible.
+ </para>
+ </listitem>
+ <listitem>
+ <para>Une liste.</para>
+ <para>Cette forme permet de passer des paramètres autres que t, y, ydot à la fonction.
+ Elle doit se présenter comme suit :
+ </para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+list(res, x1, x2, ...)
+ ]]></programlisting>
+ <para>où la séquence d'appel de la fonction
+ <literal>res</literal> est maintenant
+ </para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+r = res(t, y, ydot, x1, x2, ...)
+ ]]></programlisting>
+ <para>
+ <literal>res</literal> retourne toujours
+ <literal>r = g(t, y, ydot)</literal> comme fonction de
+ <literal>(t, y, ydot, x1, x2, ...)</literal>.
+ </para>
+ <para>Attention : cette forme ne doit pas être utilisée s'il n'y pas
+ d'argument additionnel à passer à la fonction.
+ </para>
+ </listitem>
+ <listitem>
+ <para>Une chaîne de caractères.</para>
+ <para>Elle doit se référer au nom d'une fonction C ou une routine Fortran
+ reliée à Scilab.
+ </para>
+ <para>En C, la séquence d'appel doit être :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+void res(double *t, double y[], double yd[], double r[],
+ int *ires, double rpar[], int ipar[])
+ ]]></programlisting>
+ <para>En Fortran, elle doit être :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+subroutine res(t, y, yd, r, ires, rpar, ipar)
+double precision t, y(*), yd(*),r(*),rpar(*)
+integer ires, ipar(*)
+ ]]></programlisting>
+ <para>
+ Les tableaux <literal>rpar</literal> et <literal>ipar</literal>
+ doivent être présents mais ne peuvent pas être utilisés.
+ </para>
+ </listitem>
+ </itemizedlist>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>jac</term>
+ <listitem>
+ <para>
+ <link linkend="external">external</link> (fonction, liste ou chaîne de caractères).
+ Calcule la valeur de <literal>dg/dy + cj*dg/dydot</literal> pour une valeur donnée du paramètre
+ <literal>cj</literal>.
+ </para>
+ <itemizedlist>
+ <listitem>
+ <para>Une fonction Scilab.</para>
+ <para>Sa séquence d'appel doit être
+ <literal>r = jac(t, y, ydot, cj)</literal> et doit retourner
+ <literal>r = dg(t, y, ydot)/dy + cj*dg(t, y, ydot)/dydot</literal> où
+ <literal>cj</literal> est un scalaire réel.
+ </para>
+ </listitem>
+ <listitem>
+ <para>Une liste.</para>
+ <para>Elle doit se présenter comme suit :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+list(jac, x1, x2, ...)
+ ]]></programlisting>
+ <para>où la séquence d'appel de la fonction
+ <literal>jac</literal> est désormais
+ </para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+r=jac(t, y, ydot, cj, x1, x2,...)
+ ]]></programlisting>
+ <para>
+ <literal>jac</literal> retourne toujours
+ <literal>dg/dy + cj*dg/dydot</literal> comme fonction de
+ <literal>(t, y, ydot, cj, x1, x2,...)</literal>.
+ </para>
+ </listitem>
+ <listitem>
+ <para>Une chaîne de caractères.</para>
+ <para>Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab.
+ </para>
+ <para>En C, la séquence d'appel doit être :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+void jac(double *t, double y[], double yd[], double pd[],
+ double *cj, double rpar[], int ipar[])
+ ]]></programlisting>
+ <para>En Fortran, elle doit être :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+subroutine jac(t, y, yd, pd, cj, rpar, ipar)
+double precision t, y(*), yd(*), pd(*), cj, rpar(*)
+integer ipar(*)
+ ]]></programlisting>
+ </listitem>
+ </itemizedlist>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>surf</term>
+ <listitem>
+ <para>
+ <link linkend="external">external</link> (fonction, liste ou chaîne de caractères).
+ Calcule la valeur du vecteur colonne <literal>surf(t, y)</literal> à
+ <literal>ng</literal> composantes. Chaque composante représente une surface.
+ Elle doit être définie comme suit :
+ </para>
+ <itemizedlist>
+ <listitem>
+ <para>Une fonction Scilab.</para>
+ <para>Sa séquence d'appel doit être
+ <literal>surf(t, y)</literal>
+ </para>
+ </listitem>
+ <listitem>
+ <para>Une liste.</para>
+ <para>Elle doit se présenter comme suit :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+list(surf, x1, x2, ...)
+ ]]></programlisting>
+ <para>où la séquence d'appel de la fonction
+ <literal>surf</literal> est maintenant
+ </para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+r = surf(t, y, x1, x2, ...)
+ ]]></programlisting>
+ </listitem>
+ <listitem>
+ <para>Une chaîne de caractères.</para>
+ <para>Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab.
+ </para>
+ <para>En C, la séquence d'appel doit être :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+void surf(int *ny, double *t, double y[], int *ng, double gout[])
+ ]]></programlisting>
+ <para>En Fortran, elle doit être :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+subroutine surf(ny, t, y, ng, gout)
+double precision t, y(*), gout(*)
+integer ny, ng
+ ]]></programlisting>
+ </listitem>
+ </itemizedlist>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info</term>
+ <listitem>
+ <para>
+ liste contenant <literal>14</literal> éléments. La valeur par défaut est
+ <literal>list([], 0, [], [], [], 0, [], 0, [], 0, [], [], [], 1)</literal>.
+ </para>
+ <variablelist>
+ <varlistentry>
+ <term>info(1)</term>
+ <listitem>
+ <para>réel scalaire donnant le temps maximal pour lequel
+ <literal>g</literal> peut être évalué ou une matrice vide
+ <literal>[]</literal> si aucune limite de temps n'est imposée.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(2)</term>
+ <listitem>
+ <para>
+ drapeau indiquant si <literal>daskr</literal> retourne
+ ses valeurs intermédiaires calculées (<literal>= 1</literal>)
+ ou seulement les temps indiqués par l'utilisateur
+ (<literal>= 0</literal>).
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(3)</term>
+ <listitem>
+ <para>
+ vecteur de deux éléments donnant la définition
+ <literal>[ml,mu]</literal> de la matrice bande calculeé par
+ <literal>jac</literal>; <literal>r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)"</literal>.
+ Si <literal>jac</literal> retourne une matrice pleine, fixer
+ <literal>info(3)=[]</literal>.
+ Inutile si
+ <literal>info(8)=1</literal>.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(4)</term>
+ <listitem>
+ <para>réel scalaire donnant la taille maximale du pas. Fixer
+ <literal>info(4)=[]</literal> si illimité.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(5)</term>
+ <listitem>
+ <para>réel scalaire donnant le pas initial. Fixer
+ <literal>info(5)=[]</literal> si non spécifié.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(6)</term>
+ <listitem>
+ <para>
+ fixer <literal>info(6)=1</literal> si la solution est
+ non-négative, sinon fixer
+ <literal>info(6)=0</literal>.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(7)</term>
+ <listitem>
+ <para>
+ si ydot0 est fixé tel que
+ <literal>g(t0, y0, ydot0) = 0</literal>, alors fixer
+ <literal>info(7)=[]</literal>. Sinon, fixer
+ <literal>info(7)=[+-1, ..., +-1]</literal>, avec
+ <literal>info(7)(i) = 1</literal> si y(i) est une variable différentielle et
+ <literal>info(7)(i) = -1</literal> si y(i) est une variable algébrique
+ (si ses dérivées n'apparaissent pas explicitement dans la fonction g(t, y, ydot)).
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(8)</term>
+ <listitem>
+ <para>
+ méthode directe / Krylov. Fixer <literal>info(8)=1</literal> et founrnir une routine <literal>psol</literal>
+ si vous souhaitez que le solveur utilise des itérations de Krylov, sinon (méthode directe) fixer
+ <literal>info(8)=0</literal>.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(9)</term>
+ <listitem>
+ <para>
+ paramètres de Krylov. Inutile si vous avez fixé
+ <literal>info(8)=0</literal>. Sinon, fixer
+ <literal>info(9)=[]</literal> ou
+ <literal>info(9)=[maxl kmp nrmax epli]</literal>, où :
+ </para>
+ <para>
+ - maxl = nombre maximal d'itérations de l'algorithme GMRes (par défaut
+ <literal>min(5, neq)</literal>),
+ </para>
+ <para>
+ - kmp = nombre de vecteurs sur lesquels l'orthogonalisation est faite dans GMRes
+ (par défaut maxl),
+ </para>
+ <para>
+ - nrmax = nombre maximal de redémarrages de GMRes par intération non-linéaire
+ (par défaut <literal>5</literal>),
+ </para>
+ <para>
+ - epli = constante du test de convergence de GMRes (par défaut <literal>0.05</literal>).
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(10)</term>
+ <listitem>
+ <para>
+ conditions initiales. A ignorer si
+ <literal>info(7)=[]</literal>. Fixer
+ <literal>info(10)=1</literal> si le solveur doit s'arrêter après
+ le calcul des valeurs initiales, sinon fixer
+ <literal>info(10)=0</literal>.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(11)</term>
+ <listitem>
+ <para>
+ routine pour le calcul et la factorisation LU du préconditionneur pour <literal>psol</literal>.
+ Inutile si <literal>info(8)=0</literal>. Fixer
+ <literal>info(11)=1</literal> et fournir une routine <literal>pjac</literal> si l'<link linkend="external">external</link>
+ <literal>psol</literal> doit utiliser une routine spécifique, sinon fixer
+ <literal>info(11)=0</literal>.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(12)</term>
+ <listitem>
+ <para>
+ si vous souhaitez contrôler l'erreur localement sur toutes les variables, fixez
+ <literal>info(12)=[]</literal>. Sinon, fixez
+ <literal>info(12)=[+-1, ..., +-1]</literal>, avec
+ <literal>info(12)(i) = 1</literal> si y(i) est une variable différentielle et
+ <literal>info(12)(i) = -1</literal> si y(i) est une variable algébrique
+ (si ses dérivées n'apparaissent pas explicitement dans la fonction g(t, y, ydot)).
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(13)</term>
+ <listitem>
+ <para>
+ paramètres heuristiques. Ignorer si
+ <literal>info(7)=[]</literal>. Sinon, fixer
+ <literal>info(13)=[]</literal> ou
+ <literal>info(13)=[mxnit mxnj mxnh lsoff stptol epinit]</literal>, où :
+ </para>
+ <para>
+ - mxnit = nombre maximal d'itérations de Newton par évaluation du préconditionneur (par défaut
+ <literal>5</literal> si <literal>info(8)=0</literal>, <literal>15</literal> sinon),
+ </para>
+ <para>
+ - mxnj = nombre maximal d'évaluations du préconditionneur (par défaut
+ <literal>6</literal> si <literal>info(8)=0</literal>, <literal>2</literal> sinon),
+ </para>
+ <para>
+ - mxnh = nombre maximal de valeurs artificielles du pas <literal>h</literal> à tenter si info(7) ≠ [] (par défaut
+ <literal>5</literal>),
+ </para>
+ <para>
+ - lsoff = drapeau pour désactiver l'algorithme de recherche linéaire (lsoff = 0 pour activer, lsoff = 1 pour désactiver)
+ (par défaut <literal>0</literal>),
+ </para>
+ <para>
+ - stptol = pas minimal (dimmensionné) dans l'algorithme de recherche linéaire (par défaut <literal>(unit roundoff)^(2/3)</literal>),
+ </para>
+ <para>
+ - epinit = facteur déterminant dans le test de convergence de l'itération Newton (par défaut <literal>0.01</literal>).
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>info(14)</term>
+ <listitem>
+ <para>
+ verbosité. Fixer <literal>info(14)=1</literal> pour une information minimale,
+ <literal>info(14)=2</literal> pour une information maximale, sinon fixer
+ <literal>info(14)=0</literal>.
+ </para>
+ </listitem>
+ </varlistentry>
+ </variablelist>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>psol</term>
+ <listitem>
+ <para>
+ <link linkend="external">external</link> (fonction, liste ou chaîne de caractères).
+ Résout un système linéraire <literal>P*x = b</literal>,
+ où P est le préconditionneur LU-factorisé que <literal>pjac</literal>
+ a calculé auparavant et stocké dans <literal>wp</literal> et <literal>iwp</literal>.
+ </para>
+ <itemizedlist>
+ <listitem>
+ <para>Une fonction Scilab.</para>
+ <para>Sa séquence d'appel doit être
+ <literal>[r, ier] = psol(wp, iwp, b)</literal> et doit retourner la solution du système dans
+ <literal>r</literal> et un drapeau d'erreur <literal>ier</literal>.
+ </para>
+ </listitem>
+ <listitem>
+ <para>Une liste.</para>
+ <para>Elle doit se présenter comme suit :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+list(psol, x1, x2, ...)
+ ]]></programlisting>
+ <para>
+ où la séquence d'appel de <literal>psol</literal> est
+ </para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+psol(wp, iwp, b, x1, x2, ...)
+ ]]></programlisting>
+ <para>
+ <literal>psol</literal> retourne toujours la solution dans <literal>r</literal>.
+ </para>
+ </listitem>
+ <listitem>
+ <para>Une chaîne de caractères.</para>
+ <para>Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab
+ </para>
+ <para>En C, la séquence d'appel doit être :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+void psol (int*neq, double*t, double*y, double*ydot, double*savr,
+ double*wk, double*cj, double*wght, double*wp, int*iwp, double*b, double*eplin, int*ier, double*rpar, int*ipar)
+ ]]></programlisting>
+ où les vecteurs <literal>wp</literal> et <literal>iwp</literal> contiennent le préconditionneur LU-factorisé
+ <literal>P</literal>, <literal>wp</literal> représentant les valeurs et
+ <literal>iwp</literal> les pivots utilisés dans la factorisation.
+ <para>En Fortran, elle doit être :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+subroutine psol (neq, t, y, ydot, savr, wk, cj, wght,
+ wp, iwp, b, eplin, ier, rpar, ipar)
+double precision t,y(*), ydot(*), savr(*), wk(*), cj, wght(*), wp(*),
+ b(*), eplin, rpar(*)
+integer neq, iwp(*), ier, ipar(*)
+ ]]></programlisting>
+ </listitem>
+ </itemizedlist>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>pjac</term>
+ <listitem>
+ <para>
+ <link linkend="external">external</link> (fonction, liste ou chaîne de caractères). Calcule la valeur de
+ <literal>dg/dy + cj*dg/dydot</literal> pour une valeur donnée du paramètre
+ <literal>cj</literal> et la LU-factorise en deux vecteurs, réel et entier.
+ </para>
+ <itemizedlist>
+ <listitem>
+ <para>Une fonction Scilab.</para>
+ <para>Sa séquence d'appel doit être
+ <literal>[wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)</literal> et en retour,
+ les vecteurs <literal>wp</literal> et <literal>iwp</literal>
+ doivent contenir toutes les informations du préconditionneur factorisé.
+ </para>
+ </listitem>
+ <listitem>
+ <para>Une liste.</para>
+ <para>Elle doit se présenter comme suit :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+list(pjac, x1, x2, ...)
+ ]]></programlisting>
+ <para>
+ où la séquence d'appel de <literal>pjac</literal> est
+ </para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+pjac(neq, t, y, ydot, h, cj, rewt, savr, x1, x2,...)
+ ]]></programlisting>
+ <para>
+ <literal>pjac</literal> retourne toujours
+ <literal>dg/dy + cj*dg/dydot</literal> comme fonction de
+ <literal>(neq, t, y, ydot, h, cj, rewt, savr, x1, x2, ...)</literal>.
+ </para>
+ </listitem>
+ <listitem>
+ <para>Une chaîne de caractères.</para>
+ <para>Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab
+ </para>
+ <para>En C, la séquence d'appel doit être :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+void pjac (double*res, int*ires, int*neq, double*t, double*y, double*ydot, double*rewt, double*savr,
+double*wk, double*h, double*cj, double*wp, int*iwp, int*ier, double*rpar, int*ipar)
+ ]]></programlisting>
+ <para>En Fortran, elle doit être :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+subroutine pjac (res, ires, neq, t, y, ydot, rewt, savr,
+ wk, h, cj, wp, iwp, ier, rpar, ipar)
+double precision res(*), t, y(*), ydot(*), rewt(*), savr(*),
+ wk(*), h, cj, wp(*), rpar(*)
+integer ires, neq, iwp(*), ier, ipar(*)
+ ]]></programlisting>
+ </listitem>
+ </itemizedlist>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>hd</term>
+ <listitem>
+ <para>
+ vecteur réel servant à stocker le contexte de
+ <literal>daskr</literal> et reprendre l'intégration.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>r</term>
+ <listitem>
+ <para>
+ matrice réelle. Chaque colonne est le vecteur <literal>[t; x(t); xdot(t)]</literal> où
+ <literal>t</literal> est l'indice en temps aulequel la solution a été calculée,
+ <literal>x(t)</literal> est la valeur de la solution calculée,
+ <literal>xdot(t)</literal> est la dérivée de la solution calculée.
+ </para>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>nn</term>
+ <listitem>
+ <para>
+ vecteur à deux entrées <literal>[times num]</literal>,
+ <literal>times</literal> est la valeur du temps auquel la surface est traversée,
+ <literal>num</literal> est le nombre de surfaces traversées.
+ </para>
+ </listitem>
+ </varlistentry>
+ </variablelist>
+ </refsection>
+ <refsection>
+ <title>Description</title>
+ <para>Solution de l'équation différentielle implicite :</para>
+ <programlisting role="no-scilab-exec"><![CDATA[
+g(t, y, ydot) = 0
+y(t0) = y0 et ydot(t0) = ydot0
+ ]]></programlisting>
+ <para>Retourne les temps de traversée de surface et le nombre de surfaces traversées
+ dans <literal>nn</literal>.
+ </para>
+ <para>Des exemples détaillés se trouvent dans SCI/modules/differential_equations/tests/unit_tests/daskr.tst</para>
+ </refsection>
+ <refsection>
+ <title>Exemples</title>
+ <programlisting role="example"><![CDATA[
+// 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;
+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] = daskr([y0, y0d], t0, t, atol, rtol, res1, ng, gr1);
+nn
+
+// Retourne nn = [2.4698972 2]
+ ]]></programlisting>
+ </refsection>
+ <refsection role="see also">
+ <title>Voir aussi</title>
+ <simplelist type="inline">
+ <member>
+ <link linkend="ode">ode</link>
+ </member>
+ <member>
+ <link linkend="dasrt">dasrt</link>
+ </member>
+ <member>
+ <link linkend="dassl">dassl</link>
+ </member>
+ <member>
+ <link linkend="impl">impl</link>
+ </member>
+ <member>
+ <link linkend="fort">fort</link>
+ </member>
+ <member>
+ <link linkend="link">link</link>
+ </member>
+ <member>
+ <link linkend="external">external</link>
+ </member>
+ </simplelist>
+ </refsection>
+ <refsection>
+ <title>Historique</title>
+ <revhistory>
+ <revision>
+ <revnumber>5.5.0</revnumber>
+ <revdescription>daskr ajouté</revdescription>
+ </revision>
+ </revhistory>
+ </refsection>
+</refentry>
#include "machine.h"
#include "core_math.h"
-extern void C2F(dgesv)(int*, int*, double*, int*, int*, double*, int*, int*);
-extern double C2F(dlamch)(const char*);
-typedef void (*resfunc)(double*, double*, double*, double*, int*, double*, int*);
+extern void C2F(dgefa) (double *A, int *lead_dim_A, int *n, int *ipivots, int *info);
+extern void C2F(dgesl) (double *A, int *lead_dim_A, int *n, int *ipivots, double *B, int *job);
+extern double C2F(dlamch) (const char*);
+typedef void (*resfunc) (double*, double*, double*, double*, int*, double*, int*);
-void C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *actual, double *actualP,
- double *rewt, double *savr, double *wk, double *h, double *cj, double *wp, int *iwp,
- int *ier, double *rpar, int *ipar)
+void C2F(pjac1) (resfunc res, int *ires, int *neq, double *tOld, double *actual, double *actualP,
+ double *rewt, double *savr, double *wk, double *h, double *cj, double *wp, int *iwp,
+ int *ier, double *rpar, int *ipar);
+void C2F(psol1) (int *neq, double *tOld, double *actual, double *actualP,
+ double *savr, double *wk, double *cj, double *wght, double *wp,
+ int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2);
+
+void C2F(pjac1) (resfunc res, int *ires, int *neq, double *tOld, double *actual, double *actualP,
+ double *rewt, double *savr, double *wk, double *h, double *cj, double *wp, int *iwp,
+ int *ier, double *rpar, int *ipar)
{
int i = 0;
int j = 0;
+ int info = 0;
int nrow = 0;
double tx = 0;
double del = 0;
double delinv = 0;
double ysave = 0;
double ypsave = 0;
- double* e = NULL;
+ double * e = NULL;
- int neq = *nequations;
double SQuround = sqrt(C2F(dlamch)("P"));
tx = *tOld;
- e = (double*) calloc(neq, sizeof(double));
+ e = (double *) calloc(*neq, sizeof(double));
- for (i = 0 ; i < neq ; ++i)
+ for (i = 0 ; i < *neq ; ++i)
{
- del = Max(SQuround * Max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
+ del = Max(SQuround * Max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
del *= (*h * actualP[i] >= 0) ? 1 : -1;
- ysave = actual[i];
- ypsave = actualP[i];
+ del = (actual[i] + del) - actual[i];
+ ysave = actual[i];
+ ypsave = actualP[i];
actual[i] += del;
actualP[i] += *cj * del;
res(&tx, actual, actualP, e, ires, rpar, ipar);
}
delinv = 1. / del;
- for (j = 0 ; j < neq ; j++)
+ for (j = 0 ; j < *neq ; j++)
{
wp[nrow + j] = (e[j] - savr[j]) * delinv;
/* NaN test */
if (ISNAN(wp[nrow + j]))
{
*ier = -1;
- free (e);
+ free(e);
return;
}
- iwp[nrow + j] = i + 1;
- iwp[nrow + j + neq * neq] = j + 1;
}
- nrow += neq;
+ nrow += *neq;
actual[i] = ysave;
actualP[i] = ypsave;
}
+ /* Proceed to LU factorization of P. */
+ C2F(dgefa) (wp, neq, neq, iwp, &info);
+ if (info != 0)
+ {
+ *ier = -1;
+ }
+
free(e);
}
-void C2F(psol1)(int *nequations, double *tOld, double *actual, double *actualP,
- double *savr, double *wk, double *cj, double *wght, double *wp,
- int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2)
+void C2F(psol1) (int *neq, double *tOld, double *actual, double *actualP,
+ double *savr, double *wk, double *cj, double *wght, double *wp,
+ int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2)
{
- int i = 0;
- int nColB = 1;
- int info = 0;
- int *ipiv = NULL;
+ int i = 0, job = 0;
- ipiv = (int *) malloc(*nequations * sizeof(int));
+ C2F(dgesl) (wp, neq, neq, iwp, b, &job);
- C2F(dgesv) (nequations, &nColB, wp, nequations, ipiv, b, nequations, &info);
-
- if (info == 0)
+ /* NaN test */
+ for (i = 0; i < *neq; ++i)
{
- /* NaN test */
- for (i = 0; i < *nequations; ++i)
+ if (ISNAN(b[i]))
{
- if (ISNAN(b[i]))
- {
- *ier = -1;
- free (ipiv);
- return;
- }
+ /* Indicate a recoverable error, meaning that the step will be retried with the same step size,
+ but with a call to 'jacpsol' to update necessary data, unless the Jacobian data is current,
+ if (b[i] - b[i] != 0) in which case the step will be retried with a smaller step size. */
+ *ier = -1;
+ return;
}
}
- else
- {
- *ier = -1;
- }
-
- free (ipiv);
}
ypsave = ydot(i);
y(i) = y(i) + del;
ydot(i) = ydot(i) + cj*del;
- [e ires]=res1(tx, y, ydot);
+ [e ires] = res1(tx, y, ydot);
if ires < 0 then
ires = -1;
return;
ypsave = ydot(i);
y(i) = y(i) + del;
ydot(i) = ydot(i) + cj*del;
- [e ires]=res1(tx, y, ydot);
+ [e ires] = res1(tx, y, ydot);
if ires < 0 then
ires = -1;
return;
projects / scilab.git / commitdiff