<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Hi all,</p>
<p>As discussed in this thread last month<br>
</p>
<p><a class="moz-txt-link-freetext"
href="https://www.mail-archive.com/users@lists.scilab.org/msg10679.html">https://www.mail-archive.com/users@lists.scilab.org/msg10679.html</a></p>
<p>I am glad to annouce that a first version of the sci-sundials
toolbox (maybe part of Scilab in the future) is available on Atoms
in the "Differential Equations" category (refresh the package list
if you don't see it). To have an idea of its features I pasted
below the content of the README.md file on the gitlab project of
sci-sundials (<a class="moz-txt-link-freetext" href="https://gitlab.com/mottelet/sci-sundials/">https://gitlab.com/mottelet/sci-sundials/</a>).<br>
</p>
<p>If you appreciate the work and want to help, doc, demos, or even
more if you know how to code in C and C++ you are welcome !</p>
<p>S.<br>
</p>
<h2 dir="auto" data-sourcepos="5:1-5:44">What has been done and what
is to be done</h2>
<p dir="auto" data-sourcepos="7:1-7:477">Until now only CVODE has
been interfaced but many common features have been developped in
the OdeManager class, hence interfacing IDA will be quite easy.
ARKODE (various modern explicit, implicit and mixed
explicit/implicit Runge-Kutta solvers) does not exist in the old
2.4.0 version of Sundials which is used in Scilab hence won't be
interfaced unless an upgrade is done (Sundials is now at version
5.7.0). The support of Sparse Jacobians is also missing for the
same reason.</p>
<h2 dir="auto" data-sourcepos="9:1-9:11">
Features</h2>
<p dir="auto" data-sourcepos="11:1-11:137">The CVODE gateway
cvode_solve() implements the following features which were missing
by the legacy LSODE/LSODA/LSODAR/... ode() gateway :</p>
<ul dir="auto" class="task-list" data-sourcepos="13:1-13:140">
<li class="task-list-item" data-sourcepos="13:1-13:140">
full a posteriori access to solver continuous extension of
solution at arbitrary points via a MList output when only one
lhs is given:</li>
</ul>
<pre class="code highlight js-syntax-highlight language-plaintext white" data-sourcepos="14:1-18:3" lang="plaintext"><code><span class="line" id="LC1" lang="plaintext">sol = cvode_solve(f, [t0 tf], y0)</span>
<span class="line" id="LC2" lang="plaintext">t = linspace(t0, tf, 1000)</span>
<span class="line" id="LC3" lang="plaintext">plot(t, sol(t)) </span></code></pre>
<p dir="auto" data-sourcepos="19:1-19:217">The <code>sol</code>
MList fields gathers all information related to the obtained
solution (time steps, solution at time steps, events, ...). The
solver can be restarted by giving the MLlist as first argument to
<code>cvode_extend</code> :</p>
<pre class="code highlight js-syntax-highlight language-plaintext white" data-sourcepos="20:1-22:3" lang="plaintext"><code><span class="line" id="LC1" lang="plaintext">sol2 = cvode_extend(sol, tx, ...)</span></code></pre>
<p dir="auto" data-sourcepos="23:1-23:179">where <code>tx</code> is
the time point to which solution has to be extended. The options
of the call that yielded <code>sol</code> are used and can be
changed as optional named parameters after <code>tx</code>.</p>
<ul dir="auto" class="task-list" data-sourcepos="24:1-24:102">
<li class="task-list-item" data-sourcepos="24:1-24:102">
a user-friendly access to solver options via optional named
parameters in cvode_solve call, i.e.</li>
</ul>
<pre class="code highlight js-syntax-highlight language-plaintext white" data-sourcepos="25:1-27:3" lang="plaintext"><code><span class="line" id="LC1" lang="plaintext">cvode_solve(f, tspan, y0, h0=0.01, rtol=1e-3)</span></code></pre>
<ul dir="auto" class="task-list" data-sourcepos="28:1-28:158">
<li class="task-list-item" data-sourcepos="28:1-28:158">
a really simpler way to give time span of integration allowing
to choose between error driven solver internal time steps and
user fixed time steps, i.e.</li>
</ul>
<pre class="code highlight js-syntax-highlight language-plaintext white" data-sourcepos="29:1-33:3" lang="plaintext"><code><span class="line" id="LC1" lang="plaintext">[t,y] = cvode_solve(f, [t0 tf], y0)</span>
<span class="line" id="LC2" lang="plaintext">[t,y] = cvode_solve(f, [t0 t1 ... tf], y0)</span>
<span class="line" id="LC3" lang="plaintext">[t,t] = cvode_solve(f, tspan, y0, t0=0)</span></code></pre>
<p dir="auto" data-sourcepos="34:1-34:216">the latter style being
the closest to actual <code>ode()</code> behavior where solution
is by default given at user time steps and <code>t0</code> is not
necessarily equal to <code>tspan(1)</code>, which is the default
in the two former calls above.</p>
<ul dir="auto" class="task-list" data-sourcepos="35:1-35:303">
<li class="task-list-item" data-sourcepos="35:1-35:303">
a better and user-friendly specification of events via a
variable number of outputs event function (giving value of event
equations, wheter to stop integration for a given event and
event direction selection), minimal style beeing a single
output. Information about events is also simpler to get:</li>
</ul>
<pre class="code highlight js-syntax-highlight language-plaintext white" data-sourcepos="36:1-44:3" lang="plaintext"><code><span class="line" id="LC1" lang="plaintext">function [eq,term,dir] = evfun(t,y)</span>
<span class="line" id="LC2" lang="plaintext"> eq = y(1)-1.7;</span>
<span class="line" id="LC3" lang="plaintext"> term = %f;</span>
<span class="line" id="LC4" lang="plaintext"> dir = 1;</span>
<span class="line" id="LC5" lang="plaintext">end</span>
<span class="line" id="LC6" lang="plaintext">[t,y,te,ye,ie] = cvode_solve(f, tspan, y0, events = evfun)</span>
<span class="line" id="LC7" lang="plaintext">sol = cvode_solve(f, tspan, y0, events = evfun)</span></code></pre>
<p dir="auto" data-sourcepos="45:1-45:83">in the latter call
information about events is recovered in <code>sol.te,sol.ye,sol.ie</code>.</p>
<ul dir="auto" class="task-list" data-sourcepos="46:1-46:90">
<li class="task-list-item" data-sourcepos="46:1-46:90">
support of complex solution with detection of complexity in <code>y0</code>
or <code>f(t0,y0)</code>, e.g.</li>
</ul>
<pre class="code highlight js-syntax-highlight language-plaintext white" data-sourcepos="47:1-53:3" lang="plaintext"><code><span class="line" id="LC1" lang="plaintext">function out = crhs(t,y)</span>
<span class="line" id="LC2" lang="plaintext"> out = 10*exp(2*%i*%pi*t)*y;</span>
<span class="line" id="LC3" lang="plaintext">end</span>
<span class="line" id="LC4" lang="plaintext">[t,y] = cvode_solve(crhs, [0,5], 1)</span>
<span class="line" id="LC5" lang="plaintext">plot(t,real(y),t,imag(y))</span></code></pre>
<ul dir="auto" class="task-list" data-sourcepos="54:1-54:189">
<li class="task-list-item" data-sourcepos="54:1-54:189">
Support of a callback function called after each successfull
step, giving access to current solver statistics and allowing to
stop integration (e.g. by a "stop" button on a GUI), e.g.</li>
</ul>
<pre class="code highlight js-syntax-highlight language-plaintext white" data-sourcepos="55:1-63:3" lang="plaintext"><code><span class="line" id="LC1" lang="plaintext">function stop = scicallback(t,y,flag,stats)</span>
<span class="line" id="LC2" lang="plaintext"> stop = %f</span>
<span class="line" id="LC3" lang="plaintext"> if flag == "step"</span>
<span class="line" id="LC4" lang="plaintext"> mprintf("%s : hlast=%g\n", flag, stats.hlast)</span>
<span class="line" id="LC5" lang="plaintext"> end</span>
<span class="line" id="LC6" lang="plaintext">end</span>
<span class="line" id="LC7" lang="plaintext">[t,y] = cvode_solve(f, tspan, y0, intcb=scicallback);</span></code></pre>
<ul dir="auto" class="task-list" data-sourcepos="64:1-65:0">
<li class="task-list-item" data-sourcepos="64:1-65:0">
support of an arbitrary number of dimensions of ode state</li>
</ul>
<h2 dir="auto" data-sourcepos="66:1-66:14">
Performance</h2>
<p dir="auto" data-sourcepos="68:1-68:326"><code>cvode_solve()</code>
is roughly two times faster than <code>ode()</code> for both
fixed methods (Adams and BDF). As <code>ode()</code> already did,
compiled and dynamically linked C,C++ or Fortran externals are
supported by <code>cvode_solve()</code>. When using such
externals instead of Scilab functions <code>cvode_solve()</code>
is generally an order of magnitude faster.</p>
<p> </p>
<pre class="moz-signature" cols="72">--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
<a class="moz-txt-link-freetext" href="http://www.utc.fr/~mottelet">http://www.utc.fr/~mottelet</a>
</pre>
</body>
</html>