Kshitij Singh ( Staff ) answered . 01/01/2020

You mean a function that evaluates sinx/x sinxx for x≠0 and 1 1for x=0 x=0?

Obviously it would be easy to just program something with a statement like

function result = mySinc(x)

% My implementation of the sinc function.

 

if x == 0

    result = 1;

else

    result = sin(x)/x;

end

end

But that could in fact be less than ideal. Sure, it avoids the singularity and has the correct mathematical structure, but we of course have to remember that the values in x aren't numbers: they're patterns of bits represented by switches in hardware and only approximate numbers. They work pretty well most of the time, but there are circumstances where things just don't work right. Such as when evaluating a + log(1+x) for small x, or when evaluating a-b when a and b have very similar values. You can lose a lot of precision if you evaluate the latter naïvely (although sometimes you may have no choice, alas). Division by small numbers can cause problems too, although luckily both sinx sinx and x x are comparable when |x|≪1 |x|≪1. Maybe it wouldn't be a problem, but when dealing with floating point arithmetic it's rarely good practice to have code like if x == foo.

What to do in this case? Consider the Taylor series of sinx/x sinxx centred at x=0 x=0. I won't write it down, but you surely know how to compute it.

Then replace the testif x == 0 by if abs(x) < sqrt(eps),where in Matlab eps returns the machine epsilon (comparable to 10^-16 10-16 on most systems these days), and replace the line result = 1; by result = ...;, where the ... is going to be your calculation of the Taylor series truncated after an appropriate order (usually about 2 should do, but you should check for yourself the scale of the highest order term you pick for x~10^-8 x~10^-88 and whether the result will be as accurate as you need or not).

The nice thing here is that what you lose in speed running this code you'll likely gain in accuracy, especially when evaluating the sine-integralSi(x)=∫x0f(t)dt Si(x)=∫x0f(t)dt, where f(t) f(t)is the continuous extension of the sinc sinc function to all reals as its domain.