A function calculating the value of expression (sin(x) - x)*x^(-3).

write a function (exact specification below) calculating as accurately and as quickly as possible the value of the expression f(x) = (sin(x) - x)*x^(-3).

the maximum absolute error of the results calculated using the sent function should not exceed ≈100ε,

[James Tursa]

What have you done so far? What specific problems are you having with your code? Is this an exercise in controlling the cancellation error for small x? Or ...?

[Mikl Nik]

function Y = c4(X) 

 
n = size(X,1) * size(X,2); 
 
 
Y = zeros(size(X,1), size(X,2)); 
for i=1:n 
    Y(i) = -1/6; 
    k = 1; 
    while k<200 
   Y(i) = Y(i) + (-1)^(k+1)*X(i)^(2*k)/factorial(2*k+3);
   k = k + 1; 
    end 
end 

Well, I use Taylor series to solve this problem and the values of expressions are pretty close, but they are still exceeding 100*eps.

[Mikl Nik]

And yes, it's an example for controlling the cancellation error for small x.

And as I say, I used Taylor series:

f(x) = -1/6 + x^2/120 - x^4/5040 + x^6/362880 - x^8/39916800 + ...

So with help of while loop I was trying to write Taylor series for N = 200.

Answer 1

[David Hill]

x=.5;%whatever
n=200;
f=sum((ones(1,n+1)*x).^(0:2:2*n).*(-1).^(1:n+1)./factorial(2*(0:n)+3));
g=(sin(x)-x)*x^-3;
h=abs(f-g)<10*eps;

Answer 2

[James Tursa]

So, for small values of x your code looks correct to me (although the N=200 is WAY WAY OVERKILL ... you don't need nearly this many terms). And when you compare it to extended precision from Symbolic Toolbox using vpa( ) this seems to confirm that you nailed it:

>> digits 100  % <-- use 100 digits extended precision
>> c = @(x)(sin(x)-x)/x^3
c =
  function_handle with value:
    @(x)(sin(x)-x)/x^3
>> c(1e-5)
ans =
  -0.166667284899440  % <-- lots of cancellation error in this
>> c4(1e-5)
ans =
  -0.166666666665833
>> c4(1e-5) - double(c(vpa(1e-5)))
ans =
     0   % <-- compared to extended precision you nailed it down to the last bit
>> 
>> c(1e-10)
ans =
     0  % <-- MASSIVE cancellation error here ... lost everything
>> c4(1e-10)
ans =
  -0.166666666666667
>> c4(1e-10) - double(c(vpa(1e-10)))
ans =
     0   % <-- compared to extended precision you nailed it down to the last bit

So, maybe adjust your while loop to end sooner, but your Taylor Series code is correct. E.g., as soon as your terms stop contributing to the sum because they are too small, you can stop the while loop. Consider x = 1e-2. Each term is going to go down by (1e-2)^2/factorial(something) compared to the last term. It only takes a few of these terms before you are beyond 1e-16 or so and they stop contributing to the sum in double precision.

Was there a particular range of inputs that you are interested in? For a generic function you would probably set a threshold of input magnitude for when to calculate the function outright and when to use the Taylor Series.

When you thought you were greater than 100*eps, I think you were incorrectly comparing to the original function as written using a double input instead of comparing to an extended precision routine. As the results show above, not only is your Taylor Series code under 100*eps, it is matching the extended precision calculation down to the last bit. Can't ask for better than that!

[Mikl Nik]

Wow, I have not expected it, haha. :)

It was the range from -8*pi to 8*pi. And it’s like a task, so I have test from teacher for my function that probably show me the accuracy and it always shows 0%, and now I think it’s because of the fact that for x larger than 1 or even pi, the Taylor’s series is not the best solution.

[James Tursa]

Correct. You should have an input magnitude check. Something like this for each input value

if( abs(x) < something )
    % use Taylor Series code
else
    % use original function code
end

You may have to experiment a bit to find the sweet spot for the threshold. E.g., try this for several values of x to determine when each of c(x) or c4(x) diverges from the vpa answer. Then use that to pick a threshold:

[c(x)-double(c(vpa(x)));c4(x)-double(c(vpa(x)))]

I've done this to figure it out, but I will let you work with this and discover it on your own.

[James Tursa]

P.S., this line

Y = zeros(size(X,1), size(X,2)); 

would be better as this

Y = zeros(size(X));

And this line

n = size(X,1) * size(X,2); 

would be better as

n = numel(X); 

[Mikl Nik]

function Y = c4(X) 
 
n = numel(X); 
 
Y = zeros(size(X));
 
for i=1:n 
    if(abs(X(i)) < 0.301)  
        A = X(i)^2; 
        Y(i) = -1/6 + A/120 * (1 - A/42 * (1 - A/72 * (1 - ...
        A/110))); 
    else      
   Y(i) = (sin(X(i)) - X(i)) * X(i)^(-3);  
    end 
    
end 

That's my code. Test shows 99% accuracy and 77% efficiency. And I want to do it better. Maybe you have an idea?

[David Hill]

function Y = c4(X)
Y(X<.301)= -1/6 + X(X<.301)/120 .* (1 - (X(X<.301)/42).* (1 - (X(X<.301)/72).*(1 - X(X<.301)/110)));
Y(X>=.301)=(sin(X(X>=.301))-X(X>=.301)).*X(X>=.301).^-3;
end

[James Tursa]

What does 77% efficient mean? Where does this number come from?

Same question for the 99% accuracy. How did that get determined?

[Mikl Nik]

Test function shows there results. I've got this test function from my teacher and it divides the run time of my program by the run time of the most effective program.

[David Hill]

Did you try the above code?