Finding the zeros of a function

I need to find where y=0 within 0<x<100

y=5*sin(1.9*x)+2.1*sin(9.1*x)

Answer 1

[John D'Errico]

This is an old post, so there is no real need to add to it, but for someone who may want to solve a similar problem, I'll chime in.

I am always pleased with the capabilities of chebfun (as found on the file exchange.) It seems to have had no problem in finding all 85 roots of that function over the interval of interest.

fun = @(x) 5*sin(1.9*x)+2.1*sin(9.1*x);
roots(chebfun(fun,[0,100]))
ans =
            0
       1.7019
       3.4031
       5.1027
       6.3858
       6.5009
       6.7989
       8.0633
       8.3036
       8.4871
       9.7565
       11.455
       13.155
       14.857
       16.559
       18.261
       19.961
       21.659
       22.929
       23.112
       23.353
       24.617
       24.915
        25.03
       26.313
       28.013
       29.714
       31.416
       33.118
       34.819
       36.519
       37.802
       37.917
       38.215
       39.479
        39.72
       39.903
       41.172
       42.871
       44.571
       46.273
       47.975
       49.677
       51.377
       53.075
       54.345
       54.528
       54.769
       56.033
       56.331
       56.446
       57.729
       59.429
        61.13
       62.832
       64.534
       66.235
       67.935
       69.218
       69.333
       69.631
       70.895
       71.135
       71.319
       72.588
       74.287
       75.987
       77.689
       79.391
       81.092
       82.793
       84.491
       85.761
       85.944
       86.184
       87.449
       87.747
       87.862
       89.145
       90.845
       92.546
       94.248
        95.95
       97.651
        99.35

As an alternative, one could simply sample the function at a fine interval, looking for any sign changes. Then call fzero on each interval where a sign change was found.

xs = linspace(0,100,1001);
ys = fun(xs);
scinter = find(diff(sign(ys)));

See that there were 85 intervals found where a sign change occurred. I carefully chose code such that the first interval would be found, so fzero will find the zero at 0.

ninter = numel(scinter)
ninter =
     85

xroots = NaN(1,ninter);
for i = 1:ninter
  xroots(i) = fzero(fun,xs(scinter(i) + [0 1]));
end

Of course, the code above will miss roots that were too close together compared to the discretization interval of 0.1 that I chose. (That is, 1001 steps on the interval [0,100].)

[Gaetan Foisy]

Thank you so much, this saved me on a project I have.

[Lavorizia Vaughn]

coukd you please help me. i have the code:

function p = findmanyzeros(f,a,b,n,tol)

x=linspace(a,b,n+1);

y=f(x);

scinter=find(diff(sign(y)));

i=numel(scinter);

p=NaN(1,i);

for k=1:i

p(k)=findzero(f,x(scinter(k)),x(scinter(k)+1),tol);

end

the final kine isnt working. it is supposed to findzeros

Answer 2

[Shashank Prasanna]

Did you try this?

http://www.mathworks.com/help/matlab/ref/fminbnd.html

[Tristan]

I tried but it just gives me one answer

>> y = @(x)abs(5*sin(1.9*x)+2.1*sin(9.1*x));

>> x = fminbnd(f, 0, 100)

x =

   67.5442

Answer 3

[Alan Weiss]

Or maybe some combination of ezplot and fzero.

Alan Weiss

MATLAB mathematical toolbox documentation

Answer 4

[Walter Roberson]

There are 11 zeros in that range. You can eyeball the plot to approximate the ranges and then fzero() over the distinct ranges.

You are not going to be able to find closed form solutions for the zeros; it involves roots of two 90th degree polynomials.

[Alan Weiss]

Actually, I get 85 distinct roots (84 if we don't count 0). Some are quite close together.

fun = @(x)5*sin(1.9*x)+2.1*sin(9.1*x);
t = 0:.01:100;
for ii = 1:length(t)
yy(ii) = fzero(fun,t(ii));
end
zz = unique(yy);
zz = sort(zz);
z = round(1e9*zz)*1e-9; % Get rid of spurious differences
z = unique(z);
q = find(z>100);
z(q) = [];

Alan Weiss

MATLAB mathematical toolbox documentation

[Walter Roberson]

I was not able to figure out later why I thought there were only 11!

I can say, though, that this question led me to filing a bug report against a different manufacturer's software ;-)

Answer 5

[Dennis Yeh]

EDIT: I suggested this solution in the past, but found that it does not work well.

% Find all intersections of y1 and y2, where x is the independent variable:
x = linspace(0,10,10000); y1 = sin(x); y2 = cos(x);

coeff = polyfit(x,(y2-y1),100); % Perform a degree 100 polynomial curve fit on the difference between y1 and y2. 

possible_zeros = sort(unique(abs(roots(coeff)))); % Roots of the polynomial curve fit - the absolute value is to convert complex roots, and the unique() function turns each complex pair to only 1 point

error = polyval(coeff,possible_zeros); % Evaluate the polynomial curve fit at every root - answers close to zero correspond to an intersection

tolerance = 1e-6; % Change this value to be larger if intersections are skipped, and smaller if fake intersections are returned

actual_zeros_indices = find(abs(error) < tolerance); % Find the indices of the vector "possible_zeros" that correspond to actual zeros

actual_zeros = possible_zeros(actual_zeros_indices); % Plug the indices into possible_zeros to find the real zeros

f1 = figure(1);
subplot(2,1,1);
plot(x,y1,x,y2); % Plot the original functions

subplot(2,1,2);
plot(x,(y2-y1),'b',x,polyval(coeff,x),'g',actual_zeros,polyval(coeff,actual_zeros),'r*') % Plot the difference between y1 and y2, the curve fit, and the zeros

[John D'Errico]

Please don't tell people to do 100 degree polynomial fits with polyfit, or with ANY tool in double precision! The coefficients that arise will be essentially numerical garbage. I'm sorry, but this solution is terrible advice.

Answer 6

[Mauricio Redaelli]

Sorry, it´s an old post. Hope you don´t mind. I solved using this:

 % code
s = sign(y); 
s(isnan(s))=[];
rr = find(sign(y)==0); %look up for machine number zeros 
r = find(([s(1) s]==[s s(end)])==0); %look up for nearby machine number zeros
if isempty(rr)
    return
else
    r((rr+1)==r)=[]; %get rid of duplicities
    r((rr-1)==r)=[];
end

[Walter Roberson]

What is your "y" here? Considering you are using find() on something produced from y, it appears that y must be a numeric vector rather than a symbolic expression. But if it is a numeric vector then how do you know how fine-grained to make the vector? I see your plot talks about step sizes but the code you posted does not; how do you know that you have probed finely enough ?

You could keep going finer until you stopped generating more crossings, but then in the general case you run into the possibility that there are two zeros very close together, such that you might need to refine down several decimal places to tell them apart.

[Mauricio Redaelli]

You are right, y was genereted according to the expression provided. But this code compares the sign of y to the signs of the one-element-translated y. They will be all equal, except right on the zeros. It doens´t matter how fine-grained is the vector, because the sign will have the same step, and the transition from 1 to -1 occurs only in on point, right? I´m really curious to verify where this could be wrong, as I really am a novice.

This is my first post and my first answer, so pardon my english flaws. Thank you.

[Walter Roberson]

f = @(x) 5*sin(1.9*x)+2.1*sin(9.1*x);  %from above

y = f(linspace(0,100,2));

How many of the zeros do you detect in y?

y = f(linspace(0,100,10));

How many of the zeros do you detect in this y?

x = linspace(0,100,500);
y = f(x);

There, 500 points sampled for what we happen to know are 85 roots. If it does not matter how fine-grained the vector y is, then 500 points should be enough to find all of the roots. But your code returns 78 roots.

plot(x,y,'b-', x(r), y(r), 'r*');

Now try with x = linspace(0,100,790); and observe that your code counts 84 roots (your code is discarding the root at 0 which is the 85th root.) But then try with x = linspace(0,100,800); and observe that your code counts only 76 roots. How can increasing how close together the samples are result in fewer roots being found? How is someone to know that it is safe to use 790 samples but that 800 will not find all of the roots?

Answer 7

[Mauricio Redaelli]

Walter, this is very interesting. It seems to me that x=0 will never be found because the sign(y) will never see a transition, as we have none negative values. It´s a bug and should be fixed. Another bug is that the roots I found is not always the closest to the zero line - they always lay after the transition. I must easily fix it adding a new r = find(...) line, but moving the sign vector for the left, and then comparing those two r vectors. Another interesting point of discussion is the number of roots as a function of the x step. Well, for me this is an issue that should be on the concern of who build the y vector. If the function is not suficiently fine-grained to solve all the roots, so be it. Check the figure. Thank you again for the discussion!

[Walter Roberson]

But the original poster didn't even know how many zeros there were in the range (I wouldn't either, not just be looking at that function!). How would they know whether they had used enough points in the vector?

Remember, the original poster did not give a vector of samples: the original poster gave an expression, and you choose to interpret it in terms of sampling over a vector of locations, so you as the proposer have a responsibility to say how to choose that vector in order to be sure of getting all of the roots.

[Mauricio Redaelli]

I agree with you, but I can't provide a general rule of sampling to be sure we will get all of the roots. Maybe my algorithm is not suficient to find roots af any function in any condition, but sampling is an analytical task one must do in every case, not only to find roots. Anyway, thank you a lot for the discussion!