Please help me with the optimization

Hello Dear Experts,

I wanted to ask how can I optimize the following solution:

 % Step 1: f(x) = (15/8)*(1-x^2)^2 ; 0 < x < 1.
  % Step 2: F(x) = 1.875*x - 1.25*x^3 + 0.375*x^5
  % Step 3: Let u ~ U[0,1] => X = F^(-1)(u) => u = F(x) 
  % => u = 1.875*x - 1.25*x^3 + 0.375*x^5
  % Step 4: Find the x that satisfy the equation.
  n = 100000;
  u = rand(n,1);
  y = zeros(n,1);
for i = 1:n
    y(i,1) = fzero(@(x) 0.375*x^5 - 1.25*x^3 + 1.875*x - u(i,1),0.5);
end
I = find(y > 0 & y < 1);
hist(y(I),1000);
It runs like 5-10 minutes.

Answer 1

[Roger Stafford]

A couple of possibilities for reducing the execution time:

1) Use the 'roots' function instead of 'fzero' and select the single root of the five whose imaginary part is zero for this polynomial.

2) Do a sort on u and use the y solution for one sorted u value as an initial estimate of y for the next higher u. Because successive u values will necessarily be close together, these estimates should already be fairly accurate. Also it might help to use the Newton-Raphson method instead of 'fzero' in combination with this - you already have an expression for the derivative of F(x), namely f(x). Since you are only doing a histogram of y it shouldn't matter that these are also now in sorted order, but if necessary their original order could easily be restored.

Note: The 'find' operation is unnecessary here since all y values will lie between 0 and 1, which is true because F(0)=0 and F(1)=1 and F is monotone increasing in between.