Error using symengine - Division by zero.

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

Answer 1

[Walter Roberson]

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