[Scilab-users] Command integrate() is too slow under certain expressions

[Federico Miyara]

Stéphane,

Thanks for your solution.

I've found another solution that is even slightly faster: I've modified 
the function intsplin(), changing
both instances of sum by cumsum (and some formal details), so the same 
interpolator is used for the whole set of data.

Regards,

Federico Miyara



On 11/02/2020 04:23, Stéphane Mottelet wrote:
>
> Sorry the ode call should be (y0 first, then x0)
>
> y  =  ode(-1,x0,x1,list(f,x1,y1,d1))
>
> The elapsed time seems OK:
>
> --> tic;y = ode(-1,x0,x1,list(f,x1,y1,d1));toc
>  ans  =
>
>    0.002235
>
> S.
>
> Le 11/02/2020 à 08:17, Stéphane Mottelet a écrit :
>>
>> Hello,
>>
>> The usual way to compute a primitive would be to use ode, like this:
>>
>> function  out=f(x, y, x1, y1, d1)
>>      out  =  interp(x,x1,y1,d1)
>> endfunction
>>
>> x0  =  0
>> x1  =  0:0.01:2*%pi;
>> y1  =  sin(x1);
>>
>> d1  =  splin(x1,y1);
>>
>> y  =  ode(x0,-1,x1,list(f,x1,y1,d1))
>> Your proposition is very slow because you are recomputing the spline 
>> many times.
>>
>> S.
>>
>> Le 10/02/2020 à 20:52, Federico Miyara a écrit :
>>>
>>> Dear all,
>>>
>>> The function integrate() is very handy to obtain a numerical 
>>> primtive of a function, particularly useful when there is no closed 
>>> form for the primitive or it is too involved. According to the 
>>> documentation
>>>
>>>     x = integrate(expr, v, x0, x1 [, atol [, rtol]])
>>>
>>>     expr: a character string defining a Scilab expression.
>>>     v:    a character string, the integration variable name
>>>
>>>
>>> However there is a case that would be interesting to handle, and it 
>>> is when one has a set of experimental (or simulated) data xk, yk. 
>>> Here there is no expression defining a function.
>>>
>>> One way to circumvent this is using as the expression some sort of 
>>> interpolator such as
>>>
>>> x = integrate("interp1(xk, yk, x, ''spline'')", "x", x0, x1)
>>>
>>> This works, but for some reason I don't quite understand it is very 
>>> slow.
>>>
>>> For instance,
>>>
>>> x0 = 0
>>> x1 = 0:0.01:2*%pi;
>>> y1 = sin(x1);
>>>
>>> tic
>>> X = integrate("interp1(x1, y1, x, ''spline'')", "x", x0, x1);
>>> toc
>>>
>>> Requires 27 s to execute. In the meantime, control is seemingly 
>>> returned to the console, one can enter instructions, but then the 
>>> program freezes until the integrate comand finishes.
>>>
>>> Changing "spline" to "linear" even worsens it rising to 33 s.
>>>
>>> Has anybody an idea of what can be happening?
>>>
>>> Maybe it computes the full interpolator for each single point? Even 
>>> if so, I have only 629 values of the independent variable.
>>>
>>> Regards,
>>>
>>> Federico Miyara
>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> users mailing list
>>> users at lists.scilab.org
>>> https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users
>> -- 
>> 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
>> http://www.utc.fr/~mottelet
>>
>> _______________________________________________
>> users mailing list
>> users at lists.scilab.org
>> https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users
> -- 
> 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
> http://www.utc.fr/~mottelet
>
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20200211/4552b1a9/attachment.htm>
-------------- next part --------------
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA
// Copyright (C) 2017 - Samuel GOUGEON
//
// Copyright (C) 2012 - 2016 - Scilab Enterprises
//
// This file is hereby licensed under the terms of the GNU GPL v2.0,
// pursuant to article 5.3.4 of the CeCILL v.2.1.
// This file was originally licensed under the terms of the CeCILL v2.1,
// and continues to be available under such terms.
// For more information, see the COPYING file which you should have received
// along with this program.

function v = intsplin_cum(x, s)
    //splin  numerical integration.
    //v = intsplin(x,s) computes the integral of y with respect to x using
    //splin interpolation and integration.
    //x and y must be vectors of the same dimension
    //
    //v = intsplin(s) computes the integral of y assuming unit
    //spacing between the data points.

    x = x(:);
    if argn(2) < 2 then
        s = x;
        x = (1:size(s,"*"))';
    else
        s = s(:);
        if size(x,"*")<>size(s,"*") then
            msg = _("%s: Wrong size for input arguments: Same size expected.\n")
            error(msprintf(msg, "intsplin"));
        end
    end
    if ~isreal(x) then
        if isreal(x,%eps)
            x = real(x);
        else
            msg = _("%s: Argument #%d: Real number(s) expected.\n");
            error(msprintf(msg, "intsplin",1))
        end
    end
    h = x(2:$) - x(1:$-1);
    y = real(s);
    d = splin(x, y);
    v = cumsum((h.*(d(1:$-1)-d(2:$))/12 + (y(1:$-1)+y(2:$))/2).*h);
    if ~and(imag(s)==0) then
        y = imag(s);
        d = splin(x, y);
        v = v + %i*cumsum((h.*(d(1:$-1)-d(2:$))/12 + (y(1:$-1)+y(2:$))/2).*h);
    end
    v = [0; v];
endfunction