Quartic solver using Embedded matlab function

Hello,

I have used an embedded matlab function in simulink to solve a general quartic. Following is the code that I have written in this block

function t = my_quartic (C1,C2,C3,C4,C5)
%eml.extrinsic('p');

p=[C1 C2 C3 C4 C5];
t=roots(p);
%function t=roots(p)

I have used constant blocks to provide the values of C1=3,C2=6,C3=-123, C4=-126 AND C5=1080. I have used the "to file" block as a sink to store the values of t.

Now for the values of coeficients shon above, the roots of t should come out to be 5,3,-4,-6. But when I run the simulation, I get the following error: Data 't' (#23) is inferred as a variable size matrix, while its specified type is something else.

I have tried changing the properties of t in the editor, but its not working. Can someone help me on this? I will be greatly thankful for any comments.

Best Regards, Amir

Answer 1

[Mike Hosea]

The problem here is that ROOTS always returns a variable-size result. That's because it ignores leading zeros (or extremely small entries, relative to the other coefficients). So you will find that both codegen and MATLAB versions of ROOTS behave this way. I can think of three easy options here. One is to use variable-size signals in Simulink. You can do this by navigating to the MATLAB Function Block (or Embedded MATLAB Block) in the Model Explorer, clicking on the output signal, checking the "Variable size" check-box on the Size line of the General tab and entering [4,1] for the size.

Another alternative, is to stick with fixed-size calculations by using eig(compan(p)) instead of roots(p), since you probably do not care about dealing efficiently with leading and trailing zeros and such.

Another alternative is

function t = my_quartic (C1,C2,C3,C4,C5)
r = roots([C1 C2 C3 C4 C5]);
t = complex(zeros(4,1));
t(1:numel(r)) = r(:);

That uses variable-size arrays internally, but the output signal is fixed-size.

Do not consider "extrinsic" in this matter. It won't help. -- Mike

[Amir Khan]

Thanks a lot Mike! its working like a charm!

[Amir Khan]

Hello Mike,

Can you please guide me on this. Beased on the suggestions, I have changed the properties in the model explorer and slight changes to the codes are below. The thing is that I want only real roots and want to remove the NaNs in the final matrix. But Im getting the following error when I run this code:

function r = my_quartic (C1,C2,C3,C4,C5)

R=roots([C1 C2 C3 C4 C5]);

r=R(R==real(R));

r(any(isnan(r),2),:) = [];

the error is "Removing elements from a matrix is not supported."

I have to remove the rows containg Nans so that I can use the final data obtained after this for further processing. Can you suggest a way around? Thanks very much for your help!

[Mike Hosea]

I don't have time to test this right now, but try something like

function r = my_quartic (C1,C2,C3,C4,C5)

R=roots([C1 C2 C3 C4 C5]);

r = zeros(4,1);

n = 0;

thresh = 2*eps;

for k = 1:4

if ~isnan(R(k)) && abs(imag(R(k))) < thresh

n = n + 1;

r(n) = real(R(k));

end

end

r = r(1:n);

[Mike Hosea]

Just tested that, and it seems to work. This also works:

function r = my_quartic (C1,C2,C3,C4,C5)

R = roots([C1 C2 C3 C4 C5]);

thresh = 2*eps;

r = real(R(~isnan(R) & abs(imag(R)) < thresh));

[Mike Hosea]

The numerical analyst in me wants to make one more comment. Although I couldn't bear to compare equality of imag(R) to zero, I did resist getting too fancy with "thresh". However, you may need to account for very large or very small roots, so you might want to compare abs(imag(R)) to some multiple of eps(real(R)) instead of to a constant. Point is that comparing to zero is naive, but how fancy you want to get in diagnosing "real" roots is up to you. Normally we would argue that r=1e100+1e80i is numerically indistinguishable from r=1e100 in IEEE double precision. Likewise r=1e-100+1e-80i might be regarded as indistinguishable from r=1e-80i rather than r=1e-100.

[Harshal]

Thanks Mike!! This works perfect.

Answer 2

[Walter Roberson]

Before the t = roots(p) call, add the line

t = zeros(1,4);

Yes, this will assign a value to t and then overwrite it on the very next line, but doing this will give enough information to the embedded coder to know the size to expect for t.

By the way: why are you specifying p to be extrinsic, as it is not an external function? I could understand if you had said that 'roots' was extrinsic ?

[Amir Khan]

Thanks!

I have changed the code to following:

function t = my_quartic (C1,C2,C3,C4,C5)

p=[C1 C2 C3 C4 C5];

t = zeros(1,4);

t=roots(p);

Still Im getting the same error i.e.

"Data 't' (#23) is inferred as a variable size matrix, while its specified type is something else."

[Walter Roberson]

Ummm, I'm not sure, but as an experiment, try changing the last line to

t(1:4) = roots(p);

[Walter Roberson]

By the way, I am not certain that you need to construct p explicitly. I think you _might_ be able to use

function t = my_quartic (C1,C2,C3,C4,C5)

t = zeros(1,4);

t(1:4) = roots([C1,C2,C3,C4,C5]);

Though it would not surprise me if it turned out that you needed to declare roots as extrinsic at some point.