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

3 views (last 30 days)
Lukas
Lukas on 17 May 2018
Commented: Lukas on 18 May 2018
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?

Accepted Answer

Torsten
Torsten 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.
  7 Comments

Sign in to comment.

More Answers (2)

Majid Farzaneh
Majid Farzaneh on 17 May 2018
Hi, You can easily add an epsilon to x like this:
sin(y*atan(x+eps))/(x+eps)
  1 Comment
Lukas
Lukas 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.

Sign in to comment.


Ameer Hamza
Ameer Hamza on 17 May 2018
How about defining it like this
f = @(x,y) y.*(x==0) + sin(y.*atan(x))./(x+(x==0)*1);

Products


Release

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by