Find real positive roots of 9th order polynomial

3 次查看（过去 30 天）
T S Singh
T S Singh 2020-12-14
评论： T S Singh 2020-12-16
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.
  6 个评论
显示 4更早的评论
Mathieu NOE
Mathieu NOE 2020-12-15
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 : mathieu.noe@hutchinson.com

请先登录，再进行评论。

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

回答（1 个）

Matt J
Matt J 2020-12-15
编辑：Matt J 2020-12-15
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;
  6 个评论
显示 4更早的评论
Matt J
Matt J 2020-12-15
编辑：Matt J 2020-12-15
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
T S Singh 2020-12-16
As of now I don't have any other ways to generate the parameters.

请先登录，再进行评论。

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by