Error using symengine - Division by zero.

6 次查看（过去 30 天）
WONG Weng Tuck
WONG Weng Tuck 2021-9-24
回答： Walter Roberson 2021-9-24
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

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2021-9-24
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 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by