# [Scilab-users] Numderivative and integrate

8 messages
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
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
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.htmlmode(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
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
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
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