Finding the zeros of a function

175 次查看（过去 30 天）
Tristan
Tristan 2013-10-8
评论： Lavorizia Vaughn 2021-9-30
I need to find where y=0 within 0<x<100
y=5*sin(1.9*x)+2.1*sin(9.1*x)

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2014-12-6
编辑：John D'Errico 2014-12-6
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 个评论
Gaetan Foisy
Gaetan Foisy 2019-4-7
Thank you so much, this saved me on a project I have.
Lavorizia Vaughn
Lavorizia Vaughn 2021-9-30
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

请先登录，再进行评论。

更多回答（6 个）

Shashank Prasanna
Shashank Prasanna 2013-10-8
  1 个评论
Tristan
Tristan 2013-10-8
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

请先登录，再进行评论。

Alan Weiss
Alan Weiss 2013-10-8
Or maybe some combination of ezplot and fzero.
Alan Weiss
MATLAB mathematical toolbox documentation
Walter Roberson
Walter Roberson 2013-10-8
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 个评论
Alan Weiss
Alan Weiss 2013-10-9
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 2015-6-25
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 ;-)

请先登录，再进行评论。

Dennis Yeh
Dennis Yeh 2014-12-6
编辑：Walter Roberson 2015-5-8
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 个评论
John D'Errico
John D'Errico 2014-12-6
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.

请先登录，再进行评论。

Mauricio Redaelli
Mauricio Redaelli 2015-6-25
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 个评论
显示 1更早的评论
Mauricio Redaelli
Mauricio Redaelli 2015-6-29
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 2015-6-29
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?

请先登录，再进行评论。

Mauricio Redaelli
Mauricio Redaelli 2015-6-29
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 个评论
Walter Roberson
Walter Roberson 2015-6-29
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 2015-6-30
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!

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by