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

Question

Lukas 2018-5-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/401190-numerically-stable-implementation-of-sin-y-atan-x-x

评论： Lukas 2018-5-18

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2018-5-17

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/401190-numerically-stable-implementation-of-sin-y-atan-x-x#answer_320605

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

James Tursa 2018-5-17

编辑：James Tursa 2018-5-17

在 MATLAB Online 中打开

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.

Torsten 2018-5-18

在 MATLAB Online 中打开

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 2018-5-18

Thats a nice way to do it, thanks.

请先登录，再进行评论。

Answer 2

Majid Farzaneh 2018-5-17

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/401190-numerically-stable-implementation-of-sin-y-atan-x-x#answer_320592

在 MATLAB Online 中打开

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

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Lukas 2018-5-17

在 MATLAB Online 中打开

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.

请先登录，再进行评论。