MATLAB: How to find polynomial roots in Simulink

polynomialrootssimulinksolve

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

Best Answer

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];
Related Question