How to handle a variable inside matrix without using syms toolbox?

I have a application to have a variable inside a 2x2 matrix.

The matrix multiplication needs to be done with number 2x2 matrices with variables.

At the end, the final product matrix will have polynomial equatoins of the variable as its elements.

The roots of the polynomial equation at element (1,2), need to be evaluated.

The application need to be done without using syms toolbox.

Answer 1

[Alan Stevens]

Are you looking for somethig like this;

f = @(w) [1,-0.1409*w^2;0,1]*[1,0;1/286.48,1]*[1,-0.05493*w^2;0,1]...
    *[1,0;1/1793.55,1]*[1,-28.99*w^2;0,1];
% 
% The sixth order polynomial of f(1,2) can be written as 
%  a*w^6 + b*w^5 +...+ f*w + g;
%
% Let w take on values 0, 1, 2..6
% Then you would have 7 equations in the 7 unknowns, a, b etc.
%
% Now you have a system of linear equations that can be solved
% for the unknowns using simple matrix division.
%
% The resulting values are the coefficients of the polynomial
% the roots of which can be found using the roots function.


for w = 0:6
   P=f(w);
   V(w+1,1)=P(1,2);
   M(w+1,:) = [w^6,w^5,w^4,w^3,w^2,w,1]; 
end
coeffs = M\V;
r = roots(coeffs);
disp('Coefficients')

Coefficients

disp(coeffs')

   -0.0000    0.0000    0.0175    0.0000  -29.1858    0.0000         0

disp('Roots')

Roots

disp(r)

         0
 -195.4814
  195.4814
  -41.8216
   41.8216
    0.0000

% Check that the polynomial is, in fact, zero at the calculated values
disp('Check')

Check

for i = 1:7
   F = f(r(w));
   disp(F(1,2))
end

  -1.7339e-27

  -1.7339e-27

  -1.7339e-27

  -1.7339e-27

  -1.7339e-27

  -1.7339e-27

  -1.7339e-27

[Vinothkumar Sethurasu]

Thank you so much.

Great efforts.

[Vinothkumar Sethurasu]

Hello Alan,

Thanks for your method to find the coefficients.

It is working fine for the above shown example case.

But, i have some the cases providing polynomial equation after matrix multiplication which having variable in denominator.

The equation at P(1,2) is shown below,

The final equation need to be,

The method, we are working now is not producing the same co-efficients.

Do you have any suggestion on this topic?

Thank you

Vinothkumar S

[Alan Stevens]

Multiply every term through by the term in the denominator. Expand terms in brackets and collect like terms. Looks like you will get an eighth order polynomial. Use the same technique as before, but modified for the higher order.

[Vinothkumar Sethurasu]

Actually we don't know the polynomial equation (both numerator and denominator). It is an output from matrix multiplication.

I am sharing you the new function which providing the co-efficients,

f=@(w) [1 -w^2*0.1409;0 1]*[1 0;0.003490659 1]*[1 (-0.06734*w^2+2.26E-6*w^4)/(-4.053E-5*w^2+1);0 1]*[1 0;0.000557552 1]*[1 -w^2*28.99959505;0 1];
for w = 0:8
   P=f(w);
   V(w+1,1)=P(1,2);
   M(w+1,:) = [w^8 w^7,w^6,w^5,w^4 w^3,w^2,w,1]; 
end
coeffs = M\V;
r = roots(coeffs);

[Alan Stevens]

That seems to work ok. If you are only interested in the real roots, you could add

c = 0;
for i = 1:numel(r)
    if abs(imag(r(i)))<10^-10
        c = c+1;
        r1(c) = r(i);
    end
end
disp('Real roots')
disp(r1)

at the end.

[Vinothkumar Sethurasu]

No, the script shared by me providing different equation co-efficients than actual co-efficient from manual derivation.

The actual co-efficients need to provided from script also.

So, the method we used in the script is not providing correct co-efficients for the cases like equation P(1,2) having variables in denominator.

[Alan Stevens]

Try the following:

% You need to multiply through by the term in the denominator
% then set the result to zero (and hope the denominator term
% doesn't equal zero!). The coefficients I obtained didn't match 
% those you described.

% The results match those I obtained using a combination of
% symbolic maths expansion with a final stage of numerical
% calcuation.

k = @(w)-4.053E-5*w.^2+1;
f1 = @(w) [1 -w.^2*0.1409;0 1]*[1 0;0.003490659 1];
f2 = @(w) [k(w) (-0.06734*w.^2+2.26E-6*w.^4);0 k(w)];
f3 = @(w) [1 0;0.000557552 1]*[1 -w.^2*28.99959505;0 1];
f=@(w) f1(w)*f2(w)*f3(w); % Note: don't divide through by k(w)

for w = 0:8
   P=f(w);
   V(w+1,1)=P(1,2);
   M(w+1,:) = [w^8,w^7,w^6,w^5,w^4,w^3,w^2,w,1]; 
end
coeffs = M\V;
r = roots(coeffs);
c = 0;

format shortG
disp('Coefficients')

Coefficients

disp(coeffs')

   1.7972e-11   1.1295e-14  -1.2436e-06   9.3013e-13     0.018846   7.0012e-12      -29.208   2.8784e-12            0

format short
disp('Roots')

Roots

disp(r)

         0
 -219.3031
  219.3021
 -139.1257
  139.1262
   41.7830
  -41.7830
    0.0000

[Vinothkumar Sethurasu]

The co-efficient are same, but divided by a constant.

You have made it.!!

Awesome idea. Great efforts.

Thank you.

[Vinothkumar Sethurasu]

In the same application,

Z = @(w) [1 -w.^2*0.05493;0 1]*[1 0;0.003266128 1]*[1 -w.^2*0.012409;0 1];
a = @(w) Z(1,2);
b = @(w) Z(2,2);
f2 = @(w) [b(w) a(w);0 b(w)];

Can we use the polynomial at Z(1,2) & Z(2,2) like mentioned in script?

I am facing some issue in this script, Can you please guide me?

[Alan Stevens]

LIke this?

% Since your functions contain powers of w^2 I suggest you first
% set u = w^2 and solve for u.  Then the w solutions will be +/-u.

Z = @(u) [1 -u*0.05493;0 1]*[1 0;0.003266128 1]*[1 -u*0.012409;0 1];
% Note that Z(1,2) is a quadratic in u and Z(2,2) is linear in u
% so

for u = 0:2
   P = Z(u);
   if u<2
       V22(u+1,1) = P(2,2);
       M22(u+1,:) = [u, 1];
   end
   V12(u+1,1) = P(1,2);
   M12(u+1,:) = [u^2, u, 1]; 
end
coeffs12 = M12\V12;
r12 = roots(coeffs12);
w12 = sqrt(r12);
coeffs22 = M22\V22;
r22 = roots(coeffs22);
w22 = sqrt(r22);

disp('Roots12')

Roots12

disp([w12, -w12])

         0         0
  173.9176 -173.9176

disp('')
disp('Roots22')

Roots22

disp([w22 -w22])

  157.0779 -157.0779

[Vinothkumar Sethurasu]

Yes, ofcourse!

Here, i am not finding roots of the equation at Z(1,2) & Z(2,2).

I am going to feed it on function handle f2. The co-efficients are enough, i think.

Thank you.

[Vinothkumar Sethurasu]

Hello Alan,

I have, the application to replace the varibale in the function handle into value,

j=0.05493;
Z = @(u) [1 -u*j;0 1]; % J - used as a varible

Here, the j value keeps on changing in the loop. So, inside the function handle the vaibale 'J' need to be evaluvated & replaced with the current value.

the function need to be,

Z = @(u) [1 -u*0.05493;0 1]; % J - evaluated & Replced with the value

[Alan Stevens]

You can make Z a function of j as well:

Z = @(u,j) [1 -u*j;0 1]*[1 0;0.003266128 1]*[1 -u*0.012409;0 1];

However, I'm puzzled by the statement "the j value keeps on changing in the loop". To which loop do you refer? The only loop I see is that for u, and it doesn't make any sort of sense to change j within that loop!

[Vinothkumar Sethurasu]

To make things clear, an example loop shown here,

j=[0.1409 0.003490659 0.05493];
for i=1:3
Z{i} = @(u) [1 -u*j(i);0 1];
end

For each 'i' value, 'j' need to be evaluated & placed in the Z{i} function handle.

The results from the above loop will be,

Z{1}= @(u) [1 -u*j(i);0 1];
Z{2}= @(u) [1 -u*j(i);0 1];
Z{3}= @(u) [1 -u*j(i);0 1];

But, the results i need to be,

Z{1}= @(u) [1 -u*0.1409;0 1];
Z{2}= @(u) [1 -u*0.003490659;0 1];
Z{3}= @(u) [1 -u*0.05493;0 1];

[Alan Stevens]

You can still define Z as

Z = @(u,j) [1 -u*j;0 1]*[1 0;0.003266128 1]*[1 -u*0.012409;0 1];

Then put the whole thing in a loop (indexed by k, say) and when you call Z, use Z(u, j(k)).

Time for you to take over!

[Vinothkumar Sethurasu]

Hello Alan,

Thanks for your wornderful support.

I have facing some results deviation between using syms & function handle.

I am sharing the script below.

>> Function handle

f = @(u) [1 -u*0.112;0 1]*[1 0;(1/1057.68) 1]*[1 -u*0.0358;0 1]*[1 0;(1/9.358) 1]*[1 -u*0.0578;0 1]*[1 0;(1/2113.281) 1]*[1 -u*0.2436;0 1]*[1 0;(1/10015.548) 1]*[1 -u*33.607;0 1];
for u = 0:5
   P=f(u);
   V(u+1,1)=P(1,2);
   M(u+1,:) = [u^5,u^4,u^3,u^2,u,1]; 
end
coeffs = M\V;
r = roots(coeffs);
r=abs(r);
r=sort(sqrt(r));

>> With syms

syms u

f = [1 -u*0.112;0 1]*[1 0;(1/1057.68) 1]*[1 -u*0.0358;0 1]*[1 0;(1/9.358) 1]*[1 -u*0.0578;0 1]*[1 0;(1/2113.281) 1]*[1 -u*0.2436;0 1]*[1 0;(1/10015.548) 1]*[1 -u*33.607;0 1];
P=f(1,2);
R=vpasolve(P);
R=double(R);
R=sqrt(R);

The resulst from syms looks correct.

kindly give me a suggestion to avoid the resuls deviation in function handle.

[Walter Roberson]

You already asked that at https://www.mathworks.com/matlabcentral/answers/888942-getting-different-results-from-function-handle-syms-for-a-same-equation-how-to-avoid-it?s_tid=srchtitle and you were answered.