Accuracy of Integral and Integral2 functions

I am using integral functions but I am told the integral functions do have some inaccuracy when used. That's why my answers are not same with the paper I am trying to validate. let me share my code:

format long g
clear all
clc
syms lambda
syms x
syms i
rhoEpoxy = 1200;
Em = 3e+09;
R=5;
E_L=Em;
E_r=Em;
rho_L=rhoEpoxy;
rho_r=rhoEpoxy;
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);
        H1_K(mi,ri)=integral2(fun_h1,0,1,0,@(ksei2) ksei2);
        H1_M(mi,ri)=integral2(fun_h1_lambda,0,1,0,@(ksei3) ksei3);
        H2_M(mi,ri)=integral2(fun_h2_lambda,0,1,0,1);
        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

1.0e+00 * 3.51601526867379 22.034800684001 61.7194163744804 128.398091475146 226.847068659209 1094.88143674641

this is my code why I run this code i get following results;

3.516

22.035

61.719

128.4

but the correct results are;

2.4256

18.6031

55.1791

109.5748

[Torsten]

I don't believe that such a big difference between the results can be due to inaccuracies of integral2. There must be some substantial difference between the two computations.

[Syed Abdul Rafay]

But the same code give very close results to other paper I am validating the results with

Answer 1

[Walter Roberson]

Given that code, the eigenvalues come out the 3.516 and so on. It isn't a matter of low accuracy: I switched to high accuracy here.

format long g
clear all
clc
syms lambda
syms x
syms i
Q = @(v) sym(v);

rhoEpoxy = Q(1200);
Em = Q(3)*Q(10)^9;
R = 5;
E_L=Em;
E_r=Em;
rho_L=rhoEpoxy;
rho_r=rhoEpoxy;
nu = Q(3);
A_K=zeros(R+1,R+1,'sym');
A_M=zeros(R+1,R+1,'sym');
G=zeros(R+1,R+1,'sym');
H1_K=zeros(R+1,R+1,'sym');
H1_M=zeros(R+1,R+1,'sym');
H2_M=zeros(R+1,R+1,'sym');
syms ksei1 ksei2 ksei3 ksei4 s2 s3 s4
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) = int(fun1(ksei1),ksei1,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);
        %H1_K(mi,ri) = integral2(fun_h1,0,1,0,@(ksei2) ksei2);
        H1_K(mi,ri) = int(int(fun_h1(ksei2,s2), s2, 0, ksei2), ksei2, 0, 1);
        %H1_M(mi,ri)=integral2(fun_h1_lambda,0,1,0,@(ksei3) ksei3);
        H1_M(mi,ri) = int(int(fun_h1_lambda(ksei3,s3), s3, 0, ksei3), ksei3, 0, 1);
        %H2_M(mi,ri)=integral2(fun_h2_lambda,0,1,0,1);
        H2_M(mi,ri) = int(int(fun_h2_lambda(ksei4,s4), ksei4, 0, 1), s4, 0, 1);
        
        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
dA_K = double(A_K);
dA_M = double(A_M);
sort(sqrt(eig(dA_K,-dA_M)))

ans = 6×1

1.0e+00 * 3.51601526878612 22.0347977911742 61.7162927762873 128.389334759042 223.551343192217 1006.01205007505

[Syed Abdul Rafay]

Look they are the same. I tried with the changing the accuracy of the results but I got the same value like you.

[Torsten]

@Walter Roberson wanted to show you that the results are not different with higher accuracy, as was your supposition.

[Syed Abdul Rafay]

I tried and now my first 3 modes are correct but my 4 th and 5 th mode are not correct. If i increase R number of mode shapes increases as well but only first 3 modes are correct rest have too much discrepency. Can you look why. Like 4 th mode is 121 and my 4th mode is 128. The 5th mode is 247 and mine is 563. This only start from 4 th mode first 3 are always correct.