Having symbolic in a Matrix
조회 수: 2 (최근 30일)
이전 댓글 표시
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
댓글 수: 0
채택된 답변
Paul
2022년 11월 26일
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
Paul
2022년 11월 28일
Hi John,
I agree with your recommendation in principle, but was curious enough to test.
timeit(@() (zeros(500,500,'sym')))
timeit(@() sym(zeros(500,500)))
tic
for ii = 1:5
X = sym(zeros(500,500));
end
toc
tic
for ii = 1:5
X = zeros(500,500,'sym');
end
toc
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
2022년 11월 28일
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))]
plot([t1, t2]);
legend({'sym(zeros)', 'zeros(sym)'})
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Assumptions에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!