[Scilab-users] Numderivative and integrate

classic Classic list List threaded Threaded
8 messages Options
Carrico, Paul Carrico, Paul
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] Numderivative and integrate

Dear All,

 

I’m wondering if there’s a way to speed up the following code using both “numderivative” and “integrate” ?

 

I tried to replace the loops by vectors, but it fails …

 

Of course this basic example has been put in place in order to test the code – I’m working with test results with more than 200 000 measurements.

 

Thanks for any help and/or suggestions.

 

Regards

 

Paul

 

 

###############################################################

mode(0)

 

clear

 

PATH = get_absolute_file_path("test_signal.sce");

 

// here we generate a sin signal just for testing

Gamma_ampl = 10; // in [g]

Frequency = 100; // in [Hz]

Omega = 2 * %pi * Frequency;

Number_of_periods = 10;

 

t = [0 : (1/(20 * Frequency)) : (Number_of_periods / Frequency)];

Acceleration = Gamma_ampl *sin(Omega * t);

Signal = [t ; Acceleration];

 

scf()

plot(t,Acceleration,'g');

xtitle('Ref curve - accelerogram','t','$\gamma \mbox{[m/s$^{2}$]}$');

xs2png(gcf(),'Accelerogram.png');

 

// signal first derivation : y' = (Omega * Gamma_ampl)*cos(Omega * t)

Acceleration_prime = (Omega * Gamma_ampl) * cos(Omega * t);

 

// calculation of the jerk -> numerical derivative

function y=Gamma_function(x)

    y = interpln(Signal,x)

endfunction

 

n = size(t,'*');

Numl_Gamma_derrivative = zeros(2,n);

Numl_Gamma_derrivative(1,:) = t;

 

tic();

for i = 1 : n

    Numl_Gamma_derrivative(2,i) = numderivative(Gamma_function, Numl_Gamma_derrivative(1,i),1e-3,order=4);

end

Numederivative_duration = toc();

 

scf()

plot(Numl_Gamma_derrivative(1,:),Numl_Gamma_derrivative(2,:),'b')

plot(t,Acceleration_prime,'ro')

xtitle('Jerk - numerical derivation vs derivative formula','t','$\mbox{J [m/s}^{3}\mbox{]}$');

legend('Numerical','Reference',opt=1);

xs2png(gcf(),'Jerk.png');

 

// calculation of the speed

Speed_ref = (Gamma_ampl / Omega) * (1 - cos(Omega * t));

 

Speed_num_integration = zeros(2,n);

Speed_num_integration(1,:) = t;

tic();

Speed_num_integration(2,:) = integrate('Gamma_function','x',0,t);

First_integrate_duration = toc();

 

scf()

plot(Speed_num_integration(1,:),Speed_num_integration(2,:),'m')

plot(t,Speed_ref,'cd')

xtitle('Speed - numerical integration vs antiderivative formula','t','v [m/s]');

legend('Numerical','Reference',opt=1);

xs2png(gcf(),'Speed.png');

 

 

// calculation of the displacement

Displ_ref = (Gamma_ampl / Omega) * (t - (sin(Omega * t)/Omega));

 

function y=Speed_function(x)

    y = interpln(Speed_num_integration,x)

endfunction

 

Displ_num_integration = zeros(2,n);

Displ_num_integration(1,:) = t;

tic();

Displ_num_integration(2,:) = integrate('Speed_function','x',0,t);

Second_integrate_duration = toc();

 

scf()

plot(Displ_num_integration(1,:),Displ_num_integration(2,:),'r')

plot(t,Displ_ref,'b>');

xtitle('Displacement - numerical integration vs antiderivative formula','t','$\delta \mbox{ [m]}$');

legend('Numerical','Reference',opt=4);

xs2png(gcf(),'Displacements.png');

 

// information's

printf("The duration for Numederivative = %g\n",Numederivative_duration);

printf("The duration for 1rst integration = %g\n",First_integrate_duration);

printf("The duration for 2nd integration = %g\n",Second_integrate_duration);

EXPORT CONTROL :
Cet email ne contient pas de données techniques
This email does not contain technical data

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
ol.bond ol.bond
Reply | Threaded
Open this post in threaded view
|

Re: Numderivative and integrate

This post was updated on .
Hi,

As for differentiation, since you use the experimental data spaced at
regular intervals, you can use the concept of differentiation matrix:

function D =diff_matrix(t)

    k =size(t,'*');

    dt = (t(2)-t(1));

    m = matrix(ones(5,1)*[3:k-2],1,5*(k-4))';

    n = matrix([1;2;3;4;5]*ones(1,(k-4))+ones(5,1)*[0:k-5],1,5*(k-4))';

    d = ones((k-4),1)*[1 -8 0 8 -1];

    sp = sparse([m,n],d');

    sp([1 2],1:5) = [-25 48 -36 16 -3; -3 -10 18 -6 1];

    sp([k-1 k],(k-4):k) = [-1 6 -18 10 3; 3 -16 36 -48 25];

    D = 1/(12*dt)*sp;

endfunction

It is based on the classical five-point centred, fourth-order finite
difference approximations.

and then, the differential can be computed very fast as follows:


n = size(t,'*');

Numl_Gamma_derrivative = zeros(n,2);

Numl_Gamma_derrivative(:,1) = t';

 D =diff_matrix(t);

tic();

Numl_Gamma_derrivative(:,2) =D*Signal(2,:)';

Numederivative_duration = toc();

scf();

plot(Numl_Gamma_derrivative(:,1),Numl_Gamma_derrivative(:,2),'b')

plot(t,Acceleration_prime,'ro')

xtitle('Jerk - numerical derivation vs derivative formula','t','$\mbox{J
[m/s}^{3}\mbox{]}$');

legend('Numerical','Reference',opt=1);

For the integration, you should probably also rewrite the routine to
process in matrix form.

Yours.



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users
Hermes Hermes
Reply | Threaded
Open this post in threaded view
|

Re: Numderivative and integrate

In reply to this post by Carrico, Paul
He,
here I put a solution for differentiation
  in your case we obtain The duration for Numederivative = 1.18248
and for this other: The duration for Numederivative = 0.00125101
the FourMatDiff function could be improved.
Is there any function in Scilab to create band Matrices?

//http://mailinglists.scilab.org/Scilab-users-Numderivative-and-integrate-td4037333.html
mode(0)
 
clear
 
//PATH = get_absolute_file_path("test_signal.sce");
 
// here we generate a sin signal just for testing
Gamma_ampl = 10; // in [g]
Frequency = 100; // in [Hz]
Omega = 2 * %pi * Frequency;
Number_of_periods = 10;
 
t = [0 : (1/(20 * Frequency)) : (Number_of_periods / Frequency)];
Acceleration = Gamma_ampl *sin(Omega * t);
Signal = [t ; Acceleration];
 
scf()
plot(t,Acceleration,'g');
xtitle('Ref curve - accelerogram','t','$\gamma \mbox{[m/s$^{2}$]}$');
xs2png(gcf(),'Accelerogram.png');
 
// signal first derivation : y' = (Omega * Gamma_ampl)*cos(Omega * t)
Acceleration_prime = (Omega * Gamma_ampl) * cos(Omega * t);
 
// calculation of the jerk -> numerical derivative
function y=Gamma_function(x)
    y = interpln(Signal,x)
endfunction
 
n = size(t,'*');
Numl_Gamma_derrivative = zeros(2,n);
Numl_Gamma_derrivative(1,:) = t;

Y=Gamma_function(Numl_Gamma_derrivative(1,:));
//
function M2=FourMatDiff2(n)
 for j=1:n
    for i=1:n
    //M(i,j) = 0; // main diagonal
    if pmodulo(i-j,n)==1 then
        M2(i,j)=-2/3
          elseif pmodulo(i-j,n)==2 then
        M2(i,j)=1/12
          elseif pmodulo(i-j,n)==n-1 then
        M2(i,j)=2/3
          elseif pmodulo(i-j,n)==n-2 then
        M2(i,j)=-1/12
          elseif i-j==0 then
        M2(i,j)=0
        end
    end
end
//disp( full(M2) );  
endfunction;
N=n;
block=8;
sel=4;
h= t(3)-t(2);
X=0:h:N;
MD=FourMatDiff2(block);
tic();
//for i = 1 : n
//   Numl_Gamma_derrivative(2,i) = numderivative(Gamma_function,
Numl_Gamma_derrivative(1,i),1e-3,order=4);
//end
for i=1:N-block+1
 yd2($+1)=[MD/h*Y(i:i+block-1)'](sel);
end


Numederivative_duration = toc();
 
scf()
plot(Numl_Gamma_derrivative(1,4:197)',yd2,'b');
plot(t,Acceleration_prime,'ro');
xtitle('Jerk - numerical derivation vs derivative formula','t','$\mbox{J
[m/s}^{3}\mbox{]}$');
legend('Numerical','Reference',opt=1);
xs2png(gcf(),'Jerk.png');

// information's
printf("The duration for Numederivative = %g\n",Numederivative_duration);

Gracias



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Pinçon Bruno Pinçon Bruno
Reply | Threaded
Open this post in threaded view
|

Re: Numderivative and integrate

Le 15/12/2017 à 11:37, Hermes a écrit :
> He,
> here I put a solution for differentiation
>    in your case we obtain The duration for Numederivative = 1.18248
> and for this other: The duration for Numederivative = 0.00125101
> the FourMatDiff function could be improved.
> Is there any function in Scilab to create band Matrices?

     Hi,

   you could use sparse matrix. If (as it seems the case) your matrix
   can be defined by its diagonals maybe the following function could
   be useful for you. It lets you to build fastly a sparse matrix by
given the
   different diagonals and if the diagonal is constant only the constant
   could be provided, for instance if you want to build the classical
   matrix with 2 on the diag and -1 on the first (upper and lower) diagonals

   n = 6;  // a really small example
   A = my_spdiags([n,n],0,2,-1,-1,1,-1)


   If one diagonal of the matrix is not a constant you should provide
   the entire vector of the diagonal values (instead of the constant).

   hth
  Bruno

  PS : I coded this function a long time ago it's possible that there is
  such an internal function in scilab now and I have not tested it on
  scilab-6.x

===================================================

   function A = my_spdiags(dims,varargin)
    // a function to build easily a sparse matrix by filling some given
    // diagonals.
    //
    // 1/ dims is a vector with matrix dimensions [number of rows,
number of columns]
    // 2/ the special varargin keyword stands for the other input
    //    arguments. Any number of (valid) diagonals is allowed and
    //    each diagonal is introduced with (first) its
    //    diagonal number then (second) by the vector of values of the
diagonal.
    //    When the diagonal is to be filled by a same value it is allowed to
    //    replace the vector of values by the scalar.
    //    The diagonal number is defined (as usual) by k = j-i
    //
    // So the function usage is:
    // A = my_spdiags([m,n], num_diag1, val_diag1, num_diag2,val_diag2,...)
    //
    // AUTHOR : Bruno Pinçon [hidden email]

    if ~(typeof(dims) == "constant" & length(dims)==2 &
and(floor(dims)==dims)...
     & min(dims) >= 1 ) then
       error("bad first argument")
    end
    m = dims(1); n = dims(2);

    M = length(varargin)
    if modulo(M,2) ~= 0 then, error("bad entries..."), end
    nd = M/2;

    ij = []; val = [];

    for i = 1:2:M-1
       k = varargin(i);
       values =  varargin(i+1);
       if k < -(m-1) | k > (n-1) then
      error(msprintf("argument %d is not a good diag number",i+1))
       end
       if k >= 0 then
      nbelem = min(m,n-k);
      i = (1:nbelem)'; j = i+k;
       else
      nbelem = min(m+k,n);
      j = (1:nbelem)'; i = j-k;
       end
       nv = length(values)
       if nv==1 then       // a scalar was given, build the
corresponding vector
      values = values*ones(nbelem,1);
       elseif nv == nbelem then
      values = values(:)
       else                // the vector has not the same length than
the diagonal
      error(msprintf("argument %d has not the good size",i+2));
       end
       ij = [ij;[i,j]]; val = [val;values];
    end
    A = sparse(ij,val,dims)
endfunction


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Pinçon Bruno Pinçon Bruno
Reply | Threaded
Open this post in threaded view
|

Re: Numderivative and integrate

In reply to this post by Hermes

      HI Paul,

   I have not looked deeply in your code but why not
  use a spline representation of the Signal and use the
  spline representation (see the splin and interp functions)
  to both differentiate and integrate. Differentiation is
  provided by the interp function so the remaining work
  should be to integrate a spline but it is not a big deal.

   Just my 2 cents...
  Bruno

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Numderivative and integrate

In reply to this post by Hermes

Hi Hermes,

 

Regarding your question "the FourMatDiff function could be improved. Is there any function in Scilab to create band Matrices?",

please note that following code:

 

M = toeplitz([0  -2/3  1/12  zeros(1,n-5)  -1/12  2/3]);

M = tril(M) - triu(M);

 

produces the same results for n >=5 as your function:  FourMatDiff2(n)

 

Regards,

Rafael

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Numderivative and integrate

Hi Hermes,

 

Regarding your question "the FourMatDiff function could be improved. Is there any function in Scilab to create band Matrices?", the simpler code:

 

r = [0  -2/3  1/12  zeros(1,n-5)  -1/12  2/3];

M = toeplitz(r,-r);

 

produces the same results for n >=5 as your function FourMatDiff2(n)

 

Regards,

Rafael

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Samuel GOUGEON Samuel GOUGEON
Reply | Threaded
Open this post in threaded view
|

Re: Numderivative and integrate

In reply to this post by Carrico, Paul
Hello Paul,

The loop

Le 11/12/2017 à 11:53, Carrico, Paul a écrit :
...
for
i = 1 : n

    Numl_Gamma_derrivative(2,i) = numderivative(Gamma_function, Numl_Gamma_derrivative(1,i),1e-3,order=4);

end


could be avoided by vectorizing numderivative().
There was a previous recent thread about this, started by Hermes.

If you deem that vectorizing numderivative() would be preferable,
could you please report this wish on Bugzilla?
Then, it could be implemented for Scilab 6.1, provided that there will be
a 6.0.2 release for obsoleting the usage of the x input as a row:
Please see http://bugzilla.scilab.org/show_bug.cgi?id=14440#c7

Regards
Samuel




_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users