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
[email protected]
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
[email protected]
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
[email protected]
http://lists.scilab.org/mailman/listinfo/users
// 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
_______________________________________________
users mailing list
[email protected]
http://lists.scilab.org/mailman/listinfo/users
