How to find polynomial roots in Simulink?

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

Answer 1

[Walter Roberson]

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

[ozan eren]

Dear Walter,

Thank you so much for your detailed answer and effort. I tried the method and works perfect (minor round-off and sign issues). Could you please provide the name of this method if it has any or is it just a optimzed version of roots function?

Once again thanks a lot for your elaborated reply and time. Best regards.

[Walter Roberson]

I used roots() and then matlabFunction writing to file with optimization on. It took a fair amount of time to optimize, perhaps a few hours.

[Walter Roberson]

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;

[ozan eren]

Dear Walter,

I already tested and started to built the previous solution in Simulink, but this one is much more compact than earlier one as you mentioned and I will switch to this solution. I was stuck in this point almost a month, this is a huge help.

Once more thanks a lot for your effort and best regards.

Answer 2

[Sai Teja Paidimarri]

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.

[ozan eren]

Hi,

Thanks a lot for your kind answer. But I could not understand how to find r. If I understand your explanation correctly, first I should find r value then x,y and z values. Are the roots of the equation r, x, y and z? It will be great, if you could explain a bit further.

Thanks in advance.

[Sai Teja Paidimarri]

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.

[ozan eren]

Dear Sai Teja,

Thanks a lot for your time and kind answers. But I think the above solution will be better for my application. In any case thank you for your time and best regards.