MATLAB Answers

Lukas
0

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

Asked by Lukas
on 17 May 2018
Latest activity Commented on by 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?

  0 Comments

Sign in to comment.

Products


Release

R2017b

3 Answers

Answer by Torsten
on 17 May 2018
 Accepted Answer

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 0.

Best wishes

Torsten.

  7 Comments

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.

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.

Thats a nice way to do it, thanks.

Sign in to comment.


Answer by 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

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.


Answer by 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);

  0 Comments

Sign in to comment.