How do you get first n number of solutions while solving numerically ?
3 views (last 30 days)
Show older comments
Nuri Efe TATLI
on 11 Jun 2023
Answered: John D'Errico
on 11 Jun 2023
I have a function as follows:
f = beta^3*((sin(beta*L1)*coth(beta*L1))-cos(beta*L1))-(2*k*sin(beta*L1)) == 0;
L1 and k variables are constants with values.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1408154/image.png)
I want to solve this equation for Beta numerically. But I want first 5 solutions starting from zero and going up.
Basicly first five values of beta that makes this equation zero.
vpasolve only gives 1 solution and if you want more solutions you can only get solutions at random locations.
How may I achieve this ?
4 Comments
Dyuman Joshi
on 11 Jun 2023
Edited: Dyuman Joshi
on 11 Jun 2023
Please specify the values of L1 and k. The idea will be to use to them to plot the graph of f vs beta and see where it crosses the function crosses the x axis.
The only other typo is that, 1 is missing from the last L in the sin() term.
Accepted Answer
Dyuman Joshi
on 11 Jun 2023
k = 1.3333e+06;
L1 = 0.5;
f = @(beta) beta.^3.*((sin(beta*L1).*coth(beta*L1))-cos(beta*L1))-(2*k*sin(beta*L1));
%plotting the function
fplot(f,[0 50])
hold on
%x-axis
yline(0)
xticks(0:5:50)
%Initial value guess
val=[6 12 18 24 30];
for k = 1:numel(val)
%specify a value to fzero to find a solution close to that value
out(k) = fzero(f,val(k));
end
out
More Answers (1)
John D'Errico
on 11 Jun 2023
Too late, since you have an answer that solves your problem to your satisfaction. I want to point out something you probably need to know though.
Finding the first n solutions using a numerical solver is actually an impossible task to solve. This is because you cannot insure that you have found all solutions over some interval. You may lose solutions, and that means you got close to the first n solutions, but you missed a few in there. And sometimes, a numerical solver can find solutions that are not in fact solutions. The problem is fraught with issues that touch on the weak points of numercal solvers, and if you throw a difficult problem at them, well, OOPS!
Ugh. Why am I saying this? First, how does a tool like fzero work? I'll suggest fzero is your best choice for a 1-d numerical rootfinder, since it is more robust than a tool like fsolve. Fzero prefers a bracket, where the function is known to cross zero. So it wants an interval [a,b], such that f(a) and f(b) have different signs. (Or f(x) could be exactly zero at one of the endpoints. That is ok too.) If you don't give fzero a bracket, but only ONE point (a), then fzero will search for a second point where f(x) has a different sign than f(a). And then once it has a bracket, it is happy again.
But once fzero has a bracket, then it uses a hybrid set of methods that insure efficient convergence on good problems, but also insure the convergence cannot go too slowly on some classically known bad problems that can arise. We don't need to worry about those tweaks, but only accept that a tool like fsolve does not use them, and so can have issues of its own. (A good course on numerical analysis would be nice here, but I won't be spending the time to teach it. Sorry.)
Let me suggest a simple problem to solve, then I'll modify it in just a tiny way. Find the root of the function
f1 = @(x) 1 - 0.01*x;
Yes, I know the root lies at x==100. But the computer does not. Numerical tools do NOT actually look at the function itself, but only as if the function is a black box they cannot look inside. So they send in x values, and then look at what comes out of the magical black box.
opts = optimset('fzero');
opts.Display = 'iter';
x0 = 1;
[x,fval,exitflag] = fzero(f1,x0,opts)
So fzero does eventialluy look out as far as x==100. once it gets there, it zooms right in on the solution.
But if we give fzero a bracket, then it is supremely happy, even if the bracket is very wide.
x0 = [0,1000];
[x,fval,exitflag] = fzero(f1,x0,opts)
But now, suppose we give fzero a slightly different function?
f2 = @(x) 1 - 0.01*x - exp(-100*(x - 23).^2);
f1 and f2 are virtually identical everywhere, except for a tiny blip, that occurs only around x==23.
fplot(f1,[10,40])
hold on
fplot(f2,[10,40])
hold off
So f2 actually has a zero crossing, TWO of tem, in fact, very close to x==23. And everywhere ealse, f1 and f2 are pretty much identical. I could have made the blip an even sharper one too. Visualize an approximate Dirac delta function.
The problem is, no numerical tool can find that pair of roots, especially if I make the blip a very sharp one. If I give fzero the function f2, it will still find only the root at x==100.
x0 = 1;
[x,fval,exitflag] = fzero(f2,x0,opts)
In fact, it evaluated the function at x=21.48 and at x=29.9631. It completely missed the pair of roots near x=23.
And, fzero can also be tricked, if you give it a function with a zero crossing that is not actually a root.
x0 = [1,2];
f3 = @tan;
[x,fval,exitflag] = fzero(f3,x0,opts)
So fzero finds a root at pi/2, but then it returns an exitfkag of -5, to say that it thinks this is not actlually a root. An exitflag of -5 means it thinks this is a point of singularity, but not a root.
help fzero
These are general problems of such a tool. We can always trick a solver, as long as we understand how the tool works. So insuring you have found the first n roots, this is actually an impossible task, as I said before.
It gets worse of course. Consider the function f4.
f4 = @(x) sin(1./x);
fplot(f4,[0,3])
For the function f3(x), there are infinitely many roots in the interval [0,b], for ANY positive value of b, no matter how small you make b. So again, finding the first n positive roots of f3 is a mathematical impossibility.
There are things you can do to try to make the problem semi-tractable. I posted a function called allcrossings on the file exchange some time ago, but allcrossings is not perfect. allcrossings uses a scheme where it evaluates the function at many points on an interval, looking for any sign change. Then on any sub-interval, it calls fzero, because it has a bracket that should contain a root. But allcrossings does not find the first n roots. It tries to find ALL roots, within limits on its abilities.
The point is? Finding the first n roots is a difficult problem for a numerical solver, not one for which there will ever be an easy solution.
0 Comments
See Also
Categories
Find more on Loops and Conditional Statements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!