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

조회 수: 5 (최근 30일)
Lukas
Lukas 2018년 5월 17일
댓글: 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?

채택된 답변

Torsten
Torsten 2018년 5월 17일
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
Torsten
Torsten 2018년 5월 18일
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 2018년 5월 18일
Thats a nice way to do it, thanks.

댓글을 달려면 로그인하십시오.

추가 답변 (2개)

Majid Farzaneh
Majid Farzaneh 2018년 5월 17일
Hi, You can easily add an epsilon to x like this:
sin(y*atan(x+eps))/(x+eps)
  댓글 수: 1
Lukas
Lukas 2018년 5월 17일
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.

댓글을 달려면 로그인하십시오.


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

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품


릴리스

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by