Find real positive roots of 9th order polynomial

Hi,

I am trying to find the real positive roots of a 9th order polynomial equation using roots . My program is given below

clear all
close all
clc
%%%% Parameters %%%%
gamma=0.0783;
zai=0.025;
fo=0.25;
f1=0.01;
om_data=.01:.001:0.9;
n=length(om_data);
cnt=1;
cnt1=1;
for ii=1:1:n
    A=[];
    sol=[];
    om=om_data(ii);
    P=[25*gamma^3 0 -20*gamma^2*om^2 -15*gamma^2*fo +4*gamma*om^2*(om^2+4*zai^2) ...
    16*gamma*fo*om^2 3*gamma*(2*f1^2-3*fo^2) -4*fo^2*om^2*(om^2+4*zai^2) ...
    4*fo^2*om^2 -fo^3];                 % Polynomial Coefficients 
    sol=roots(P);                       % Roots of polynomial 
    A=sol(imag(sol)==0 & sol>=0);       % Select real and positive roots
    no_A=length(A);
    if no_A==1
        plot(om,A,'.');
        hold on
    elseif no_A==3
        om_same = om*ones(1,3);
        plot(om_same,A,'.')
        hold on
    elseif no_A==5
        om_same2 = om*ones(1,5);
        plot(om_same2,A,'.')
        hold on
    end
end

The plot which I obtained is different from the actual answer. The correct plot is (I. Kovacicetal et. al. Int. J. Non-linear Mech. 43, 2008)

Q.1 Is there any limitation with the roots command or is there any problem with my code.

Q.2 For fo=0.2 and f1=0.1 the plot should look like this

That is, there are some parametric region where we obtain multiple roots. But my program is giving a differen plot.

Any help please

Thank You.

[Mathieu NOE]

hello

until now, I have not found anything in your code that is incorrect

are you sure about your polynomial coefficients ? (I don't have the refernce publication)

[T S Singh]

@ Mathieu Yes, I have cross checked the polynomial coeffcients (I also have check its mathematical derivations). I can share the paper if it is ok with you.

[Mathieu NOE]

yes please

[T S Singh]

@ Mathieu I have not shared any research papers before in this forum, so can you please guide me :-) I don't know wether it is Ok to share the paper here or shall I have to send you personally through mail?

[Mathieu NOE]

hi

I suppose both options work as soon as this publication is not covered by any NDA or other restrictions

if you prefer , you can send it to : [e-mail]

[T S Singh]

@ Mathieu Thanks. Please check your mail

Answer 1

[Matt J]

The roots of a polynomial become increasingly sensitive to small variations in the polynomial coefficients as the order of the polynomial gets large, and order 9 is a bit high. Are you sure gamma and the other parameters affecting the coefficients are accurately given, to the maximum precision possible, compared to the original paper? In particular, does gamma really only have 3 significant digits?

Also, instead of,

om=om_data(ii);

it would be more accurate to have,

om=0.01+ii*0.001;

[T S Singh]

@ Matt As far as the parameters are concerned they are case for a non-dimenionalized equation of motion. So the only way I can implement in my code is choose the value of parameters given in the paper. Am I missing something?

[Matt J]

Also, you cannot general expect that imag(sol)==0 will be ideally satisfied. You need to use tolerances, such as

 A=sol(abs(imag(sol))<=abs(real(sol))*1e-10 & real(sol)>=0);

[T S Singh]

@ Matt Just to cross check the solutions I did tried with pseudo arclength continuation method coded in MATLAB I did found a different plot (although not the correct one :-(). I think the possible reason behind this is the tolerance which you are mentioning here. But I did cross check the solutions and it was giving the reals postive roots only!

[T S Singh]

@ Matt even if I put the tolerance the roots obtained are still the same.

In another program using numerical technique (`pseudo arclength continuation `) I obtained the following plots.

fo=0.25 and f1=0.01

fo=0.2 and f1=0.1

The first plot is almost matching with the paper plot(Fig. 2(a)) but in the second plot all the roots are not obtained in some parameter range (there might be a constraint in the code which I need to figure it out). Although, in the second plot it is not showing the whole solutions but the range of Ao is same with that of the paper (Fig. 3(a)). This observation is not obtained when I am using roots that's why I was currious whether there is some limitation in this command.

[Matt J]

My best guess is that the paper only published their parameters (e.g., gamma) to 4 decimal places, whereas more precision is needed to obtain ttheir actual plots. Do you have any way of regenerating the input parameters from first principles? Otherwise, maybe contact the authors.

[T S Singh]

As of now I don't have any other ways to generate the parameters.