MATLAB Answers


Numerically stable implementation of sin(y*atan(x))/x

Asked by Lukas
on 17 May 2018 at 10:21
Latest activity Commented on by Lukas
on 18 May 2018 at 10:33

I am trying to implement a modified version of the Magic Tyre Formula. The simplified version of my problem ist that i need to calculate this function:


especially at and around x = 0. I know that this function is defined for x = 0 because:

sin(y*atan(x))/x = sin(y*atan(x))/(y*atan(x)) * (y*atan(x))/x = sin(c)/c * atan(x)/x * y with c = y*atan(x) 


sin(c)/c and atan(x)/x

are defined for x = 0 and c = 0.

I would like to use built in functions to solve my problem, because im not that good at numerics.

What I have tried until now is:

1) I can calculate sin(c)/c by using the built in sinc function, but then I still have to calculate atan(x)/x which i have found no solution for by now.

2) I know that

sin(atan(x)) = x/(sqrt(1+x^2))

But i havent found a way to rewrite this equation using


Does anyone have an idea how to solve my problem?


Sign in to comment.




3 Answers

Answer by Torsten
on 17 May 2018 at 12:09
 Accepted Answer

By L'Hospital, lim (x->0) sin(y*atan(x))/x = y.

Thus define your function to be y if x=0 and sin(y*atan(x))/x if x 0.

Best wishes



Thanks a lot for this detailed answer.

So the reason this all works is basically that in both cases, the derivatives of the nominator and the denominator are the same at zero and that dividing x/x near zero ist numerically stable. Thats definitely something I did not know before.

The case with a very small y won't happen in my calculation, because its a fixed tyre parameter, which is normally somewhere in between 1 and 2.

So if i understood everything correctly it should suffice to just exclude x==0 and calculate sin(y*atan(x))/x in all other cases.

Again thanks a lot for all the answers. For me, everything seems clear now.

I suggest

function z = your_function(x,y)
z = y.*ones(size(x));
i = find(x);
z(i) = sin(atan(x(i)).*z(i))./x(i);

Best wishes


Thats a nice way to do it, thanks.

Sign in to comment.

Answer by Majid Farzaneh on 17 May 2018 at 10:40

Hi, You can easily add an epsilon to x like this:


  1 Comment

Thank you for the idea, unfortunately thats exactly what i am trying to avoid, because I want to use this formula in a simulation and i cant guarantee that for example x does not equal -eps.

I even thought about using for example the power series of atan(x) and divide it by x:

atan(x) = sum((-1)^k * (x^(2k+1))/(2k+1),k = 0..inf)
atan(x)/x = sum((-1)^k * (x^(2k))/(2k+1),k = 0..inf)

But then i still dont know how many iterations i need, to use the full range of double precision and I am still not sure if this would be an efficient implementation.

Sign in to comment.

Answer by Ameer Hamza
on 17 May 2018 at 14:34

How about defining it like this

f = @(x,y) y.*(x==0) + sin(y.*atan(x))./(x+(x==0)*1);


Sign in to comment.

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today