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

### Lukas (view profile)

on 17 May 2018
Latest activity Commented on by Lukas

on 18 May 2018

### Torsten (view profile)

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:
sin(y*atan(x))/x
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)
Both
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
sin(y*atan(x))
Does anyone have an idea how to solve my problem?

R2017b

### Torsten (view profile)

on 17 May 2018

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 href = ""</a> 0.
Best wishes
Torsten.

Lukas

### Lukas (view profile)

on 18 May 2018
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.
Torsten

### Torsten (view profile)

on 18 May 2018
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
Torsten.
Lukas

### Lukas (view profile)

on 18 May 2018
Thats a nice way to do it, thanks.

### Majid Farzaneh (view profile)

on 17 May 2018

Hi, You can easily add an epsilon to x like this:
sin(y*atan(x+eps))/(x+eps)

Lukas

### Lukas (view profile)

on 17 May 2018
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.