Eigen Function for modes of natural frequency

clear all
clc
syms lambda
syms x
syms i
rhoEpoxy = 1200;
Em = 3e+09;
R=5;
E_L=70e+9;
E_r=200e+9;
rho_L=2702;
rho_r=5700;
nu=3;
A_K=zeros(R+1,R+1);
A_M=zeros(R+1,R+1);
G=zeros(R+1,R+1);
H1_K=zeros(R+1,R+1);
H1_M=zeros(R+1,R+1);
H2_M=zeros(R+1,R+1);
for ri=1:R+1
    r = ri-1;
    for mi=1:R+1
        m = mi-1;
        fun1 = @(ksei1) (ksei1.^(r+m)).*((1-(exp(nu*ksei1)-1)./(exp(nu)-1))+E_r/E_L*(exp(nu*ksei1)-1)./(exp(nu)-1));
        %G(mi,ri)=integral(fun1,0,1);
        fun_h1=@(ksei2,s2) (-2*nu/(exp(nu)-1)*(E_r/E_L-1)*exp(nu*s2)).*(s2.^r).*(ksei2.^m)+...
            (nu^2/(exp(nu)-1)*(E_r/E_L-1)*exp(nu*s2)).*(s2.^r).*(ksei2.^(m+1))-...
            (nu^2/(exp(nu)-1)*(E_r/E_L-1)*exp(nu*s2)).*(s2.^(r+1)).*(ksei2.^m);
        fun_h1_lambda=@(ksei3,s3) (-1/6*((ksei3-s3).^3).*((1-(exp(nu*s3)-1)./(exp(nu)-1))+...
            rho_r/rho_L*(exp(nu*s3)-1)./(exp(nu)-1))).*(s3.^r).*(ksei3.^m);
        fun_h2_lambda=@(ksei4,s4) (1/6*(ksei4.^3).*((1-(exp(nu*s4)-1)./(exp(nu)-1))+...
            rho_r/rho_L*(exp(nu*s4)-1)./(exp(nu)-1))-1/2*(ksei4.^2).*s4.*((1-(exp(nu*s4)-1)./(exp(nu)-1))+...
            rho_r/rho_L*(exp(nu*s4)-1)./(exp(nu)-1))).*(s4.^r).*(ksei4.^m);
        
        options = {'RelTol', 1e-22, 'AbsTol', 1e-24};
        G(mi,ri) = integral(fun1, 0, 1, options{:});
        H1_K(mi,ri) = integral2(fun_h1, 0, 1, 0, @(ksei2) ksei2, options{:});
        H1_M(mi,ri) = integral2(fun_h1_lambda, 0, 1, 0, @(ksei3) ksei3, options{:});
        H2_M(mi,ri) = integral2(fun_h2_lambda, 0, 1, 0, 1, options{:});
        A_K(mi,ri)=G(mi,ri)+H1_K(mi,ri);
        A_M(mi,ri)=H1_M(mi,ri)+H2_M(mi,ri);
    end
end
sort(sqrt(eig(A_K,-A_M)))

ans = 6×1

2.8545 21.4957 63.7919 135.2966 242.8566 957.7791

For some reason my 4th, 5th 6th mode values are not correct they don't match the results. My first 3 mode values for vibration are correct adn they match the research paper. after 3rd mode no matter what boyndary condition I am using they don't match after 3rd mode. Can someone tell me why.

[Syed Abdul Rafay]

these are correct values

2.854

21.496

63.679

126.60

and these are my values

2.854

21.496

63.679

135

[Torsten]

Seriously: How should anybody not involved in the physical background of your problem be able to answer this ?

And your RelTol and AbsTol values are much too low as to any solver could be able to satisfy these error tolerances.

[Syed Abdul Rafay]

I have already tried different tolerance values but result is same. I have tried more than 100 different combinations. The hand calculations give correct 4th 5 th values but when transferred to code they change only first 3 values are correct rest diverge. The problem here is not the equations it has to do something with matlab implementations that's why I am asking question here.

Answer 1

[Rupesh]

Hi Syed Abdul Rafay,

It looks like the problem you're having with finding the higher mode values might be because of some tricky parts in the math calculations MATLAB does. To address the discrepancy in the 4th mode value in your vibration analysis problem, you can focus on two main modifications: refining the numerical integration technique and enhancing the accuracy of the eigenvalue problem solver. Below are some thoughts and potential adjustments that might guide you towards resolving the issue:

1] Numerical Integration Precision

Here the goal is to improve the precision for numerical inetgration, especially for integrands that exhibit complex behavior over the integration range. A useful strategy is to segment the integration range into parts where the integrand behaves more uniformly. This approach allows for more accurate integration over each segment, which is particularly beneficial for capturing the nuances of higher modes.

 
% Assuming 'fun1' is your integrand and it behaves differently across the range 
splitPoint = 0.5; % This is an illustrative value; adjust based on your integrand's behavior 
options = {'RelTol', 1e-6, 'AbsTol', 1e-8}; % Fine-tune these tolerances as needed 
integralValue = integral(fun1, 0, splitPoint, options{:}) + integral(fun1, splitPoint, 1, options{:});  

2] Eigenvalue Solver Accuracy

The choice and configuration of the eigenvalue solver are crucial for accurately determining the vibration modes. When dealing with large matrices or seeking a subset of eigenvalues, the eigs function can offer a more efficient and potentially more accurate approach than the standard eig function, especially for higher modes.

% Assuming A_K and A_M are your stiffness and mass matrices, respectively 
numModes = 6; % Number of modes you're interested in 
opts.tol = 1e-6; % Adjust the tolerance to improve accuracy 
[eigenvectors, eigenvalues] = eigs(A_K, A_M, numModes, 'smallestabs', opts); 
sorted_eigenvalues = sort(sqrt(diag(eigenvalues)));  

3] Matrix Conditioning

The condition of the matrices used in the eigenvalue problem significantly affects the accuracy of the results. Ensuring that the stiffness and mass matrices (A_K and A_M) are well-conditioned .This might involve revisiting how these matrices are constructed and ensuring they are symmetric when theoretically expected to be so.

% Check for symmetry and condition number of A_K and A_M 
if ~isequal(A_K, A_K.') || ~isequal(A_M, A_M.') 
    warning('Matrices A_K or A_M are not symmetric. Check their construction.'); 
end 
if cond(A_K) > 1e10 || cond(A_M) > 1e10 
    warning('One of the matrices is poorly conditioned. Consider revising its construction.'); 
end  

You can also refer to below document for the usage of above-mentioned functions .

https://in.mathworks.com/help/matlab/ref/eigs.html

Hope this helps!