Having symbolic in a Matrix

1 view (last 30 days)
Felis Flora
Felis Flora on 26 Nov 2022
Commented: Walter Roberson on 28 Nov 2022
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

Accepted Answer

Paul
Paul on 26 Nov 2022
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));
  6 Comments
Walter Roberson
Walter Roberson on 28 Nov 2022
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)'})

Sign in to comment.

More Answers (0)

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!