Finding the zeros of a function

263 views (last 30 days)
I need to find where y=0 within 0<x<100
y=5*sin(1.9*x)+2.1*sin(9.1*x)

Accepted Answer

John D'Errico
John D'Errico on 6 Dec 2014
Edited: John D'Errico on 6 Dec 2014
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].)
  2 Comments
Gaetan Foisy
Gaetan Foisy on 7 Apr 2019
Thank you so much, this saved me on a project I have.
Lavorizia Vaughn
Lavorizia Vaughn on 30 Sep 2021
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

Sign in to comment.

More Answers (6)

Shashank Prasanna
Shashank Prasanna on 8 Oct 2013
  1 Comment
Tristan
Tristan on 8 Oct 2013
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

Sign in to comment.


Alan Weiss
Alan Weiss on 8 Oct 2013
Or maybe some combination of ezplot and fzero.
Alan Weiss
MATLAB mathematical toolbox documentation

Walter Roberson
Walter Roberson on 8 Oct 2013
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.
  2 Comments
Alan Weiss
Alan Weiss on 9 Oct 2013
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
Walter Roberson on 25 Jun 2015
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 ;-)

Sign in to comment.


Dennis Yeh
Dennis Yeh on 6 Dec 2014
Edited: Walter Roberson on 8 May 2015
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
  1 Comment
John D'Errico
John D'Errico on 6 Dec 2014
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.

Sign in to comment.


Mauricio Redaelli
Mauricio Redaelli on 25 Jun 2015
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
  3 Comments
Mauricio Redaelli
Mauricio Redaelli on 29 Jun 2015
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
Walter Roberson on 29 Jun 2015
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?

Sign in to comment.


Mauricio Redaelli
Mauricio Redaelli on 29 Jun 2015
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!
  2 Comments
Walter Roberson
Walter Roberson on 29 Jun 2015
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
Mauricio Redaelli on 30 Jun 2015
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!

Sign in to comment.

Categories

Find more on Introduction to Installation and Licensing in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!