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

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?

Answer 1

[Torsten]

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]

Thanks a lot for this this Idea. I have implemented a simple test case for your suggestion, and it seems to work correctly but i dont understand why. This is the function:

function [solution,if_used] = SinAtan(x,y)
    if(x == 0)
        solution = y;
        if_used = 1;
    else
        solution = sin(y*atan(x))/x;
        if_used = 0;
    end
end

And thats the test case:

clear all;
close;
clc;

x = linspace(-1e-30,1e-30,1e4+1);
y = 2*ones(size(x));
solution = NaN(size(x));
if_used = NaN(size(x));

for i=1:size(x,2)
    [solution(i),if_used(i)] = SinAtan(x(i),y(i));
end

figure
plot(x,solution)

figure
plot(x,if_used)

What i don't understand is that even if i get ridiculously close to zero the values seem to be okay. Honestly my expectation was that the values get worse if i get really close to 0 and that they are only okay if i am exactly at zero or far enuogh away.

Is it possible to understand why this seems to work, without knowing the numerics?

[Torsten]

Why do you think that the values get worse when you approach 0 ?

Your function is analytic if you set it to y at x=0:

https://www.wolframalpha.com/input/?i=Series%5BSin%5By*Atan%5Bx%5D%5D%2Fx,%7Bx,0,10%7D%5D

Thus you can evaluate it around 0 without any problems (like sin(x)/x).

Best wishes

Torsten.

[Lukas]

Maybe thats a misconception i always had about numerics. I thought dividing somithing colse to zero by something close to zero numerically is generally a bad idea.

Honestly i am not even really sure why I am allowed to do it with sin(x)/x ...

Anyways thanks a lot for your help, i think you answered all my questions. Maybe I need to take a closer look at how numeric division works :D

[James Tursa]

Consider the Taylor series for these functions:

atan(x) = x - x^3/3 + x^5/5 - ...

sin(x) = x - x^3/3! + x^5/5! - ...

The closer x gets to 0, the closer these series will be to simply x. So when x is small enough for that x^3/3 term to become smaller than eps(x) in double precision, the atan(x) calculation will evaluate to x exactly, and you will have the equivalent of

sin(y*x)/x

Then as long as y is not too large, at some point as x gets smaller the (y*x)^3/3! in the sin calculation will become smaller than eps(y*x) and you will get exactly y*x as the result. So, when x is very small you will numerically be doing the following calculation:

y*x/x

This is numerically stable even for very small values of x (scaling and un-scaling by the same amount) as long as y is also not small, in which case you could eventually have underflow/denormalized problems. E.g.,

>> y = 100
y =
   100
>> x = 1e-200
x =
    1.000000000000000e-200
>> y*x/x
ans =
   100         <-- recovered y nicely
>> y = 1e-200
y =
    1.000000000000000e-200
>> y*x/x
ans =
     0     <-- underflow caused y*x to be 0, so did not recover y

If you were concerned about this calculation for even these potentially tiny values of y and x, you would need to test for these types of conditions as part of your if-test and simply return y for these cases. E.g., you might test to see if x is small and x*y is smaller than the largest denormalized number, then just return y. Or something similar. Stuff like that would have to be in your if-test in addition to testing for x == 0 exactly.

[Lukas]

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]

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]

Thats a nice way to do it, thanks.

Answer 2

[Majid Farzaneh]

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

sin(y*atan(x+eps))/(x+eps)

[Lukas]

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.

Answer 3

[Ameer Hamza]

How about defining it like this

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