Error using symengine - Division by zero.
38 views (last 30 days)
Show older comments
WONG Weng Tuck
on 24 Sep 2021
Answered: Walter Roberson
on 24 Sep 2021
I am trying to solve a vibration by using the Rayleigh-Ritz method with an admissable function.
The error was as follows -
Error using symengine
Division by zero.
Error in sym/subs>mupadsubs (line 160)
G = mupadmex('symobj::fullsubs',F.s,X2,Y2);
Error in sym/subs (line 145)
G = mupadsubs(F,X,Y);
Error in ca2draft3 (line 54)
MassMatrix(rows, columns) = subs(Mij, [i,j], [rows, columns])
Tried using eps and still got the error. Could anyone please assist?
Code is as follows -
clear all
close all
clc
YoungMod = 2e12; %young's modulus
WingDensity = 300; %density of aircraft wing
syms i j y
%length of wing 44 m (span)
i(i == 0) = eps;
j(j == 0) = eps;
y(y == 0) = eps;
AreaWing = 0.3 - (y/440) + ((y^2)/9680); %area dependence on y (outboard distance) from wing, expanded out
SecMomAreaWing = 0.001 - (y/132000) + ((y^2)/3872000) - ((y^3)/85184000); %varies as above, expanded out
phi_i_y = ((y^2)/(44^2))*cos((i*pi*y)/44); %define assumed admissable function, scaled appropriately to make it dimensionless
phi_j_y = ((y^2)/(44^2))*cos((j*pi*y)/44);
phi_doubleprime_i_y = diff(phi_i_y,y,2); %differentiate admissable function twice
phi_doubleprime_j_y = diff(phi_j_y,y,2);
phi_i_inbdEng = ((8.8^2)/44^2)*cos(0.2*i*pi); %inboard engine at 8.8 m from centerline
phi_i_midEng = ((17.6^2)/44^2)*cos(0.4*i*pi); %mid engine, 17.6 m from centerline
phi_i_outerEng = ((26.4^2)/44^2)*cos(0.6*i*pi); %outboard engine, 26.4 m from centerline
phi_j_inbdEng = ((8.8^2)/44^2)*cos(0.2*j*pi);
phi_j_midEng = ((17.6^2)/44^2)*cos(0.4*j*pi);
phi_j_outerEng = ((26.4^2)/44^2)*cos(0.6*j*pi);
%elements of the mas matrix is the integral of the beam plus the three
%engines
MassinbdEng = 4100*(phi_i_inbdEng)*(phi_j_inbdEng); %mass of engines
MassmidEng = 4100*(phi_i_midEng)*(phi_j_midEng);
MassouterEng = 4100*(phi_i_outerEng)*(phi_j_outerEng);
%stiffness of the wing (ignoring the engines)
StiffnessWing = 2e12*(SecMomAreaWing)*(phi_doubleprime_i_y)*(phi_doubleprime_j_y);
MassWing = 300*(AreaWing)*(phi_i_y)*(phi_j_y);
%elements of mass matrix is Mij and stiffness matrix, Kij
Mij = int(MassWing,y,0,44) + MassinbdEng + MassmidEng + MassouterEng;
Kij = int(StiffnessWing,y,0,44);
ii = 4; %change this parameter here
MassMatrix = zeros(ii);
StiffMatrix = zeros(ii);
tic
for rows = 1:ii;
for columns = 1:ii;
MassMatrix(rows, columns) = subs(Mij, [i,j], [rows, columns])
StiffMatrix(rows, columns) = subs(Kij, [i,j], [rows, columns])
end
end
Systemmatrix = inv(MassMatrix)*StiffMatrix;
[V,D] = eig(Systemmatrix);
[D_sorted, ind] = sort(diag(D),'ascend');
V_sorted = V(:,ind);
nat_freq_one = sqrt(D_sorted(1))
nat_freq_two = sqrt(D_sorted(2))
mode_shape_one = V_sorted(:,1)
mode_shape_two = V_sorted(:,2)
% un-% mode_shape_1, mode_shape_2, nat_freq_3 and nat_freq_4 for ii exceeding 3!
nat_freq_three = sqrt(D_sorted(3))
nat_freq_four = sqrt(D_sorted(4))
mode_shape_three = V_sorted(:,3)
mode_shape_four = V_sorted(:,4)
toc
0 Comments
Accepted Answer
Walter Roberson
on 24 Sep 2021
int(MassWing,y,0,44)
That integral ends up with a division by (i-j) . When rows == cols that is a division by 0.
You can get around that by subs() for i, and then taking the limit to j
syms i j y
i(i == 0) = eps;
j(j == 0) = eps;
y(y == 0) = eps;
i, j, y are indefinite variables at that point. They do not "equal" anything at that point.
Furthermore, you later substitute in row and column numbers for i and j, so i and j are always positive.
Here is the new loop:
for rows = 1:ii
for columns = 1:ii
MassMatrix(rows, columns) = simplify(limit(subs(Mij, i, rows), j, columns));
StiffMatrix(rows, columns) = simplify(limit(subs(Kij, i, rows), j, columns));
end
end
0 Comments
More Answers (0)
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!