Multiple roots with fzero, or an alternative to find multiple root

Dear MATLAB users, I want to get several roots using "fzero" function. I also tried "chebfun" but I could not run the program. So I would highly appreciate if anyone can give me any advice:

% clear all
clc
N = 100 ; 
d = 6 ; 
kappa = 1 ; 
lambda = linspace(0.05,0.175, N ); 
nu_p = linspace(0.1,0.35, N ); 
[Lm,NP] = meshgrid(lambda, nu_p ); 
cosfi = 1;
f = @(x,lambda,nu_p) ((lambda - nu_p)+(lambda*x)-cosfi.*(((sqrt((1+x)./2)).*exp(-(d*(sqrt(2*lambda*(1+x))))))).*(kappa*(sqrt((1-x)./(1+x)))-(sqrt((1+x)./(1-x ))))); 
for k1 = 1:length(lambda ) 
for k2 = 1:length(nu_p ) 
z(k1,k2) = fzero(@(x) f(x,lambda(k1),nu_p(k2)), 0.000000000001 ); 
t1 = sqrt(1+(z(k1,k2)));
t2 = sqrt(1-(z(k1,k2)));
q = sqrt(.5).*t1.*exp(-d.*sqrt(2.*lambda(k1)).*t1);
J11 = 0;
J12(k1,k2) = -q.*((1 + kappa)/2).*sqrt(1-(z(k1,k2)).^2);
J21(k1,k2) = lambda(k1) + (1/2).*q.*[(1./(2.*(1 + (z(k1,k2)))))-d.*sqrt(lambda(k1)./(2.*(z(k1,k2))))].*[(1./((1+(z(k1,k2))).^2)).*(t1./t2).*kappa + (1./(1-(z(k1,k2)).^2)).*(t2./t1)];
J22 = 0;
% Eigenvalues
Eig1 = sqrt(J12.*J21);
Eig2 = -sqrt(J12.*J21);
EigReal = real(Eig1);
EigIm = imag(Eig1);
end 
end 
% DeltaE = Lambda - nu_p;
Z=abs(z);
figure(1 ) 
mesh(nu_p,lambda, EigReal)
view(2)
% title('kd=6');
 xlabel('nu_p');
 ylabel('lambda');
 zlabel('Z');

When I change "d" 9 or 12 there will be several roots which this code would just find the first one, and I am interested in finding Eigenvalue for those other roots also.

Thanks in advance

[Walter Roberson]

With d = 9 your equations cause fzero to fail. You should change your fzero to

fzero(@(x) f(x,lambda(k1),nu_p(k2)), [-1+eps 1-eps] )

[Walter Roberson]

You could try using https://www.mathworks.com/matlabcentral/fileexchange/55206-findroots

[Argumanh]

Thanks Walter,

I have two questions now:

1)Now with the alteration you suggest, for which root we are plotting the eigenvalues?

2)How may I use the link to plot eigenvalues for every root?

Thanks in advance

[Walter Roberson]

For (1): which-ever root fzero() happens to find.

In the runs I checked, your functions have singularities at +/- 1. fzero() and fsolve() do not have the ability to put in an upper and lower bound that is independent of the initial guess: they are restricted to putting in either an initial guess or an upper and lower bound. vpasolve() has the same restriction, of being able to specify a starting point or a bounds but not both.

I am still working on how best to integrate FindRoots. But in the meantime, I suggest that you rewrite your code under the assumption that an fzero-like call can return multiple values, something like

these_roots = fzero_like_call(@(x) f(x,lambda(k1),nu_p(k2)), [-1+eps 1-eps] );

Now rewrite your code under the assumption that these_roots is non-scalar. (Test with a couple of random values.)

[Walter Roberson]

When I investigate for d = 9, I see that it is not possible for there to be multiple roots in those ranges of lambda and nu_p .

[Argumanh]

Thanks Walter,

I could not understand how to use these_roots, should I just put that sentence in my code and try? Also about Findroots, Are you going to fix this using Findroots?

I really do appreciate your help

[Argumanh]

Dear Walter, I also tried to use chebfun, can you figure where the problem is in this code?

z = linspace(-1+1.e-5,1-1.e-3,101);
d = 6;
Lambda= linspace(0.01,0.11,101);
DeltaE = linspace(0.01,0.11,101);
kappa = 1;
[Lm,DE] = meshgrid(Lambda, DeltaE ); 
% equations
 t1 = sqrt(1+z);
 t2 = sqrt(1-z);
 q = sqrt(.5).*t1.*exp(-d.*sqrt(2.*Lambda).*t1);
cosfi = 1;
for k1 = 1:length(Lambda) 
for k2 = 1:length(DeltaE) 
myfun(k1,k2) = chebfun(@(z) (DeltaE(k2)+Lambda(k1)*z - .5*cosfi*sqrt(.5*(1-z))...
     .*(kappa-(1+z)./(1-z)).*exp(-d*sqrt(2*Lambda(k1)*(1+z))))/pi,'splitting','on');
% r(k,:) = roots(myfun);
r(k1,k2) = roots(myfun);
if length(r) > 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% jacobbian EigenValue:
t1 = sqrt(1+r);
t2 = sqrt(1-r);
q = sqrt(.5).*t1.*exp(-d.*sqrt(2.*Lambda).*t1);
J11 = 0;
J12 = -q.*((1 + kappa)/2).*sqrt(1-r.^2);
J21 = Lambda + (1/2).*q.*[(1./(2.*(1 + r)))-d.*sqrt(Lambda./(2.*r))].*[(1./((1+r).^2)).*(t1./t2).*kappa + (1./(1-r.^2)).*(t2./t1)];
J22 = 0;
% Eigenvalues
Eig1 = sqrt(J12.*J21);
Eig2 = -sqrt(J12.*J21);
EigReal = real(Eig1);
EigIm = imag(Eig1);
% plot
mesh(DeltaE,Lambda,EigReal);
end
end
end
%axis([-1 1 -0.02 0.02]);
axis auto;

[Walter Roberson]

If you are using https://www.mathworks.com/matlabcentral/fileexchange/47023-chebfun-current-version then the object that gets returned from chebfun claims to of size infinity by something. That cannot be stored in myfun(k1,k2) as that expects a scalar object, not something infinite in size. You need to store in myfun{k1,k2}

Then you try to take roots() of all of the myfun that have been created to this point, instead of taking just the root of the current item.

When you obtain the roots, you attempt to store them in r(k1,k2) -- but you are using chefun in hopes that the you will get multiple roots, so you cannot expect to store the result in the single scalar location r(k1,k2) . You need to store into a cell array again.

On the next line you test if length(r) > 0 but you are testing all of the roots accumulated so far, not just the roots that you just calculated.

The calculations that happen after that are probably going to fail if you produced multiple roots.

The calculations that happen after that involve the 1 x 101 vector Lambda, so EigReal is going to come out as a 1 x 101 vector. That is a problem for the mesh() call, as the mesh() call needs a 101 x 101 object there.

If you manage to create a 101x101 array to mesh() there, then it is not going to show up on screen. You are doing that mesh() plot every iteration, but MATLAB does not display anything until you use drawnow() or pause() or figure() or return to the keyboard.

I wonder if perhaps the calculations should be using Lambda indexed at k1 and building a single value to store in a 101 x 101 array that you mesh() once at the end. But if so then I still do not understand what output you want if there are multiple roots at a location.

[Argumanh]

Dear Walter, thank you very much for your comments. In the case of just one root, I want to calculate the eigenvalue for that specific root. And in the case of several roots(let's say three roots) for every root I want to plot the mesh plot, Every root is a fixed point in this case, I want to calculate the Eigenvalue for every fixed point found for a range of variables we defined as "lambda" and "DeltaE". Again, I really do appreciate your valuable explanations on my query.

[Walter Roberson]

To confirm: if the maximum number of roots found anywhere in the area is R, then you want to draw R separate mesh plots, with the values at the single-solution points being the same but the value for the particular location that has the multiple roots varying over the different plots?

If so then suppose at one point there are two different roots and at a different point there are three different roots, then you might have something like

A1 0
0 B1

and

A2 0
0  B2

as your first two plots, but should the third plot be

A1 0
0  B3

or should it be

A2 0
0  B3

?

[Argumanh]

Dear Walter, I think the first two answers are OK, but I do not get how we get the third plot like that,you say that we will not have A3 0 0 B3 right? So then lets try with A1 0 0 B2 and see what happens.

Thank you very much, Armin

[Walter Roberson]

There is at least one point where the function only has a single root. The points at which it does not have a single root: can you guarantee that the number of roots will be the same between all of those points? If not, then the upper left corner might have two roots, A1 and A2, but the lower right might have three roots, B1, B2, B3. If you want three plots because some point on the grid has three roots, then you need to decide what content you want for the places that have different numbers of roots: when you are presenting the third plot, the one that shows the B3 root, then what should be in the upper left corner, which A root? The minimum among all of the A roots? The maximum?

So far I see no evidence that any of your points has anything other than 1 root. I will attach the script after I do another test.

[Argumanh]

Dear Walter,

thank you very much for your attention. The ponit is that, there is a constant named "kappa" which I put it as "1", I am going to use it with the values of 0.2 0.4 0.6 0.8 and 1 . So there may be some points which we get multiple roots. Anyway, I would highly appreciate your help.

Answer 1

[Walter Roberson]

Script enclosed. You will want to bump up M and N to your 101 each. I did not test that high with this version of the script because it is not fast.

I see no evidence at all that any point has more than one root.

When I was testing the fzero() version, I was able to show that the function under consideration then did not have more than one root for fixed parameters. What I did not do is go back to check to see if the function you are passing to chebfun is the same as the one I investigated earlier.

[Argumanh]

Thank you Walter. You mean I should check whether fzero() and chebfun() give the same answer or not? but for other values of "d" and "kappa" there should be more than one root. For example d=12 must have multiple roots.

I really appreciate your kind attention.

Armin

[Walter Roberson]

What is a combination of parameters that you know should have multiple roots? Not just "d=12": specific lambda and deltaE and kappa as well ?

[Walter Roberson]

Your z is exactly the same size as your Lambda, and has to be the same size because you have

 t1 = sqrt(1+z);
 t2 = sqrt(1-z);
 q = sqrt(.5).*t1.*exp(-d.*sqrt(2.*Lambda).*t1);

so t1 and t2 are the same size as z. You are using .* between an expression in Lambda and t1, so they have to be compatible sizes. The exp() returns that size and you t1 .* that result, so t1 and the exp() have to be compatible sizes.

But I have to wonder whether it is really true that you want Lambda and z to be in one-to-one correspondence with each other, or if instead you want to be calculating over all combinations of z and Lambda ?

[Argumanh]

Dear Walter,

To answer your first question, I have to say if you use cosfi=-1, d=9, kappa=1 you will get three roots.

To answer your second question, I got a little confused, I think I should calculate for all combinations of z and Lambda. Does the code do this?

[Walter Roberson]

You only ever used cosfi = 1 in the above :(

The current code matches z and lambda, one z for each lambda.

[Argumanh]

It is OK Dear Walter, let me try to see the results. Then I will open another question to calculate for all combinations of z and Lambda. Thank you very much for your time.

[Walter Roberson]

In the eigroot2 code that I posted, if you change cosfi to -1 then it will show three plots.

Plotting with all the combinations of z and lambda will take it up to like 6 dimensions, with all of these being independent:

1. Lambda
2. DeltaE
3. z
4. kappa
5. d

Even if you make kappa and d temporary constants, you have three independents, plus the "z" variable of your anonymous function that is confusingly named the same as your z of (-1,+1) interval, with you expecting that multiple values will give zeros. Anyhow, you can't really mesh plot anymore, too many dimensions.

[Argumanh]

Dear Walter,

I am also interested to see which roots give me the eigenValues, is there any way to get every roots? I think I can only see in "r" the third root.

Best, Armin

[Walter Roberson]

Please confirm which variables you want to hold constant for the purposes of determining whether there are multiple roots? I would like to be sure we are still talking about the code in https://www.mathworks.com/matlabcentral/answers/399986-multiple-roots-with-fzero-or-an-alternative-to-find-multiple-root#comment_566505 and that what you are looking for is all of the z values within line

myfun(k1,k2) = chebfun(@(z) (DeltaE(k2)+Lambda(k1)*z - .5*cosfi*sqrt(.5*(1-z))...
   .*(kappa-(1+z)./(1-z)).*exp(-d*sqrt(2*Lambda(k1)*(1+z))))/pi,'splitting','on');

So that would be the set of values, z, such that

(DeltaE(k2)+Lambda(k1)*z - .5*cosfi*sqrt(.5*(1-z))...
   .*(kappa-(1+z)./(1-z)).*exp(-d*sqrt(2*Lambda(k1)*(1+z)))

is 0 ?

I am not sure it is possible to find an analytic form in terms of the other variables.

[Walter Roberson]

In the meantime here is some code for exploring the number of roots under various conditions.

The code executes with constant cosfi.

It loops over lambda values and d values. For each, it creates a 3D grid of deltaE, kappa, and z (the same z as mentioned above), and evaluates the function on that grid. It draws an isosurface in 3D of the deltaE / kappa / z values that make the function 0. Below that it draws a heat diagram of sorts of the number of zeros detected in deltaE vs kappa -- so the number of z values in that mix where there was a root (detected by way of a sign change.)

The looping creates one figure for each distinct lambda value used. Within each figure, subplots are created for the different d values used.

Unfortunately the isosurfaces are monochrome because they are all graphing the same constant, 0. (I did not check to see if that could be enhanced.)

With cosfi = -1, the maximum roots I detected was 3.

[Walter Roberson]

Experiments hint that for cosfi >= 0, there are some areas with no roots and otherwise the maximum is 1 root. For the negative cosfi I have tested, I have seen 0, 1, 2, and 3 roots at various locations. In the cases with multiple roots, the blob at high z value appears to be mostly flat-ish, but the blob at low z values is mostly oriented at a different angle and is generally reverse-C shape. That leads to 1 zero where the high z blob does not overlap with the reverse-C at all, and leads to 2 zeros where reverse-C shape has no high z blob above it or where the angle is such that only one of the two branches of the reverse-C is overhung by part of the high blob; and leads to 3 zeros where the high blob and two branches of the reverse-C shape are all present.

[Argumanh]

Dear Walter,

Thank you very much for your time and attention. The parameters that I want to be constant are "d" "cosfi" and "kappa" and the variables would be "Lambda" and "DeltaE". In addition to that, I want to know the Value or the magnitude of the roots. For example we have Eigenvalues for first root plotted in the 3D(eigroot2.m) I want to know the magnitude of the root which gave me that 3D plot.

Best,

Armin

[Walter Roberson]

In the eigroot3 code I posted yesterday, I do not attempt to find exact roots: the pair of lines

        T = sign(gFd);
        zc = sum(T(:,:,1:end-1) ~= T(:,:,2:end),3);

looks for sign changes along the z dimension.

How precise do you need to z values? At present the code does a grid search with numz samples taken over -1 to +1 . Is that good enough, or do you need the root to a number of decimal places? If you do then it would probably make the most sense to use the approximate zeros (detected by sign change) to create upper and lower bounds for a numeric root finder.