How to find polynomial roots in Simulink?
3 views (last 30 days)
Show older comments
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 Comments
Accepted Answer
Walter Roberson
on 22 Jun 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 Comments
Walter Roberson
on 23 Jun 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;
More Answers (1)
Sai Teja Paidimarri
on 18 Jun 2020
Edited: Walter Roberson
on 22 Jun 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 Comments
Sai Teja Paidimarri
on 23 Jun 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.
See Also
Categories
Find more on Loops and Conditional Statements in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!