Accuracy of Integral and Integral2 functions
1 回表示 (過去 30 日間)
古いコメントを表示
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)))
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
2 件のコメント
Torsten
2024 年 2 月 9 日
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.
回答 (1 件)
Walter Roberson
2024 年 2 月 10 日
編集済み: Walter Roberson
2024 年 2 月 10 日
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)))
3 件のコメント
Torsten
2024 年 2 月 12 日
@Walter Roberson wanted to show you that the results are not different with higher accuracy, as was your supposition.
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!