Re: Numerical Differentiation and Integration of Array Data

[Carlo de Falco]

On 4 Dec 2011, at 10:22, c. wrote:

> 
> On 3 Dec 2011, at 22:14, syberraith wrote:
> 
>> Actually, I'm unsure if I could apply those simplifacations to the actual
>> functions I am using.  It works well enough for a function simple as sin(a *
>> t), however my function include a dynamic correction factor as well.
>> 
>> Perhaps you like to look over the paper I'm writing. for publication rather
>> than homework, that has the full explanation of what my function are trying
>> to express mathematically?
>> 
>> I'd be glad to give you a credit.
> 
> Although it might be more complicated to implement in case of a 
> different function the algorithm should work identically.
> 
> Find below a generaliziotion of your procedure.
> 
> c.
> 

And here is the application to your problem as described in "gyro.c"
as you can see the numeriacl integration step appears to be useless
and (either with or without it) the computation is very fast:

>> tic, gyro, toc
err_hor =  2.9681e-15
err_tan = -6.2235e-19
err_ver = -1.0915e-14
The tangential sum is:  0.000000e+00
The horizontal sum is:  7.604169e-16
The vertical sum is:    0.000000e+00
Elapsed time is 0.008413 seconds.

c.


##----------------------------------------------------------------
## file gyro.m
##----------------------------------------------------------------
## algorithm parameters
n = 1000;


## function parameters
gyr_rpm = 1200;
p = 60 / gyr_rpm;
gyr_rad = 2.50e-2;
gyr_vel = 2 * pi * 1200 / 60;
xrm_rad = 1.50e-1;
xrm_vel = 2 * pi * 31.4 / 60;

## if you need to apply the same procedure
## to different functions you'll want to 
## make it into a function

function [f, s, s2] = process (fun, t)
  
  h = t(2) - t(1); 

  f  = diff (fun (t)) / h;
  s2 = h * sum (f);
  s  = fun (t (end)) - fun (t (1));

endfunction

## define here the function(s) you want to process
function  res = raw_vel (x, gyr_rad, gyr_vel, xrm_vel)
  res = sqrt ((-gyr_rad * gyr_vel * sin (gyr_vel * x) + gyr_rad * xrm_vel) .^ 2 
+ 
              (gyr_rad * gyr_vel * cos (gyr_vel * x)) .^ 2);
endfunction

function res = del_the (x, gyr_rad, xrm_rad, gyr_vel, xrm_vel)
  f_ang = atan2 (-gyr_rad * gyr_vel* sin (gyr_vel * x) + xrm_rad * xrm_vel, 
                 gyr_rad * gyr_vel * cos (gyr_vel * x)); 

  r_ang = atan2 (sin (gyr_vel * x), 
                 cos (gyr_vel * x));

  res = f_ang + r_ang;
  res(res > pi) -=  2 * pi;
endfunction

tan_vel = @(t) cos (del_the (t, gyr_rad, xrm_rad, gyr_vel, xrm_vel)) .* ...
    raw_vel(t, gyr_rad, gyr_vel, xrm_vel);

hor_vel = @(t) -sin (gyr_vel * t) .* ...
    cos (del_the (t, gyr_rad, xrm_rad, gyr_vel, xrm_vel)) .* ...
    raw_vel (t, gyr_rad, gyr_vel, xrm_vel);

ver_vel = @(t) cos (gyr_vel * t) .* ...
    cos (del_the (t, gyr_rad, xrm_rad, gyr_vel, xrm_vel)) .* ...
    raw_vel(t, gyr_rad, gyr_vel, xrm_vel);

## time grid
t = linspace (0, p , n);

## process the first function
[f, hor_sum, sf] = process (hor_vel, t);
# consisency error
err_hor = hor_sum - sf

## process the second function
[f1, tan_sum, sf1] = process (tan_vel, t);
# consisency error
err_tan = tan_sum - sf1

## process the third function
[f2, ver_sum, sf2] = process (ver_vel, t);
# consisency error
err_ver = ver_sum - sf2

printf ("The tangential sum is:  %1.6e\n", tan_sum );
printf ("The horizontal sum is:  %1.6e\n", hor_sum );
printf ("The vertical sum is:    %1.6e\n", ver_sum );