How to find polynomial roots in Simulink?
13 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
ozan eren
le 8 Mai 2020
Commenté : ozan eren
le 24 Juin 2020
Hi,
I am trying to solve a 4th order polynomial equation in Simulink. I need to solve the equation by using Simulink blocks. The coefficients are calculated in Simulink blocks as well and I need to find the roots of this equations for each iteration.
To be more clear, I am using a for iterator block in Simulink and iterate 1000 times, in this for loop certain coefficients (C1,C2, C6 etc) are calculated from math expressions in each iteration. And I neet to find the roots of the polynomial whose coefficients are those calculated C values and I need find the roots for each iteration. I do not want to use script or built-in functions, but solve this problem with Simulink blocks if possible. Below I share a piece of m script which does exactly what I need to do in Simulink.
Any help will be appriciated, thanks in advance.
r = roots([C5 (-C1+2*C4-C6) (3*C2-3*C5) (C1-2*C3+C6) -C2]);
r = r(r==conj(r)); % discard complex-conjugate roots
r = r(r>0); % discard negative roots
0 commentaires
Réponse acceptée
Walter Roberson
le 22 Juin 2020
By using roots() on symbolic variables, you can get four closed form expressions for the roots. They occur in pairs, A+/-B and P+/-Q where B and Q are sqrt(), so by detecting whether the sqrt() involve imaginary quantities you can eliminate conjugate pairs as you wanted.
The closed form expressions are long, but if you are willing to run some code overnight, you can matlabFunction() writing to a 'File' with 'optimize' turned on. That produces a serquence of simple MATLAB code to calculate the expression.
You indicated that you do not want to use scripts, just Simulink blocks. Each of the statements from the optimized code is simple enough to make it practical to implement as a small number of Math blocks (or similar).
It would be a bit tedious to create these expressions manually, but it it only about 100 of them. And you would get exact solutions without iteration.
t2 = C2.*3.0;
t3 = C3.*2.0;
t4 = C4.*2.0;
t5 = C5.*3.0;
t9 = 1.0./C5;
t13 = sqrt(3.0);
t14 = sqrt(6.0);
t6 = -t3;
t7 = -t4;
t8 = -t5;
t10 = t9.^2;
t11 = t9.^3;
t15 = C2.*t9;
t12 = t10.^2;
t16 = C1+C6+t6;
t17 = C1+C6+t7;
t18 = t15.*1.2e+1;
t19 = t2+t8;
t20 = -t18;
t21 = t17.^2;
t22 = t17.^3;
t24 = t9.*t16;
t25 = t9.*t19;
t26 = (t9.*t17)./4.0;
t32 = t10.*t16.*t17.*3.0;
t35 = (t10.*t16.*t17)./4.0;
t37 = (t10.*t17.*t19)./2.0;
t23 = t21.^2;
t27 = t10.*t21.*(3.0./8.0);
t28 = (t11.*t22)./8.0;
t36 = -t35;
t38 = t11.*t19.*t21.*(3.0./4.0);
t39 = (t11.*t19.*t21)./1.6e+1;
t29 = -t27;
t30 = -t28;
t31 = t12.*t23.*(3.0./2.56e+2);
t33 = t12.*t23.*(9.0./6.4e+1);
t40 = -t39;
t34 = -t33;
t41 = t25+t29;
t47 = t24+t30+t37;
t53 = t15+t31+t36+t40;
t42 = t41.^2;
t43 = t41.^3;
t48 = t47.^2;
t54 = t53.^2;
t55 = t53.^3;
t58 = t41.*t53.*(4.0./3.0);
t59 = t41.*t53.*7.2e+1;
t44 = t42.^2;
t45 = t43.*2.0;
t46 = t43./2.7e+1;
t49 = t48.^2;
t50 = t48.*2.7e+1;
t52 = t48./2.0;
t56 = t55.*2.56e+2;
t57 = t43.*t48.*4.0;
t61 = t42.*t54.*1.28e+2;
t62 = t41.*t48.*t53.*1.44e+2;
t51 = t49.*2.7e+1;
t60 = t44.*t53.*1.6e+1;
t63 = t51+t56+t57+t60+t61+t62;
t64 = sqrt(t63);
t65 = t13.*t64.*3.0;
t66 = (t13.*t64)./1.8e+1;
t67 = t45+t50+t59+t65;
t69 = t46+t52+t58+t66;
t68 = sqrt(t67);
t70 = t69.^(1.0./3.0);
t72 = 1.0./t69.^(1.0./6.0);
t71 = t70.^2;
t74 = t41.*t70.*6.0;
t76 = t14.*t47.*t68.*3.0;
t73 = t71.*9.0;
t75 = -t74;
t77 = -t76;
t78 = t20+t32+t34+t38+t42+t73+t75;
t79 = sqrt(t78);
t80 = 1.0./t78.^(1.0./4.0);
t81 = t42.*t79;
t83 = t53.*t79.*1.2e+1;
t84 = t73.*t79;
t85 = t71.*t79.*-9.0;
t86 = (t72.*t79)./6.0;
t88 = t41.*t70.*t79.*1.2e+1;
t82 = -t81;
t87 = -t86;
t89 = -t88;
t90 = t76+t82+t83+t85+t89;
t91 = t77+t82+t83+t85+t89;
t92 = sqrt(t90);
t93 = sqrt(t91);
t94 = (t72.*t80.*t92)./6.0;
t95 = (t72.*t80.*t93)./6.0;
GG = [t26+t87-t94; t26+t87+t94; t26+t86-t95; t26+t86+t95];
4 commentaires
Walter Roberson
le 23 Juin 2020
Maple finds a more compact sequence:
t2 = C1-2*C4+C6;
t3 = 1/C5;
t5 = 1/4*t2*t3;
t6 = C2-C5;
t8 = t2^2;
t9 = C5^2;
t10 = 1/t9;
t13 = 3*t6*t3-3/8*t8*t10;
t14 = t13^2;
t15 = C2*t3;
t17 = 3^(1/2);
t18 = t14*t13;
t20 = C1-2*C3+C6;
t24 = 1/t9/C5;
t30 = t20*t3-1/8*t8*t2*t24+3/2*t6*t2*t10;
t31 = t30^2;
t34 = t8^2;
t35 = t9^2;
t37 = t34/t35;
t40 = t20*t2*t10;
t43 = 3*t6*t8*t24;
t45 = t15+3/256*t37-1/4*t40-1/16*t43;
t46 = t45^2;
t49 = t31^2;
t53 = t14^2;
t60 = (144*t13*t31*t45+128*t14*t46+4*t18*t31+256*t45*t46+16*t45*t53+27*t49)^(1/2);
t61 = t17*t60;
t65 = t13*t45;
t67 = 1/18*t61+1/27*t18+1/2*t31+4/3*t65;
t68 = t67^(1/3);
t69 = t13*t68;
t71 = t68^2;
t76 = t14-12*t15-6*t69+9*t71-9/64*t37+3*t40+3/4*t43;
t77 = t76^(1/2);
t78 = t67^(1/6);
t79 = 1/t78;
t81 = 1/6*t77*t79;
t83 = 12*t45*t77;
t85 = 9*t71*t77;
t86 = t14*t77;
t88 = 12*t69*t77;
t89 = 6^(1/2);
t96 = (3*t61+2*t18+27*t31+72*t65)^(1/2);
t98 = 3*t89*t30*t96;
t100 = (t83-t85-t86-t88+t98)^(1/2);
t102 = t76^(1/4);
t103 = 1/t102;
t105 = 1/6*t100*t79*t103;
t109 = (t83-t85-t86-t88-t98)^(1/2);
t112 = 1/6*t109*t79*t103;
F(1) = t5-t81-t105;
F(2) = t5-t81+t105;
F(3) = t81+t5-t112;
F(4) = t81+t5+t112;
Plus de réponses (1)
Sai Teja Paidimarri
le 18 Juin 2020
Modifié(e) : Walter Roberson
le 22 Juin 2020
Hi Ozan Eren,
You can find roots in Simulink and Matlab as well.
For a polynomial c5 x^4 + c4 x^3 +c3 x^2 +c2 x + c1, roots r can be found out using below steps
x = c4 + r*c5;
y = c3 + r*x;
z = c2 + r*y;
r value for which c1+r*z is zero is root of equation.
Above equations are from the concept of synthetic division method.
Counter, logical operator and constants can be used in case of Simulink implementation.
For or while loop can be used in case of Matlab Implementation.
3 commentaires
Sai Teja Paidimarri
le 23 Juin 2020
Hi ozan eren
First you need to find x,y,z because z depends on y , y depends upon x .Then condition associated with z.After that you can use r (using counter in simulink and any loop in matlab) to check the condition assocaited with z. if condition satisfies you can move r value in roots (declare roots as an array).
You can refer synthetic division method for understanding the above formulae.
Voir également
Catégories
En savoir plus sur Loops and Conditional Statements dans Help Center et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!