Having symbolic in a Matrix

Hi guys!

In this code, I am constructing two matrices via a loop and also asks for some user input so the matrices can change according to the userinputs.

In matrix Z2, I want to have a symbolic 'p' inside the matrix because I want to calculate what p value would make the determinant to zero(I know how to do this). But it is giving an error msg:

Error using symengine

Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.

Which means I cannot put symbolic as variable in matrices. Is there a way that I can have this 'p' inside the matrix Z2? Thank you!

% define variables
L = 1; % in metre
w = 10 * 10^-3; % width in m
h = 3.5 * 10^-3; % height in m
I = (1/12) * w * (h^3) % area moment of inertia about bending axies


% ask user for # of element
prompt = "What is the element #? ";
n = input(prompt);

% coifficient matrix for stiffness
L = 1 / n;
cons1 = [12, 6*L, -12, 6*L; 
        6*L, 4*L^2, -6*L, 2*L^2;
        -12, -6*L, 12, -6*L;
        6*L, 2*L^2, -6*L, 4*L^2]; 
    
cons2 = [36, 3*L, -36, 3*L; 
        3*L, 4*L^2, -3*L, -L^2;
        -36, -3*L, 36, -3*L;
        3*L, -L^2, -3*L, 4*L^2];
    
% initialize stiffness matrix with 0  
Z1 = zeros(2*n+2,2*n+2);
Z2 = zeros(2*n+2,2*n+2);

% ask user for E values to compute EI/L^3
syms p
K1 = zeros (1, n);
K2 = -p / (30*L);
for i = 1:n
    prompt = "Enter E value in GPa (10^9) of " + i + "th element: ";
    K1(1, i) = 10^9 * (input(prompt)) * I / (L^3);
end    
    
% stiffness matrix Z1
for i=1:n
    %1st row elements
    Z1(2*i-1, 2*i-1) = Z1(2*i-1, 2*i-1) + K1(1, i) * cons1(1,1);
    Z1(2*i-1, 2*i) = Z1(2*i-1, 2*i) + K1(1, i) * cons1(1,2);
    Z1(2*i-1, 2*i+1) = Z1(2*i-1, 2*i+1) + K1(1, i) * cons1(1,3);
    Z1(2*i-1, 2*i+2) = Z1(2*i-1, 2*i+2) + K1(1, i) * cons1(1,4);
    
    %2nd row elements
    Z1(2*i, 2*i-1) = Z1(2*i, 2*i-1) + K1(1, i) * cons1(2,1);
    Z1(2*i, 2*i) = Z1(2*i, 2*i) + K1(1, i) * cons1(2,2);
    Z1(2*i, 2*i+1) = Z1(2*i, 2*i+1) + K1(1, i)* cons1(2,3);
    Z1(2*i, 2*i+2) = Z1(2*i, 2*i+2) + K1(1, i) * cons1(2,4);
    
    %3rd row elements
    Z1(2*i+1, 2*i-1) = Z1(2*i+1, 2*i-1) + K1(1, i) * cons1(3,1);
    Z1(2*i+1, 2*i) = Z1(2*i+1, 2*i) + K1(1, i) * cons1(3,2);
    Z1(2*i+1, 2*i+1) = Z1(2*i+1, 2*i+1) + K1(1, i) * cons1(3,3);
    Z1(2*i+1, 2*i+2) = Z1(2*i+1, 2*i+2) + K1(1, i) * cons1(3,4);
    
    %4th row elements
    Z1(2*i+2, 2*i-1) = Z1(2*i+2, 2*i-1) + K1(1, i) * cons1(4,1);
    Z1(2*i+2, 2*i) = Z1(2*i+2, 2*i) + K1(1, i) * cons1(4,2);
    Z1(2*i+2, 2*i+1) = Z1(2*i+2, 2*i+1) + K1(1, i) * cons1(4,3);
    Z1(2*i+2, 2*i+2) = Z1(2*i+2, 2*i+2) + K1(1, i) * cons1(4,4);
end

% stiffness matrix Z2
for i=1:n
    %1st row elements
    Z2(2*i-1, 2*i-1) = Z2(2*i-1, 2*i-1) + K2 * cons2(1,1);
    Z2(2*i-1, 2*i) = Z2(2*i-1, 2*i) + K2 * cons2(1,2);
    Z2(2*i-1, 2*i+1) = Z2(2*i-1, 2*i+1) + K2 * cons2(1,3);
    Z2(2*i-1, 2*i+2) = Z2(2*i-1, 2*i+2) + K2 * cons2(1,4);
    
    %2nd row elements
    Z2(2*i, 2*i-1) = Z2(2*i, 2*i-1) + K2 * cons2(2,1);
    Z2(2*i, 2*i) = Z2(2*i, 2*i) + K2 * cons2(2,2);
    Z2(2*i, 2*i+1) = Z2(2*i, 2*i+1) + K2 * cons2(2,3);
    Z2(2*i, 2*i+2) = Z2(2*i, 2*i+2) + K2 * cons2(2,4);
    
    %3rd row elements
    Z2(2*i+1, 2*i-1) = Z2(2*i+1, 2*i-1) + K2 * cons2(3,1);
    Z2(2*i+1, 2*i) = Z2(2*i+1, 2*i) + K2 * cons2(3,2);
    Z2(2*i+1, 2*i+1) = Z2(2*i+1, 2*i+1) + K2 * cons2(3,3);
    Z2(2*i+1, 2*i+2) = Z2(2*i+1, 2*i+2) + K2 * cons2(3,4);
    
    %4th row elements
    Z2(2*i+2, 2*i-1) = Z2(2*i+2, 2*i-1) + K2 * cons2(4,1);
    Z2(2*i+2, 2*i) = Z2(2*i+2, 2*i) + K2 * cons2(4,2);
    Z2(2*i+2, 2*i+1) = Z2(2*i+2, 2*i+1) + K2 * cons2(4,3);
    Z2(2*i+2, 2*i+2) = Z2(2*i+2, 2*i+2) + K2 * cons2(4,4);
end

Answer 1

[Paul]

Hi Felis,

I didn't run the code because I don't know what inputs to provide. Does changing Z2 to a sym object solve the problem?

Z2 = sym(zeros(2*n+2,2*n+2));

[John D'Errico]

Better is this use of zeros:

Z2 = zeros(2*n+2,2*n+2,'sym');

where you can specify the class directly. As such it is surely more efficient than creating an array of many zeros, and then forcing syms to convert them all into syms. Alternatively, this would probably also be efficient:

Z2 = repmat(sym(0),2*n+2,2*n+2);

though I've not tested them to see.

However, my guess is if n is at all large (5 is probably large enough) then the problem will have no analytical solution, since the result will be a matrix that is a function of p, where p appears in lots of places in the array. And that will mean the determinant will be a polynomial in p of order at least 5, and so no analytical solution will exist.

rootfinders could still be used of course, so vpasolve would suffice.

Regardless, without knowing the various inputs, I could do no more either.

[Felis Flora]

Hi Paul and John,

I tried "Z2 = repmat(sym(0),2*n+2,2*n+2);" and it worked for having the Z2 matrix! Have not thought of changing the class or object.

The vpa(solve()) will be fine to solve for a near zero, have not try it yet but I will update when I can work on it tomorrow.

This is a code for Finite element analysis for buckling, so n is something around 20 - 100 :D.

I need to compare different 'p' values based on E inputs to select the best material combination for largest 'p' so I do not need an exact answer, numerical approximation will be okay.

Thank you guys so much! happy holidays

[Walter Roberson]

Be careful about the fact that there might be multiple solutions. If it does become a polynomial then vpasolve will find them all, but for non polynomial it will only find one .

[Felis Flora]

thanks for the declearation. I think it will have some complex number and real number solutions and i ll try to make physical sense from these answers

[Paul]

Hi John,

I agree with your recommendation in principle, but was curious enough to test.

timeit(@() (zeros(500,500,'sym')))

ans = 0.9729

timeit(@() sym(zeros(500,500)))

ans = 1.0486

tic
for ii = 1:5
    X = sym(zeros(500,500));
end
toc

Elapsed time is 5.274078 seconds.

tic
for ii = 1:5
    X = zeros(500,500,'sym');
end
toc

Elapsed time is 5.417881 seconds.

Maybe the parser is smart enough to realize that sym(zeros(...)) doesn't require creation of double and then convert to sym? Or maybe that's in the noise and most of the time is spent on some other aspect of the operation.

[Walter Roberson]

sym(zeros) is slightly slower in my tests, but the margin of error is enough that if you exclude the first point (which could be out due to time spent compiling) then you pretty much cannot tell the two apart reliably.

25 iterations of each test was right at the boundary of all that I could fit into 55 seconds

N = 25;
t1 = zeros(N,1);
t2 = zeros(N,1);
for K = 1 : N; start = tic; sym(zeros(500,500)); stop = toc(start); t1(K) = stop; end
for K = 1 : N; start = tic; zeros(500,500,'sym'); stop = toc(start); t2(K) = stop; end

[mean(t1(2:end)), mean(t2(2:end))]

ans = 1×2

1.0335 1.0353

plot([t1, t2]);
legend({'sym(zeros)', 'zeros(sym)'})