Which is the best method to find zero ? and choose the max value and plot a graph comproving the zeros.

For a given parameters

ss=0.6;
T=0.1;
g=1;
b=0.5;

 f1 =  @(r)  1 - (b./r).^2 - (g^-2)*((2/15)*(ss)^9 *(1./(r - 1).^9 - 1./(r + 1).^9 - 9./(8*r).*(1./(r - 1).^8 - 1./(r + 1).^8)) -(ss)^3 *(1./(r-1).^3 - 1./(r+1).^3 - 3./(2*r).*(1./(r-1).^2 - 1./(r+1).^2)));

I want to find the r value when

f1 == 0

So, which is the best method ? findzero, fminbnd, newtonraphson ? fsolve ?

i try

a1=  fzero(f1,1.00000001)
assume(r, 'real');
t=vpasolve(1 - (b/r)^2 - (g^-2)*(2/15*[(s/R)]^9 *(1/(r - 1)^9 - 1/(r + 1)^9 - 9/(8* r) *(1/(r - 1)^8 - 1/(r + 1)^8)) - [(s/R)]^3 *(1/(r -1)^3 - 1/(r + 1)^3 - 3/(2* r)* (1/(r - 1)^2 - 1/(r + 1)^2))) == 0,r,[-Inf,Inf]);
t

options=optimset('MaxIter',1e3,'TolFun',1e-10);
t1=fminbnd(func,0.1,15,options);
t1

but i dont know if it is correct. I want to plot a graph r x f1 to see the zeros.

Any help how i can find zeros exactly ?

Answer 1

[dpb]

>> fsolve(f1,sqrt(2))

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

ans =
    1.4121
>> f1(ans)
ans =
  -2.4100e-12
>> r=linspace(1.2, 1.6);
>> figure
>> plot(r,f1(r))
>> hold on
>> line([r(1) r(end)],[0 0])
>> ylim([-20 20])
>>

I started by plotting f1(r) for r=linspace(1+eps,100,20); and exploring from there to get a reasonable initial guess for fzero Looking at the algebra first would probably lead to a way to guess that as a closer starting point rather than just brute force, but it works, too, often... :)

[Lucas Pollito]

Hi, thans a lot ! but, i want all roots of f1, and use vpasolve to look for it. COrrect ?

But, in all roots, i want to choose the max real value, in this case 1.4421, or i want to choose the min value. f1 has many roots. How can i choose the max ou min value of roots ? how can i plot f1 x r that show me all the roots ??

Answer 2

[Walter Roberson]

You cannot get exact values for r. The exact solutions involve roots of a 10th degree polynomial that does not factor. There is no closed form solutions for such a polynomial.

[Lucas Pollito]

but

   vpasolve does not find exactly all roots ?

Answer 3

[Star Strider]

It’s not easy, but you can find many of them using fsolve and a loop. I plotted these first to get the approximate range (-10,10):

ss=0.6;
T=0.1;
g=1;
b=0.5;

f1 =  @(r)  1 - (b./r).^2 - (g^-2)*((2/15)*(ss)^9 *(1./(r - 1).^9 - 1./(r + 1).^9 - 9./(8*r).*(1./(r - 1).^8 - 1./(r + 1).^8)) -(ss)^3 *(1./(r-1).^3 - 1./(r+1).^3 - 3./(2*r).*(1./(r-1).^2 - 1./(r+1).^2)));

r = linspace(-10, 10, 1000);

r_e = -50:1.2:50;
for k1 = 1:length(r_e)
    re_r_sol(k1) = fsolve(f1, r_e(k1));
end
re_r = uniquetol(re_r_sol, 0.05);
for k1 = 1:length(r_e)
    cx_r_sol(k1) = fsolve(f1, 1i*r_e(k1));
end
[cx_ru,ia,ic] = uniquetol(abs(cx_r_sol), 0.05);
cx_r = cx_r_sol(ia);
fprintf(1,'\nReal roots:\n')
fprintf(1,'\t%+8.3f\n', re_r)
fprintf(1,'\nComplex roots:\n')
fprintf(1,'\t%+8.3f %+8.3fi\n', [real(cx_r); imag(cx_r)])

Real roots:
    -2.528
    -1.412
    -0.576
    +0.576
    +1.412
    +2.528

Complex roots:
    +0.335   +0.377i
    +0.000  +79.550i
    +0.000  -84.050i

I am not entirely confident of these results, because I would expect the complex roots to be complex conjugates of each other. Using a symbolic solver might be more productive.